<figure><img src="image_1.png" /></figure>

HAL
open science

Numerical study of oil-water emulsion formation in
stirred vessels: effect of impeller speed

Fuyue Liang, Lyes Kahouadji, Juan Pablo Valdes, Seungwon Shin, Jalel

Chergui, Damir Juric, Omar Matar

► To cite this version:

Fuyue Liang, Lyes Kahouadji, Juan Pablo Valdes, Seungwon Shin, Jalel Chergui, et al.. Numerical
study of oil-water emulsion formation in stirred vessels: effect of impeller speed.
FLOW, 2022, 2,
pp.E34. 10.1017/flo.2022.1 . hal-03799387

HAL Id: hal-03799387

https://hal.science/hal-03799387

Submitted on 5 Oct 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Numerical study of oil-water emulsion formation in stirred vessels: effect of impeller speed

Fuyue Liang<fnref n="1" /> Ⓓ ID, Lyes Kahouadji<fnref n="1" />* Ⓓ ID, Juan Pablo Valdes<fnref n="1" /> Ⓓ ID, Seungwon Shin<fnref n="2" /> Ⓓ ID, Jalel Chergui<fnref n="3" /> Ⓓ ID, Damir Juric<fnref n="3" />,<fnref n="4" /> Ⓓ ID and Omar K Matar<fnref n="1" /> Ⓓ ID.

**Received:** XX 2022; **Revised:** XX XX 2022; **Accepted:** XX XX 2022

**Keywords:** Oil-water Emulsion, Multiphase mixing, LES

## Abstract

The mixing of immiscible oil and water by a pitched blade turbine in a cylindrical vessel is studied numerically. Three-dimensional simulations combined with a hybrid front-tracking/level-set method are employed to capture the complex flow and interfacial dynamics. A large eddy simulation approach, with a Lilly-Smagorinsky model, is employed to simulate the turbulent two-phase dynamics at large Re = 5663 – 56632. The numerical predictions are validated against previous experimental work involving single-drop breakup in a stirred vessel. For small Re, the interface is deformed but does not reach the impeller hub, assuming instead the shape of a Newton’s Bucket. As the rotating speed increases, the deforming interface attaches to the impeller hub which leads to the formation of long ligaments that subsequently break up into small droplets. For the largest Re studied, the system dynamics becomes extremely complex wherein the creation of ligaments, their breakup, and the coalescence of drops occur simultaneously. The simulation outcomes are presented in terms of spatio-temporal evolution of the interface shape and vortical structures. The results of a drop size analysis in terms of the evolution of the number of drops, and their size distribution, is also presented as a parametric function of Re.

**Impact Statement** Emulsions are of primary interest in modern engineering research given their heavy usage both in our daily life (e.g., food and cosmetics) and in high-end industrial applications (e.g., drug development). An insufficient understanding of mixing dynamics has commonly led to over-design practices to guarantee a certain threshold of quality and performance, which can be easily overlooked causing severe economical and environmental consequences. Attention has been focused previously on the effect of operational parameters (e.g., impeller speed) on mixing performance, which, however, has led to limited insight due to the difficulty in obtaining clear visualisation of the liquid-liquid mixing process. The results from this study are intended to provide a detailed picture of the underlying physics and interfacial dynamics within mixing systems as a function of impeller speed, attained through our high-fidelity three-dimensional numerical simulations. Based on this, substantial improvement to optimal operation at the design stage can be achieved.

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

<footnote n="1">¹Department of Chemical Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, UK</footnote>

<footnote n="2">²Department of Mechanical and System Design Engineering, Hongik University, Seoul 04066, Republic of Korea</footnote>

<footnote n="3">³Université Paris Saclay, Centre National de la Recherche Scientifique (CNRS), Laboratoire Interdisciplinaire des Sciences du Numérique (LISN), 91400 Orsay, France</footnote>

<footnote n="4">⁴Department of Applied Mathematics and Theoretical Physics, University of Cambridge, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 0WA, UK</footnote>

<footnote>*Corresponding author. E-mail: l.kahouadji@imperial.ac.uk</footnote>
# 1. Introduction

The applications of the mixing of two immiscible liquids by a stirrer range from the production of daily necessities (e.g., food and cosmetics) to large-scale industrial processing (e.g., petrochemicals and pharmaceuticals). Mixing in a stirred vessel is the most common operation unit in this field where a rotating impeller controls the continuous-phase flow dynamics and exerts an influence on the dispersed-phase behaviour. The flow dynamics within such a system has attracted research interest for decades since it determines the interfacial area and, consequently, the product quality in addition to the efficiency of heat and mass transfer between the two liquid phases. At sufficiently large impeller rotational speeds, turbulent eddies are introduced by the rotating impeller, which subsequently interact with the dispersed phase giving rise to dispersed-phase deformation, ligament production, and, finally, drop breakup; the competition between drop breakup and coalescence determines the dispersed drop numbers, their size distribution, and therefore the interfacial area. From an industrial standpoint, the inherent complexity of these processes hamper the understanding and control of the flow dynamics and interface behaviours within a mixing vessel which are important for product design and scale-up.

Numerous researchers have conducted studies on the characteristic drop size distribution during the emulsification process (with a focus on $d_{max}$, the maximum stable drop diameter, and $d_{32}$, the Saunter mean drop diameter, typically) and investigated the effect of various parameters (e.g., the physical properties of the two phases, operating conditions, and the design of the vessel and impeller) on drop breakup in stirred vessels; this has led to several reviews of this field (see, for instance, the work of Afshar Ghotli et al. (2013) and Hasan (2017) and references therein). There have been many experiments on drop size distribution in liquid-liquid dispersions exemplified by the work of Wang & Calabrese (1986); Lovick et al. (2005); El-Hamouz et al. (2009); Boxall et al. (2010); Becker et al. (2014); De Hert & Rodgers (2017); Naeeni & Pakzad (2019b). On the other hand, Hinze (1955) carried out pioneering work in this field proposing two dimensionless groups, the Weber number, *We*, and the viscosity group, *Vi*, to account for the contributions of the external and the dispersed-phase viscous forces to the drop breakup. The analysis of the correlation for $d_{32}$ and $d_{max}$ started with Sprow (1967) in which it was verified through his experimental data that $d_{32} = Cd_{max}$ where *C* is a proportionality constant. Later, Zhou & Kresta (1998a) reported that *C* varies with operating conditions and fluid properties. Though Pacek et al. (1998) doubted the relationship between $d_{32}$ and $d_{max}$, several investigators adopted the assumption and modified it by considering the effect of dispersed phase viscosity and volume fraction (Calabrese et al., 1986; Kraume et al., 2004). More recently, the combination of population balance models with CFD (Ramkrishna, 2000) has become widely used to predict the drop size distribution in stirred vessels (Roudsari et al., 2012; Qin et al., 2016; Naeeni & Pakzad, 2019a).

The above mentioned determination of drop size distribution, however, does not provide clear insights into the interfacial phenomena during the breakup of drops and their re-coalescence in complex flow fields. Previous studies have reported that a dispersed drop experiences various stages before reaching an equilibrium size which include random deformation, stretching, and elongation (Andersson & Andersson, 2006; Sanjuan-Galindo et al., 2015; Nachtigall et al., 2016). Although numerous experiments have been carried out to study drop breakup in stirred vessels (Stamatoudis & Tavlarides, 1985; Calabrese et al., 1986; Wang & Calabrese, 1986; Giapos et al., 2005; Maaß et al., 2010; Zhou & Kresta, 1998b), a number of contradictory conclusions remain in the literature. For example, Sathyagal et al. (1996) have found that the breakup rate decreases considerably with increasing dispersed-phase viscosity and interfacial tension, with a stronger dependence on the latter, while Kraume et al. (2004) suggested that the viscosity exerts a larger influence on the drop breakup. These contradictions may be ascribed to differences in the experimental setups used to capture the breakup process.

