VOF SIMULATION OF SINGLE RISING DROPS IN THREE LIQUID-LIQUID EXTRACTION SYSTEMS USING CSF AND CSS INTERFACIAL FORCE MODELS

* Corresponding author: capepub@cape.iust.ac.ir Abstract In liquid-liquid extraction contractors, mass transfer and stage efficiency are closely related to drop hydrodynamics. In the present study, hydrodynamic simulation of three standard liquid-liquid extraction systems recommended by the EFCE (European Federation of Chemical Engineering) has been investigated. Toluene/water, n-butyl acetate/water, and n-butanol/water with different drop diameters were considered in the simulations, representing systems with high, medium, and low interfacial tension respectively. In the current research, for the first time simulations have been carried out using the VOF-PLIC (Volume of Fluid Piecewise Linear Interface Calculation) model, implementing two surface tension force models of CSS (Continuum Surface Stress) and CSF (Continuum Surface Force) as a source term in the momentum equation. Simulations have been carried out in an axisymmetric geometry with a moving droplet in the static zone. The stages of droplet acceleration, deformation, and stability in terms of shape and velocity have been captured through simulations. Simulation results show that the average relative error reduces by using the CSS model and the most enhanced effect is observed in the toluene/water system, followed by the n-butyl acetate/water and n-butanol/water systems, respectively. This is due to higher parasitic current effects in the highest surface tension system (toluene/water). The onset of oscillations in the toluene/water system was correctly predicted by the CSS model, while the CSF model could not. Droplet shapes, aspect ratio, terminal and transient velocity and streamlines were also reported in the two surface tension models and compared.


INTRODUCTION
The velocity of free rising or falling droplets has great importance in liquid-liquid extraction systems.Especially, it will help to predict the residence time and drop size in single drop systems, which are the main parameters that influence the mass transfer coefficient and stage efficiencies.Although a swarm of droplets rather than single droplets exists in large industrial systems, deep understanding of this phenomenon will help in the design of large liquid-liquid extraction columns (Asadollahzadeh et al., 2011;Huang et al., 2016, andTorab-Mostaedi et al., 2011), with the new developed methods in order to enhance mass transfer (Goodarzi and Esfahany, 2016;Moghadam et al., 2017;Rahbar-Kelishami et al., 2015;Saien and Daliri, 2012).Due to the presence of contaminants and surface-active agents in the experimental studies, there are limitations on obtaining more precise results and, in some cases, requires very clean systems (Bäumler et al., 2011).Nowadays, numerical methods are widely used in different industrial applications (Ferreira et al., 2016;Oliveira et al., 2017, Pelissari et al., 2016), as well as in the new technologies of multiphase flow (Akbari et al., 2016;Goodarzi et al., 2014;Safaei et al., 2014Safaei et al., 2016).In the field of liquid-liquid extraction, numerical simulations are confidently used in order to attain a better hydrodynamic insight.This work focuses on rising single drops using the volume of fluid (VOF) method.
The most well-known methods to simulate moving boundaries are front tracking and front capturing methods.In the front tracking method, the location of the interface is represented by additional elements following the interface movement, while in the interface capturing method an additional indicator function is defined to implicitly find the location of the interface.The alternative approaches are the phase field or diffuse interface method and the Lattice-Boltzmann technique (Bertakis et al., 2010).Bertakis et al. (2010) did a hydrodynamic study, including both simulation and experimental assessments of n-butanol droplets with the finiteelement and extended finite element method.The level set function was used for capturing the interface movement in the non-stationary full 3D computational model.Bäumler et al. (2011) simulated toluene/ water, n-butyl acetate/water, and n-butanol/water systems with the mesh moving method as a front tracking approach in the moving reference frame and axisymmetric setting.The simulation results included terminal and transient velocities and aspect ratio variations with respect to time.Their study also included experimental evaluations of n-butyl acetate drops.Eiswirth et al. (2011) presented a numerical and experimental study of toluene drops using the level set method and they also visualized their results in a transformed reference frame.New fitting parameters were considered based on the Henschke correlation.Engberg and Kenig (2014) performed a numerical study for prediction of the hydrodynamic behavior of three standard test systems for liquid-liquid extraction with the level set front capturing method.Two surface tension models were used, including the continuum surface force (CSF) and ghost fluid model (GFM).Engberg and Kenig (2015) performed simulations reproducing the appearance of two different rise velocities for the same droplet diameter by the level set front capturing method.The different rise velocities had been previously shown by Wegener et al. (2010).Komrakova et al. (2013) studied the rising of n-butanol drops of 1-4 mm drop size range with both stationary and moving reference frame.Deviations of 5% for terminal velocities of smaller drops and up to 20% for large oscillating drops were reported using the Lattice-Boltzmann technique.
On the application of the VOF method in bubble hydrodynamics, one study was carried out by Bothe and Fleckenstein (2013).They used the VOF method for hydrodynamics and mass transfer simulation of oxygen transfer from a single bubble at a moderate Re number.Good agreements were reported with known Sherwood correlation.Several studies were also reported using the VOF method for the hydrodynamics and mass transfer of the rising bubbles, none including liquid-liquid extraction droplets (Alke et al., 2009;Bao et al., 2015;Fleckenstein and Bothe, 2013;Weiner and Bothe, 2017).
Some previous studies have considered the VOF method as a front capturing method for bubble hydrodynamic studies.To the best of the author's knowledge, there is no computational study using the VOF method for the liquid-liquid extraction system for a wide spectrum of droplet diameters.Additionally, the effect of two different surface tension models was evaluated in the simulations.Drop aspect ratio, transient and terminal velocity and oscillations were obtained from simulations.Continuum surface stress (CSS) and CSF models' results were compared with each other and with the terminal velocities of existing experimental data and simulation results of the previous authors with level set-CSF&GFM, Lattice-Boltzmann, and dynamic mesh models.The most important focus of the present study is to establish the precise VOF model accurately validated to be applied for different LLE systems hydrodynamics.The implementation of mass transfer effects in the presence of Marangoni effects is the other purpose of the present study in the current VOF model.The two mentioned aspects are the most important steps for the design and optimization of large-scale industrial LLE systems.

