**1.4. Issues of numerical simulation**

422 Numerical Simulation – From Theory to Industry

and space, expressed as

across the interface.

equation (13)

interface only.

method is based on the conservation of the volume fraction function *F* with respect to time

In VOF method, the computational grid is kept fixed and the interface between the two fluids is tracked within each cell through which it passes. In a computational grid cell, the interface can be effectively represented by line of the slope. To reconstruct the interface, the piecewise linear interpolation calculation (PLIC) method developed by Youngs (1982) is used in the computation. The interfacial surface forces are incorporated as body forces per unit volume in the Navier-Stokes equations; hence no extra boundary condition is required

The basic working principle of front tracking method is based on the marker and cell (MAC) formulation (Harlow and Welch, 1966; Daly, 1967). The interface is represented discretely by Lagrangian markers connected to form a front which lies within and moves through a stationary Eulerian mesh. As the front moves and deforms, interface points are added, deleted, and reconnected as necessary. Further details of the method may be found in

Level set methods have been developed by Osher and Sethian (1988). This method can compute the geometrical properties of highly complicated interface without explicitly tracking the interface. The basic principle of the level set method is to embed the propagating interface Γ(t) as the zero level set of a higher dimensional function *φ* , defined as *φ* (*x*, *t* =0) = ±*d*, where *d* is the distance from *x* to *Γ* (*t* = 0). The function is chosen to be positive (negative) if *x* is outside (inside) the initial position of the interface *Γ* (*t* = 0) = *φ* (*x*, *t* =0)=0. Afterwards a dynamical equation for *φ* (*x*, *t* =0) that contains the embedded motion for *Γ*(*t*) as the level set *φ* =0 can be derived similarly as in the volume of fluid conservation

In phase field method, interfacial forces are modeled as continuum forces by smoothing interface discontinuities and forces over thin but numerically resolvable layers. This smoothing allows conventional numerical approximations of interface kinematics on fixed grids. The method has been used for investigating the problems governed by Navier-Stokes

Boundary integral methods are designed to track the interface explicitly, as in front tracking methods, although the flow solution in the entire domain is deduced solely from information possessed by discrete points along the interface. The advantage of these methods is the reduction of the flow problem by one dimension involving quantities of the

Particle-based methods use discrete "particles" to represent macroscopic fluid parcels. Here, Lagrangian coordinates are used to solve the Navier-Stokes equations on "particles" having properties such as mass, momentum, and energy. The nonlinear convection term is modeled simply as particle motion and by knowing the identity and position of each particle,

�� <sup>+</sup> (�� �)��� (13)

��

(Glimm et al., 1985; Churn et al., 1986; Tryggvason et al., 1998).

equations (Antanovskii, 1995; Jacqumin, 1996).

Numerical investigation of active droplet formation is reasonably complex as it is potentially a multi-physics problem governed by a large number of partial differential equations (PDEs). Numerical simulation of droplet dynamics in membrane emulsion process requires accurate capturing of the evolving liquid–liquid interfaces. This gives sufficient challenge since the boundary between the two phases is not known priori and it is a part of the solution. In addition, the solution technique has to deal with different properties of both the phases such as density, viscosity, and velocity ratios. Numerical investigation of membrane emulsification process requires two main issues to be dealt with: one is permeation of the disperse phase through the membrane pores; the other is the mechanism of droplet detachment. Several earlier investigations have treated these issues separately; however both of them simultaneously contribute to droplet evolution and should be considered in the same framework. Both these issues have been considered in the present work. In the present work, numerical simulation of droplet dynamics in a membrane emulsification has been carried out considering multipore membrane. The hydrodynamic effects due to multipores and the effect of different operating parameters on the droplet dynamics have been investigated.

## **2. Problem description**