The results of single-drop breakup experiments have also been reported in the literature. Konno et al. (1983) carried out the first single-drop experiment in three stirred vessels of different diameters where the authors tracked the drop breakup trajectory and reported the path length along with the distribution of breakup positions. Maaß et al. (2007) used a more sophisticated breakup cell device in which the continuous phase fluid flows through a fixed single blade in a rectangular channel to produce a

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51
comparable flow field to a stirred liquid-liquid system, and a single drop is injected to mimic the breakup phenomena; this cell facilitated the recording of the breakup details. Subsequent studies adopted this setup and reported data on breakup event location (Maaß et al., 2009), rate (Maaß & Kraume, 2012), mechanisms (Nachtigall et al., 2012), and deformation dynamics (Nachtigall et al., 2016).

In terms of numerical predictions of flows in stirred vessels, work in this area must overcome the significant challenges associated with simulating interfacial deformation in a drop-laden turbulent flow driven by a rotating impeller. Nevertheless, various approaches have been used including the volume-of-fluid (Hirt & Nichols, 1981), phase field (Anderson et al., 1998; Chiu & Lin, 2011), front tracking (Unverdi & Tryggvason, 1992; Tryggvason et al., 2001), and level set methods (Sethian & Smereka, 2003). On the other hand, a large number of modelling techniques for stirred vessels have been developed such as the black-box (Harvey et al., 1982), the multiple reference frame (MRF) (Luo & Gosman, 1994), and sliding mesh (Bakker et al., 1997) methods. The Lattice Boltzmann Method (LBM) has also been used to simulate flows in stirred vessels (Eggels, 1996; Derksen & Van den Akker, 1999; Derksen, 2003; Tyagi et al., 2007; Derksen, 2011, 2012; Gillissen & Van den Akker, 2012; Shu & Yang, 2018); these studies have yielded flow field information only rather than the detailed interfacial dynamics. Scarbolo & Soldati (2013); Scarbolo et al. (2015) used a phase field approach to study drops deforming in a fully-developed turbulent channel flow with a focus on the effect of the turbulent field on the droplet interface. Scarbolo et al. (2016) built on this previous work to study turbulence modification by the dispersed phase; subsequently, Roccon et al. (2017) extended this study to show that an increase in viscosity or interfacial tension decreases the breakup rate. Albernaz et al. (2017) used a hybrid LBM to study the deformation of a single drop in stationary isotropic turbulence while Komrakova (2019) applied a diffuse interface free energy LBM to examine single-drop breakup in homogeneous isotropic turbulence. These high-resolution simulations confirmed the presence of previously reported breakup mechanisms and elucidated an additional, so-called “burst breakup” mechanism.

As the foregoing review indicates, there remains a gap in our understanding of all the stages of the emulsification process, from the onset of drop breakup through to the attainment of a dynamic steady state. Bridging this gap requires a transparent connection between the interfacial deformation, the prevailing flow conditions (parameterised by the impeller rotational speed) and their associated vortical structures. Although previous studies have suggested that increasing the impeller speed produces smaller droplets (Zhou & Kresta, 1998a; El-Hamouz et al., 2009; Naeeni & Pakzad, 2019a), such a connection has not been established. In an effort to address this in the present study, we will deploy massively parallel, three-dimensional, large eddy simulations of oil and water emulsification in a cylindrical vessel stirred by a pitched blade turbine, combined with a hybrid front-tracking and level-set interface-capturing algorithm. This approach will provide detailed, high-fidelity visualisations of the intricate interfacial dynamics coupled to the turbulent flow fields, and will elevate our fundamental understanding, and, ultimately, our control of the emulsification process. Allied to this, an analysis of the drop size distribution as a function of the impeller speed will also be carried out; this analysis demonstrates that the number of dispersed drops is maximised for an intermediate value of the impeller rotational speed.

The rest of this article is organised as follows: Section 2 describes the configurations, governing equations, and the computational methods employed. Section 3 provides a comprehensive discussion of the vortical structures and interfacial dynamics inside the stirred vessel as a function of the impeller speed, followed by an analysis of the drop size distribution. Finally, conclusions are summarised and ideas for future work are outlined in Section 4.

## 2. Simulation Configuration and Methods

The configuration considered in this work is shown in Figure 1-(a). It is composed of a cylindrical vessel of diameter $D_T = 8.5$ cm filled with water in the lower half and oil in the upper half (volume fraction of oil, $\alpha = 0.5$). The impeller employed is a pitched blade turbine (PBT), which consists of four blades of 2.5 cm length, 1 cm height, and 0.2 cm thickness. The PBT is immersed in the water phase with a clearance of $C = 1$ cm from the bottom of vessel, and rotates with frequency $f$. The impeller

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
diameter is $D_i = 4.25$ cm which corresponds to an impeller-to-vessel diameter ratio of $D_i/D_T = 0.5$. A snapshot from a typical simulation is displayed in Figure 1-(b), which provides an example of the level of complexity associated with the flows under consideration in the present work.

Under the assumptions of incompressible and immiscible viscous fluids, the mass and momentum conservation equations are solved in a Cartesian domain $\mathbf{x} = (x, y, z) \in [0, 8.5]^2 \times [0, 12.75]$ cm:

$$ \nabla \cdot \bar{\mathbf{u}} = 0, \qquad (1) $$

$$ \rho \left( \frac{\partial \bar{\mathbf{u}}}{\partial t} + \bar{\mathbf{u}} \cdot \nabla \bar{\mathbf{u}} \right) = -\nabla \bar{\rho} + \nabla \cdot \left[ \left( \mu + \rho C_s^2 \Delta^2 |\bar{S}| \right) \left( \nabla \bar{\mathbf{u}} + \nabla \bar{\mathbf{u}}^T \right) \right] + \rho \mathbf{g} + \mathbf{F} + \mathbf{F}_{\text{fsi}}, \qquad (2) $$

where $\bar{\mathbf{u}}$ and $\bar{\rho}$ are the ensemble-averaged fluid velocity and pressure, respectively, $t$ denotes time, and $\mathbf{g}$ is the gravitational acceleration. In Eq. (2), we use a single-field formulation for the density $\rho$ and viscosity $\mu$ expressed by:

$$ \rho(\mathbf{x}, t) = \rho_o + (\rho_w - \rho_o) \mathcal{H}(\mathbf{x}, t), \qquad (3) $$

$$ \mu(\mathbf{x}, t) = \mu_o + (\mu_w - \mu_o) \mathcal{H}(\mathbf{x}, t), \qquad (4) $$

where the subscripts $o$ and $w$ designate the oil and water phases, respectively. The indicator function $\mathcal{H}$, is essentially a numerical Heaviside function, zero in the oil phase and unity in the water phase. Here, $\mathcal{H}$ is resolved with a sharp but smooth transition across 3-4 grid cells, and is generated using a vector distance function $\varphi(\mathbf{x})$, positive for the water phase and negative for the oil phase, computed directly from the tracked interface (Shin & Juric, 2009a; Shin et al., 2017). The density and viscosity of the water and oil phases are 998 kg/m$^3$ and $10^{-3}$ Pa.s, and 824 kg/m$^3$ and $5.4 \times 10^{-3}$ Pa.s, respectively; the interfacial tension is $\sigma = 0.035$ N/m. The oil phase corresponds to a type of silicone-oil used in previous work (Ibarra, 2017; Constante-Amores, 2021; Constante-Amores et al., 2021a).