GOVERNING EQUATIONS
The VOF method proposed by Hirt and Nichols (1981), and level set method proposed by Osher and Sethian (1988) are amongst the most popular methods in front capturing techniques for moving boundaries.Recently published papers have used the level set method for the hydrodynamic study of liquid-liquid extraction systems.Hence, in the present study, the performance of the VOF method has been evaluated for the three-standard liquid-liquid extraction systems.

Volume Fraction Equation
The VOF method is used when the shape and flow of the process occurring near interface are of interest (Ranade, 2001).The important aspect of the model is the interface tracking between dispersed and continuous phases in order to capture accurate terminal velocity and surface.The VOF approach proposed by Hirt and Nichols (1981) meets this requirement with the definition of a variable called the volume fraction (a) with: α = 0 : When the cell is empty of drop phase α = 1 : When the cell is full of drop phase 0 < α < 1 : When the cell contains the interface of two phases In order to capture interface movement in the velocity field of U, the a equation has been solved for the drop phase for incompressible flow as (Ranade, 2001): Eq.1 was used in order to obtain of the a values of Eqs.3-7 in every computational cell.The surface tension of the water/dispersed phase makes the pressure jump across the interface with a gradient equaling the surface tension force appearing as a source term (F SF ) in the momentum equation (Eq.5).The surface tension force model has been utilized in the interface area, where the two phases exist.This term was neglected at the other sides of the computational domain.The CSF model proposed by Brackbill et al. (1992) and the CSS model presented by Lafaurie et al. (1994) have been used in the current research.

Interfacial Force Term of the Momentum Equations
In the present study, two surface tension model forces have been used and their effects on the hydrodynamics were investigated.
Continuum surface force model -The CSF interfacial force model replaces surface force by a smoothly varying volumetric force acting on all the fluid elements in the interface transition region.In multiphase flow, the interfacial tension force is taken into consideration through the below equation (Brackbill et al., 1992): The above equation can be discretized by standard methods of discretization.In the current study, the finite volume method was used by explicit formulation as below: where n+1 and n are the indices for new and previous time steps, respectively.a f n is the face value of dispersed phase volume fraction, V f n is volume flux through the face based on normal velocity, and V is the cell volume.The a equation was not solved for the continuous phase and the volume fraction of continuous phase was calculated from:

Fluid Dynamics
In the VOF approach, single momentum and continuity equations are solved for multiphase flow.So, the pressure and velocity component values are shared between the phases.The averaged mass and momentum conservation equations for the mixture of dispersed and continuous phases (toluene/water, n-butyl acetate/water, n-butanol / water) are indicated below (Pozzetti and Peters, 2017;Ranade, 2001) The averaged density and viscosity are defined as (Pozzetti and Peters, 2017;Ranade, 2001): (1 ) r = ar + − a r d c (1 ) In the present study, only two phases exist; thus κ j = -κ i and ∇αj = -∇αi, therefore, Eq.8 will be: In this equation, σ is interfacial tension, κ is local surface curvature, n is surface unit normal, and r is volume fraction averaged density, which was computed from Eq.6, and local surface curvature κ is defined as: and the surface normal vector n is defined as the gradient of phase volume fraction as below: ( ) Hence, κ will be: velocity component.The Green-Gauss Cell-Based gradient evaluation method has been applied in order to solve the gradient of the convective and diffusive terms in the flow conservation equations.In order to compute the gradient of the scalar f at the cell center C 0 , the following discrete form has been written: Constant surface stress model -An alternative way of modeling surface tension is the CSS model.Unlike the formulation of the CSF, CSS avoids the explicit calculation of curvature and could be considered as an anisotropic variant of modeling capillary forces based on the surface stresses.The CSS model converts the surface tension force to the stress form of T as below (Lafaurie et al., 1994): ) This is tangential to the surface.Then the surface tension force by the CSS model is written by: . .
The surface normal vector, n, normal unit vector n, curvature κ, and, a have been calculated at the cell centers.