The flow configuration of membrane emulsification considered in the present work has been shown in Fig. 3. Here the membrane emulsification process with two pores has been considered. A uniform pore arrangement in the membrane surface has been considered and the simulation has been mode for a row of pores consisting of two pores. The dispersed phase liquid has been injected through cylindrical pores of 10 *μ*m diameter and a length of 100 *μ*m. The distance between the pores in cross-flow direction has been considered as 100 *μ*m. The height of the rectangular channel through which the continuous phase flows has been considered as 150 *μ*m (*y*-direction). The width of the computational domain has been considered as 150 *μ*m (*z*-direction) and a length of 500 *μ*m (*x*-direction) has been considered. Since the present computational domain is a small element of whole membrane emulsification process, symmetrical boundary conditions have been considered in both the sides of computational domain in *z*-direction. The fluid properties used in the present simulation are within the range of properties of o/w emulsion system. Simulations have been made for different values of the non-dimensional numbers and other flow properties as shown in Table 1.

Numerical Simulation of Droplet Dynamics in Membrane Emulsification Systems 425

*f* is the surface force, which includes surface tension force of the interface of the

Parameters Ranges Diameter of pore 10 *μ*m

Viscosity of continuous phase 0.001 Pa.s

Density of continuous phase 1000 kg/m3 Density of dispersed phase 827 kg/m3

**Table 1.** Values of physical properties used in the simulation

) are defined as:

distribution of volume fraction is solved from its transport equation.

**2.2. Interface tracking** 