<figure><img src="image_2.png" /><figcaption>Figure 1: (a) Schematic illustration of the computational domain for the oil-water mixing system, which corresponds to a stirred vessel filled with oil and water in the upper and lower halves, respectively, with a pitched blade turbine immersed in the water phase. The domain is of size $8.6 \times 8.6 \times 12.75$ cm$^3$ and it is divided into $4 \times 4 \times 6$ subdomains. The Cartesian structured grid per subdomain is 64$^3$, which gives a total grid number of $256 \times 256 \times 384$. (b) An illustrative snapshot of the interface at $t = 25 \times T$ coloured by velocity field for $f = 7$ Hz; here, a quarter of the interface is hidden to provide a clear view of the impeller and the interface in its vicinity.</figcaption></figure>

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
A range of impeller rotation frequencies, $f = 1-10$ Hz, is studied which corresponds to $Re = 5663-56632$ and $We = 22-2160$ where the Reynolds and Weber numbers are respectively given by

$$Re = \frac{\rho_w \pi f D_i^2}{\mu_w}, \quad We = \frac{\rho_w \pi^2 f^2 D_i^3}{\sigma}, \qquad (5)$$

which provide a measure of the relative significance of inertial to viscous and interfacial tension forces. Here, $D_i$ and $\pi f D_i$ are the characteristic length and velocity scales, respectively. The large Reynolds numbers in the upper end of the Re range motivate the use of a large eddy simulation approach; here, we employ a Smagorinsky turbulence model (Smagorinsky, 1963) which is included in Eq. (2) wherein $C_s$ is the Smagorinsky-Lilly coefficient, $\Delta = V^{1/3}$ where $V$ is the volume of a grid cell, $V = \Delta x \Delta y \Delta z$ and $|\bar{S}| = \sqrt{2\bar{S}_{ij}\bar{S}_{ij}}$ with $\bar{S}_{ij}$ being the strain rate tensor

$$\bar{S}_{ij} = \frac{1}{2} \left( \frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i} \right) \qquad (6)$$

(Pope & Pope, 2000; Pope, 2004; Meyers & Sagaut, 2006). In the current work, $C_s = 0.2$ which is an intermediate value in the range $0.1-0.3$ quoted in the literature (Lilly, 1966, 1967; Deardorff, 1970; McMillan & Ferziger, 1979; Pope & Pope, 2000); setting $C_s = 0$ eliminates the turbulence model and Eq. (2) reduces to the Navier-Stokes equations.

The forces $\mathbf{F}$ and $\mathbf{F}_{\text{fsi}}$ in Eq. (2) denote the local surface tension and the solid-fluid interaction forces, respectively, where $\mathbf{F}$ is defined using a hybrid formulation (Shin et al., 2017; 2018):

$$\mathbf{F} = \sigma \kappa_H \nabla \mathcal{H}(\mathbf{x}, t), \qquad (7)$$

in which $\sigma$, the surface tension, is assumed constant. In Eq. (7), $\kappa_H$ is twice the mean interface curvature field calculated on an Eulerian grid using:

$$\kappa_H = \frac{\mathbf{F}_L \cdot \mathbf{G}}{\sigma \mathbf{G} \cdot \mathbf{G}}, \qquad (8)$$

in which $\mathbf{F}_L$ and $\mathbf{G}$ are given by:

$$\mathbf{F}_L = \int_{\Gamma(t)} \sigma \kappa_f \mathbf{n}_f \delta_f (\mathbf{x} - \mathbf{x}_f) \, ds, \qquad (9)$$

$$\mathbf{G} = \int_{\Gamma(t)} \delta_f (\mathbf{x} - \mathbf{x}_f) \, ds. \qquad (10)$$

In these formulae, $\mathbf{n}_f$ corresponds to the unit normal vector to the interface, and $ds$ is the length of the interface element; $\kappa_f$ is twice the mean interface curvature obtained from the Lagrangian interface. The geometric information corresponding to the unit normal $\mathbf{n}_f$ and the length of the interface element $ds$ in $\mathbf{G}$ are computed directly from the Lagrangian interface and then distributed on an Eulerian grid using the discrete delta function and following an immersed boundary approach (Peskin, 1977). A detailed computation of the force and construction of the function field $\mathbf{G}$ can be found in Shin (2007); Shin & Juric (2009a,b); Shin et al. (2017).

The Lagrangian interface is advected by integrating

$$\frac{d\mathbf{x}_f}{dt} = \mathbf{V}, \qquad (11)$$

with a second-order Runge-Kutta method where the interface velocity $\mathbf{V}$ is interpolated from the Eulerian velocity. Incorporating the complex geometry of the impeller and its rotation requires the implementation of the so-called Direct Forcing Method (Mohd-Yusof, 1997; Fadlun et al., 2000); this, in
turn, is implemented by incorporating a fluid-solid interaction force $\mathbf{F}_{\text{fsi}}$ into Eq. (2) defined numerically using the latest step of the temporal integration of Eq. (2):

$$ \rho \frac{\bar{\mathbf{u}}^{n+1} - \bar{\mathbf{u}}^n}{\Delta t} = \text{local}^n + \mathbf{F}_{\text{fsi}}^n. \qquad (12) $$

Here, “local” represents the right-hand-side terms of Eq. (2) that comprise the convective, pressure gradient, viscous, turbulent, gravitational, and surface tension force terms; the superscripts denote the discrete temporal step in the computation.

In the solid part of the domain corresponding to the impeller, the rotational motion $\mathbf{V}^{n+1}$ is enforced:

$$ \bar{\mathbf{u}}^{n+1} = \mathbf{V}^{n+1} = 2\pi f ((y - y_0), -(x - x_0)), \qquad (13) $$

where $(x_0, y_0) = (4.25, 4.25)$ cm are the coordinates of the impeller axis. Hence, $\mathbf{F}_{\text{fsi}}^n$ is given by

$$ \mathbf{F}_{\text{fsi}}^n = \rho \frac{\mathbf{V}^{n+1} - \bar{\mathbf{u}}^n}{\Delta t} - \text{local}^n. \qquad (14) $$

The no-slip condition is imposed on the velocity and the interface at the edge of the impeller parts.

The computational domain (see Figure 1-(a)) is a rectangular parallelepiped discretized by a uniform fixed three-dimensional finite-difference mesh and has a standard staggered MAC cell arrangement (Harlow & Welch, 1965). The velocity component $\bar{u}$, $\bar{v}$, and $\bar{w}$ are defined on the corresponding cell faces while the scalar variables (pressure $\bar{p}$ and the distance function $\psi$) are located at the cell centres. All spatial derivatives are approximated by second-order centred differences. The velocity field is solved by a parallel generalized minimal residual (GMRES) method (Saad, 2003) and the pressure field by a modified parallel 3D Vcycle multigrid solver based on the work of Kwak & Lee (2004) and described in Shin et al. (2017). Parallelisation is achieved using an algebraic domain decomposition where communication across processes is handled by Message Passing Interface (MPI) protocols.

The chosen pitched blade turbine shown in Figure 1-(a) and described in the beginning of this section is built using a combination of primitive geometric objects (planes, cylinders, and rectangular blocks) where each object is defined by a static distance function $\psi(x, y, z)$, positive in the fluid and negative in the solid. The resulting shape in Figure 1-(a) corresponds to the iso-value $\psi(x, y, z) = 0$. Details on how to construct similar complex objects are described in Kahouadji et al. (2018).

The temporal scheme is based on a second-order Gear method (Tucker, 2014) with implicit solution of the viscous terms of the velocity components. The time-step $\Delta t$ is chosen at each temporal iteration in order to satisfy a criterion based on:

$$ \Delta t = \min \{ \Delta t_{\text{cap}}, \Delta t_{\text{vis}}, \Delta t_{\text{CFL}}, \Delta t_{\text{int}} \} \qquad (15) $$