SOLUTION PROCEDURE
Eqs.1, 4, and 5 were solved using the finite volume method.The momentum equation was solved using the Quick discretization method of Leonard and Mokhtari (1990).Discretization of Eq. 1 through the control volume faces with the finite volume method requires fluxes of a through all cell faces.Computational errors in the evaluation of convective fluxes lead to inaccuracy in the interface representation.In order to reduce the errors, the two popular methods of Simple Line Interface Calculation (SLIC), proposed by Hirt and Nichols (1981), and Piecewise Linear Interface Calculation (PLIC), proposed by Youngs (1982) were used.In the current research, the PLIC method was used due to its accuracy, strict mass conservation, and significant improvement of interface representation (Gerlach et al., 2006).In the PLIC method, the interface line within the computational cell is approximated with a straight-line segment with its slope obtained from the interface normal.The computational cell is divided by the line segment where the drop phase volume fraction is a.The resulting polygon shape of the drop phase is then used for calculation of drop fluxes along four faces of the computational cell with the available ( ) where, f f is the value of f at the cell face centroid, and V and A f , are computational face surface and cell volume respectively.The summation is over all faces enclosing the cell.f f is the arithmetic average of the values at the neighboring cell centers and is defined as: where, f c0 and f c1 are the value of f in the neighbor cells.
Unsteady terms in Eqs.1, 4 and 5 were discretized using a first order implicit method.The pressureimplicit with the splitting of operators algorithm (PISO), developed by Issa (1986), was used in order to relate the pressure and velocity coupling.The Pressure-Implicit with Splitting of Operators (PISO) coupling scheme, as a part of the SIMPLE family of algorithms, is recommended for transient calculations in the VOF method with large time steps.It improves the efficiency of SIMPLE and SIMPLEC algorithms through satisfying the momentum balance equation after the solution of the pressure correction equation.The PRESTO! interpolation scheme is applied to compute the face values of pressure from the cell values in the VOF multiphase flow simulations (Patankar, 1980).In the PISO algorithm, when the momentum equation was solved, pressure corrections were computed by neighbor correction, skewness correction, and skewness -neighbor coupling processes.Based on the recently updated values of pressure, face mass fluxes and velocities were updated.The variable time stepping method was used through simulation with a minimum time step of 0.01 millisecond.The maximum Courant number was considered to be 0.15 and assumed constant through all of the simulations.Calculation reached convergence when the residuals of all equations fell below 10 -5 .The codes have been developed in Visual Studio 2010.

PROBLEM DETAILS
The binary material properties of standard liquidliquid extraction systems are given in Table 1.N-butanol/water and toluene/water binary systems have been used through simulations as the lowest and the highest surface tension systems, respectively.The schematic diagram of the computational domain with boundary condition and mesh structure is shown in Fig. 1.The axis symmetric boundary condition was considered in this research, effectively reducing the computational effort (Bäumler et al., 2011;Engberg and Kenig, 2014).No slip boundary condition was assumed at the bottom & left, while the pressure boundary condition was assumed at the upper side of the computational domain (Fig. 1A).The displayed pressure on the figure is the modified pressure based on Rusche (2003), to simplify the specification of boundary conditions and is defined as: (Engberg and Kenig, 2014;Engberg et al., 2014b andEngberg et al., 2014a) and different domain heights with respect to the drop type and size.Due to the important effect of the domain height on the transient and terminal velocity results, different domain height sizes were considered for every type of the droplets of a specified diameter.For each droplet having reached a certain height, the results did not show considerable new variations with increasing domain heights.This was regarded as the final height of the droplet.Toluene droplets have the highest computational domain in comparison to the other two systems of the same droplet size.The domain heights were changed from 10 to 55cm from small to large drops.

