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

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.

Key words:
Single drop; VOF; Liquid-liquid extraction; Hydrodynamics; CSF; CSS

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, and Torab-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., 2014, Safaei 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 finite-element 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):

$$\frac{\partial \alpha }{\partial t}+\left(U.\nabla \alpha \right)=0$$ (1)

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:

$$\frac{{\alpha }^{n+1}-{\alpha }^n}{∆t}V+\sum _{ƒ}{V}_{ƒ}^n{\alpha }_{ƒ}^n=0$$ (2)

where n+1 and n are the indices for new and previous time steps, respectively. α_fⁿ is the face value of dispersed phase volume fraction, V_fⁿ is volume flux through the face based on normal velocity, and V is the cell volume. The α equation was not solved for the continuous phase and the volume fraction of continuous phase was calculated from:

$${\alpha }_d+{\alpha }_c=1$$ (3)

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

$$\frac{\partial \rho }{\partial t}+\nabla .\rho U=0$$ (4)

$$\frac{\partial \left(\rho U\right)}{\partial t}+\nabla .\rho UU=-\nabla P+\nabla .\mu \left(\nabla U+\nabla U^T\right)+\rho g+F_{SF}$$ (5)

The averaged density and viscosity are defined as (Pozzetti and Peters, 2017; Ranade, 2001):

$$\rho =\alpha {\rho }_d+\left(1-\alpha \right){\rho }_c$$ (6)

$$\mu =\alpha {\mu }_d+\left(1-\alpha \right){\mu }_c$$ (7)

Eq.1 was used in order to obtain of the α 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):

$$F_{SF}=\sum _{pairsij,i<j}{\sigma }_{ij}\frac{{\alpha }_i{\rho }_i{\kappa }_j\nabla {\alpha }_j+{\alpha }_j{\rho }_j{\kappa }_j\nabla {\alpha }_i}{{\displaystyle \frac{1}{2}}\left({\rho }_i+{\rho }_j\right)}$$ (8)

In the present study, only two phases exist; thus κ_j = - κ_i and ∇αj = - ∇αi, therefore, Eq.8 will be:

$$F_{SF}=\sigma \frac{\rho \kappa \nabla \alpha }{\left({\rho }_i+{\rho }_j\right)/2}$$ (9)

In this equation, σ is interfacial tension, κ is local surface curvature, n is surface unit normal, and ρ is volume fraction averaged density, which was computed from Eq.6, and local surface curvature κ is defined as:

$$\kappa \left(X,t\right)=\nabla .\alpha$$ (10)

in which;

$$\hat{n}=\frac{n}{\left|n\right|}$$ (11)

and the surface normal vector n is defined as the gradient of phase volume fraction as below:

$$n\left(X,t\right)=\nabla \alpha$$ (12)

Hence, κ will be:

$$\kappa \left(X,t\right)=\frac{1}{\left|n\right|}\left[\left(\frac{n}{\left|n\right|}·\nabla \right)\left|n\right|-\left(\nabla ·n\right)\right]$$ (13)

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

$$T=\sigma \left|n\right|\left(I-\hat{n}\otimes \hat{n}\right)=\sigma \left(\left|\nabla \alpha \right|I-\frac{\nabla \alpha \otimes \nabla \alpha }{\left|\nabla \alpha \right|}\right)$$ (14)

This is tangential to the surface. Then the surface tension force by the CSS model is written by:

$$F_{SF}=\nabla ·T=\nabla ·\left[\sigma \left(\left|\nabla \alpha \right|I-\frac{\nabla \alpha \otimes \nabla \alpha }{\left|\nabla \alpha \right|}\right)\right]$$ (15)

The surface normal vector, n, normal unit vector n, curvature κ, and, α 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 α 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 α. 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 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 ϕ at the cell center C₀, the following discrete form has been written:

$${\left(\nabla \varphi \right)}_{c_0}=\frac{1}{V}\sum _{ƒ}{\overline{\varphi }}_{ƒ}{A}_{ƒ}$$ (16)

where, ϕ_f is the value of ϕ 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 is the arithmetic average of the values at the neighboring cell centers and is defined as:

$${\overline{\varphi }}_{ƒ}=\frac{{\varphi }_{c0}+{\varphi }_{c1}}{2}$$ (17)