where $\Delta t_{\text{cap}}$, $\Delta t_{\text{vis}}$, $\Delta t_{\text{CFL}}$, and $\Delta t_{\text{int}}$ represent the capillary, viscous, Courant-Friedrichs-Lewy (CFL), and interfacial CFL time-steps, respectively, defined by:

$$ \Delta t_{\text{cap}} \equiv \frac{1}{2} \left( \frac{\rho_o + \rho_w}{\pi \sigma} \frac{\Delta x_{\text{min}}^3}{\Delta t} \right)^{1/2}, \quad \Delta t_{\text{vis}} \equiv \min \left( \frac{\rho_w}{\mu_w}, \frac{\rho_o}{\mu_o} \right) \frac{\Delta x_{\text{min}}^2}{6}, \qquad (16) $$

$$ \Delta t_{\text{CFL}} \equiv \min_j \left( \min_{\text{domain}} \left( \frac{\Delta x_j}{u_j} \right) \right), \quad \Delta t_{\text{int}} \equiv \min_j \left( \min_{\Gamma(t)} \left( \frac{\Delta x_j}{\|V\|} \right) \right), \qquad (17) $$

where $\Delta x_{\text{min}} = \min_j (\Delta x_j)$ refers to the minimum size $x$ at a given cell $j$, $u_j$, and $V$ are the maximum fluid and interface velocities, respectively.

Mass is conserved to within ±0.02‰ in all of our simulations. Our numerical procedure (Shin et al., 2018) has also been employed to simulate aeration in gas-liquid mixing (Kahouadji et al., 2022), vortex-interface dynamics and dispersion formation in turbulent water jets injected into an oil

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53
<figure><img src="image_3.png" /><figcaption>Figure 2: Spatio-temporal evolution of the vortical structures for $Re = 5663$ and $We = 22$ from an initial static state (a) until $t = 20 \times T$ (l). In this and subsequent figures, the results are displayed in the $(x, z)$ plane at $y = 4.25$, and cyan and yellow are used to designate the water and oil phases, respectively, while the red solid line represents the interface location. Enlarged snapshots of the vortical structures with corresponding annotations are shown for $t = 6 \times T$ and $t = 9.5 \times T$ in (m) and (n), respectively.</figcaption></figure>

medium (Constante-Amores et al., 2020b; Constante-Amores, 2021; Constante-Amores et al., 2021a), the breakup of ligaments (Constante-Amores et al., 2020a), the coalescence of drops with interfaces (Constante-Amores et al., 2021b), the bursting of bubbles through interfaces (Constante-Amores et al., 2021c), and the dynamics of falling films (Batchkarov et al., 2020). We have also carried out a validation study against the experimental work of Hançil & Rod (1988); the results are in the Supplementary Information section.

### 3. Results and Discussion

In this section, we present the spatio-temporal evolution of oil-water emulsion formation in the stirred vessel with a 4-pitch-blade turbine. First, a discussion of the vortical flow structures, and associated interfacial dynamics, is provided for the full range of $Re$ and $We$ studied starting with the flows generated for the lowest rotational speeds. Then, the results of an examination of the drop size distribution evolution and its dependence on $Re$ and $We$ are presented.

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52
## 3.1. Flow structures

The evolution of the flows for the two phases in the (x, z) plane at y = 4.25 shows the existence of rich dynamics as demonstrated by the example depicted in Figure 2 for the lowest Reynolds and Weber numbers studied, Re = 5663 and We = 22. A large primary vortex is created at early times (t = 0.5 × T, see Figure 2-(b)), which starts at the blade tips toward the bottom of the vessel and then develops upwards near the vessel wall toward the interface. A small secondary vortex is also generated near the bottom of the impeller, which rotates radially in the opposite direction compared to the primary vortex. A suite of vortical structures accompany the development of the primary vortex which correspond to Moffatt, blade-tip, Kelvin-Helmholtz, and wall-end vortices (see Figure 2-(c) – (e)). As the velocity is continuous across the interface, and the flow in the oil phase is dragged radially inward with the meridional flow in the water phase, a large vortex is also generated in the oil phase which has the opposite circulation to that in the water. At t = 6 × T, a swirling flow is observed near the interface in the water, which centrifuges the upper oil adjacent to it radially outward generating a small vortex of opposite circulation to the large one (see Figure 2-(f)). At t = 9.5 × T, a small vortex breakdown appears at the interface in the oil (see Figure 2-(g)), which grows squeezing the large meridional flow to the fixed sidewall until it gives way to two large circulation cells with a direction of radially outward (Figure 2-(h) – (k)); we will discuss vortex breakdown at higher Re and We below. The oil flow then remains in a configuration where the original meridional flow is confined to a thin layer in the vicinity of interface while the counter one above dominates in the oil phase (see Figure 2-(l)). On the other hand, flow separation (Figure 2-(g)) occurs in the water close to the edge of the rotating impeller (as described in the work of Piva & Meiburg (2005)). The water flow is characterised by a high level of

<figure><img src="image_4.png" /><figcaption>Figure 3: Spatio-temporal evolution of interface for Re = 22653 and We = 346 from a state where the interface starts to deform (a) until $t = 20 \times T$ (l).</figcaption></figure>

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52
turbulence involving the development of a multitude of vortices of varying scales. Similar behaviour in the water phase was demonstrated by Kahouadji et al. (2022) who carried out simulations of air-water mixing in stirred vessels.

We now present the spatio-temporal evolution of the interface shape for $Re = 22653$ and $We = 346$ in Figure 3. As the impeller rotates, the azimuthal velocity generated by the impeller exerts a torque on the interface deforming it and subsequently drags the oil in the centre downwards (see Figure 3-(a) – (b)). The impeller rotation is not sufficiently rapid so as to cause the interface to reach the impeller; instead, the interface deforms into a helical shape, as shown in Figure 3-(c). Afterwards, a long ligament detaches from the deforming interface, which subsequently breaks into a large drop and several smaller ones (Figure 3-(d) – (e)). When $t > 14 \times T$, the dispersed large drop retracts to the deforming interface spirally stretching in the meantime to become longer and thinner before attaching to the impeller (see Figure 3-(f) – (i)). The ligament detachment, breakup, and retraction cycle is repeated until the system reaches a dynamic steady state where small drops remain dispersed in the water phase together with the deforming interface in the shape of a Newton’s Bucket (see Figure 3-(j) – (l)). The existence of both ligament retraction, due to the action of interfacial tension, and elongation, due to the interaction with the turbulent flow in the water phase, highlights the presence of the delicate interplay amongst the various competing forces which results in the complexity associated with the mixing process.

An enlarged view of the breakup process is provided in Figure 4, which presents the concurrence of multiple drop pinch-offs along one single ligament. Close inspection reveals that the ligament initially takes time to elongate and thin giving rise to thread-like regions between nodules, which continue to thin until they break up into individual droplets via a Rayleigh-Plateau mechanism. To explain the phenomenon observed within our mixing system, we examine closely the vorticity profiles around the ligament. Figure 5 presents four horizontal planes that depict the flow coloured according to the magnitude of the vorticity. As shown in these plots, which also illustrate the position of the ligament relative to the vorticity distribution in the flow, the ligament lies in the high vorticity region generated by the impeller rotation. The interaction of the vorticity with the interface contributes to the instability of the ligament and this is particularly pronounced in the near-impeller regions in which the magnitude of the vorticity is significantly larger than those close to the vessel walls.