Grid Convergence
Center of mass velocity of the droplet has been calculated from Eq.19 for the dispersed phase as below (Jeon et al., 2009): * .= − r P P g X More details have been presented in Rusche (2003).The initial position of the droplet was considered to be slightly higher from the bottom wall of the domain (Fig. 1B).The droplet was initiated as a sphere and quiescent at the beginning.It can be seen from Fig. 1B that the mesh type was structured and the most refined within the vicinity of the droplet.In terms of initial drop diameter D 0 , the computational domain has 4D 0 width  1 , where V cm,y is the y component of the center of mass velocity, a i is the droplet volume fraction of the i th cell, v i is the volume of the i th cell and v i is the y direction velocity component of the drop phase at the center of the computational cell.Droplet center of mass velocity was calculated from Eq. 19 and results are provided in Table 2 using the CSF and CSS models for an initial drop diameter of 2.5 mm for the three systems types.
It can be seen that the toluene droplet results did not change considerably after D0/80, while this value is D0/60 for n-butanol and n-butyl acetate droplets where D 0 is the initial drop diameter.Increasing the interfacial tension and an increase in the aspect ratio required a more refined grid size to evaluate the droplet curvature and consequently to reach grid independence of the results. (

Drop Shapes, Axial Velocity Contours, and Streamlines
Axial velocity contours are presented in Fig. 2 for two liquid-liquid extraction systems.It can be seen that the drop shapes and axial velocity profiles changed with the variations in drop diameters and types.The maximum value of the axial velocity exists at the center of the droplet and the velocity values decrease close to the interface due to the shear stress of the continuous phase.
For the three systems, axial velocities and pressure profiles are demonstrated in Figs.3A & B. It can be seen that the axial velocity values are maximum at the center of the droplet and gradually decrease close to the interface, as was observed in Fig. 2. As expected, toluene has the largest axial velocity (Fig. 3A).Because of the larger surface tension of toluene drops, the pressure jump is larger for the toluene/ water system (Fig. 3B).Any symptom of instabilities close to the drop interface is not observable; therefore, the results of this section showed that no numerical instabilities exist close to the drop interface.
CSS and CSF simulation result comparisons in terms of axial velocities are presented in Figs.3C&D for the toluene/water system.It can be seen from Fig. 3C that, for drop diameters ranging from 1 to 2 mm, there is a considerable difference between the CSS and CSF simulation results.As can be seen from Fig. 3C, after 2 mm drops, the differences between CSS and CSF simulation results are reduced.It can be seen from Fig. 3D that, for droplets larger than 3mm, the axial velocity evaluated by the CSF method is slightly larger compared to the CSS method.forces with respect to the surface tension force acting as a resisting force.Unlike the toluene and n-butyl acetate drops (Figs. 5 and 6), n-butanol drops (Fig. 4) showed strong shape deformations in the diameters of 3.48mm and 4mm due to their lower surface tension force.The 5mm n-butanol drop shapes prior to breakup were simulated with the CSS model (Fig. 7).It is observable that the shape changes are from spherical to crescent and finally to toroidal rim prior to breakup.This shape is no longer stable, since N f and Eö numbers are relatively higher.Consequently, the buoyant force is dominant over the viscous and surface tension forces.This force pushes the drop from the rear during the transition process and finally breakup occurs.Similar drop shapes were previously reported by several authors (Engberg and Kenig, 2014;Bertakis et al., 2010;Liu et al., 2013) for the most deforming and breaking drops.
During simulations, the droplet initial mass was compared with the final value for all types and sizes of droplets.No mass losses were observed and mass conservation was found.The aspect ratio values (Tables 3-5) agreed well with previously reported data and confirmed the drop shapes obtained from simulations.circulating, and oscillating regimes.As the drop diameter increases, there is a tendency for shape change from spherical to oblate, as is more evident at the very large diameter of every test system.This is due to the increase in Eö and We numbers and the increase of buoyant and inertial forces as deforming as the ratio of the drop vertical to horizontal diameter and calculated from Eq. 20 for the axis-symmetric conditions (Engberg and Kenig, 2014).Streamlines for different drop diameters and types are shown in Fig. 8.As can be seen from Fig. 8A, the vortex appears in the wake of n-butanol droplets from the drop diameter of 2.26mm.The recirculation zone inside of the wake enlarges as the drop diameter increases (Fig. 8E).With the increase of drop diameter and longitudinal extension of drop shape, the vortex zone size increases.Toluene and n-butyl acetate oblate drops do not have vortexes inside of the wakes for drop diameters similar to n-butanol drops.This was due to their larger aspect ratio and smaller longitudinal extension.Wakes start to form for 6 mm n-butyl acetate drops from 0.4s and become larger with time (Figs.8G & H).

Aspect Ratio
The Aspect ratio is the exact comparison parameter for drop shapes obtained from simulations.It is defined Here x min , x max , y min and y max were computed from the simulations as the maximum and minimum values of the droplet elemental positions in the x and y directions, where x min is zero in the current axisymmetric conditions.The aspect ratio values of the CSF and CSS models, as well as the comparison with the level set-CSF and GFM models of Engberg and Kenig (2014), are provided in Tables 3-5.Close agreement between VOF results of the present work and the level set simulation results of Engberg and Kenig ( 2014) is observable.It can be seen that, with the increase of drop diameter, the aspect ratio values reduce.This was due to the increase of the buoyant and inertial forces as the deforming forces with respect to the surface tension force acting as the resisting force.Therefore, drop shapes tend to an oblate shape with lower aspect ratio values.
One of the important parameters to distinguish the regime change of rising droplets is the aspect ratio.The spherical regime ends when the aspect ratio of the droplets falls below 0.95 and the oscillation regime starts when the oscillations in the aspect ratio are more than 2 percent (Bäumler et al., 2011).The results in Table 3 show that only 1mm n-butanol drops were within the spherical regime.The drop diameter for the oscillating regime startup was obtained for the droplets of 3.48 and 4 mm diameter by the VOF (20) method for n-butanol droplets.Regarding n-butyl acetate droplets (Table 4), it can be seen from Table 4 that, for drop diameters between 1.5 and 2 mm, the aspect ratio falls below 0.95, so the spherical regime endpoint in the VOF method was located between 1.5 and 2 mm drops.The oscillations were also observed at drop diameters about 4mm and above.Simulation results reported in Table 5 show that the aspect ratio of toluene droplets falls below 0.95 for the droplet diameters between 2 and 2.5 mm.Therefore, the end of the spherical regime was located between 2 and 2.5 mm diameter.The oscillations were observed initially for 4.4 mm diameter toluene drops using the CSS model.Therefore, the start of the oscillating regime for toluene droplets is 4.4 mm, which has been verified by experimental investigations of Wegener et al. (2010).Good agreement has been achieved between experimental results and the current simulation predictions.
For the toluene/water system, simulation results showed that the aspect ratio values predicted by the VOF-CSS model were larger in the spherical regime with respect to the VOF-CSF.The prediction of correct curvature by CSS model, especially for small toluene drops, causes accurate prediction of transient and terminal velocities.

Droplets Transient and Terminal Velocities
Transient and terminal velocities of three different liquid-liquid extraction systems for a wide range of droplets and relative deviations with respect to existing experimental data are discussed in this section.
N-butanol drops -N-butanol drop transient velocities obtained from the VOF-CSS model are shown in Fig. 9. Droplets within a diameter range of 1 to 3mm demonstrate an increasing trend of transient velocities.Above 3mm, oscillations and a decreasing trend in velocity values were observed for 3.48 mm and 4mm diameters (Fig. 9A).In the spherical regime, with the increase of drop diameter, the velocity increases due to the increase of the buoyant force up to its maximum value in the transition regime.An additional increase in drop diameter causes a reduction in terminal velocity due to an increase of interfacial area and consequently an increase of drag coefficient.Transient velocities of n-butanol oscillating drops are shown in Fig. 9B.It is observable that a n-butanol 4 mm droplet has almost the same transient profile with both CSF and CSS models (Fig. 9B).
The terminal velocity of n-butanol drops with both the CSF and CSS models is shown in Fig. 10A.

N-butyl acetate droplets -
The transient velocities of n-butyl acetate droplets are demonstrated in Figs.11A & B. An increasing trend of transient/terminal velocity is observable until the drop diameter of 3mm, above which the velocity decreases.Due to the need for long rise times to understand oscillation behaviour, oscillating droplets are given in a separate plot.N-butyl acetate large droplets of 4 and 6 mm were simulated using the CSS model (Fig. 11B).As expected, with the increase of drop diameter, the terminal velocity decreases due to the larger interfacial area and larger drag coefficient, while the amplitude of oscillations increases (Fig. 11B).It can be seen from Fig. 11B that, changing the surface tension model for a n-butyl acetate 4 mm droplet did not change the terminal velocity significantly, but caused differences in the amplitude and frequency predictions.
The terminal velocity of n-butyl acetate droplets is provided in Fig. 12A using the CSF and CSS models.It can be seen that the 3 mm droplets have maximum velocities in the CSS and CSF models.Relatively small deviations are observable for VOF-CSS results with respect to the VOF-CSF model due to the reduction of parasitic current effects of CSS with respect to CSF.This concept will be discussed below in more detail.Excellent agreement is also observable between the current simulation results and existing experimental data.Relative errors of current simulations and the Engberg and Kenig (2014) simulation results are given  Average relative errors are reported in Table 7 with respect to experimental data of Bäumler et al. (2011).In this table, U char is the characteristic mean velocity and is computed by dividing the travel distance by the time needed to pass.The travel distance starts after the acceleration stage is finished and ends when the upper end of the column is reached (Baumler et al. 2011).It can be observed that the average relative error was reduced from 8.68 to 5.96 in terms of experimental U max velocity when the interfacial force is changed from CSF to CSS.The average relative error of the dynamic mesh method of Bäumler et al. (2011) is smallest within the four models compared in terms of U max .
Toluene droplets -The transient velocity of toluene droplets is presented in Figs.13A & B. As stated in Figs.9A and 11A, transient/terminal velocity increases to the maximum value when the drop diameter increases.After that, a small increase in the drop diameter causes a reduction in terminal /transient velocity.It can be seen from Fig. 13A, toluene 3.5 mm drops have the maximum velocity among the simulated drops.Oscillating drops are presented in Fig. 13B.It can be seen that the oscillating regime started from 4.4 mm drops for the CSS model while the CSF model could not predict oscillations for 4.4mm drops (Fig. 13B).The experimental investigation of Wegener et al. (2010) showed that the startup of the oscillating regime for toluene droplets was 4.4 mm.The CSS model successfully predicted the oscillation regime for toluene drops.VOF-CSS simulations continued for the large droplets of 6 and 7mm (Fig. 13B).The simulation results showed that the increase in drop diameter from 4.4 to 7 mm reduced the terminal velocity due to the large surface area and large drag coefficient, but the amplitude and frequencies increase.The simulated transient velocities of toluene drops were compared with experimental data of (Wegener et al., 2010) and good agreement was observed (Fig. 13C).It can be seen from Fig. 13B that the amplitude of the oscillations was 5 mm/s, while the periods were 0.07 s for 4.4 mm drops.With the increase of drop diameter to 5 mm, the amplitude and periods were 42 mm/s and 0.09 s respectively.In order to verify the simulated periods and amplitudes, experiments should have been carried out at different frame rates in which amplitudes and frequencies were constant, while the frame rate increased up to a sufficient level.Then, it would be confirmable whether periods and amplitudes were real or arose from numerical instabilities.For 5 mm toluene drops, Engberg and Kenig (2015) reported that discrepancies in amplitudes and frequencies were due to the axisymmetric assumption and initial drop shape.But for 4.4 mm drops, an 0.07 s period is uncertain under the present conditions as there are insufficient data on oscillations with different frame rates.Similar discrepancies were also reported between simulated amplitudes and frequencies and experimental data in the previously published literature (Engberg et al., 2014b;Bäumler et al., 2011).

Method
The terminal velocity of toluene droplets is presented in Fig. 14A.For droplet diameters ranging from 1 to 3.5mm, the droplet velocity increased as droplet diameter raised, after which a further increase in droplet diameter caused a reduction in terminal velocity.It was also observed that large toluene droplets in the diameter range of 5 to 7 mm were successfully simulated with the current VOF-CSS model for the first time.Referring to Figs. 14B and C, the CSS model considerably reduced the VOF method relative errors.The results were in agreement with the simulation results shown in Fig. 3C&D.Referring to Fig. 3C&D, it was concluded that changing the surface tension force model from CSF to CSS reduces the parasitic current effects and gives a larger axial velocity.Therefore, more precise results are obtainable compared to experimental data.Similar observations were also possible in Fig. 14A.Referring to Fig. 14A, the CSS method also provides better results for the 3.5 and 4 mm drops.The average relative error reduced from 9.88 to 3.63 percent with respect to U max upon changing surface tension force from CSF to CSS.Comparison of average relative errors for the four models are presented in Table 8.It can be seen that the level set-GFM method gives the best results among the different models.
The results of the present study showed that the better performance of the level set method with respect to VOF is most evident in the toluene/water system.A study by Abadie et al. (2015) showed that spurious currents of the level set-CSF method were less than in the VOF-CSF method.Hence, application of the level set method reduced the parasitic currents and gave more precise terminal velocities.The CSS model average maximum parasitic velocity was about 0.056m/s, while for the CSF model this value was about 0.064 m/s for 2 mm toluene drops in different grid sizes.The average relative error for pressure jump was about 6.1 percent for the CSS model, while for the CSF model this value was about of 7.6 percent.Hence, it can be seen that the CSS model reduces parasitic current effects and consequently reduces the average relative errors with respect to the experimental data.Similar findings of the parasitic current reduction when the model was changed from VOF-CSF to VOF-CSS and better performance of the VOF-CSS method have also been reported previously (Ahmadi Nadooshan and Shirani, 2008;Albert et al., 2012;Seifollahi et al., 2008).The results of the present study agree well with the previous findings of better performance of the CSS model with respect to the CSF.
Among the three systems, the surface tension value differences are larger (ranging from 1.63 to 35 mN/m) with respect to density and viscosity value changes.Consequently, the increase of surface tension coefficient increases the surface tension force compared to other effective forces such as inertia and viscous forces and reduces the Capillary and We numbers.As stated by Harvie et al. (2006), large parasitic currents are observable at higher Re and lower We and Ca numbers.Referring to the dimensionless values in the Appendix, it was predictable from the above discussion that the toluene/ water system, with a larger Re number and lower Capillary and We numbers, has the larger parasitic current value, followed by the n-butyl acetate/water and n-butanol/water systems, respectively.This is the reason why the maximum deviation reduction from 9.88 to 3.63 occurred in the toluene system,when the surface tension model was changed from CSF to CSS, followed by the n-butyl acetate and n-butanol systems, respectively.Hence, the surface tension model change shows its enhanced effect most evidently in the system with larger parasitic current effects (toluene/water).

CONCLUSIONS
VOF simulation of liquid-liquid extraction droplets with two interfacial forces of the CSF and CSS models were considered in the current research for the first time.Droplet shapes, aspect ratios, axial velocities, transient and terminal velocities, and streamlines were reported for the threedifferent liquid-liquid extraction systems of toluene/ water, n-butyl acetate/water, and n-butanol/water.Additionally, large droplets, including 5, 6, and 7 mm toluene drops, 6 mm n-butyl acetate drops, and 5 mm n-butanol drops, were simulated in the breakup regime for the first time.The following results were achieved for the three systems: 1. Good agreement was achieved between VOF results and previous experimental data (Bertakis et al., 2010;Wegener et al., 2010;Bäumler et al., 2011) for the three test systems.The results included: terminal velocity, drops with maximum velocity, regime change from spherical to circulating and from circulating to oscillating, and shape variations with time.Amplitude and periods of oscillations had differences in comparison with experimental data due to the axisymmetric assumption and initial drop shape.
2. The simulated transient velocities of 3.5 and 4 mm toluene drops were compared with experimental data and good agreement was achieved.
3. Changing the surface tension force from CSF to the CSS model reduced the average relative error with respect to the existing experimental data.This reduction was emphasized in the toluene/water system, with an average relative error reduction from 9.88 to 3.63 percent with respect to the U max experimental data of Wegener et al. (2010).The average relative error of the n-butyl acetate/water system reduced from 8.68 to 5.96 percent with respect to the U max experimental data of Bäumler et al. (2011).The average relative error of the CSS model was 2.48 percent for the n-butanol/ water system, while this value was 3.01 percent for the CSF model with respect to the existing experimental data of Bertakis et al. (2010).
4. For the n-butanol/water system, the VOF-CSS method provided the best performance, with a minimum average relative error.The main part of the deviation reduction was observed in the drop size range from 3.06 to 3.48 mm.For the n-butyl acetate/ water and toluene/water systems, the dynamic mesh method of Bäumler et al. (2011) and level set-GFM method of Engberg and Kenig (2014) provided the minimum average relative error, respectively, when compared to the U max experimental data.
5. The VOF-CSS method could predict the oscillating regime for 4.4 mm toluene drops, which was in agreement with experimental findings.The startup of the oscillating regime was 5mm in the VOF-CSF model.
6. Simulation results showed that, above 2.26 mm n-butanol droplets, a vortex was observed in the wake of n-butanol droplets and enlarged as the drop diameter became larger.
7. 5mm n-butanol drops were in the breakup regime.But due to the need for very fine mesh refinement in the interfacial area, further investigations are needed.

Figure 1 .
Figure 1.Schematic diagram of the computational domain (A) boundary conditions and (B) mesh structure.
Drop shapes obtained from simulations are shown in Figs.4-6 for three liquid-liquid extraction systems.The shapes in the figures are in the spherical,

Figure 2 .
Figure 2. Axial velocity contours for different drop diameters of liquid-liquid extraction systems.

Figure 3 .
Figure 3. A) Axial velocity, B) Static pressure plots in three liquid-liquid extraction systems, surface tension force model effect on the axial velocity profiles for different drop sizes C: A) 1 to 2mm, and D) 2.5 to 4 mm.