where, ϕ_c0 and ϕ_c1 are the value of ϕ in the neighbor cells.

Unsteady terms in Eqs.1, 4 and 5 were discretized using a first order implicit method. The pressure-implicit 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 liquid-liquid 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:

$$P^{*}=P-\rho g·X$$ (18)

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

	Density (kg/m3)	Viscosity (mPa.s)	Surface Tension (mN/m)	Data Source for Validation	Viscosity Ratio	Density Ratio
water (continuous phase)	986.51	1.39	1.63	Bertakis et al. (2010)	2.36	0.86
n-butanol (dispersed phase)	845.44	3.28
water (continuous phase)	996.6	0.904	14	Bäumler et al. (2011)	0.91	0.88
n-butyl acetate (dispersed phase)	877.19	0.82
water (continuous phase)	997.02	0.8903	35	Wegener et al. (2010)	0.62	0.86
toluene (dispersed phase)	862.3	0.552

Figure 1
Schematic diagram of the computational domain (A) boundary conditions and (B) mesh structure.

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₀, the computational domain has 4D₀ width (Engberg and Kenig, 2014; Engberg et al., 2014b and Engberg 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.

RESULTS AND DISCUSSION

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

$$V_{cm,y}=\frac{\sum _{i=1}^n{\alpha }_i{V}_i{v}_i}{\sum _{i=1}^n{\alpha }_i{V}_i}$$ (19)

where V_cm,y is the y component of the center of mass velocity, α_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₀ 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.

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

	d/20	d/40	d/60	d/80	d/100	d/120
toluene/Water -CSS	131.02	142.08	144.73	146.37	146.32	146.35
toluene/Water -CSF	130.73	141.56	142.018	142.76	142.56	142.59
n-butyl acetate/Water-CSS	106.2	111.43	114.01	114.27	114.055	114.28
n-butyl acetate/Water-CSF	105.53	110.51	111.92	112.0	111.93	112.334
n-butanol/water-CSS	58.79	60.41	60.73	60.85	60.89	60.94
n-butanol/water-CSF	59.1	60.42	60.65	60.78	60.81	60.85

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.

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

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.

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.

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.

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

Figure 4
Droplet shapes obtained from CSS simulations for n-butanol drops with the initial diameters of 2, 2.48, 3, 3.48 and 4 mm.

Figure 5
Droplet shapes obtained from CSS simulations for n-butyl acetate drops with the initial diameters of 1.5, 2.0, 2.5, 3.0, 3.5 and 4.0 mm.

Figure 6
Droplet shapes obtained from CSS simulations for toluene drops with the initial diameters of 1.5, 2.0, 2.5, 3.0, 3.5, 4.0 and 4.4 mm.

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.

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.

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

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 of Engberg and Kenig (2014).

Diameter (mm)	ARVOF-CSF	ARVOF-CSS	AR level set-CSF	AR level set-GFM
1	0.9595	0.9591	0.95	0.95
1.5	0.8340	0.8354	0.84	0.84
2	0.6567	0.6653	0.66	0.66
2.5	0.4876	0.4786	0.48	0.47
3	0.3591	0.3587	0.35	0.36

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 of Engberg and Kenig (2014).

Diameter (mm)	ARVOF-CSF	ARVOF-CSS	AR level set-CSF	AR level set-GFM
1	1.00	1.00	0.99	0.99
1.5	0.9830	0.9821	0.97	0.97
2	0.9316	0.9289	0.92	0.91
2.5	0.8140	0.7996	0.79	0.79
3	0.6618	0.6429	0.65	0.64
3.5	0.5509	0.5410	0.55	0.56

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 of Engberg and Kenig (2014).

Diameter (mm)	ARVOF-CSF	ARVOF-CSS	AR level set-CSF	AR level set-GFM
1	1.00	1.00	1	1
1.5	0.9930	1.00	0.99	0.99
2	0.9767	0.9779	0.96	0.95
2.5	0.8916	0.8969	0.89	0.88
3	0.7818	0.7693	0.77	0.76
3.5	0.6684	0.6350	0.65	0.65
4	0.5820	0.5470	0.58	-

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
Streamline of the different drop diameters in three liquid-liquid extraction systems using the CSS model.

Aspect Ratio

The Aspect ratio is the exact comparison parameter for drop shapes obtained from simulations. It is defined as the ratio of the drop vertical to horizontal diameter and calculated from Eq. 20 for the axis-symmetric conditions (Engberg and Kenig, 2014).

$$AR=\frac{y_{max}-y_{min}}{2\left(x_{max}-x_{min}\right)}$$ (20)

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

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.

The terminal velocity of n-butanol drops with both the CSF and CSS models is shown in Fig. 10A. It can be 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.

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

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.

Table 6
The average relative error of the present VOF method (VOF-CSS and VOF-CSF), level set model (level set-CSF and level set-GFM) of Engberg and Kenig (2014), dynamic mesh model of Bäumler et al. (2011), and Lattice-Boltzmann model of Komrakova et al. (2013) with respect to the experimental data of Bertakis et al.(2010) for n-butanol drops.

Method	% Average Relative Error
VOF-CSS	2.485
VOF-CSF	3.01
level set-CSF	3.059
level set-GFM	3.016
dynamic mesh	2.578
Lattice-Boltzmann	5.375

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.

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.

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

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.

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.

Method	% Average Relative Error
	Umax	Uchar
VOF-CSS	5.964	4.597
VOF-CSF	8.682	4.586
level set-CSF	5.127	2.911
level set-GFM	5.106	4.026
dynamic mesh	5.089	4.026

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

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.

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

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.

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

Method	% Average Relative Error
	Umax	Uchar
VOF-CSS	3.628	12.162
VOF-CSF	9.885	12.531
level set-CSF	2.463	11.661
level set-GFM	2.125	13.419
dynamic mesh	2.437	13.733

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 three-different 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.

NOMENCLATURE

• A_f  computational face surface in Eq. 16, m²
• AR  aspect ratio
• F_SF  the interfacial tension force, N/m³
• g  gravitational acceleration, m/s²
• I  identity matrix
• N  the normal vector of the boundary in Fig.1, m
• n  the normal vector of the interface in Eq. 12, 1/m
• $$\hat{n}$$  normal unit vector
• P  pressure, kg/ms²
• P*  modified pressure (Eq. 18), kg/ms²
• T  stress form of surface tension (Eq. 14), kg/ms²
• t  time, s
• U  velocity vector, m/s
• u  velocity component in the x direction in Fig.1, m/s
• v  velocity component in the y direction in Fig.1, m/s
• V  cell volume in Eq.2, m³
• V_f  volume flux through the face in Eq.2, m³/s
• v_i  velocity component in the y direction in the i^h cell (Eq. 19), m/s
• V_i  volume of the i^th cell (Eq. 19), m³
• V_cm,y  y component of the center of mass velocity (Eq.19), m/s
• X  position vector in Eq. 18, m
• x  x-axis of domain, m
• y  y-axis of domain, m

Dimensionless numbers

• Ca  Capillary number, $$Ca=\frac{{\mu }_c{u}_t}{\sigma }$$
• Eö  Eötvös number, $$Eö=\frac{g\Delta \rho d_p^2}{\sigma }$$
• N_f  Inverse viscosity number, $$N_f=\frac{\sqrt{{\rho }_c\left({\rho }_c-{\rho }_d\right)g{d}^3}}{{\mu }_c}$$
• Re  Reynolds number, $$Re=\frac{u_t{d}_p{\rho }_c}{{\mu }_c}$$
• We  Weber number, $$We=\frac{u_t^2{d}_p{\rho }_c}{\sigma }$$

Greek letters

• α  volume fraction
• μ  dynamic viscosity, Pa.s
• σ  interfacial tension, N/m
• ρ  density, kg/m³
• k  curvature in Eq.10, 1/m

Subscripts

• c  continuous phase
• d  dispersed phase, droplet
• max  maximum
• min  minimum

REFERENCES

• Abadie, T., Aubin, J., Legendre, D., On the combined effects of surface tension force calculation and interface advection on spurious currents within volume of fluid and level set frameworks, Journal of Computational Physics, 297, 611-636 (2015).
• Ahmadi Nadooshan, A., Shirani, E., Interface pressure model for surface tension force for VOF-based methods in interfacial flows, Engineering Applications of Computational Fluid Mechanics, 2, 496-513 (2008).
• Akbari, O.A., Safaei, M.R., Goodarzi, M., Akbar, N.S., Zarringhalam, M., Shabani, G.A.S., Dahari, M.A., Modified two-phase mixture model of nanofluid flow and heat transfer in a 3-D curved microtube, Advanced Powder Technology, 27, 2175-2185 (2016).
• Albert, C., Raach, H., Bothe, D., Influence of surface tension models on the hydrodynamics of wavy laminar falling films in Volume of Fluid-simulations, International Journal of Multiphase Flow, 43, 66-71 (2012).
• Alke, A., Bothe, D., Kröger, M., Warnecke, H., VOF-based simulation of conjugate mass transfer from freely moving fluid particles, Computational Methods in Multiphase Flow V, 157-168 (2009).
• Asadollahzadeh, M., Haghighi-Asl, A., Safdari, J., Torab-Mostaedi, M., Flooding characteristics in pulsed packed extraction columns, Brazilian Journal of Chemical Engineering, 28, 625-638 (2011).
• Bao, D., Zhang, X., Dong, H., Ouyang, Z., Zhang, X., Zhang, S., Numerical simulations of bubble behavior and mass transfer in CO2 capture system with ionic liquids, Chemical Engineering Science, 135, 76-88 (2015).
• Bäumler, K., Wegener, M., Paschedag, A., Bänsch, E., Drop rise velocities and fluid dynamic behavior in standard test systems for liquid/liquid extraction—experimental and numerical investigations, Chemical Engineering Science, 66, 426-439 (2011).
• Bertakis, E., Groß, S., Grande, J., Fortmeier, O., Reusken, A., Pfennig, A., Validated simulation of droplet sedimentation with finite-element and level-set methods, Chemical Engineering Science 65, 2037-2051 (2010).
• Bothe, D., Fleckenstein, S., A volume-of-fluid-based method for mass transfer processes at fluid particles, Chemical Engineering Science, 101, 283-302 (2013).
• Brackbill, J., Kothe, D.B., Zemach, C., A continuum method for modeling surface tension. Journal of Computational Physics, 100, 335-354 (1992).
• Eiswirth, R., Bart, H.-J., Atmakidis, T., Kenig, E., Experimental and numerical investigation of a free rising droplet. Chemical Engineering and Processing: Process Intensification, 50, 718-727 (2011).
• Engberg, R.F., Kenig, E.Y., Numerical simulation of rising droplets in liquid–liquid systems: A comparison of continuous and sharp interfacial force models, International Journal of Heat and Fluid Flow, 50, 16-26 (2014).
• Engberg, R.F., Kenig, E.Y., An investigation of the influence of initial deformation on fluid dynamics of toluene droplets in water, International Journal of Multiphase Flow, 76, 144-157(2015).
• Engberg, R.F., Wegener, M., Kenig, E.Y., The impact of Marangoni convection on fluid dynamics and mass transfer at deformable single rising droplets–A numerical study, Chemical Engineering Science, 116, 208-222 (2014a).
• Engberg, R.F., Wegener, M., Kenig, E.Y., The influence of Marangoni convection on fluid dynamics of oscillating single rising droplets, Chemical Engineering Science, 117, 114-124 (2014b).
• Ferreira, D., Arnao, J., Moreira, B., Rangel, L., Park, S., The impact of radiation in the gas combustion model for sugarcane bagasse grate boiler, Brazilian Journal of Chemical Engineering, 33, 617-626 (2016).
• Fleckenstein, S., Bothe, D., Simplified modeling of the influence of surfactants on the rise of bubbles in VOF-simulations, Chemical Engineering Science, 102, 514-523 (2013).
• Gerlach, D., Tomar, G., Biswas, G., Durst, F., Comparison of volume-of-fluid methods for surface tension-dominant two-phase flows, International Journal of Heat and Mass Transfer, 49, 740-754 (2006).
• Goodarzi, H.H., Esfahany, M.N., Experimental investigation of the effects of the hydrophilic silica nanoparticles on mass transfer and hydrodynamics of single drop extraction, Separation and Purification Technology, 170, 130-137 (2016).
• Goodarzi, M., Safaei, M., Vafai, K., Ahmadi, G., Dahari, M., Kazi, S., Jomhari, N., Investigation of nanofluid mixed convection in a shallow cavity using a two-phase mixture model, International Journal of Thermal Sciences, 75, 204-220 (2014).
• Hirt, C.W., Nichols, B.D., Volume of fluid (VOF) method for the dynamics of free boundaries, Journal of Computational Physics, 39, 201-225 (1981).
• Huang, Z., Ye, C., Li, L., Zhang, X., Qiu, T., Measurement and correlation of the mass transfer coefficient for a liquid-liquid system with high density difference, Brazilian Journal of Chemical Engineering, 33, 897-906 (2016).
• Issa, R.I., Solution of the implicitly discretised fluid flow equations by operator-splitting, Journal of Computational Physics, 62, 40-65 (1986).
• Jeon, S.-S., Kim, S.-J., Park, G.-C., CFD simulation of condensing vapor bubble using VOF model, World Academi of Science, Engineering and Technology, 60, 209-215 (2009).
• Komrakova, A., Eskin, D., Derksen, J., Lattice Boltzmann simulations of a single n-butanol drop rising in water, Physics of Fluids (1994-present), 25, 042102 (2013).
• Lafaurie, B., Nardone, C., Scardovelli, R., Zaleski, S., Zanetti, G., Modelling merging and fragmentation in multiphase flows with SURFER, Journal of Computational Physics, 113, 134-147 (1994).
• Leonard, B., Mokhtari, S., ULTRA-SHARP nonoscillatory convection schemes for high-speed steady multidimensional flow, NASA Technical Memorandum, 102568 (1990).
• Liu, L., Tang, H., Quan, S., Shapes and terminal velocities of a drop rising in stagnant liquids, Computers & Fluids, 81, 17-25 (2013).
• Moghadam, E.H., Bahmanyar, H., Heshmatifar, F., Ziaei-Azad, H., The investigation of mass transfer coefficients in a pulsed regular packed column applying SiO 2 nanoparticles, Separation and Purification Technology, 176, 15-22 (2017).
• Oliveira, T.C.L. de, Pereira Neto, A.T., Alves, J.J.N., Morais Júnior, A.A. de, CFD Simulation of an Industrial Reactor for Thermal Cracking of 1, 2-Dichloroethane, Brazilian Journal of Chemical Engineering, 34, 541-555 (2017)
• Osher, S., Sethian, J.A., Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics, 79, 12-49 (1988).
• Patankar, S. Numerical heat transfer and fluid flow. CRC press(1980).
• Pelissari, D., Alvarez-Castro, H., Mori, M., Martignoni, W., Study of feedstock injection to improve catalyst homogenization in the riser of a fcc, Brazilian Journal of Chemical Engineering, 33, 559-566 (2016).
• Pozzetti, G., Peters, B., A multiscale DEM-VOF method for the simulation of three-phase flows, International Journal of Multiphase Flow, 99, 186-204 (2017).
• Rahbar-Kelishami, A., Ashrafizadeh, S.N., Rahnamaee, M., The Effect of Type and Concentration of Nano-Particles on the Mass Transfer Coefficients: Experimental and Sherwood Number Correlating, Separation Science and Technology, 50, 1776-1784 (2015).
• Ranade, V.V., Computational flow modeling for chemical reactor engineering, Academic press(2001).
• Rusche, H., Computational fluid dynamics of dispersed two-phase flows at high phase fractions, Ph.D. Thesis, Imperial College London (University of London) (2003).
• Safaei, M., Togun, H., Vafai, K., Kazi, S., Badarudin, A., Investigation of heat transfer enhancement in a forward-facing contracting channel using FMWCNT nanofluids, Numerical Heat Transfer, Part A: Applications, 66, 1321-1340 (2014).
• Safaei, M.R., Jahanbin, A., Kianifar, A., Gharehkhani, S., Kherbeet, A.S., Goodarzi, M., Dahari, M., Mathematical modeling for nanofluids simulation: a review of the latest works. In Modeling and Simulation in Engineering Sciences. InTech, 189-220 (2016).
• Saien, J., Daliri, S., Mass Transfer from single drops and the influence of temperature. Industrial & Engineering Chemistry Research, 51, 7364-7372 (2012).
• Seifollahi, M., Shirani, E., Ashgriz, N., An improved method for calculation of interface pressure force in PLIC-VOF methods, European Journal of Mechanics-B/Fluids, 27, 1-23 (2008).
• Torab-Mostaedi, M., Ghaemi, A., Asadollahzadeh, M., Pejmanzad, P., Mass transfer performance in pulsed disc and doughnut extraction columns, Brazilian Journal of Chemical Engineering, 28, 447-456 (2011).
• Wegener, M., Kraume, M., Paschedag, A.R., Terminal and transient drop rise velocity of single toluene droplets in water, AIChE Journal, 56, 2-10 (2010).
• Weiner, A., Bothe, D., Advanced subgrid-scale modeling for convection-dominated species transport at fluid interfaces with application to mass transfer from rising bubbles, Journal of Computational Physics, 347, 261-289 (2017).
• Youngs, D.L., Time-dependent multi-material flow with large fluid distortion, Numerical Methods for Fluid Dynamics, 24, 273-285 (1982).

APPENDIX

Table A1
n-butanol drop terminal velocities and dimensionless numbers.

drop diameter (mm)	VOF-CSF mm/s	VOF-CSS mm/s	REcss	WEcss	Eö	Nf	Ca
1	27.4	27.61	19.60	0.46	0.61	26.57	0.024
1.5	45.69	45.72	48.67	1.90	1.36	48.81	0.039
1.56	45.71	45.95	50.87	1.99	1.47	51.77	0.039
1.79	51.88	52.19	66.30	2.95	1.94	63.63	0.045
2	56.15	56.53	80.24	3.87	2.42	75.15	0.048
2.12	56.28	58.2	87.57	4.35	2.72	82.01	0.05
2.26	59.51	59.62	95.63	4.86	3.09	90.27	0.051
2.48	60.57	60.84	107.08	5.56	3.72	103.76	0.052
2.5	60.65	60.73	107.75	5.58	3.78	105.02	0.052
3	61.72	61.75	131.48	6.92	5.45	138.05	0.053
3.06	61.66	61.6	133.78	7.03	5.67	142.22	0.053
3.48	62.55	62.19	153.60	8.15	7.33	172.48	0.053
4	59.28	59.24	168.18	8.50	9.68	212.55	0.051

Table A2
n-butyl acetate drop terminal velocities and dimensionless numbers.

drop diameter (mm)	VOF-CSF mm/s	VOF-CSS mm/s	REcss	WEcss	Eö	Nf	Ca
1	28.45	35.14	38.74	0.09	0.07	37.78	0.002
1.5	55.36	59.76	98.82	0.38	0.16	69.40	0.004
2	84.69	87.61	193.17	1.09	0.28	106.85	0.006
2.5	111.92	114.01	314.22	2.31	0.44	149.33	0.007
3	122.17	124.0	410.11	3.28	0.64	196.29	0.008
3.5	122.06	123.0	474.60	3.77	0.87	247.36	0.008
4	119.0	119.64	527.58	4.08	1.14	302.21	0.008
4.4		118.61	575.34	4.41	1.38	348.66	0.008
6		111.032	734.43	5.27	2.56	555.20	0.007

Table A3
Toluene drop terminal velocities and dimensionless numbers.

drop diameter (mm)	VOF-CSF mm/s	VOF-CSS mm/s	REcss	WEcss	Eö	Nf	Ca
1	32.0	42.45	47.54	0.05	0.03	40.75	0.0011
1.5	59.45	71.92	120.81	0.22	0.06	74.87	0.0018
2	94.34	110.25	246.93	0.69	0.11	115.26	0.0028
2.5	142.76	146.37	409.79	1.53	0.18	161.08	0.0037
3	167.602	174.72	586.99	2.61	0.26	211.75	0.0044
3.5	175.85	178.11	698.11	3.16	0.35	266.84	0.0045
4	174.11	174.72	782.65	3.48	0.46	326.01	0.0044
4.4	170.84	170.08	838.06	3.63	0.55	376.12	0.0043
5		163.29	914.32	3.80	0.71	455.62	0.0042
6		155.71	1046.25	4.14	1.03	598.92	0.0040
7		149.96	1175.55	4.48	1.40	754.73	0.0038

Publication Dates

• Publication in this collection
  Dec 2018

History

• Received
  27 Nov 2017
• Reviewed
  03 Feb 2018
• Accepted
  20 Feb 2018