Next, we return to the occurrence of vortex breakdown in the oil phase which we had discussed above and note that it is also observed in the range of $Re = 5663 – 28316$ and $We = 22 – 540$ where the interface remains in a dynamic steady state assuming a Newton’s Bucket-like shape, as shown in Figure 6. The top panels of Figure 6 highlight the vortical structures when the flow has reached a dynamic steady state where the vortex breakdown has developed into the dominant circulation cells along with the counter-rotating ones confined to a thin layer near the interface. For the higher $Re$ and $We$ cases shown in Figure 6(a) – (e), the vortex breakdown appears earlier, for example, at around $t = 7 \times T$ for ($Re, We$) = (28316, 540) compared to $t = 9.5 \times T$ for ($Re, We$) = (5663, 22) as shown in Figure 2. The bottom panels of Figure 6 also show that with increased $Re$ and $We$ the depth of Newton’s Bucket increases approaching the impeller with the interface assuming a tear-drop shape. It may be possible to use the onset of vortex breakdown as an indicator of, and a surrogate for, the fact that the interface will remain sufficiently far from the impeller, and therefore unattached to it, at steady state. A methodology that focuses on the detection of vortex breakdown (if it occurs) for a given set of parameters, and then tracks the dependence of the breakdown over a range of such parameters, can lead to a map in ($Re, We$) space in which the region of limited mixing is demarcated; such an approach would obviate the need for lengthy computations to resolve the intricate and rapidly evolving interfacial dynamics.

We now consider a range of $Re$ and $We$ for which no retraction of the deforming interface from the impeller is observed. In Figure 7, we show the spatio-temporal evolution of the vortical structures for $Re = 39642$ and $We = 1059$ from an early-time state (at $t = 2 \times T$, see Figure 7-(a)) to a stage of the flow at which the deforming interface has reached the impeller (at $t = 9 \times T$, see Figure 7-(f)). During the early stages of the flow, the vortical structures observed are similar to those depicted in Figure 2 discussed above. Instead of the appearance of vortex breakdown seen at lower $Re$ and $We$, the streamline configuration at higher $Re$ and $We$ indicates a strong tendency of the interface to be deflected

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
<figure><img src="image_5.png" /><figcaption>Figure 4: Enlarged views of ligament breakup for $Re = 28316$ and $We = 540$ at $16.375T$ and $16.625T$.</figcaption></figure>

towards the impeller (see Figure 7-(e)). Although the water phase is dominated by the primary vortex, there are a number of secondary vortices close to the vessel sidewall in the vicinity of the interface which act to pump the water upward along the wall (see Figure 7-(c)-(f)). In Figure 8, we illustrate the spatio-temporal evolution of interface shape for ($Re, We$) = (33979, 778) and (39642,1059) shown for $t = 8 - 25 \times T$. After the interface reaches the impeller (see Figure 8-(b)), it is further stretched radially and long oil ligaments are created when the impeller blades cut through the rotating curtains (see Figure 8-(c)). Long ligaments subsequently stretch to become thinner until pinch-off occurs giving rise to multitudes of smaller drops (see Figure 8-(d)). Then, the deforming interface evolves differently depending on the value of $Re$ and $We$. Whereas for $Re = 28316$ and $We = 540$, the interface eventually detaches from the impeller and no further ligaments are generated, as shown in Figure 4, the formation of new ligaments remains in place at the higher $Re$ and $We$ (see Figure 8-(e)). The small drops created from ligament breakup either travel upwards under the action of buoyancy, eventually merging with the oil phase, or remain dispersed in water until they coalesce with other drops. As a result, although the number of dispersed drops decreases with time for $Re = 28316$ and $We = 540$, this is not the case at higher $Re$ and $We$ as the drops ‘lost’ due to coalescence with the oil-water interface are replenished via

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52
continual ligament formation and subsequent breakup (see Figure 8-(f)). A discussion of the evolution of the dispersed drop number and size distribution will be provided in subsection 3.2.

At larger *Re* and *We*, the vortical structures in the oil and water phases during the early and intermediate stages of the dynamics resemble those discussed above; this is exemplified by Figure 9-(a)–(f) in which we show the evolution of these structures for *Re* = 56631 and *We* = 2160 between *t* = 2 × *T* and 9 × *T* which culminate in the interface reaching the impeller. Inspection of Figure 9-(g)–(l), however, in which we illustrate the evolving interfacial dynamics, reveals that drop coalescence dominates over breakup. In contrast to the lower *Re* and *We* cases discussed above, a larger torque is generated in the present case which results in a deforming interface of a helical shape with four rotating curtains (see Figure 9-(g)). Following their formation, the long ligaments merge to form a relatively large mass of oil that collects around the bottom of the vessel before breaking up into smaller drops as they are rapidly and continually created from the impeller tips (see Figure 9-(h)). As can be seen in Figure 9-(i), the dispersed phase in this high *Re* case is composed mainly of stretched oil ligaments rather than the individual drops observed to accompany the flow at lower *Re*. With increasing time, an increasing amount of oil is entrained into the water phase; as a result, the latter is displaced upwards along the vessel sidewalls lifting the interface as can be seen in Figure 9-(j) – (k). At *t* = 25 × *T*, the dynamics are characterised by complex interfacial topology. It is notable that the number of small droplets has decreased due to coalescence giving way to large dispersed phase structures resulting from the merging of ligaments and drops, as shown in Figure 9-(l). The drop size distribution as a function of the impeller speed is discussed next.

## 3.2. Drop size distribution

From an industrial point of view, it is often important to track the evolution of the dispersed phase statistics and metrics in terms of number of drops, their size distribution, the interfacial area, and hold up. Figure 10 presents the temporal evolution of the number of dispersed drops ($t = 5 - 25 \times T$) for

<figure><img src="image_6.png" /></figure>

Figure 5: Vorticity in the ligament surroundings for *Re* = 28316 and *We* = 540. Four horizontal planes generated for *z* = 0.82, 0.70, 0.65, 0.58, corresponding to the red-, blue-, green- and yellow-labelled planes, are coloured according to the vorticity magnitude. The side- and top-views of the ligament shown in the middle and parts of the ligament are highlighted by white solid lines in the vorticity panels.
<figure><img src="image_7.png" /><figcaption>Figure 6: Vortical flow structures, top panels, and their corresponding interfacial shape, bottom panels, for (Re, We) combinations of (5663,22), (11326,86), (16989,194), (22653,346), and (28316,540) shown in (a) – (e), respectively, at $t = 25 \times T$ when the flow has reached a dynamic steady state.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Spatio-temporal evolution of the vortical structures for Re = 39642 and We = 1059 at $t = 2 - 9 \times T$ shown in (a) – (f), respectively.</figcaption></figure>

$f = 5 - 10$ Hz (Re = 28316 – 56632 and We = 540 – 2160), along with the corresponding interface shapes at $t = 25 \times T$. For all displayed frequencies, no dispersed drop appears until $t \approx 8 \times T$, and the appearance starts earlier with increasing frequency. Figure 10 shows that there is a steep increase in drop number from zero, caused by the above-mentioned Rayleigh-Plateau-type breakup from long

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
<figure><img src="image_9.png" /><figcaption>Figure 8: Spatio-temporal evolution of the interface toward the development of an oil-in-water emulsion for (Re, We) = (33979, 778) and (39642, 1059) shown in the top and bottom panels, respectively, for $t = 8 \times T$ between $25 \times T$ in (a) – (f), respectively.</figcaption></figure>

In addition, Figure 8-(c) provides clear evidence that when rotation frequency increases from $f = 5$ Hz to $f = 7$ Hz, a larger part of the deforming interface contacts the impeller and therefore more ligaments are produced at the beginning which corresponds to the steeper increase, as shown in Figure 10 for $f = 5 - 7$ Hz. After that, the number of dispersed oil drops for $f = 5$ Hz falls off gently since the deforming interface no longer makes contact with the impeller while dispersed drops tend to coalesce with the oil overlying the water phase. With increasing frequency, the dispersed oil drop number continues to increase until the breakup rate balances that of coalescence. At $f = 8$ Hz, the number increases rapidly reaching a peak before decreasing towards a dynamic steady state.