Figure 7 .
Figure 7. Transition process of the drop shape changes for a 5mm n-butanol drop from spherical to breakup at times of: 0, 0.1, 0.2, 0.3, 0.4, 0.5 and 0.54s using the CSS model.Hence, the VOF-CSF & CSS methods combined with PLIC were successfully used in the simulations.Regarding the above-mentioned points, no concern remained about mass losses during the simulations.Streamlines for different drop diameters and types are shown in Fig.8.As can be seen from Fig.8A, the vortex appears in the wake of n-butanol droplets from the drop diameter of 2.26mm.The recirculation zone inside of the wake enlarges as the drop diameter increases (Fig.8E).With the increase of drop diameter and longitudinal extension of drop shape, the vortex zone size increases.Toluene and n-butyl acetate oblate drops do not have vortexes inside of the wakes for drop diameters similar to n-butanol drops.This was due to their larger aspect ratio and smaller longitudinal extension.Wakes start to form for 6 mm n-butyl acetate drops from 0.4s and become larger with time (Figs.8G & H).

Figure 8 .
Figure 8. Streamline of the different drop diameters in three liquid-liquid extraction systems using the CSS model.

Figure 9 .
Figure 9. A) n-butanol drop transient velocities using the CSS model, B) transient velocity of oscillating n-butanol drops with the CSF and CSS models.
seen that both the CSS and CSF models show nearly the same results and good agreement was achieved.The deviations from Bertakis et al. (2010) experimental data are demonstrated for the two models of VOF-CSS&CSF in Fig.10 B. The Level set-CSF & GFM deviations are also presented in Fig. 10 B from the Engberg et al. (2014) simulation results.The CSS model gives less deviation compared to the CSF in almost all diameter ranges except 1.56 mm.It can be seen that at the beginning of the circulating and oscillating regimes, the deviations of the VOF method are less than the level set.It has been previously reported that the simulation rise time of the study of Bertakis et al. (2010), was too short to reach terminal velocity (Engberg and Kenig, 2014; Komrakova et al., 2013).They concluded that their large deviations from the Bertakis et al. (2010) terminal velocity data were due to not considering a sufficiently long rise time in Bertakis et al. (2010) for the oscillating large drops.Using long rise times for n-butanol oscillating droplets (Fig. 9B) resulted in a good agreement with the Bertakis et al. (2010) simulation and experimental data in comparison with previously published data (Figs.10A & B).Average relative errors are provided in Table 6 for the current CSS and CSF models and existing simulation results of the literature.It can be seen that the VOF-CSS model gives the best result among the other methods.

