Document Type : Research Paper
Authors
 Amir Hosein Tahershamsi ^{1}
 Ahad Ghaemi ^{} ^{} ^{2}
 Mansour Shirvani ^{3}
^{1} M.S. Student, Department of Chemical Engineering, Iran University of Science and Technology (IUST), Tehran, Iran
^{2} Assistant Professor of Department of Chemical Engineering, Iran University of Science and Technology (IUST), Tehran, Iran
^{3} Associate Professor of Department of Chemical Engineering, Iran University of Science and Technology (IUST), Tehran, Iran
Abstract
In this study, liquidliquid extraction process in a Kuhni extraction column was modeled and simulated. A nonequilibrium dynamic model was developed for modeling liquidliquid extraction processes based on a ratebased model. The model equations are inclusive of partial and ordinary differential equations which were discretized in column height direction. The population balance model was used for the calculation of droplet size distribution in the dispersed phase and the column hydrodynamic parameters. The equations were solved simultaneously through the finite difference method and the numerical method of lines. Experimental data on a bench scale Kuhni extraction column was used to evaluate the simulation results. The average correlation coefficient error of the mean diameter of the dispersed phase and mass transfer in various operating conditions are less than 2% 4 % respectively. A comparison between the experimental data and the simulation results proves the better productivity of the presented nonequilibrium dynamic model.
Keywords
Solvent extraction is a process, in which a particle of a fluid solvent transmits to another fluid solvent; this phenomenon is called distribution. The component distributed between the two phases is considered as the distributed or transmitted component. These two fluids are immiscible. One of them is usually water, and the other is a lowpolarity organic fluid, which plays the extractor role. The extractor is usually dissolved in a suitable fluid called diluent, which plays a role in the extractor physical property improvement (Thakur, 2000). The range of available extractor types and variables are complex to perform a certain vast responsibility; thus, choosing the suitable extractor is a complex decision and must be performed with precaution and in the best mental condition in the absence of the previous experience with the same system on pilot scale. A useful primary evaluation of extractors has been conducted by Morello and Poffenberger (Morello and Poffenberger, 1953). A classification of most important extractors and their applications has been introduced by Godfrey and Slater (1994). Nowadays, stroke or finned extraction columns are being used due to their high power operation and high separation efficiency. Kuhni extraction column is one of the finned column in which many researches have been conducted. This column is covered with worn blades on a central axis and fixed perforated plates through which levels of this column are separated from each other. The blades have an eddy current pattern that is somewhat similar to creating a rotating disk column. Also, perforated panels are used between the columns not only to provide an axial possibility, but also to separate the levels more completely than OldshueRushton columns. These decreases reverse mixing, especially for the scattered phase. Kumar et al. studied changes in the dispersed phase profile and the outburst points in this column, and they evaluated the effect of mixing speed on drop size for wateroxylene system changes in various heights of this column (Kumar et al., 1986). They witnessed that at a fixed mixing speed, changes in height do not really affect droplet sizes, while increasing the mixing speed causes recovering droplet breakage and consequently decreases the droplet size. Hufnagl et al. presented a differential model for Kuhni column for tolueneacetonewater mass transmission from continues to diffusion phase (Hufnagl et al., 1991). They used the selected data for fatigue, droplet diameter, current volume, and components of the continuous phase percentage to acquire experimental parameters.
Weinstein et al. applied a differential model to Kuhni column with tolueneacetonewater system and 40% chaos in input sporadic phase (Weinstein et al. 1998). Their numerical solution included a general balance equation in the sporadic phase, which was considered as dynamic effects. Mansur et al. evaluated zinc concentration changes in a chemical reaction including D2EHPA solution (Mansur et al. 2003). Hasseine et al. studied the application of population balance model for the prediction of the dispersed phase and droplet diameter in Kuhni extraction columns and the rotating disk (2005). They evaluated the mean droplet size and the amount of dispersed phase through the two columns of water and toluene experimentally and using modeling and population balance equations. Gomes et al. modeled an industrial Kuhni column using population balance equations (Gomes et al., 2006). They found that increasing column height raised the amount of dispersed phase in the column due to an increase in shear forces. Oliveira et al. evaluated the effect of operating parameters consisting of mixing speed, dispersed and continuous phases of the droplet size distribution debit, and dispersed phase supply along a Kuhni extraction column for water and Exxsol D80 system (Oliveira et al., 2008). They found out that with increasing mixing speed the reduction of droplet size was witnessed, and, in conclusion, droplet size distribution became more symmetric and narrower, and it was shifted toward the left side of the diagram; they also found out that due to the shear forces resulting by blade rotation and the inner components of the column, droplet size distribution at the top of the column is narrower and more symmetric than the one at the bottom of the column. Rodea et al. investigated the undesirable phenomenon of outburst in Kuhni columns (Rodea et al., 2013). Attarakiha et al. used the OPOSPM method using droplet population balance and studied it in a small Kuhni column (Attarakiha et al., 2013). The results of their modeling were very close to experimental data. Artur and Marcelo achieved results coherent to the laboratory results with dynamic modeling of Kuhni extraction column and mixing it with droplet population balance equation (Attarakiha et al., 2013). According to the conducted studies in this work, a nonequilibrium dynamic model based on ratebased model was developed for modeling extraction process in a Kuhni column. Population balance equations were used for the hydrodynamic modeling of Kuhni column considering breakage effects and droplet coalition for two chemical systems, including nbutyl acetatewateracetone with medium surface tension and toluenewateracetone with high surface tension simultaneously. The population balance method is recommended due to its accuracy, its low cost, and its extent of equations for determining various behaviors of chemical processes. The model results were evaluated using a bench scale Kuhni experimental data.
2. Population balance model in Kuhni extraction column
Researchers in different sciences evaluate molecules, particles, and cells with regards to existing rules in their science in order to achieve a certain distribution of material proportions and its variation with time and location. The distribution must be explainable for all materials regarding the general laws specific to a certain science, and this general law, which is included in a general mutual concept, is called population balance equation(Attarakih and Bart, 2014). In the liquidliquid extraction process, the population balance model is used to explain hydrodynamic interactions and mass transmissions which are explainable based on the macroscopic nature of the dispersed and continuous phases. Droplet population balance is a powerful method to forecast extraction column behavior which has been used in many works during recent years. When this method is solved, a great number of differential equations are acquired, and droplet size distribution is obtained by solving these equations. For hydrodynamic simulation, the column is divided into elements with a certain height and droplet classification, and a certain diameter is created due to mass displacement, breakage, and droplet coalition. In this model, the basic variable, droplet size distribution, is shown by p(z,d,t), and location fatigue can be obtained by Equation 1 (Casamatta and Vogelpohl, 1985).
(1) 
The basic population balance in a dynamic mode is explained by Equation 2 (Modes et al., 1999).
(2) 
In Equation 2, p(z,d) is production and consumption statement at the same time, which includes breakage effects, coalition, nucleation, and droplet growth. Equation 2 demonstrates the p(z,d) phrase.
(3) 
In which PB^{+}(z,d) and PB^{}(z,d) respectively show the increase and decrease in droplet numbers due to droplet breakage; also, the phrase PC^{+}(z,d) and PC^{}(z,d) respectively shows the increase and decrease in droplet number due to droplet coalition (Godfrey and Slater, 1994).
(4) 

