Document Type : Research Paper


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


In this study, liquid-liquid extraction process in a Kuhni extraction column was modeled and simulated. A non-equilibrium dynamic model was developed for modeling liquid-liquid extraction processes based on a rate-based 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 non-equilibrium dynamic model.


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 low-polarity 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 Oldshue-Rushton 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 water-o-xylene 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 toluene-acetone-water 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 toluene-acetone-water 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 D-80 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 non-equilibrium dynamic model based on rate-based 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 n-butyl acetate-water-acetone with medium surface tension and toluene-water-acetone 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 liquid-liquid 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).



The basic population balance in a dynamic mode is explained by Equation 2 (Modes et al., 1999).



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.



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





Equation 5 is multiplied by 0.5 because particle collision with (d1,d2) size is not different from particle collision with (d1,d2) 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).









The coalition rate has been obtained using Equation 10 (Wang et al., 2003).



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, d0 (Schmidt, 2006).



Here, the mean of daughter droplet is calculated using the following equation.



Critical diameter for every drop beginning to break out is obtained using the following equation.



Also, Modes et al. presented these two parameters, i.e. term destruction and failure of the process as follows (Modes et al., 1999):





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:



The final form of the developed model can be obtained by substitution of Equations 4, 5, 14, and 15 in Equation 16 as follows:




Number of the droplets in volume unit is as defined by:



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:



Now we can rewrite p(z,d,t) as given by:



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 Log-Normal distribution function (Oliveira et al., 2008).



According to Equation 21, the amounts of  and  were corrected by using experimental data by minimizing sum of squared errors for this extraction column:






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
















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):





Considering that different droplet classes influence mass balances, the following equation is used in the continuous phase mass balance (Artur et al., 2013):



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:





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:



where, p is the number of process variables, and (RExp.)i and (Rsim.)i are correlation coefficient for experimental and simulation process variable i respectively. Rave is average error percentage for the process.


Building a network for droplets

Building a network for time and space

For each droplet set:

1-Calculating hydrodynamic parameters through population balance

2-Calculating mass transfer coefficients

3-Considering droplet breakage and coalition

Solving the resulting differential equations for the concentration distribution simultaneously

For each stage of the column:

1-Calculating density in the dispersed and continuous phases

2-A 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 semi-industrial 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 liquid-liquid systems are presented in Tables 1 and 2.

Table 1

Geometric properties of the column.



Column component



Column height



Column diameter



Height of each stage



Number of stages



Number of holes per tray

40 degrees


Hole angle comparing to others



Plate thickness



Active volume of the column

Table 2

Physical properties of liquid-liquid systems.

n-Butyl acetate-water



Physical properties



Kg m-3

Continuous phase density



Kg m-3

Dispersed phase density




Continuous phase viscosity




Dispersed phase viscosity



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 n-butyl acetate-water and toluene-water systems. Considering that the operating condition is constant, mean droplet diameter size is bigger in toluene-water than the one in n-butyl acetate-water. The main reason for this occurrence is related to surface tension between the phases. Since the interfacial tension in toluene-water is almost twice as large as interfacial tension between n-butyl acetate-water, 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 toluene-water system (system with a high interface tension) is more than that in n-butyl acetate-water 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 n-butyl acetate-water and toluene-water 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 toluene-water system is more than that of n-butyl acetate-water 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.

n-Butyl acetate-water

N (rpm)


Qd (lit/hr.)

Experimental d32 (mm)

Simulation d32 (mm)





























































































Figure 6

Volume distribution of droplets versus the droplet diameter and the length of the column: a) n-butyl acetate-water-system and b) toluene-water 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 toluene-water 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 n-butyl acetate-water and toluene-water 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.


Column type

Deviation percentage average



Structure packing


Morales et al.


fixed-bed extractor


Lorena et al.


RDC extractor


Schmidt et al.


Kuhni extractor


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 n-butyl acetate-water-acetone system (with a low interface tension) is less than that in the toluene-water-acetone 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 non-spherical droplets is proposed.



: Interfacial area (m2/m3)


: Sauter mean drop diameter (m)


: Axial mixing coefficient (m2/s)


: Coalition rate


: Overall mass transfer coefficient (m/s)


: Rotor speed (1/s)


: Mass transfer flux (1/s)


: Breakage possibility


: Flow rate of the continuous or dispersed phase (m3/s)


: Time (s)


: Velocity (m/s)


: Slip velocity (m/s)


: Terminal velocity of droplet (m/s)


: Weber number


: Weber number critical


: Droplet diameter daughter

Greek letters


: Beta distribution function


: Density (kg/m3)


: Viscosity (Pa.s)


: Interfacial tension (N/m)



: Continuous phase


: Dispersed phase


: x-phase (continuous phase in present case)


: y-phase (dispersed phase in present case)



: n-butyl acetate-acetone-water


: Toluene-acetone-water


: Equilibrium value


: Inlet to column