Figure 10 .
Figure 10.A) n-butanol drop terminal velocities obtained from the current CSS and CSF models in comparison with the experimental data of Bertakis et al. (2010), CSF and GFM models of Engberg and Kenig (2014), Lattice-Boltzmann model of Komrakova et al. (2013), and dynamic mesh model of Bäumler et al. (2011) B) Relative errors of the CSF, CSS models of the present work, and CSF and GFM results of Engberg and Kenig (2014) with respect to the experimental data of Bertakis et al. (2010).

Figure 11 .
Figure 11.A) n-butyl acetate drops transient velocities using the CSS model, B) Transient velocity of oscillating n-butyl acetate drops with the CSF and CSS models.
A)n-butyl acetate drop terminal velocities obtained from the current CSS and CSF models in comparison with experimental data of Bäumler et al. (2011) , CSF and GFM models of Engberg and Kenig (2014), and dynamic mesh model of Bäumler et al. (2011) B) Relative errors of CSF, CSS models of the present work, and CSF, GFM models of Engberg and Kenig (2014) with respect to U char experimental data of Bäumler et al. (2011) and C) relative errors of CSF, CSS models of the present work, and CSF , GFM models of Engberg and Kenig (2014) with respect to U char experimental data of Bäumler et al. (2011). in Figs.12B & C. It is observed from Fig. 12C that the relative error of the CSS model is less compared to CSF model for all diameters except 1mm drops.