(5) 
Equation 5 is multiplied by 0.5 because particle collision with (d_{1},d_{2}) size is not different from particle collision with (d_{1},d_{2}) size. Considering the fact that p(d) is a breakage possibility, r(z,d) is called collision frequency obtained by an equation achieved by the provided case. Breakage frequency is also dependent on droplet existence time. P(d) is breakage possibility of a droplet with diameter d, and r(z,d) is breakage frequency obtained by Equation 7 (Steiner et al., 1988).
(6) 

(7) 

(8) 

(9) 
The coalition rate has been obtained using Equation 10 (Wang et al., 2003).
(10) 
To calculate child droplet size distribution, beta equation is usually used although beta equation is very complex in practice. The following equation is based on the hypothesis that if any droplet breakage does not occur, the impact energy will be absorbed in the droplet as much as possible. Then, this energy turns into surface energy in order to make the best interface area possible. This means that the creation possibility of interface area is more suitable for the child droplets with a bigger droplet size distribution. Daughter droplet size distribution is acquired by beta distribution based on the mother droplet diameter, d_{0} (Schmidt, 2006).
(11) 
Here, the mean of daughter droplet is calculated using the following equation.
(12) 
Critical diameter for every drop beginning to break out is obtained using the following equation.
(13) 
Also, Modes et al. presented these two parameters, i.e. term destruction and failure of the process as follows (Modes et al., 1999):
(14) 

(15) 
In order to describe the hydrodynamic problem, the population balance model considering the failure of the coalition drops were used simultaneously. Also, the effect of axial mixing in the dispersed phase was discarded. In conclusion, the main population balance equation will change as follows:
(16) 
The final form of the developed model can be obtained by substitution of Equations 4, 5, 14, and 15 in Equation 16 as follows:
(17) 