It is evident from Figure 10 that the steady number of dispersed drops has a non-monotonic dependence on the impeller speed while its peak increases with increasing rotation frequencies. On the other hand, our results also indicate that the interfacial area of the dispersed phase at steady state ($t = 25 \times T$) shows no such dependence as shown in ??-(a). This interfacial area is calculated as $A_d/A_{cap}$, where $A_d$ refers to the surface area of all dispersed drops and $A_{cap}$ is the sphere surface area with the capillary length as its diameter. The interfacial area appears to be proportional to the increased impeller speed (up to $f = 10$ Hz) even though the number of dispersed drops reaches its peak at $f = 7$ Hz. The peak in the number of drops at this intermediate frequency is due to the fact that, at larger frequencies, long ligaments, with larger interfacial area, dominate the dispersed phase morphology (see Figure 9). From Figure 11, the optimal impeller speed can be selected depending on the relative importance of the interfacial area and dispersed drop number to the particular application in question. For all frequencies, the drop size distribution up to $t = 25 \times T$ can be found in the Supplementary Information, where the smallest droplet size captured is $r_{min} \approx 0.0325$ mm.
<figure><img src="image_10.png" /><figcaption>Figure 9: Spatio-temporal evolution of the vortical structure between $t = 2 \times T$ and $9 \times T$, $(a) - (f)$, and of the interface between $t = 9 \times T$ and $25 \times T$, $(g) - (l)$, for $Re = 56632$ and $We = 2160$.</figcaption></figure>

4. Conclusion and Perspective

In this paper, using our high-fidelity three-dimensional simulations, we study the effect of impeller rotational speed on the complex interfacial dynamics during oil-water mixing in a stirred vessel with a 4-pitch-blade turbine; we explore rotational speeds in the range $f = 1 - 10$ Hz, which correspond to Reynolds and Weber number ranges of $Re = 5663 - 56632$ and $We = 22 - 2160$, respectively. We have shown that the interface shape and vortical structures within such a system are strongly dependent on the impeller speed. In particular, we have demonstrated the possibility of using vortex breakdown as an indicator that the interface undergoes little breakup due to very little interaction with the impeller. In addition, we have also presented a drop size analysis over the full range of frequencies studied in terms of temporal evolution of the dispersed drop number and its size distribution at steady state, and shown that an optimal impeller speed can be selected depending on the needs of the application in question. In future research, it can be interesting to explore other effects on the mixing system, which include impeller geometry, fluid phase rheology, and the presence of surfactants; the latter are of particular interest since surfactants may be present in the system either as a contaminant or by design.

**Funding Statement.** This work is supported by the Engineering & Physical Sciences Research Council, UK, through the MEMPHIS (EP/K003976/1) and PREMIERE (EP/T000414/1) Programme Grants, and by computing time at HPC facilities provided by the Research Computing Service of Imperial College London. DJ and JC acknowledge support through computing time at the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Equipement National de Calcul Intensif) Grant 2021 A0102B06721. Simulations were performed using code BLUE (Shin et al., 2017) and the visualisations were generated using Paraview.

**Declaration of Interests.** The authors declare no conflict of interest.

01
02
03
04
05
06
07
08
09
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
<figure><img src="image_11.png" /></figure>

Figure 10: Temporal evolution of the number of dispersed drops for frequencies $f = 5, 6, 7, 8, 9, 10$ Hz, which correspond to $Re = 28316 – 56632$ and $We = 540 – 2160$. Snapshots of the interface shapes at $t = 25 \times T$ are also presented for each frequency.

<figure><img src="image_12.png" /></figure>

Figure 11: Number of drops (blue line, inverted triangle markers) and interfacial area (red line, cross-shaped markers) of the dispersed oil phase.

**Author Contributions.** O.K.M. secured the funding; F.L., L.K., J.P.V., D. J. and O.K.M. contributed to the design of the research plan and manuscript editing; F.L., L.K., S.S., J.C., and D.J. contributed to the computations.

**Data Availability Statement.** Raw data are available from the corresponding author L. K.

**Ethical Standards.** The research meets all ethical guidelines, including adherence to the legal requirements of the study country.
References

AFSHAR GHOTLI, REZA, RAMAN, ABDUL AA, IBRAHIM, SHALIZA & BAROUTIAN, SAEID 2013 Liquid-liquid mixing in stirred vessels: a review. *Chemical Engineering Communications* **200**(5), 595-627.

ALBERNAZ, DANIEL L., DO-QUANG, MINH, HERMANSON, JIM C & AMBERG, GUSTAV 2017 Droplet deformation and heat transfer in isotropic turbulence. *Journal of Fluid Mechanics* **820**, 61–85.

ANDERSON, DANIEL M, MCFADDEN, GEOFFREY B & WHEELER, ADAM A 1998 Diffuse-interface methods in fluid mechanics. *Annual review of fluid mechanics* **30**(1), 139–165.

ANDERSSON, RONNIE & ANDERSSON, BENGT 2006 On the breakup of fluid particles in turbulent flows. *AIChE Journal* **52**(6), 2020–2030.

BAKKER, ANDRÉ, LAROCHE, RICHARD D, WANG, MIN-HUA & CALABRESE, RICHARD V 1997 Sliding mesh simulation of laminar flow in stirred reactors. *Chemical Engineering Research and Design* **75**(1), 42–44.

BATCHVAROV, A., KAHOUADJI, L., MAGNINI, M., CONSTANTE-AMORES, C. R., CRASTER, R. V., SHIN, S., CHERGUI, J., JURIC, D. & MATAR, O. K. 2020 Effect of surfactant on elongated bubbles in capillary tubes at high reynolds number. *Phys. Rev. Fluids* **5**, 093605.

BECKER, PER JULIAN, PUEL, FRANÇOIS, CHEVALIER, YVES & SHEIBAT-OTHMAN, NIDA 2014 Monitoring silicone oil droplets during emulsification in stirred vessel: Effect of dispersed phase concentration and viscosity. *The Canadian Journal of Chemical Engineering* **92**(2), 296–306.

BOXALL, JOHN A, KOH, CAROLYN A, SLOAN, E DENDY, SUM, AMADEU K & WU, DAVID T 2010 Measurement and calibration of droplet size distributions in water-in-oil emulsions by particle video microscope and a focused beam reflectance method. *Industrial & Engineering Chemistry Research* **49**(3), 1412–1418.

CALABRESE, RICHARD V, CHANG, TPK & DANG, PT 1986 Drop breakup in turbulent stirred-tank contactors. part i: Effect of dispersed-phase viscosity. *AIChE Journal* **32**(4), 657–666.

CHIU, PAO-HSIUNG & LIN, YAN-TING 2011 A conservative phase field method for solving incompressible two-phase flows. *Journal of Computational Physics* **230**(1), 185–204.

CONSTANTE-AMORES, C. R., KAHOUADJI, L., BATCHVAROV, A., SHIN, S., CHERGUI, J., JURIC, D. & MATAR, O. K. 2021a Direct numerical simulations of transient turbulent jets: vortex-interface interactions. *J. Fluid Mech.* **922**, A6.

CONSTANTE-AMORES, C. R. 2021 Three-dimensional computational fluid dynamics simulations of complex multiphase flows with surfactants. *Imperial College London PhD Thesis*.

CONSTANTE-AMORES, C. R., KAHOUADJI, L., BATCHVAROV, A., S., SEUNGWON, CHERGUI, J., JURIC, D. & MATAR, O. K. 2020a Dynamics of retracting surfactant-laden ligaments at intermediate ohnesorge number. *Phys. Rev. Fluids* **5**, 084007.

CONSTANTE-AMORES, C. R., KAHOUADJI, L., BATCHVAROV, A., S., SEUNGWON, CHERGUI, J., JURIC, D. & MATAR, O. K. 2021b Role of surfactant-induced marangoni stresses in drop-interface coalescence. *J. Fluid Mech.* **925**, A15.