viscosity (

surface force *<sup>S</sup>*

and

curvature is calculated from:

where

 and 

Flow rate of continuous phase 0.27 and 0.54 liter/h Flow rate of dispersed phase 0.0014 to 0.007 liter/h

Viscosity of dispersed phase 0.0036 to 0.014 Pa.s

Surface tension 0.0008 to 0.0024 N/m Weber no (*We*) 0.0021 to 0.215 Capillary no (*Ca*) 0.0208 to 0.0625 Froude no (*Fr*) 6.37 to 637

Volume of fluid method (VOF) has been used for tracking the interface. In this method, the

. 0 *<sup>F</sup> V F*

With the inclusion of volume fraction in the calculation, the volume averaged density and

 

 

1 21

 

1 21

*f* appearing in momentum equation has been calculated using continuous

are the surface tension and surface curvature respectively. Surface

 

where subscript 1 and 2 denote the continuous and dispersed phase respectively. The

, *<sup>S</sup> f F* 

> . . *<sup>F</sup> F*

(17)

*F*( ) (18)

*F*( ), (19)

(20)

(21)

*t* 

surface force (CSF) model proposed by Brackbill et al. (1992) . It is defined as:

In above *<sup>S</sup>*

two fluids.

**Figure 3.** Schematic of computational domain

#### **2.1. Governing equations**

In the simulation of membrane emulsification system, both the phases have been considered incompressible, isothermal and laminar flow. The phase properties and the surface tension have been assumed to be constant throughout the flow domain. The conservation of mass on the whole domain (both the fluid phases and interface) leading to continuity equation is written as:

$$\frac{\partial \rho}{\partial t} + \nabla.(\rho V) = 0\tag{15}$$

The momentum equation or unsteady Navier-Stokes equation is written as:

$$
\rho \left( \frac{\partial V}{\partial t} + V .\nabla V \right) = -\nabla p + \rho \,\mathrm{g} + \nabla \cdot \left[ \mu (\nabla V + \nabla V^T) \right] + f\_S \tag{16}
$$

In above *<sup>S</sup> f* is the surface force, which includes surface tension force of the interface of the two fluids.


**Table 1.** Values of physical properties used in the simulation

#### **2.2. Interface tracking**

424 Numerical Simulation – From Theory to Industry

as shown in Table 1.

**Figure 3.** Schematic of computational domain

 **2.1. Governing equations** 

written as:

considered. A uniform pore arrangement in the membrane surface has been considered and the simulation has been mode for a row of pores consisting of two pores. The dispersed phase liquid has been injected through cylindrical pores of 10 *μ*m diameter and a length of 100 *μ*m. The distance between the pores in cross-flow direction has been considered as 100 *μ*m. The height of the rectangular channel through which the continuous phase flows has been considered as 150 *μ*m (*y*-direction). The width of the computational domain has been considered as 150 *μ*m (*z*-direction) and a length of 500 *μ*m (*x*-direction) has been considered. Since the present computational domain is a small element of whole membrane emulsification process, symmetrical boundary conditions have been considered in both the sides of computational domain in *z*-direction. The fluid properties used in the present simulation are within the range of properties of o/w emulsion system. Simulations have been made for different values of the non-dimensional numbers and other flow properties

In the simulation of membrane emulsification system, both the phases have been considered incompressible, isothermal and laminar flow. The phase properties and the surface tension have been assumed to be constant throughout the flow domain. The conservation of mass on the whole domain (both the fluid phases and interface) leading to continuity equation is

.( ) 0 *V*

. .( ) *<sup>T</sup>*

 

(16)

*<sup>V</sup> VV p g V V f <sup>t</sup>*

(15)

*S*

*t* 

The momentum equation or unsteady Navier-Stokes equation is written as:

Volume of fluid method (VOF) has been used for tracking the interface. In this method, the distribution of volume fraction is solved from its transport equation.

$$\frac{\partial F}{\partial t} + V.\nabla F = 0 \tag{17}$$

With the inclusion of volume fraction in the calculation, the volume averaged density and viscosity ( and ) are defined as:

$$
\rho = \rho\_1 + F(\rho\_2 - \rho\_1) \tag{18}
$$

$$
\mu = \mu\_1 + \mathcal{F}(\mu\_2 - \mu\_1),
\tag{19}
$$

where subscript 1 and 2 denote the continuous and dispersed phase respectively. The surface force *<sup>S</sup> f* appearing in momentum equation has been calculated using continuous surface force (CSF) model proposed by Brackbill et al. (1992) . It is defined as:

$$f\_S = -\sigma \mathbf{k} \nabla F\_{\prime} \tag{20}$$

where and are the surface tension and surface curvature respectively. Surface curvature is calculated from:

$$
\kappa = \nabla. \frac{\nabla F}{\left| \nabla F \right|}. \tag{21}
$$

### **2.3. Boundary and initial conditions**

Various types of boundary conditions have been used in the simulation of the membrane emulsification. A fully developed laminar duct flow has been considered at the inlet to the continuous phase channel and the flow has been assumed to be dominant along the stream direction thus: *cp u u* , *v* = 0, *w* = 0. At the inlet of continuous phase, the value of volume fraction has been set to zero (*F* = 0). At the inlets to the micro-pores, a velocity inlet type boundary condition has been used: *v* = *v0*. The other components of the velocity have been put as *u* = 0, *w* = 0 and the value of volume fraction has been set to *F* = 1. Outflow boundary condition has been put at the outlet of the flow channel. Symmetry boundary conditions have been put at both sides of the computational domain. The top of the channel has been treated as walls where no-slip and impermeability boundary conditions have been set. At the starting of the computation, the whole flow domain has been assumed to be filled up with the continuous phase fluid. The static contact angle between the two phases and the wall influences the droplet growth in the dispersed and continuous phase channel. It has been observed that the effect of contact angle in the droplet formation is negligible when its value is greater than 165o (Sang et al., 2009). The value of static angle used in the present simulation has been set to 170o.

To reduce the number of dependable variables, the governing equations have been expressed in dimensionless form. The diameter of the pore ( *D*<sup>0</sup> ) has been used as length scale and the average velocity ( <sup>0</sup> *v* ) of the dispersed phase has been used as velocity scale for making non-dimensional form of the equation. The non-dimensional equations are:

$$\frac{\partial \stackrel{\star}{\rho}^\*}{\partial t^\*} + \nabla. (\rho^\* V^\*) = 0 \tag{22}$$

Numerical Simulation of Droplet Dynamics in Membrane Emulsification Systems 427

) and *R* is the velocity ratio ( <sup>0</sup> / *Ru v cp* ).