Figure 13 .
Figure 13.A) Toluene drop transient velocities using the CSS model, B) Transient velocity of oscillating toluene drops with the CSF and CSS models, C) Transient velocity of the current CSS model in comparison with Wegener et al. (2010) experimental data.

Figure 14 .
Figure 14.A) Toluene drop terminal velocities obtained from the current CSS and CSF models in comparison with experimental data of Wegener et al. (2010), CSF and GFM models of Engberg and Kenig (2014), and dynamic mesh model of Bäumler et al. (2011) B) Relative errors of CSF, CSS models of the present work, and CSF, GFM models of Engberg and Kenig (2014) with respect to U char experimental data of Wegener et al. (2010), and C) Relative errors of CSF, CSS models of the present work, and CSF , GFM models of Engberg and Kenig ( 2014) with respect to U max experimental data of Wegener et al. (2010).

Table 1 .
Binary material properties for standard liquid-liquid extraction systems: n-butyl acetate/water and toluene /water systems at 25 o C and n-butanol/water system at 20 o C.

Table 2 .
Mesh independence analysis for 2.5 mm drop size for terminal velocity (mm/s) for three binary liquidliquid extraction systems using the CSF and CSS models.

Table 3 .
Simulated aspect ratio values of n-butanol drops using the VOF method (AR VOF-CSF and AR VOF-CSS ) under steady state conditions and comparison with the level set simulation results ofEngberg and Kenig (2014).

Table 5 .
Simulated aspect ratio values of toluene drops using the VOF method (AR VOF-CSF and AR VOF-CSS ) under steady state conditions and comparison with the level set simulation results ofEngberg and Kenig (2014).

Table 4 .
Simulated aspect ratio values of n-butyl acetate drops using the VOF method (AR VOF-CSF and AR VOF-CSS ) under steady state conditions and comparison with the level set simulation results ofEngberg and Kenig (2014).

Table 7 .
The average relative error of the present VOF model (VOF-CSS and VOF-CSF), level set model(level set-CSF and level set -GFM) of Engberg and Kenig (2014) and dynamic mesh model of Bäumler et al., 2011) with respect to the experimental data of Bäumler et al. (2011) for n-butyl acetate drops.

Table 8 .
Relative error average of VOF model (VOF-CSS and VOF-CSF) , level set model(level set-CSF and level set -GFM) of Engberg and Kenig (2014), and dynamic mesh model of Bäumler et al.(2011) with respect to the experimental data of Wegener et al. (2010) for toluene drops.