CONSTANTE-AMORES, C. R., KAHOUADJI, L., BATCHVAROV, A., SHIN, S., CHERGUI, J., JURIC, D. & MATAR, O. K. 2020b Rico and the jets: Direct numerical simulations of turbulent liquid jets. *Phys. Rev. Fluids* **5**, 110501.

CONSTANTE-AMORES, C. R., KAHOUADJI, L., BATCHVAROV, A., SHIN, S., CHERGUI, J., JURIC, D. & MATAR, O. K. 2021c Dynamics of a surfactant-laden bubble bursting through an interface. *J. Fluid Mech.* **911**, A57.

DE HERT, SERGIO CARRILLO & RODGERS, THOMAS L 2017 On the effect of dispersed phase viscosity and mean residence time on the droplet size distribution for high-shear mixers. *Chemical Engineering Science* **172**, 423–433.

DEARDORFF, JAMES W 1970 A numerical study of three-dimensional turbulent channel flow at large reynolds numbers. *Journal of Fluid Mechanics* **41**(2), 453–480.

DERKSEN, JJ 2003 Numerical simulation of solids suspension in a stirred tank. *AIChE Journal* **49**(11), 2700–2714.

DERKSEN, JJ 2011 Blending of miscible liquids with different densities starting from a stratified state. *Computers & fluids* **50**(1), 35–45.

DERKSEN, JJ 2012 Direct simulations of mixing of liquids with density and viscosity differences. *Industrial & engineering chemistry research* **51**(19), 6948–6957.

DERKSEN, JOS & VAN DEN AKKER, HARRY EA 1999 Large eddy simulations on the flow driven by a rushton turbine. *AIChE Journal* **45**(2), 209–221.

EGGELS, JACK GM 1996 Direct and large-eddy simulation of turbulent fluid flow using the lattice-boltzmann scheme. *International journal of heat and fluid flow* **17**(3), 307–323.

EL-HAMOUZ, AMER, COOKE, MIKE, KOWALSKI, ADAM & SHARRATT, PAUL 2009 Dispersion of silicone oil in water surfactant solution: Effect of impeller speed, oil viscosity and addition point on drop size distribution. *Chemical Engineering and Processing* **48**(2), 633–642.

GIAPOS, A, PACHATOURIDIS, C & STAMATOUDIS, M 2005 Effect of the number of impeller blades on the drop sizes in agitated dispersions. *Chemical Engineering Research and Design* **83**(12), 1425–1430.

GILLISSEN, JJJ & VAN DEN AKKER, HEA 2012 Direct numerical simulation of the turbulent flow in a baffled tank driven by a rushton turbine. *AIChE Journal* **58**(12), 3878–3890.

HANČIL, VLADISLAV & ROD, VLADIMÍR 1988 Break-up of a drop in a stirred tank. *Chemical Engineering and Processing: Process Intensification* **23**(3), 189–193.

HARLOW, FRANCIS H & WELCH, J EDDIE 1965 Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface. *The physics of fluids* **8**(12), 2182–2189.

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52
HARVEY, PS, PS, HARVEY & OTHERS 1982 Turbulent flow in an agitated vessel. ii: Numerical solution and model predictions.

HASAN, BASIM O 2017 Breakage of drops and bubbles in a stirred tank: A review of experimental studies. *Chinese journal of chemical engineering* **25** (6), 698–711.

HINZE, JULIUS O 1955 Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes. *AIChE journal* **1** (3), 289–295.

HIRT, CYRIL W & NICHOLS, BILLY D 1981 Volume of fluid (vof) method for the dynamics of free boundaries. *Journal of computational physics* **39** (1), 201–225.

IBARRA, R. 2017 Horizontal and low-inclination oil-water flow investigation using laser-based diagnostic techniques. *Imperial College London PhD Thesis*.

KAHOUADJI, LYES, LIANG, FUYUE, VALDES, JUAN, CHERGUI, JALEL, JURIC, DAMIR, SHIN, SEUNGWON, CRASTER, RICHARD V & MATAR, OMAR K 2022 The transition to aeration in turbulent two-phase mixing in stirred vessels. *Flow Submitted*.

KAHOUADJI, LYES, NOWAK, EMILIA, KOVALCHUK, NINA, CHERGUI, JALEL, JURIC, DAMIR, SHIN, SEUNGWON, SIMMONS, MARK JH, CRASTER, RICHARD V & MATAR, OMAR K 2018 Simulation of immiscible liquid–liquid flows in complex microchannel geometries using a front-tracking scheme. *Microfluidics and nanofluidics* **22** (11), 1–12.

KOMRAKOVA, ALEXANDRA E 2019 Single drop breakup in turbulent flow. *The Canadian Journal of Chemical Engineering* **97** (10), 2727–2739.

KONNO, MIKIO, AOKI, MINORU & SAITO, SHOZABURO 1983 Scale effect on breakup process in liquid–liquid agitated tanks. *Journal of chemical engineering of Japan* **16** (4), 312–319.

KRAUME, MATTHIAS, GÄBLER, ANSOR & SCHULZE, KAI 2004 Influence of physical properties on drop size distribution of stirred liquid–liquid dispersions. *Chemical Engineering & Technology: Industrial Chemistry-Plant Equipment-Process Engineering-Biotechnology* **27** (3), 330–334.

KWAK, DO Y & LEE, JUN S 2004 Multigrid algorithm for the cell-centered finite difference method ii: discontinuous coefficient case. *Numerical Methods for Partial Differential Equations: An International Journal* **20** (5), 742–764.

LILLY, DOUGLAS KEITH 1966 On the application of the eddy viscosity concept in the inertial sub-range of turbulence. *NCAR manuscript* **123**.

LILLY, DOUGLAS K 1967 The representation of small-scale turbulence in numerical simulation experiments. *IBM Form pp. 195–210*.

LOVICK, J, MOUZA, AA, PARAS, SV, LYE, GJ & ANGELI, P 2005 Drop size distribution in highly concentrated liquid–liquid dispersions using a light back scattering method. *Journal of Chemical Technology & Biotechnology: International Research in Process, Environmental & Clean Technology* **80** (5), 545–552.

LUO, JY & GOSMAN, AD 1994 Prediction of impeller-induced flow in mixing vessels using multiple frames of reference. *INSTITUTE OF CHEMICAL ENGINEERS SYMPOSIUM SERIES*.

MAASS, S, GÄBLER, A, ZACCONE, A, PASCHEDAG, AR & KRAUME, M 2007 Experimental investigations and modelling of breakage phenomena in stirred liquid/liquid systems. *Chemical Engineering Research and Design* **85** (5), 703–709.

MAASS, SEBASTIAN & KRAUME, MATTHIAS 2012 Determination of breakage rates using single drop experiments. *Chemical Engineering Science* **70**, 146–164.

MAASS, SEBASTIAN, METZ, FLORIAN, REHM, TORSTEN & KRAUME, MATTHIAS 2010 Prediction of drop sizes for liquid–liquid systems in stirred slim reactors—part i: Single stage impellers. *Chemical Engineering Journal* **162** (2), 792–801.

MAASS, S, WOLLNY, S, SPERLING, R & KRAUME, M 2009 Numerical and experimental analysis of particle strain and breakage in turbulent dispersions. *Chemical Engineering Research and Design* **87** (4), 565–572.

MCMILLAN, ODEN J & FERZIGER, JOEL H 1979 Direct testing of subgrid-scale models. *Aiaa Journal* **17** (12), 1340–1346.

MEYERS, JOHAN & SAGAUT, PIERRE 2006 On the model coefficients for the standard and the variational multi-scale smagorinsky model. *Journal of Fluid Mechanics* **569**, 287–319.