Number of the droplets in volume unit is as defined by:
(18) 
in which, V(d) is droplet volume. To solve Equation 17, it is important for the equation to be independent of the droplet diameter. Hence, droplet diameter range was classified into 9 classes with equal difference:
(19) 
Now we can rewrite p(z,d,t) as given by:
(20) 
After the independence of the equation of the droplet diameter, discretization is performed on the premises of location. By discretization on the premises of location, a series of partial differential equations appear, which will be solved. Boundary conditions at the bottom of the tower are obtained by the LogNormal distribution function (Oliveira et al., 2008).
(21) 
According to Equation 21, the amounts of and were corrected by using experimental data by minimizing sum of squared errors for this extraction column:
(22) 

(23) 
The population balance equations were solved simultaneously using the algorithm shown in Figure 1.
Entering system physical properties

Entering P(n)

Building network for location

Building network for droplets

Calculating interaction functions

Simultaneously solving normal differential equations to obtain droplet distribution function 
Next Step

Figure 1
Population balance equation solution algorithm.
3. Process modeling
The first step of designing and simulating a process is obtaining an accurate mathematical model for the system under study. Functional model accuracy is a function of the complexity of the target system and computer facilities used. Since all the data in solving the equation reduce the accuracy of the results, minimum simplifying assumptions must be used. In this section, by considering the shortcomings of the provided models, a dynamic model based on population balance to better describe mass transfer is provided to predict mass transfer behavior of liquidliquid extraction processes in Kuhni extraction column. The following assumptions are used to provide a developed model:
ü Changes in the radius of the average velocity and the density of matter are not considered;
ü Phase mass transfer rate depends on the concentration of the solute;
ü Dispersed phase droplets are considered spherical distributed in nine droplet diameter classes;
ü Physical properties of the column are considered as variables.
ü Mass transfer is considered in one direction and only from dispersed phase to continuous phase.
The schematic of the model is presented in Figure 2.
Figure 2
Model developed for mass transfer.
The leading dynamic mixing method was used to model the column mass transfer. In this method, each component of the dispersed phase follows an independent trajectory, and the effects of the appeared transmission are in the pattern of droplet size distribution. Each of the droplets of the dispersed phase velocity is determined with a specific velocity which is dependent on droplet diameter. Unlike reverse mixing method, the leading mixing method is more dependent on scattering phase. Figure 3 shows differential mass balance related to different classes of droplets.

d_{1} 
d_{i} 
d_{N} 








dz 
… 
… 
Figure 3
Differential mass balance model for the direct incorporation.
The overall mass balances for continuous and dispersed phase according to Figure 3 are as follows (Artur et al., 2013):
(24) 

(25) 
Considering that different droplet classes influence mass balances, the following equation is used in the continuous phase mass balance (Artur et al., 2013):
(26) 
Changes made in the scattering phase equation for the continuous phase and mass balances based on each droplet set are true. Mass balance for the dispersed phase includes all j droplet classes. We replace Equation 26 in the mass transfer equations:
(27) 

(28) 
4. Numerical solution of the model equations
The process model was a set of partial and ordinary differential equations. The equations were discretized in the axial direction using the method of lines and finite difference methods resulting in a system of coupled ordinary differential and algebraic equations. The algorithm of model equations solution is shown in Figure 4. All of the equations of the stages were solved simultaneously. The simulation of the process was preformed while different numbers of stages were applied. In addition, it was tried to calculate more accurate results for mass transfer by pairing the balance equations of mass transfer equations. Since mass transfer coefficients in the column is dependent on droplet size distribution, it is expected that the results will be better correlated with experimental data and have fewer error comparing to other models.
The average of correlation coefficient error was supposed to be the criterion for choosing the number of the stages which resulted in the lowest error. The equation employed to calculate the average of correlation coefficient error is given by:
(29) 
where, p is the number of process variables, and (R_{Exp.})_{i }and (R_{sim.})_{i} are correlation coefficient for experimental and simulation process variable i respectively. R_{ave} is average error percentage for the process.
Building a network for droplets 
Building a network for time and space 
For each droplet set: 1Calculating hydrodynamic parameters through population balance 2Calculating mass transfer coefficients 3Considering droplet breakage and coalition 
Solving the resulting differential equations for the concentration distribution simultaneously 
For each stage of the column: 1Calculating density in the dispersed and continuous phases 2A new calculation of the dispersed phase 
Next time step 
Entering system physical properties and the initial and boundary conditions 
Figure 4
Algorithm for solving the equations of mass transfer.
5. Kuhni column and phases properties
Kuhni extraction column used in this study is a contactor on the semiindustrial scale, and the active part of the column is made of glass to be resistant against corrosion. Kuhni column has one axis rotated by an electromotor. The geometric properties of the column and physical properties of the liquidliquid systems are presented in Tables 1 and 2.
Table 1
Geometric properties of the column.
Dimension 
Unit 
Column component 
0.700 
m 
Column height 
0.117 
m 
Column diameter 
0.055 
m 
Height of each stage 
10 
 