Attarakiha, M. and Hlawitschkac, M. M. W., A Hyperbolic Population Balance Model for Dynamic Analysis of Liquid Extraction Columns, Kuala Lumpur: PSE ASIA, 25, 2013.
Attarakih, M. and Bart H. J., Solution of the Population Balance Equation Using the Differential Maximum Entropy Method (DMaxEntM): An Application to Liquid Extraction Columns, Chemical Engineering Science, Vol. 108, p. 123, 2014.
Artur, P. and Marcelo, B. M., Transient Modeling of Zinc Extraction With D2EHPA in a Kühni Column, Chemical Engineering Research Design, Vol. 91, p. 2323, 2013.
Casamatta, G. and Vogelpohl, A., Modeling of Fluid Dynamics and Mass Transfer in Extraction Column, German Chemical Engineering, Vol. 8, p. 96, 1985.
Cerutti, M. L., Ulson de Souza A. A., Ulson de Souza S. M., Solvent Extraction of Vegetable Oils: Numerical and Experimental Study,Food and Bioproducts Processing, Vol. 90, p. 199, 2012.
Godfrey, J. C. and Slater, M. J., Liquid-liquid Extraction Equipment, John Wiley and Sons, New York, 1994.
Gomes, L., N., Regueiras, P. F. R. L., and Guimaraes, M., Simulated and Experimental Dispersed-phase Breakage and Coalescence Behavior in a Kuhni Liquid-liquid Extraction Column-steady State, Industrial Engineering Chemical Research, Vol. 45, p. 3955, 2006.
Hasseine, A., Menai, A. H., Bencheikh, L. M., and Bart, H. J. Assessment of Drop Coalescence and Breakup for Stirred Extraction Columns, Chemical Engineering Technology, Vol. 28, p. 552, 2005.
Hufnagl, H., McIntyre, M., and Blab, E., Dynamic Behavior and Simulation of a Liquid–liquid Extraction Column, Chemical Engineering Technology, Vol. 14, p. 301, 1991.
Morello, V.S. and Poffenberger, N., Commercial Extraction Equipment, Industrial Engineering Chemistry, Vol. 42, p. 5, 1953.
Modes, G., Bart, H-J., Rodriguez‐Perancho, D., and Broder, D., Simulation of Solvent Extraction Column Fluid Dynamics from Single Droplet Parameters, Chemical Engineering Technology, Vol. 22, p. 231, 1999.
Mansur, M. B., Slater, M. J., and Junior, E.C.B., Reactive Extraction of Zinc Sulfate with Bis (2-ethylhexyl) Phosphoric Acid in a Short Kuhni Column Used in Batch Mode, Industrial Engineering Chemical Research, Vol. 42, p. 4068-4081, 2003.
Morales, C., Elman, H., and Pérez, A., Modeling and Simulation of a Liquid Extraction Column with Structured Packing, Computers and Chemical Engineering, Vol. 31, p. 1694, 2007.
Oliveira, N. S., Moraes Silva, D., and Gondim, M. P., A Study of the Drop Size Distributions and Hold-up in Short Kuhni Columns, Brazilian Journal of Chemical Engineering, Vol. 25, p. 729, 2008.
Kumar, A., Steiner, L., and Hartland, S., Capacity and Hydrodynamics of an Agitated Extraction Column, Industrial Engineering Chemical Research, Vol. 25, p. 728-733, 1986.
Rodea, S., Duranda, A., Mabilleb, I., and Favrea, E., Flooding Characteristics of an Aqueous two-phase System in a Counter-current Kühni-type Column, Chemical Engineering Science, Vol. 98, p. 98, 2013.
Steiner, L., Kumar, A., and Hartland, S., Determination and Correlation of Axial-mixing Parameters in an Agitated Liquid–liquid Extraction Column, The Canadian Journal of Chemical Engineering, Vol. 66, p. 241, 1988.
Schmidt, A. S., Populations Dynamische Modellierung Geruhrter Extraktions-Kolonnen Basierend auf Einzeltropfen, Aachen, Shaker Verlag, 2006.
Schmidt, S. A., Simon, M., Attarakih, M. M., Lagar, L. G., and Bart H. J., Droplet Population Balance Modeling-hydrodynamics and Mass Transfer, Chemical Engineering Science, Vol. 61, p. 246, 2006.
Wang, T., Wang, J., and Jin, Y., A Novel Theoretical Breakup Kernel Function for Bubble/Droplets in a Turbulent Flow, Chemical Engineering Science, Vol. 58, p. 4629, 2003.
Weinstein, O., Semiat, R., and Lewin, D.R., Modeling, Simulation and Control of Liquid-liquid Extraction Columns, Chemical Engineering Science, Vol. 53, p. 325-339, 1998.
Thakur, N.V., Separation of Rare Earths by Solvent Extraction, Mineral Processing Extraction Metallorgy Review, Vol. 21, p. 277-306, 2000.