NACHTIGALL, STEPHANIE, ZEDEL, DANIEL & KRAUME, MATTHIAS 2016 Analysis of drop deformation dynamics in turbulent flow. *Chinese journal of chemical engineering* **24** (2), 264–277.

NACHTIGALL, S, ZEDEL, D, MAASS, S, WALLE, A, SCHÄFER, M & KRAUME, M 2012 Determination of drop breakage mechanisms by experimental and numerical investigations of single drop breakages. In *14th European Conference on Mixing, Poland*.

NAEENI, SEPEHR KHAJEH & PAKZAD, LEILA 2019a Droplet size distribution and mixing hydrodynamics in a liquid–liquid stirred tank by cfd modeling. *International Journal of Multiphase Flow* **120**, 103100.

NAEENI, SEPEHR KHAJEH & PAKZAD, LEILA 2019b Experimental and numerical investigation on mixing of dilute oil in water dispersions in a stirred tank. *Chemical Engineering Research and Design* **147**, 493–509.

PACEK, AW, MAN, CC & NIENOW, AW 1998 On the sauter mean diameter and size distributions in turbulent liquid/liquid dispersions in a stirred vessel. *Chemical Engineering Science* **53** (11), 2005–2011.

PESKIN, CHARLES S 1977 Numerical analysis of blood flow in the heart. *Journal of computational physics* **25** (3), 220–252.

PIVA, M & MEIBURG, E 2005 Steady axisymmetric flow in an open cylindrical container with a partially rotating bottom wall. *Physics of Fluids* **17** (6), 063603.

POPE, STEPHEN B 2004 Ten questions concerning the large-eddy simulation of turbulent flows. *New journal of Physics* **6** (1), 35.

POPE, STEPHEN B & POPE, STEPHEN B 2000 *Turbulent flows*. Cambridge university press.

QIN, CHENGPENG, CHEN, CHAO, XIAO, QI, YANG, NING, YUAN, CANSHENG, KUNKELMANN, CHRISTIAN, CETINKAYA, MURAT & MÜLHEIMS, KERSTIN 2016 Cfd-pbm simulation of droplets size distribution in rotor-stator mixing devices. *Chemical Engineering Science* **155**, 16–26.

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53
RAMKRISHNA, DORAISWAMI 2000 *Population balances: Theory and applications to particulate systems in engineering*. Elsevier.

ROCCON, ALESSIO, DE PAOLI, MARCO, ZONTA, FRANCESCO & SOLDATI, ALFREDO 2017 Viscosity-modulated breakup and coalescence of large drops in bounded turbulence. *Physical Review Fluids* **2** (8), 083603.

ROUDSARI, SHIDEH FATHI, TURCOTTE, GINETTE, DHIB, RAMDHANE & EIN-MOZAFFARI, FARHAD 2012 Cfd modeling of the mixing of water in oil emulsions. *Computers & chemical engineering* **45**, 124–136.

SAAD, YOUSEF 2003 *Iterative methods for sparse linear systems*. SIAM.

SANJUAN-GALINDO, RENÉ, SOTO, ENRIQUE, ZENIT, ROBERTO & ASCANIO, GABRIEL 2015 Viscous filament fragmentation in a turbulent flow inside a stirred tank. *Chemical Engineering Communications* **202** (9), 1251–1260.

SATHYAGAL, AN, RAMKRISHNA, D & NARIMHAN, G 1996 Droplet breakage in stirred dispersions. *Chemical Engineering Science* **51** (9), 1377–1391.

SCARBOLO, LUCA, BIANCO, FEDERICO & SOLDATI, ALFREDO 2015 Coalescence and breakup of large droplets in turbulent channel flow. *Physics of Fluids* **27** (7), 073302.

SCARBOLO, LUCA, BIANCO, FEDERICO & SOLDATI, ALFREDO 2016 Turbulence modification by dispersion of large deformable droplets. *European Journal of Mechanics-B/Fluids* **55**, 294–299.

SCARBOLO, LUCA & SOLDATI, ALFREDO 2013 Turbulence modulation across the interface of a large deformable drop. *Journal of Turbulence* **14** (11), 27–43.

SETHIAN, JAMES A & SMEREKA, PETER 2003 Level set methods for fluid interfaces. *Annual review of fluid mechanics* **35** (1), 341–372.

SHIN, SEUNGWON 2007 Computation of the curvature field in numerical simulation of multiphase flow. *Journal of Computational Physics* **222** (2), 872–878.

SHIN, SEUNGWON, CHERGUI, JALEL & JURIC, DAMIR 2017 A solver for massively parallel direct numerical simulation of three-dimensional multiphase flows. *Journal of Mechanical Science and Technology* **31** (4), 1739–1751.

SHIN, SEUNGWON, CHERGUI, JALEL, JURIC, DAMIR, KAHOUADI, LYES, MATAR, OMAR K & CRASTER, RICHARD V 2018 A hybrid interface tracking–level set technique for multiphase flow with soluble surfactant. *Journal of Computational Physics* **359**, 409–435.

SHIN, SEUNGWON & JURIC, DAMIR 2009a A hybrid interface method for three-dimensional multiphase flows based on front tracking and level set techniques. *International Journal for Numerical Methods in Fluids* **60** (7), 753–778.

SHIN, SEUNGWON & JURIC, DAMIR 2009b Simulation of droplet impact on a solid surface using the level contour reconstruction method. *Journal of mechanical science and technology* **23** (9), 2434–2443.

SHU, SHULI & YANG, NING 2018 Gpu-accelerated large eddy simulation of stirred tanks. *Chemical Engineering Science* **181**, 132–145.

SMAGORINSKY, JOSEPH 1963 General circulation experiments with the primitive equations: I. the basic experiment. *Monthly weather review* **91** (3), 99–164.

SPROW, FB 1967 Distribution of drop sizes produced in turbulent liquid—liquid dispersion. *Chemical Engineering Science* **22** (3), 435–442.

STAMATOUDIS, MICHAEL & TAYLARIDES, LAWRENCE L 1985 Effect of continuous-phase viscosity on the drop sizes of liquid-liquid dispersions in agitated vessels. *Industrial & Engineering Chemistry Process Design and Development* **24** (4), 1175–1181.

TRYGGVASON, GRÉTAR, BUNNER, BERNARD, ESMAEELI, ASGHAR, JURIC, DAMIR, AL-RAWAHI, N, TAUBER, W, HAN, J, NAS, S & JAN, Y-J 2001 A front-tracking method for the computations of multiphase flow. *Journal of computational physics* **169** (2), 708–759.

TYAGI, MAYANK, ROY, SOMNATH, HARVEY III, ALBERT D & ACHARYA, SUMANTA 2007 Simulation of laminar and turbulent impeller stirred tanks using immersed boundary method and large eddy simulation technique in multi-block curvilinear geometries. *Chemical engineering science* **62** (5), 1351–1363.

UNVERDI, SALIH OZEN & TRYGGVASON, GRÉTAR 1992 A front-tracking method for viscous, incompressible, multi-fluid flows. *Journal of computational physics* **100** (1), 25–37.

WANG, CY & CALABRESE, RICHARD V 1986 Drop breakup in turbulent stirred-tank contactors. part ii: Relative influence of viscosity and interfacial tension. *AIChE journal* **32** (4), 667–676.

ZHOU, GENWEN & KRESTA, SUZANNE M 1998a Correlation of mean drop size and minimum drop size with the turbulence energy dissipation and the flow in an agitated tank. *Chemical Engineering Science* **53** (11), 2063–2079.

ZHOU, GENWEN & KRESTA, SUZANNE M 1998b Evolution of drop size distribution in liquid–liquid dispersions for various impellers. *Chemical Engineering Science* **53** (11), 2099–2113.

01

02

03

04

05

06

07

08

09

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52