Re

 

*cp cp <sup>u</sup> RWe Ca* 

Commercial code Ansys Fluent (V12) based on finite volume method has been used in the simulation. The momentum and volume fraction transport equation have been discretized with 2nd order upwind scheme. The PISO (pressure implicit with splitting of operators) algorithm has been used for pressure correction. The VOF/CSF techniques have been used to track the fluid interface between the two immiscible fluids. A geometry reconstruction scheme has been used in the simulation to avoid the diffusion at the interface. The interface was reconstructed by the piecewise-linear interface calculation (PLIC) technique (Youngs, 1982). The unsteady term was treated with first-order implicit time stepping. Simulations were made with very small time steps (~10-7 s). The solutions have been assumed to be converged and therefore iterations have been terminated when the normalized sum of residual mass was less than 10-4 and variation of other variables in successive iteration was less than 10-2. A non-uniform grid was used in the simulation where grids were clustered near the walls and the injection portion of the dispersed phase. The channel was decomposed into 12 105 numbers of control volumes and the pore with 12 103 after a grid

In the present work, the dynamics of droplet formation in two pores of membrane emulsification has been investigated for different flow rates (velocities) of dispersed phase, continuous phase, surface tension and viscosity of the two phases. It is to be noted that In case of membrane emulsification process, the dispersed phase fluid gets more space to interact with the continuous phase fluid compared to the confined geometries in case of Tjunction emulsification. Due to this the evolution of dispersed phase becomes different than the case of T-junction emulsion and the dependence of the process on different properties of

In membrane emulsification system, the flow rate of dispersed phase controls the droplet dynamics via its inertial force competing with the drag force imparted by the continuous phase and interfacial tension force. In order to investigate the effect of dispersed phase flow rate, simulations have been made for different values of *We* number by changing the dispersed phase velocity i.e. inertial force and keeping the surface tension force fixed. Before discussing the effects of *We*, the growth and detachment of the droplet for a constant values of *We* (0.0086) and *Ca* (0.028) at different time levels have been shown in Fig. 4. It has been observed that both the droplet grow at their respective micro-pore and detached by the

where 

**2.4. Numerical method** 

independent study.

**3. Results and discussions** 

both the phases also changes.

**3.1. Growth and detachment of the droplets** 

is the viscosity ratio ( / *dp cp*

 

$$\nabla \left( \frac{\partial V^\*}{\partial t} + V^\* \nabla V^\* \right) = -\nabla p^\* + \frac{\rho^\*}{Fr} + \frac{1}{\text{Re}} \nabla \cdot \left[ \left( \nabla V^\* + \nabla V^{\*T} \right) \right] + \frac{1}{\text{We}} \cdot \nabla \cdot \frac{\nabla F}{\left| \nabla F \right|} \tag{23}$$

In above starred quantities are non-dimensional parameters. The three non-dimensional numbers appeared in the problem are: Reynolds number (*Re*) , Weber number (*We*) and Froude number (*Fr*). They are defined as:

$$Re = \frac{\rho\_{dp} v\_0 D\_0}{\mu\_{dp}} ; \text{ } \text{ } We = \frac{\rho\_{dp} v\_0^2 D\_0}{\sigma} ; \text{ } Fr = \frac{v\_0^2}{g D\_0} \text{ }$$

The Reynolds numbers (*Re*) is defined as the ratio of viscous force to the inertia force, Weber number (*We*) is defined as the ratio of inertia force relative to the surface tension force and Froude number (*Fr*) is defined as the ratio of inertia to gravity force. In the present work to incorporate the effect of the continuous phase fluid, the Capillary number has been introduced and is defined as:

$$Ca = \frac{\mu\_{cp}\mu\_{cp}}{\sigma} = \frac{RWe}{\lambda \text{Re}}$$

where is the viscosity ratio ( / *dp cp* ) and *R* is the velocity ratio ( <sup>0</sup> / *Ru v cp* ).