Number of stages 
27 
 
Number of holes per tray 
40 degrees 
 
Hole angle comparing to others 
0.002 
m 
Plate thickness 
0.00752 
m^{2} 
Active volume of the column 
Table 2
Physical properties of liquidliquid systems.
nButyl acetatewater 
Toluenewater 
Unit 
Physical properties 
998.2 
998.2 
Kg m^{3} 
Continuous phase density 
880.9 
865.2 
Kg m^{3} 
Dispersed phase density 
1.0274 
0.9630 
mPa.s 
Continuous phase viscosity 
0.734 
0.584 
mPa.s 
Dispersed phase viscosity 
14.1 
36.0 
mN s^{1} 
Interfacial tension 
6. Results and discussion
The model equations including mass balance, momentum balance, and population balance equations were solved simultaneously using numerical techniques. In the following, the hydrodynamic and mass transfer results are indicated. In Figure 5, the range of mean diameter as a function of rotor speed is demonstrated for nbutyl acetatewater and toluenewater systems. Considering that the operating condition is constant, mean droplet diameter size is bigger in toluenewater than the one in nbutyl acetatewater. The main reason for this occurrence is related to surface tension between the phases. Since the interfacial tension in toluenewater is almost twice as large as interfacial tension between nbutyl acetatewater, this causes an increase in drag force and resistance to the motion of the dispersed phase. The effect of mixing speed on the mean droplet size in the toluenewater system (system with a high interface tension) is more than that in nbutyl acetatewater system because droplet breakage in system with a medium interface tension is limited. Moreover, the diagram has a logical theme and the mean size of droplets decreases with an increase in rotor speed.
Weinstein et al. considered only average diameter to solve equations by solving their differential model, and they solved the equations using this diameter. But in this work, the population balance equations were employed. The droplets were classified in 9 classes according to the assumptions. The modeling results were consistent with the experimental results and the changes in physical properties of the system are fixed because the column length is low.
In Table 3, mean droplet diameter and their deviation with experimental data for the two chemical systems are presented. The results show that with an increase in the phase flow rate, the mean droplet diameter rise. It shows that increasing the flow rate of both phases enhances the coalescences of drops. Additionally, the results indicate that raising rotor speed drops mean droplet diameter. Flow rate rise in integrated phase has poor effects on the increase of droplet mean size, and it could lead to an increase in retention time and droplet collisions. Moreover, it could result in a decrease in the relative velocity of the droplets.
Figure 5
Mean diameter versus rotor speed for nbutyl acetatewater and toluenewater systems.
It is observed that along the column height, droplet breakage overcomes their coalition (Figure 6); the reduced surface tension causes droplet breakage and reducing the size of the larger droplets. On the other hand, in the dispersed phase, the amount of circular and rotational currents leads to more stability, increasing the time of coalition, thereby reducing the amount of droplet coalition. Thus, the peak of the diagram leans toward left, which shows higher breakage of droplets and shrinkage of the droplet mean diameter. Furthermore, according to the fact that the interface tension of the toluenewater system is more than that of nbutyl acetatewater system, higher power is needed to break the droplets in this system to achieve a class of droplets with a smaller diameter.
Table 3
Experimental and simulation results of the droplet mean diameter.
nButyl acetatewater 

N (rpm) 
Q_{c}(lit/hr.) 
Q_{d} (lit/hr.) 
Experimental d_{32} (mm) 
Simulation d_{32} (mm) 
Error 

150 
20 
24 
1.4901 
1.5017 
0.7784 

28 
24 
1.6099 
1.6468 
2.2920 

32 
24 
1.6541 
1.6214 
1.9769 

36 
24 
1.7097 
1.7020 
0.4503 

150 
24 
20 
1.4901 
1.5017 
0.7784 

24 
28 
1.6099 
1.6468 
2.2920 

24 
32 
1.6541 
1.6214 
1.9769 

24 
36 
1.7097 
1.7020 
0.4503 

Toluenewater 

180 
18 
24 
1.6286 
1.6580 
1.8052 

30 
24 
1.7765 
1.7950 
1.0413 

36 
24 
1.8013 
1.8300 
1.5932 

42 
24 
1.9868 
1.9950 
0.4127 

180 
24 
18 
1.5205 
1.5012 
1.2693 

24 
30 
1.9416 
1.9180 
1.2154 

24 
36 
2.3874 
2.3514 
1.5079 

24 
42 
2.6761 
2.6448 
1.1696 

a) 

b) 
Figure 6
Volume distribution of droplets versus the droplet diameter and the length of the column: a) nbutyl acetatewatersystem and b) toluenewater system.
Oliveira et al. found that by increasing the speed of mixing, a reduction in droplet size was seen, and droplet size distribution became narrower and more symmetrical and passed to the left of the graph. As indicated in the graph, by increasing the speed of mixing, the droplet size lessens and shifts to the left. Because of the small diameter of the column, the changes in the radial average speed and concentration can be ignored.
Figure 7 presents the influence of rotor speed on the interface tension of the two existing phases. By increasing the rotor speed, and thus increasing the shear forces, mother droplets break into child droplets and droplet coalition becomes minimum. Reducing the droplet size raises the amount of the dispersed phase in the extraction column. Moreover, due to a rise in interfacial tension, droplet coalition is enhances and consequently results in bigger droplets. In addition, relative phase velocity increases, which could lead to a reduction in the amount of the dispersed phase. Furthermore, it can be seen that higher rotor speeds influence the dispersed phase in the toluenewater system (with a higher interfacial tension) because droplet breakage in this system is more effective, thereby leading to an increase in the hold up of dispersed phase.
Figure 7
Effect of two rotor speeds on the amount of hold up for nbutyl acetatewater and toluenewater systems.
Comparison of the presented model results with the previous model results is presented in Table 4; it is obvious that the presented model is more rigorous than previous models.
Table 4
Comparison of the presented model results with other models.
Reference 
Column type 
Deviation percentage average 
Model 
(2007) 
Structure packing 
8.0 
Morales et al. 
(2012) 
fixedbed extractor 
20.0 
Lorena et al. 
(2006) 
RDC extractor 
15.0 
Schmidt et al. 

Kuhni extractor 
4.0 
This work 
7. Conclusions
In this study, the droplet breakage and the coalition and its impact on the number of droplets in the Kuhni extraction column were modeled and simulated by considering the balance equation. The most important parameters affecting the droplet breakage are rotor speed, interface tension, and dispersed and continuous phase flow rates. An increase in rotor speed causes stronger shear forces and thus results in the breakage of droplets; therefore, dispersed phase holdup increases. The holdup of the dispersed phase decreased by increasing the surface tension and droplet diameter. Increased surface tension and droplet diameter raise speed but reduces the relative phases of the dispersed phase. The results showed that droplet coalition plays a major role in population balance equation, and it was increased by increasing the flow rates of the dispersed and continues phases. The comparison of droplet mean diameter changes and the dispersed phase holdup shows that the effect of rotor speed in the nbutyl acetatewateracetone system (with a low interface tension) is less than that in the toluenewateracetone system (with a high interface tension). The results indicate that the simulation results of the process are in good agreement with the experimental data; therefore, using the presented model is recommended for modeling industrial Kuhni column. Furthermore, this model is suggested for modeling other types of extraction columns, especially RDC and structure columns. Finally, an investigation on droplet mean diameter and mass transfer coefficients based on nonspherical droplets is proposed.
Nomenclature
a 
: Interfacial area (m^{2}/m^{3}) 

d_{32} 
: Sauter mean drop diameter (m) 

E 
: Axial mixing coefficient (m^{2}/s) 

h 
: Coalition rate 

K 
: Overall mass transfer coefficient (m/s) 

N 
: Rotor speed (1/s) 

N_{A} 
: Mass transfer flux (1/s) 

P 
: Breakage possibility 

Q 
: Flow rate of the continuous or dispersed phase (m^{3}/s) 

t 
: Time (s) 

U 
: Velocity (m/s) 

V_{s} 
: Slip velocity (m/s) 

V_{t} 
: Terminal velocity of droplet (m/s) 

We 
: Weber number 

: Weber number critical 

x_{m} 
: Droplet diameter daughter 

Greek letters 

: Beta distribution function 


: Density (kg/m^{3}) 

: Viscosity (Pa.s) 

: Interfacial tension (N/m) 

Subscripts 

c 
: Continuous phase 

d 
: Dispersed phase 

x 
: xphase (continuous phase in present case) 

y 
: yphase (dispersed phase in present case) 

Superscripts 

nBuacwa 
: nbutyl acetateacetonewater 

Toacwa 
: Tolueneacetonewater 

* 
: Equilibrium value 

o 
: Inlet to column 