Droplet size distribution and mixing hydrodynamics in a liquid–liquid stirred tank by CFD modeling

International Journal of Multiphase Flow 120 (2019) 103100

Contents lists available at ScienceDirect

International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmulflow

Droplet size distribution and mixing hydrodynamics in a liquid–liquid stirred tank by CFD modeling Sepehr Khajeh Naeeni, Leila Pakzad∗ Department of Chemical Engineering, Lakehead University, 955 Oliver Road, Thunder Bay, ON P7B 5E1, Canada

a r t i c l e

i n f o

Article history: Received 28 November 2018 Revised 12 August 2019 Accepted 28 August 2019 Available online 4 September 2019 Keywords: Emulsions Computational fluid dynamics Population balance modeling Droplet size distribution Rushton impeller

a b s t r a c t Stirred tanks are the common unit operations employed for liquid–liquid dispersion in many industries. In this study, computational fluid dynamic (CFD) coupled with population balance modeling (PBM) were applied to simulate mixing of water in crude oils in a stirred tank equipped with Rushton turbine. An Eulerian multiphase model and standard k-ε turbulence model were employed to simulate the flow field inside the tank. Experimental results were then used to validate simulation results. The effect of impeller speed, volume fraction of dispersed phase, and oil viscosity on droplet size distribution were investigated in this work. An increase in agitation speed was found to decrease the mean and Sauter mean diameter while increasing the homogeneity of the system. A wider droplet size distribution profile was achieved at higher oil volume fractions with almost no significant change in droplet size. Additionally, increasing the viscosity of the oil phase resulted in a gradual shift towards smaller droplets and narrower droplet size distribution. Liquid–liquid dispersions were occurred in two different turbulent regimes depending on the viscosity of continuous phase. Related correlations for each regime were obtained and the values of maximum droplet size were fitted with the dimensionless numbers. © 2019 Published by Elsevier Ltd.

1. Introduction Liquid–liquid dispersion plays a major role in many industrial applications whether chemical, petroleum, biopharma, food, or paint. The design and scale-up of mixing tanks employed for agitation of two immiscible liquids depend on comprehensive knowledge on droplet size distribution (DSD) as well as hydrodynamic characteristics of the system. DSD has a significant effect on the hydrodynamic and the stability of a liquid–liquid dispersion; additionally, this parameter affects the interphase mass and heat transfer rates through interfacial areas (Sjoblom, 2001). Thus, a stable dispersion with the desired DSD in a large-scale industrial tank is always demanding though challenging. Several experimental studies have been conducted on DSD in liquid–liquid dispersions (Wang and Calabrese, 1986; Ohtake et al., 1987 and 1988; Angle and Hamza, 2006; El-Hamouz et al., 2009; Khajeh Naeeni and Pakzad, 2019; Mirshekari et al., 2019). The analysis of numerical correlations for droplet size began with various studies such as Sprow (1967), Brown and Pitt (1972), Calabrese et al. (1986), and Lemenand et al. (2003). Zhou and Kresta (1998) proposed empirical correlations between Sauter mean diameter, d32 , and maximum stable droplet size, dmax , in ∗

Corresponding author. E-mail address: [email protected] (L. Pakzad).

https://doi.org/10.1016/j.ijmultiphaseflow.2019.103100 0301-9322/© 2019 Published by Elsevier Ltd.

liquid–liquid dispersion in a stirred tank equipped with a Rushton impeller. They found that the distribution shape changed from unimodal to bimodal distribution with increasing the impeller speed. According to their study, the proportional relation between d32 and dmax is nonlinear. However, earlier study by Sprow (1967) shows a linear relation between d32 and dmax in the case where breakage is dominant in the system. Pacek et al. (1999) concluded that the shape of the distribution profile depends on the continuous phase; being unimodal in aqueous continuous phase and a bimodal in organic continuous ones. They observed that a bimodal distribution profile at the low impeller speeds changed to a unimodal distribution as the impeller speed increased in the aqueous continuous phase. Lovick et al. (2005) used the optical reflectance measurement (3D ORM) particle size analyser and Endoscope to obtain on-line DSD in Kerosene in water dispersion with volume fraction of 10%–60%. They concluded that the effect of dispersed phase volume fraction on droplet size is insignificant in highly concentrated dispersions. Boxall et al. (2010) focused their study on the effect of continuous phase viscosity on DSD. An empirical correlation for the mean size of droplets obtained by focused beam reflectance measurement (FBRM) and particle video microscope (PVM) was reported in their study. Similarly, Vladisavljevic et al. (2011) presented the droplet size distribution as a function of oil viscosity in an oil-in-water micro-channel.

2

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

Nomenclature Abbreviation CFD computational fluid dynamic CLD chord length distribution DSD droplet size distribution FBRM focused beam reflectance measurement MRF multiple reference frame ORM optical reflectance measurement PBM population balance modeling PBT pitched blade impeller PSD particle size distribution PVM particle video microscope RSM Reynolds Stress Model Symbols Ai B C1ε C2ε C3ε CD d dmax,viscous dmax,D D D F f g Gb Gk I Kpq k N n P P Re −→ R pq S t u u u¯  u u’ U˜ v w We X

interfacial area concentration, m2 birth rate in Eq. (14), m3 s−1 k-ε turbulence model constant, k-ε turbulence model constant, k-ε turbulence model constant, drag coefficient, droplet diameter, m maximum diameter of stable drop in turbulent viscous regime according to Eq. (31), m Maximum diameter of stable drop in turbulent inertial regime according to Eq. (30), m death rate in Eq. (13), m3 s−1 impeller diameter, m external force in continuity equation, N drag force in Eq. (7), N gravitational acceleration, ms−2 generation rate of turbulent kinetic energy due to Buoyancy, kgm−1 s−3 generation rate of turbulent kinetic energy, kgm−1 s−3 unit tensor, Nm−2 interface momentum exchange coefficient, s−1 turbulence kinetic energy, m2 s−2 impeller speed, rps number density function of particles with volume V, m−3 pressure, Nm−2 probability of collision in Eq. (24), Reynolds number, interfacial force between phases p and q, ms−2 general variable, time, s x-component of velocity, ms−1 particle velocity, ms−1 characteristic velocity of collision, ms−1 velocity vector, ms−1 fluctuating velocity, ms−1 average velocity, ms−1 y-component of velocity, ms−1 z-component of velocity, ms−1 weber number, bin fraction size, -

Greek symbols ∇ vector differential operator ρ density, kgm−3 ɛ turbulence dissipation rate, m2 s−3 τp particle relaxation time, s

τ τi τp α σ k, ɛ σ σ υ υ B (υi , υk ) ag

θ ij λ λ β ω μ μ ∅

stress tensor, Nm−2 static pressure on particle or droplet, Nm−2 particle surface stress, Nm−2 phase volume fraction, turbulence model constants, surface tension, Nm−1 standard deviation, volume of particle, m3 variance, dimensionless breakup rate of drop with diameter di in to daughter drop size diameter dk , m3 s−1 (υi , υ j ) coalescence rate of a drop with volumes of vi , and vj , m3 s−1 collision rate function, m3 s−1 dimensionless eddy size, bulk viscosity in Eq. (5), Pa.s daughter droplet size distribution, Kg−1 collision frequency, m−3 median of distribution, viscosity, Pa.s scalar variable

Subscript B B C c D g I L m max p p pq q t i j k i&j v ag

μ

breakage Buoyancy in Eqs. (3) and (6) coalescence continuous drag gravitational index number lift mixture maximum dispersed phase particle in Eq. (6) interphase continuous phase turbulent counting variable counting variable counting variable directions in Eq. (11) virtual mass aggregation viscous

Carrillo de hert & Rodgers (2017) observed a transition from monomodal to bimodal DSD as the viscosity of the dispersed phase increased. Since the experimental study on liquid–liquid dispersion system is usually costly and time-consuming (Boxall et al., 2010; Rueger and Calabrese, 2013), it has been supported by numerical ones through using the population balance models (PBM) combined with computational fluid dynamics (CFD) (Ramkrishna, 1985). One of the most challenging parts involved in PBM is the deviation of calculated droplet sizes from experimental results which stems from inherent dependency of PBM on the empirical correlations and measurements. Among several available methods for modeling of particle size distribution (PSD) in gas and liquid dispersions, discrete method has shown promising results. The discrete method, in which the population is discretized into relatively small size intervals, has the advantage of computing PSD directly (Hounslow et al., 1988). The past decade has witnessed a huge growth in the use of discrete method for gas–liquid and

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

liquid–liquid dispersions. For example, Rathore et al. (2012) modeled an aerated bioreactor equipped with an e-blade impeller. They employed the discrete method with thirteen bins to evaluate the effect of impeller speed, gas flow rate, and liquid height on mass transfer coefficient. In another study, Gelves et al. (2014) proposed the use of the discrete method to study a bioreactor equipped with a Rushton impeller. Azargoshasb et al. (2016) used CFD-PBM with the discrete method to simulate a gas–liquid dispersion in a bioreactor with eleven bin classes of arbitrary size to investigate the effect of hydrodynamic parameters on mass transfer, reaction behavior, and the bubbles Sauter mean diameter. Agahzamin and Pakzad (2019) simulated a bubble column with dense vertical internals. They have applied the Eulerian–Eulerian model coupled with discrete method (twenty bins) to evaluate the effect the internals (rods) gas holdup and local gas and liquid velocities Schütz et al. (2009) employed the mixture model and Lehr’s breakage model (Lehr, 2001) for modeling the DSD in water-diesel separation process in hydrocyclone. Their model successfully described the droplet breakage and coalescence in liquid–liquid dispersion. In another study by Håkansson et al. (2009), PBM was used to study fragmentation in emulsification process in homogenizer. Further research in this area by Srilatha et al. (2010) includes the droplet size distribution in a mixing tank for two different liquid–liquid dispersion systems (i.e. water-in-tri butyl phosphate and water-in-xylene). However; they assumed that aggregation and breakage occurred only among the droplets with equal diameter. Hermann et al. (2011) concluded that droplet breakage probability is higher at low dispersed phase viscosity. Roudsari et al. (2012) used Luo’s breakage model coupled with the Eulerian approach and k-ε turbulence model for the mixing of water in crude oils. They successfully validated the model with experimental study by Boxall et al. (2010). In another article proposed by Gelves et al. (2014), the discrete method was used to investigate the effect of turbulence, rotating flow, and breakagecoalescence in a stirred tank by using Rushton turbine or pitched blade impeller (PBT). They found that a PBT impeller can produce smaller bubbles and higher mass transfer for the cell culture process. Vonka and Soos (2015) used CFD and PBM for a liquid– liquid dispersion in a stirred tank equipped with Rushton turbine. They found that the critical Weber number is only dependent on droplet viscosity and it is independent of interfacial tension and local energy dissipation rate. The dependency of Sauter mean diameter on impeller tip speed was also predicted by their model. Sutherland et al. (2019) employed volume of fluid (VOF) multiphase method and PBM to simulate a moving-baffle oscillatory baffled column (moving-baffle OBC). The study has been marked for using the moving overset meshing to simulate agitator motion in a fluid system. Population balance results have been successfully validated with the experimental data for the inverse-suspension of non-reactive aqueous acrylamide in Isopar oil. Notwithstanding all fore-mentioned studies, more investigations with CFD modeling on hydrodynamic characteristics and droplet size distribution in the mixing of liquid–liquid dispersions are required. The CFD-PBM coupled models have been applied to study the impacts of variable volume fractions, viscosities of continuous phase, or interfacial tensions in the limited studies such as Li et al. (2017a,b). Available correlations for droplet size are limited to the assumption of dominant breakage in the dispersion system. In this study, CFD was coupled with PBM to evaluate the hydrodynamic mixing and size distribution of liquid–liquid dispersions with viscosity ranging from 3 to 262 mPa s, and volume fractions between 2.5% and 15.0%. The effect of impeller speed, oil viscosity, volume fraction, and the effect of the number of bin fractions on droplet size distribution have been investigated. The data for droplet size is fitted with available correlations in literature for two different turbulent regimes. The simulation was vali-

3

dated with the experimental droplet size distribution reported by Boxall et al. (2010). 2. CFD model development The mixing tank modeled in this study is a flat-bottom cylindrical tank with an internal diameter of 0.102 m. The liquid height in the tank is 0.098 m providing approximately 0.8l L volume of the fluid. A Rushton impeller with a diameter of 0.0508 m and blade thickness of 0.002 m was attached to a central shaft (0.0063 m) with an off-bottom clearance of 0.0515 m and four equal-spaced baffles (0.0102 m) were mounted to the wall. A schematic diagram of the mixing vessel with more details on dimensions is shown in Fig. 1. The tank was filled with a crude oil as the continuous phase and tap water as the dispersed phase. Four different crude oils have been considered in this study. The physical properties (Boxall et al., 2010) of the fluids are presented in Table 1. The volume fraction of the dispersed phase (i.e. water) was varied from 2.5% to 15.0% (i.e. 2.5%, 5.0%, 7.5%, 10.0%, and 15.0%). Three levels for impeller speed; 30 0 rpm, 40 0 rpm, and 60 0 rpm were modeled. The geometry and mesh were generated in ANSYS Design Modeler and ANSYS Mesh (ANSYS Inc., 2017), respectively. The fluid domain inside the tank was discretized with unstructured tetrahedral mesh (Pakzad et al., 2008) as shown in Fig. 1(b). The mesh density should be fine in the regions close to the impeller and wall. The mesh growth is controlled using size function. Size function allows the meshes to grow slowly as a function of the distance from the impeller blades and the tank wall. The accuracy of the results strongly depends on the quality of the grids (Nastac et al., 2012). Fig. 2 shows the average of the velocity magnitude of the continuous phase on an arbitrary line, 10 mm below the impeller for different grid size at an impeller speed of 300 rpm. Four different mesh sizes were generated as: 60,0 0 0; 150,0 0 0; 30 0,0 0 0; and 60 0,0 0 0. The discrepancy between them was calculated in terms of the root-mean-square (RMS) deviation (Pakzad et al., 2008). The RMS difference between the velocity magnitude for 30 0,0 0 0 and 60 0,0 0 0 mesh was calculated as 3.84%. In the grid independency study, the population balance modeling was not included. The results confirmed that 30 0,0 0 0 cells were sufficient for further simulation. To understand the dynamics of the system, the Reynolds num2 ber (Re = N D ρM/ where N is the impeller speed (rps), D is the

μM

impeller diameter (m), and ρ M and μM are the density (kg/m3 ) and viscosity (mPa.s) of the mixture, respectively) was calculated. The viscosity and density of the mixture were calculated through the following linear expression:

∅M = (1 − α )∅q + α ∅ p

(1)

where ∅M is the scalar variable of the fluid (viscosity or density) and α is the volume fraction of the dispersed phase. The subscript q and p refer to the continuous phase and dispersed phase, respectively. Calculated Reynolds numbers for different agitation speed confirm the turbulent regime in the tank. A complete time- and space-dependent solution of the exact Navier–Stokes equations for high-Reynolds-number turbulent flows is unlikely to be achievable. FLUENT offers two alternative methods to overcome this problem: Reynolds averaging and filtering. The Reynolds-averaged Navier– Stokes (RANS) equations are the transport equations based on the mean flow quantities. Thus, for the velocity vector of phase q as  q (u, v, w) with the mean value of U˜ (U, V, W ), the mass conu servation and momentum equations Drew, 1983) are written as Eqs. (2)–(3) (Bird et al. 2002), respectively:



   ∂ (αq ρq ) + ∇ . αq ρqU˜ = 0 ∂t

(2)

4

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

Fig. 1. Schematic diagram of the (a) mixing tank and Rushton impeller (dimensions in m), and (b) Numerical grid (unstructured tetrahedral mesh). Table 1 Physical properties of the fluids at standard condition (Boxall et al., 2010). Fluid

Density (kg/m3 )

Viscosity (mPa.s)

Interfacial Tension (mN/m)

Water Pure Conroe oil Troika oil 20% Conroe oil- 80% brightstock Petronius

1000 842 869 886 883

1.0 3.1 20.0 262.0 100.0

– 20 11 20 13

  v    t ∂ αq ρqU˜ + ∇ . αq ρqU˜ U˜ = −αq ∇ p¯ − ∇ . τq + τq ∂t + FD + Fg + FB + FL + FV

of phase q, U˜ is the average velocity of the same phase, p¯ is the v

(3)

where ρ q is the density of the phase q, α q is the volume fraction

t

average pressure, τq is the viscous momentum flux tensor, τq is the turbulent momentum flux tensor, FD is the drag force between phases, Fg is the gravitational force, FB is the Buoyancy force, FL is the lift force, and finally Fv is the virtual mass force.

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

5

Fig. 2. Velocity (m/s) magnitude of Conroe oil (continuous phase) on a line 10 mm below the impeller at 300 rpm with 15.0% water volume fraction in Conroe oil as a function of mesh size.

The drag force (FD ) is expressed as Eq. (4). The drag coefficient (CD ) proposed by Schiller and Naumann (1933) was used in this equation.

CD Re FD = 24   24 1 + 0.15Re0.687 CD = Re CD = 0.44

for Re ≤ 10 0 0 for Re > 10 0 0

(4)

where Reynolds number (Re) can be calculated based on the conq and V  p , respectively: tinuous and dispersed phase velocities, i.e. V

Re =

ρq|Vp −Vq | d p μq

(5)

where ρ q , μq , dp are the density of phase q (continuous phase) and viscosity of the phase q and the droplet diameter, respectively. The gravity force (Fg ) and buoyancy force (FB ) can be combined (Versteeg and Malalasekera, 2007) as follows:

Fg + FB =

π 6

d3p (ρ p − ρq )g

(6)

where g is the gravitational acceleration. The impact of the lift force (FL ) and the virtual mass force (Fv ) compared to the drag force are negligible in this study (Drew and Lahey, 1993) due to this fact that the density ratio between two phases are close to one and the size of the droplets are small. v The time-average viscous momentum flux tensor (τq ) presented in Eq. (3) (Bird et al., 2002) is calculated as follows: v

τq = −αq μq (∇ .U˜ + ∇ U˜ T ) + αq

2 3



μq − λq (∇ .U˜ )I

(7)

where μq and λq are the shear and bulk viscosity of phase q, ret

spectively. I is the unit tensor. The τq presented in Eq. (3) called turbulent momentum flux tensor is defined as follows: t

τq = −ρ ui uj

∂uj ∂ xi

(8)

where ui and uj are the fluctuating velocity component due to tur-

bulence in the directions i and j, respectively. The terms ui uj are the Reynolds stresses. For turbulence modeling, Reynolds-averaged Navier–stokes equations uses such models as Spalart–Allmaras model, k-ε models, k-ω models, and Reynolds stress model (RSM). The standard k-ε model Launder and Spalding, 1974) is a semi-empirical model based on model transport equations for the turbulence kinetic energy (k) and its dissipation rate (ε ) as presented in Eqs. (9) and (10), respectively. The turbulence kinetic energy (k) is derived from the exact equation (i.e. Eq. (9)), however this model underestimates the dissipation rate (ε ). In the derivation of the k-ε model, a fully turbulent flow has been assumed, and the effects of molecular viscosity have been neglected. Therefore, this model is valid only for fully turbulent flows. Among different available turbulence models, the standard k-ε model has been widely applied due to its robustness and accuracy in predicting the turbulent parameters, especially where the shear stresses are of interest (Heath and Koh, 2003; Marchisio et al., 2003; Kerdouss et al., 2006; Schütz et al., 2009; Gholamzadehdevin and Pakzad, 2018).

   ∂ (ρ k ) μt  ∇ k + Gk + Gb − ρε (9) + ∇ · ρ kU˜ = ∇ · μ + ∂t σk    ∂ (ρε ) μt  ∇ε + ∇ · ρεU˜ = ∇ · μ + ∂t σε ε ε2 + C1ε (Gk + C3ε Gb ) − C2ε ρ k

k

(10)

where U˜ is average velocity. The above model constants were retained as default values: C1ε = 1.44, C2ε = 1.92, and C3ε = tanh| υu |, where υ is the component of the flow velocity parallel to the gravitational vector and u is the component of the flow velocity perpendicular to the gravitational vector. The σk = 1.0, and σε = 1.3 indicate the k and ɛ turbulent Prandtl numbers. The Prandtl numbers σ k and σ ɛ connect the diffusivities of k and ε to the eddy viscosity μt . Gk and Gb are the source of turbulent kinetic energy, due to the mean velocity gradients and due to Buoyancy,

6

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

respectively. To compute the Reynolds stresses the familiar Boussinesq relationship (Eq. (11)) is applied as follows:

  2 τi j = −ρ ui uj = μt ∇ .U˜ + ∇ .U˜ T − ρq kI 3

1 ´2 2 (u

(11)

μt = ρqCμ

k2

(12)

ε

where Cμ is a dimensionless constant. Flow model for each phase based on Eulerian–Eulerian multiphase along with the standard kε turbulence model has been solved separately. The liquid–liquid dispersion is fundamentally coupled with droplet coalescence and breakage phenomena resulting in the droplet size distribution. The population balance equation for the dispersed phase (Hagesaether et al., 2002) is presented in Eq. (13):

∂ (ρ ndi ) + (ρvq ndi ) = ρ [BB − DB + BC − DC ]i ∂t

(13)

where BB and DB denote birth and death due to breakage, respectively, while BC and DC are representative of birth and death due to coalescence, respectively. The modeling of breakage and coalescence was performed using the discrete method Hounslow et al., 1988). Each droplet size interval in discrete method is assumed as an independent Eulerian phase (Schütz et al., 2009) since this balance is generally similar to other balances such as mass, energy, and momentum. This balance is applied to track changes in the droplet/particle size distribution due to imposed impacts of movement or reaction. The birth of a small size droplet from a larger one by breakage is usually assumed as a breakage source term. The droplet birth and death due to breakage are defined through Eq. (14)–(15), respectively (Hagesaether et al., 2002): N

BB ( i ) =

B (υK , υi ) +

K=i+1,i=N i−1

+

i



1 − xi,k B (υi , υK ),

i = 1, . . . , N

(14)

K=1,i=1

DB ( i ) =

i−1

d´ = 0.5

 5/3  ρ 7/5 q

d´

σ

 √ 

ε

19/5

exp −

2

3

(d´)

σ ρp

i = 2, . . . , N

(15)

k=1

The term B (υi , υK ) is the breakage rate of a droplet from size υi into the size of υK . Different empirical or semi-empirical models for breakage and coalescences phenomena can be found in literature (Nambier et al., 1992; Alopaeus et al., 1999, 2002; Luo and Svendsen, 1996; Lehr, 2001; Ghadiri and Zhang, 2002). In this study, breakage in droplets was modeled using Lehr’s breakage model (Lehr, 2001). Lehr (2001) proposed a breakage model for a liquid–gas dispersion system, however, the model was also applied for liquid dispersion



1

ε 6/5

In Eq. (16), the droplet breakage in liquid–liquid system happens due to turbulent eddies smaller than d´, the diameter of daughter droplet. External stresses from eddies of continuum destroy the fluid droplets; whereas, the surface and viscous stresses resist the droplet’s distortion (Revankar, 2001). The calculation of the DSD demands on daughter distribution function. Liao and Lucas (2009) have extensively discussed different daughter distribution functions. The daughter droplet size distribution proposed by Laakkonen et al. (2006) was used as follows:

β







d, d´ = 180

d2 d´3



d3 d´3

2 

d3 1− d´3

2

(17)

where d and d´ are the mother droplet size and daughter droplet size, respectively. More details about Lehr’s breakage mechanism can be found in Lehr (2001). The coalescence between droplets is defined by the priori assumption of having a sufficient rate of collisions among droplets with diameter of d and d . Coalescence depends on two major factors; collision rate and coalescence efficiency (the ratio of successful collisions to all possible collisions). Another key parameter in droplet coalescence is the droplet film thickness (Coulaloglou and Tavlarides, 1977). At the critical film thickness, the film breaks, droplets join together, and coalescence happens. In other words, the contact time between two droplets must be higher than the drainage time (Liu and Li, 1999). The droplets birth and death rate due to coalescence are defined as follows:

BC (i ) =

i−1

xi, j ag



vi , v j



j=1,i=N

+

i−1

  (1 − xi−1, j )ag vi−1 , v j ,

i = 2, . . . , N

(18)

j=1

DC (i ) =

N−1

  ag (vi , v j ) + ag vi , v j ,

i = 2, . . . , N

(19)

j=1

In this study, the Luo coalescence model (Luo, 1993) was applied. This model is a function of collision frequency and collision probability as follows:

(20)

where ωag (vi , vj ), the collision frequency, is defined by Eq. (21) and Pag (vi , vj ), the probability of the collision resulting in coalescence, is defined by Eq. (24):

  π 2 ωag vi , v j = di + d j ni n j u¯ i j

(21)

4

B (υi , υK ),

9 / 5

(16)

      ag vi , v j = ωag vi , v j Pag vi , v j

xi+1,k B (υi+1 , υK )

k=1,i=N



B

 

w´ 2 )

where k = + + is the turbulent kinetic energy per unit mass. The first term of the right hand side of Eq. (11) is analogous to Eq. (7) except for the turbulent or eddy viscosity (μt ), while the second term of the right hand side is devoted to the normal Reynolds stresses (where i = j). The main assumption of the k–ε model is that the turbulent (or eddy) viscosity, μt , is isotropic; it means that the ratio between Reynolds stress and mean rate of deformation is equal in all directions. The k–ε model (with the isotropic eddy viscosity assumption) are found to be almost insensitive to the system rotation (Aubin et al., 2004). However, in spite of such expected problem, the k–ε model has also showed reasonable results in the modeling of turbulence in the rotating flows such as studies done by Roudsari et al. (2012), Hossieni et al. (2010), and Aubin et al. (2004). The turbulent (or eddy) viscosity (μt ) is defined as:

v´ 2

in a study by Schütz et al. (2009). The equation for breakage rate according to Lehr’s model is expressed as follows:



u¯ i j = u¯ 2i + u¯ 2j

1 / 2

(22)

where u¯ i j is the characteristic velocity of collision between two droplets, ni and nj are number densities, and u¯ i is defined as follows:

u¯ i = 1.43(ε di )

1/3

(23)

The probability of the collision (Luo, 1993) is defined as follows:

Pag = exp −c1





0.75 1 + x2i j

 ρ2



1 + x3i j

1 / 2 

ρ 1 + 0.5



1/2

1 + xi j

3

W e1i j/2

(24)

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

where ρ 2 is density of dispersed phase, ρ 1 is density of contindi uous phase, and xi j = dj for droplets i and j. Thus, the dynamic pressure difference across the droplet (τ i ) against its surface stress (τ s ) can determine whether or not the droplet breaks in the system (Liao and Lucas, 2009). The ratio of the dynamic pressure to τ droplet surface stress is called Weber number, W e = τsi . This number indicates whether the kinetic energy or surface tension of the droplet is dominant. The We number presented in Eq. (24) is defined as follows:

 2

W ei j =

ρl di u¯ i j σ

7

sented by Eq. (26)) was applied in this study to present DSDs.



1 (lnX − μ ) f d (X ) = √ exp − 2σ 2 2πσ X

2



(26)

where X is the bin fraction size, σ is the standard deviation, and μ is the median of the distribution defined as follows:

 

σ=

ln 1 +

ν 2 dm

 and

μ = ln 

dm 1 + dν2

(27)

m

(25)

In order to model relative motion between the stationary zone and the rotating impeller blades, many authors have used either the multiple reference frame (MRF) (Rathore et al., 2012; Pakzad et al., 2008; Luo et al., 1994) or the sliding mesh method (Wutz et al., 2016; Kazemzadeh et al., 2016b). Since the sliding mesh approach, albeit more accurate, is time consuming, we employed the MRF approach to capture the motion in the vicinity of the impeller. The MRF approach has been also applied for mixing in multiphase flow such as Hossieni et al. (2010), Fan et al. (2010), and Roudsari et al. (2012). The standard wall functions in FLUENT were applied for the modeling of the near wall regions. The logarithmic law for mean velocity is known to be valid for 30 < Y + < 300. The Y+ values calculated in this study were within these limits (Launder and Spalding, 1974). The symmetry boundary condition has been applied frequently in modeling of mixing tanks for example in such studies as Pakzad et al. (2008), Kazemzadeh et al. (2016a) and Roudsari et al. (2012). Through this boundary condition, the normal velocity component and tangential velocity gradient were set to zero for the free surface of the liquid as a symmetry surface (i.e. there is no flow, convective, or diffusive flux across the symmetry surface). The SIMPLE algorithm was used for pressure-velocity coupling. The second-order discretization scheme was used for all equations except volume fraction. Volume fraction is generally discretized using QUICK scheme for rotating or swirling flows as it is more accurate than second order Upwind scheme (Rathore et al., 2012; Escue and Cui, 2010). The Multiple Sized-Group (MUSIG) method was used to solve the PBM equations (Yeoh et al., 2014). The MUSIG method is a discrete method that discretizes the droplet population into a finite number of size intervals (i.e. bins). The convergence was checked for the following equations: continuity equation, momentum volume fraction, turbulence equations, and water bin fractions equations. FLUENT software (version 17.2) was employed to solve both hydrodynamic and population balance equations. The time step of 0.01 s has been chosen to solve all equations for 35 s. To use the simultaneous solution method, supercomputing facilities of HPCVL was used. Each simulation was partitioned into 24 parts and each CPU was then used for each part. The simulation was considered converged when the scaled residuals for all transport equations were below 10−4 . The convergence for each simulation was achieved after 22 h. 3. Results and discussion The sufficient number of bin fractions in discrete method is a prominent key for population balance modeling (PBM). Small number of bin classes leads to inaccurate DSD while using large numbers of bin classes is computationally expensive. Among various distribution equations such as single probability density function (Zhuo and Kresta, 1998), linear normal distribution (Calabrese et al., 1986), and logarithmic normal distribution (Boxall et al., 2010), the logarithmic normal distribution (as pre-

where ν is the variance of DSDs and dm is the mean of all bin sizes. Fig. 3 shows droplet size distributions for four different number of bin fractions as 7, 15, 20, and 35 bin fractions for 5.0 vol% of water in Conroe oil agitated at 400 rpm. A nearly identical distribution for 20 and 35 bin fractions can be observed. The similarity in the profiles of droplet size distribution for 20 and 35 bin fractions resulted in choosing 20 bin fractions for further simulation. Fig. 4 shows the cumulative probability droplet size distribution of 15.0 vol% water in Conroe oil as a function of impeller speed. As the impeller speed increased, the size of the water droplets gradually shifted to the smaller droplets. As a result, the maximum droplet size of approximately 290 μm for 300 rpm were substituted by 210 and 135 μm at 400 and 600 rpm, respectively. This shows that the breakage phenomenon was more dominant at higher impeller speeds as expected. In fact, an increase in impeller speed results in higher turbulent pressure fluctuation. The turbulent kinetic energy provided by eddies with diameter of λ ≤ d break the mother droplet to the smaller droplets (Simon et al., 2003). An increased flow velocity gradient also increases the oscillations and the breakage probability of oil droplets (Stamatoudis and Tavalarides, 2007). Consequently, the fraction of oil droplets with small diameter is increased and the distributions are shifted to smaller sizes. The shift in droplet size distributions toward smaller size with higher agitation intensity was also noticed in other studies (Gerstlauer, 1999; Ramkrishna, 20 0 0; Alopaeus et al., 2002). Fig. 5 presents the logarithm number density of seven bin fractions at different impeller speeds for 15.0% water volume fraction. Seven arbitrary bin fractions were chosen in a way to have two small size droplet classes, three medium size droplet classes, and two large size droplet classes. The comparison between the bin classes at 300 rpm and 400 rpm shows an increase in number densities for most of the bins due to turbulent intensity induced by higher impeller rotation (Ruiz et al., 2002; Ok et al., 20 03; Schütz et al., 20 09). Further analysis of logarithmic number densities for the middle-sized droplets (i.e. classes 3, 4, and 5) especially at 400 rpm and 600 rpm shows higher numbers in these classes which justify higher coalescence rate among small droplets. The higher coalescence rate has been occurred due to increased characteristic velocity of collision between droplets as shown in Eq. (20). As the impeller speed increased to 600 rpm, a higher coalescence rate was induced by a higher adhesion energy and a higher collision rate among the middle size classes. This behavior has been also reported in Ruiz et al. (2002), Ok et al. (2003), and Schütz et al. (2009). A decrease in the size of droplets leads to an increase in cohesive stress and a decrease in disruptive stress, until droplet reaches a certain size, at which point the rate of breakage and coalescence will be the same (Rueger and Calabrese, 2013). Fig. 6 shows the vertical cross section contour plots for the 15.0% volume fraction of dispersed phase at different impeller speeds. This figure depicts a relative non-homogeneity in the vessel at the impeller speed of 300 rpm compared to 400 and 600 rpm. As can be seen in Fig. 6, the water volume fraction is

8

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

Fig. 3. Droplet size distribution for 5.0% water volume fraction in Conroe oil agitated at 400 rpm in different number of classes.

Cumulative probability distribution, (-)

1

EXP data-300rpm EXP data-400rpm EXP data-600rpm CFD data-300rpm CFD data-400rpm CFD data-600rpm

0.8

0.6

0.4

0.2

0 10

100 Droplet size, (µm)

1000

Fig. 4. Experimental data (Boxall et al., 2010) and CFD results for Conroe oil with 15.0% water volume fraction at different impeller speeds.

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

9

Fig. 5. Logarithmic number density of droplets for 7 bin classes (small, medium, and large size droplets) in Conroe oil with 15.0% water volume fraction.

Fig. 6. Contour plots of water volume fraction for Conroe oil with 15.0% volume fraction at impeller speed of: (a) 30 0 rpm, (b) 40 0 rpm, and (c) 60 0 rpm. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

higher in the region beneath the impeller at 300 rpm. The higher density of water, and insufficient pumping power generated by the Rushton impeller at 300 rpm are the two factors leading to the higher volume fraction of water at the bottom of the tank. The insufficient flow motion was recovered by increasing the impeller speed to 400 rpm, therefore, the green region at the tank surface becomes smaller. As can be seen, still some degree of inhomogeneity can be observed at the bottom of the tank. With an increase in impeller speed to 600 rpm, a full homogenous flow was obtained due to higher turbulent kinetic energy and improved recirculation rate induced by impeller (Maggioris et al., 1998). Since the dispersed phase is located at the bottom of the tank before mixing, and the pumping rate of impeller is insufficient at 300 and

400 rpm; a lower impeller off-bottom clearance could also help in higher circulation rate and more homogenous flow within the tank at lower impeller speeds. Fig. 7 shows the contour of turbulent kinetic energy on five horizontal planes for Conroe oil at 15.0% volume fraction and impeller speed of 300 rpm. This figure shows that the turbulent kinetic energy was more intense at the center of plane 3 in compared with other horizontal planes, whereby less turbulent kinetic energy was induced by flow circulation at the center of other planes. The distribution of turbulent kinetic energy on plane 1 and plane 2 shows lower values of turbulent kinetic energy at the center while more energy is distributed near the walls. The average turbulent kinetic energy on planes are as follows: Plane 1: 0.0108 m2 s−2 , Plane 2:

10

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

Fig. 7. Contour of turbulent kinetic energy for Conroe oil with 15.0% water volume fraction at 300 rpm, for plane 1: 0.0825 m, plane 2: 0.0670 m, plane 3: 0.0515 m, plane 4: 0.0360 m, and plane 5: 0.0205 m.

0.0153 m2 s−2 , Plane 3: 0.0457 m2 s−2 , Plane 4: 0.0176 m2 s−2 , and Plane 5: 0.0125 m2 s−2 . The average turbulent kinetic energy shows higher values for the planes above the impeller. This result can be the evidence for lower droplet breakage below the impeller (Maggioris et al., 20 0 0) which can be enhanced by decreasing the impeller off-bottom clearance. Fig. 8 shows the effect of dispersed phase volume fraction on DSD at 400 rpm impeller speed. Similar bell-shaped curves are observed at different volume fractions, however, dispersions with

higher volume fraction show larger number densities which was expected. A slightly narrower size distribution profile can be observed as volume fraction decreased which agrees with a study by Liu et al. (2013). According to Liu et al. (2013), an increase in volume fraction of the dispersed phase (i.e. ranging from 1.0% to 5.0%) results in a wider droplet size distribution. Wang et al. (2013) obtained the same result for evolution of chord length distribution (CLD) at different oil volume fractions (10.0–60.0%) in a mixing tank. However, Qi et al. (2015) reported no pronounced

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

11

Fig. 8. Droplet size distribution for Conroe oil at different water volume fractions and 400 rpm impeller speed.

Fig. 9. Droplet size distribution for different crude oils at impeller speed of 400 rpm and 5.0% water volume fraction.

effect of volume fraction on DSD for high concentrated emulsions. Boxall et al. (2010) also mentioned no pronounced difference in DSD of water-crude oil dispersion for volume fractions of 10.0%– 35.0%. The viscosity of continuous phase is another important factor as it affects the stability of the dispersion system and also droplets size (Boxall et al., 2010). Fig. 9 shows droplet size distribution for different crude oils at the impeller speed of 400 rpm

and 5.0% water volume fraction. As mentioned in Table 1, four different crude oils with viscosities ranging from 3.1 mPa s to 262 mPa s were considered as the continuous phase. The effect of interfacial tension and viscosity on droplet size as well as the breakage/coalescence rate in a liquid–liquid dispersion system is not well-defined due to the inherent complexity existing among the droplets’ interactions (Patel et al., 2015). Increased viscosity resulted in lower film drainage rate and lower coalescence rate

12

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

Fig. 10. Contour plots of: (a) breakage source term, (b) coalescence source term, (c) Sauter mean diameter, and (d) velocity vector of continuous phase (Conroe oil) with 10.0% water volume fraction and impeller speed of 400 rpm.

between liquid droplets in high-viscous fluids (Roudsari et al., 2012), thus, the distributions were shifted toward smaller sizes, and narrower size distributions were achieved. This can be easily observed between Conroe oil and Conroe oil/Brightstock and between Troika and Petronius; with different viscosities but almost same interfacial tensions. However, result for Conroe-Brightstock and Petronius confirms the effect of interfacial tension, although the viscosity of Petronius is lower, its interfacial tension is lower too. Therefore, lower interfacial tension resulted in a reduction

in the retaining force of a droplet, a higher breakage rate, and a smaller droplet size in Petronius oil. The evolution of droplet size distribution is influenced by the dynamic breakage and coalescence happening in the process until the steady state condition is reached. Fig. 10(a) shows the contour of breakage source term for 10.0% water volume fraction dispersed in Conroe oil at 300 rpm impeller speed. The birth/death source term due to breakage and coalescence are presented in Eqs. (14)– (15) and Eqs. (18)–(19), respectively. As can be seen, the breakage

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

13

Fig. 10. Continued

source term within the tank is negligible compared to the region around the impeller edges with high turbulent dissipation rate, and on the adjacent wall. Although the size of this region is small, high velocity magnitude of the flow induced by impeller in this region (Fig. 10(d)) results in a high pumping rate through the tank which has considerable impact on droplet size. A successful coalescence between droplets depends on collision efficiency and collision frequency. The contour of coalescence source term at the same operational condition is presented in Fig. 10(b). As shown in this figure, the coalescence happens everywhere in the vessel. According to

the study by Vonka and Soos (2015), higher turbulence in impeller region leads to lower coalescence rate for droplets since there is not enough time for droplets to coalesce. The coalescence source term is nonzero within the tank, although its value is smaller than the value of breakage source term. By comparing the rate of breakage and coalescence source terms, we can conclude that the droplet size distribution is mainly controlled by the breakage. Fig. 10(c) shows the contour of Sauter mean diameter at 10.0% water volume fraction agitated at 300 rpm in Conroe oil which was obtained by using Eq. (28) through discrete method in PBM as

14

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

Fig. 11. Comparison of experimental data (Boxall et al., 2010) and calculated values (using Eq. (29) and CFD inputs) for maximum stable droplet size (dmax,D ) for various continuous phase viscosities and impeller speeds at 15.0% water volume fraction.

follows (Khakpay and Abolghasemi, 2010):

n

d32 =

3 i=1 ni di n 2 i=1 ni di

(28)

where ni and di are the number and diameter of droplet in DSD, respectively. A relatively homogenous profile of Sauter mean diameter can be observed due to good mixing in the tank and strong velocity profile as shown in Fig. 10(d). Fig. 10(d) shows the velocity vector of the continuous phase. It can be noticed that two symmetrical liquid vortices were formed in the upper and lower region of the impeller as expected for the radial-flow Rushton impeller. A strong radial flow generated around the impeller results in strong flow circulation rate in breakage region. According to Kolmogorov’s theory (1949), two different turbulent regimes were defined; “turbulent inertial” and “turbulent viscous” regimes. To determine the mechanism of energy dissipation contributing to droplet break-up, Kolmogorov (1949) defined an eddy length scale (λ). Based on Eq. (29), droplets with the diameter larger than the eddy length scale fall into turbulent inertial regime where the interfacial stresses are balanced by inertial stresses and viscous forces can be neglected. In the turbulent inertial regime, the maximum droplet size (dmax ) is determined by the balance between turbulent fluctuations against the droplet capillary pressure. However, in the turbulent viscous regime the droplet diameter is smaller than the eddy length scale and dmax is determined by the balance between the viscous forces against the droplet capillary pressure (Vankova et al., 2007). Shinnar (1961) concluded that turbulent inertial regime is the dominant mechanism for droplet breakup in turbulent dispersions in stirred tank.



λ=

ν3 ε

 14

(29)

where λ is the Kolmogorov eddy length scale, ν is the kinematic viscosity, and ɛ is the local energy dissipation rate. To distinguish the regime proposed by Kolmogorov (Kolmogorov, 1949), the size of smallest eddies should be determined. The change in the regime depends not only on the size of the smallest eddies, but also on the viscosity of continuous

phase and the maximum stable droplet size (dmax ) in dispersion (Vankova et al., 2007). Among various correlations for dmax, an equation by Davies (1985) can describe the maximum stable drop size, in turbulent inertial regime, for the viscosity range of 3 and 100 mPa s which is applicable to our study:



1/3 μD ε 1/3 dmax,D dmax,D = 0.86 1 + 0.37 σ

3 / 5

σ 3/5 ρC−3/5 ε −2/5

(30)

where dmax, D is the maximum stable droplet size in inertial regime according to the study by Davies (1985) and μD is the viscosity of dispersed phase. Fig. 11 presents the comparison between experimental data (Boxall et al., 2010) and calculated values (from Eq. (30)) for the maximum stable droplet size for different viscosities and different impeller speeds at volume fraction of 15.0%. The rate of energy dissipation rate for calculation of maximum droplet size was obtained by CFD model after the steady state condition was reached. As can be seen in this figure, although the calculated data for Conroe oil and Troika oil are well correlated by Davies equation; a significant discrepancy between the experimental data and the theoretical data for viscous oils, Conroe-Brightstock and Petronius, was observed. The maximum droplet size decreased significantly as the viscosity of continuous phase changed from 20 mPa s to 100 mPa s. Therefore, Eq. (30) cannot predict the maximum stable droplet size at higher viscosity of the continuous phase (ηc > 20 mPa.s). This is a good indication of a change in turbulent liquid–liquid dispersion regime. According to Kolmogorov theory (Kolmogorov, 1949), as the viscosity of continuous phase increased, the size of the smallest turbulent eddies decreased. Consequently, droplets could become smaller than λ and the regime would change to viscous turbulent regime. The maximum stable droplet size in viscous regime can be obtained by following equation (Kolmogorov, 1949): /2 −1/2 dmax,viscous = 4ε −1/2 μ−1 ρc σ c

(31)

where dmax, viscous is the maximum stable droplet size in turbulent viscous regime, μc is the viscosity of continuous phase, and ρ c is the density of continuous phase. Table 2 shows a comparison between the calculated values by Eq. (31) and experimental data

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

15

Fig. 12. Correlation of dma x with: (a) We in turbulent inertial regime, and (b) We.Re in turbulent viscous regime, at different impeller speeds and 15.0% water volume fraction.

obtained by Boxall et al. (2010). As can be noticed in this table, the maximum droplet sizes are less than eddy length scales (calculated using Eq. (29)) indicating the viscous regime in dispersions based on Kolmogorov’s theory. As can be seen in Table 2, the calculated values of maximum drop size using the values of ɛ by CFD

model agree well with corresponding experimental data with the error less than 20.0%. The discrepancy in values of theoretical and experimental maximum drop size is related to imprecision to calculate the rate of ɛ, insufficient time to reach the steady state condition in experiments, and the theoretical constant in Eq. (31).

16

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100 Table 2 Comparison of experimental data (Boxall et al., 2010) and calculated values (using Eq. (30) and CFD inputs) for maximum stable droplet size (dmax,viscous ) for various continuous phase viscosities and impeller speeds at 15.0% volume fraction in turbulent viscous regime.

ηc (mPa.s)

Impeller speed (rpm)

dmax , Experimental (Boxall et al., 2010) (μm)

dmax ,viscous, Calculated using Eq. (30) & CFD inputs (μm)

λ ( μm )

100 100 262 262

300 400 300 600

58.8 39.8 103.9 38

64.1 44.6 92.5 33.4

78.7 76.2 118.6 62.3

The results of Fig. 11 and Table 2 show that the maximum droplet size is under control by different parameters depending on the type of turbulent regime. The identification of regimes are of interest for the researchers since much smaller droplets can be obtained in the viscous regime in compared to the inertial regime at the equivalent condition. One can identify the regime by investigation the effect of continuous phase viscosity on maximum droplet size. If the variations in droplet size are not significant (e.g. up to 2–3 times), then the turbulent regime falls in turbulent inertial (Davies, 1985). However, the change in the viscosity of continuous phase results in considerable decrease of droplet size with dmax , this is a good indication of turbulent viscous regime. A closer look at Eq. (31) reveals that viscosity of continuous phase is the most efficient factor in changing the regime due to its significant effect on eddy length size and maximum stable droplet size in viscous regime. However, viscosity of continuous phase has no effect on the maximum stable droplet size in inertial regime (Davies, 1985) as presented in Eq. (30). The first studies on liquid–liquid dispersion were limited to stirred tanks with low viscosity continuous phase or dilute emulsions (Rueger and Calabrese, 2013); however, industrial processes take advantage of different liquid–liquid dispersion system with a wide range of viscosity. Various models and theories have been established to correlate the mean and maximum droplet sizes to dimensionless parameters such as Weber number and Reynolds number. An appropriate mechanistic correlation between the size of droplets and dimensionless hydrodynamic parameters such as Weber number needs to be determined based on the droplet size relative to the Kolmogorov eddy length scale (Kolmogorov, 1949). Many authors such as Rueger and Calabrese (2013), El-Hamouz et al. (2009), and Pacek et al. (1999) obtained a correlation between the maximum stable droplet size and the Weber number in turbulent inertial regime as follows:

dmax ∝ εT−0.4 ∝ W e D

−0.6

ρc N 2 D 3 σ

−1 dmax ∝ (W e Re ) 3 D

4. Conclusions Computational fluid dynamics (CFD) in conjunction with the population balance modeling (PBM) were employed to study water-in-crude oil dispersions in a stirred tank equipped with a Rushton turbine. The Eulerian–Eulerian multiphase approach, k-ɛ turbulent model, and MRF technique, were implemented in the model. The PBM provided satisfactory results with the use of Lehr’s breakage model and Luo’s coalescence model in liquid–liquid dispersion in the vessel. The effect of agitation speed, the volume fraction of dispersed phase, and the viscosity of continuous phase were investigated. •

•

•

(33)

where ρ c is the density of continuous phase, N is the impeller speed, D is the impeller diameter, and σ is the surface tension. Fig. 12(a) shows the correlation for maximum droplet size and Weber number in turbulent inertial regime for Conroe oil and Troika oil at 15.0% volume fraction and different impeller speeds. The theoretical maximum droplet sizes are well correlated with Eq. (32) with regression value over 98%. This is in agreement with the data reported by El-Hamouz et al. (2009) for silicon oil with addition of surfactant; however, the correlation for our results has higher exponent (−0.67) which is attributed to higher coalescence rate of droplets without using surfactant in our system. Our findings show that dispersions with higher volume fractions have higher value of exponent compared to the value of −0.6. On the other hand, Shinnar (1961) presented a correlation for turbulent

(34)

where dmax is the maximum droplet size in dispersion, D is impeller diameter, We and Re are the Weber number and Reynolds number, respectively. As shown in Fig. 12(b), the dmax of ConroeBrightstock and Petronius oils with 15.0% water volume fraction at different impeller speeds vary inversely with almost cube root of the (We.Re), which is in good agreement with the literature (Shinnar, 1961). Rueger and Calabrese (2013) mentioned that in dispersion when the eddy length scale is greater than the droplet size, but not by orders of magnitude, the inertial stress can be ignored. This correlation provides satisfactory results for viscous oils in turbulent viscous regime with regression value over 99%.

(32)

where dmax is maximum stable droplet size, D is impeller diameter, ɛT is mean energy dissipation rate per unit mass, and We is Weber number. The Weber number in stirred tanks is defined as follows:

We =

viscous regime as follows:

An increased turbulent energy dissipation rate, generated at higher impeller speeds, leads to droplets of smaller size and narrower size distribution. Higher impeller speed results in predominant breakage in the system and a shift in droplet size distribution (DSD) towards smaller droplets. The effect of volume fraction on DSD was found to be negligible while wider distributions were obtained at higher volume fractions. Additionally, an increase in the viscosity of the continuous phase, as well as a decrease in the interfacial tension, lead to smaller droplet size and narrower droplet size distribution. The results showed that the dispersion could be performed in two different regimes, i.e. turbulent inertial and turbulent viscous regimes.

The maximum droplet diameter depends slightly on volume fraction and the viscosity of continuous phase in turbulent inertial regime. Transition from inertial to turbulent viscous regime was accomplished with an increase in continuous phase viscosity and smaller droplets were obtained in turbulent viscous regime. The theoretical droplet size was calculated by Davies expression using the values of energy dissipation rate from CFD model, and then were compared by experimental data in inertial regime. The maximum droplet size in turbulent viscous regime was successfully correlated with the product of Weber number and Reynolds number.

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

Acknowledgment We gratefully acknowledge the financial support for this research provided by the Natural Sciences and Engineering Research Council of Canada (NSERC-DISCOVERY: 11-50-14050129). The authors would also like to express their appreciation to the Center for Advanced Computing (formerly HPCVL) for running intense simulations for this study using their high-performance facilities. Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ijmultiphaseflow.2019. 103100. References Agahzamin, S., Pakzad, L., 2019. A comprehensive CFD study on the effect of dense vertical internals on the hydrodynamics and population balance model in bubble columns. Chem. Eng. Sci. 193, 421–435. Alopaeus, V., Koskinen, J., Keskinen, J., 1999. Simulation of the population balance for liquid–liquid system in nonideal stirred tank, part 1- description and qualitative validation of the model. Chem. Eng. Sci. 54, 5887–5899. Alopaeus, V., Koskinen, J., Keskinen, K., 2002. Simulation of the population balance for the liquid–liquid systems in noideal stirred tank. Part 2: parameter fitting and the use of multi block model for dense dispersion. Chem. Eng. Sci. 57, 1815–1825. Angle, C.W., Hamza, H.A., 2006. Predicting the sizes of toluene-diluted heavy oil emulsions in turbulent flow part 2: Hinze–Kolmogorov based model adapted for increased oil fractions and energy dissipation in a stirred tank. Chem. Eng. Sci. 61, 7325–7335. Aubin, J., Fletcher, D.F., Xuereb, C., 2004. Modeling turbulent flow in stirred tanks with CFD: the influence of the modeling approach, turbulence model and numerical scheme. Exp. Therm. Fluid Sci. 28, 431–445. Azargoshasb, H., Mousavu, S., Jamialahmadi, O., Shojaosadati, S., Mousavi, S., 2016. Experiments and a three-phase computational fluid dynamics (CFD) simulation coupled with population balance equations of a stirred tank bioreactor for high cell density cultivation. Can. J. Chem. Eng. 94, 20–32. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 1960. Transportphenomena. New York, Wiley. Boxall, J., Koh, C., Sloan, E., Sum, A., Wu, D., 2010. Measurement and calibration of droplet size distribution in water-in-oil emulsions by particle video microscope and a focused beam reflectance method. Ind. Eng. Chem. Res. 49, 1412–1418. Brown, D.E., Pitt, K., 1972. Drop size distribution of stirred non-coalescing liquid–liquid system. Chem. Eng. Sci. 27, 577–583. Calabrese, R., K., C.T.P., Dang, P., 1986. Drop breakup in turbulent stirred-tank contactors part I: effect of dispersed phase viscosity. AICHE J. 32, 657–666. Carrillo de hert, S., Rodgers, T., 2017. On the effect of dispersed phase viscosity and mean residence time on the droplet size distribution for high-shear mixers. Chem. Eng. Sci. 172, 423–433. Coulaloglou, C.A., Tavlarides, L.L., 1977. Description of interaction processes in agitated liquid-liquid dispersions. Chem. Eng. Sci. 32 (11), 1289– 1297. Davies, J.T., 1985. Drop sizes of emulsions related to turbulent energy dissipation rates. Chem. Eng. Sci. 40, 839–842. Drew, D.A., 1983. Mathematical modeling of two-phase flow. Annu. Rev. Fluid Mech. 15, 261–291. Drew, D.A., Lahey, R., 1993. In Particulate Two-Phase Flow. Butterworth-Heinemann, Boston. USA. El-Hamouz, A., Cooke, M., Kowalski, A., Sharratt, P., 2009. Dispersion of silicone oil in water surfactant solution: effect of impeller speed, oil viscosity and addition point on drop size distribution. Chem. Eng. Process. Process Intensification 48, 633–642. Escue, A., Cui, J., 2010. Comparison of turbulence models in simulating swirling pipe flows. Appl. Math. Model. 34, 2840–2849. Fan, L., Xu, N., Wang, Z., Shi, H., 2010. PDA experiments and CFD simulation of a lab-scale oxidation ditch with surface aerators. Chem. Eng. Res. Des. 88, 23–33. Gelves, R., Dietrich, A., Takors, R., 2014. Modeling of gas–liquid mass transfer in a stirred tank bioreactor agitated by a Rushton turbine or a new pitched blade impeller. Bioprocess Biosyst. Eng. 37, 365–375. Gerstlauer, A., 1999. Herleitung und reduktion populationsdynamischer modelle am beispiel der flussig-flussig- extraction (Ph.D. Thesis). University of Stuttgart, Germany. Ghadiri, M., Zhang, Z., 2002. Impact attrition of particulate solids part 1. A theoretical model of chipping. Chem. Eng. Sci. 57, 3656–3669. Gholamzadehdevin, M., Pakzad, L., 2018. The hydrodynamic characteristics of an activated sludge bubble column through computational fluid dynamics (CFD) and response surface methodology (RSM). Can. J. Chem. Eng. 97, 967–982. Hagesaether, L., Jakobsen, H., Svendsen, H., 2002. A model for turbulent binary breakup of dispersed fluid particles. Chem. Eng. Sci. 57, 3251–3267. Håkansson, A., Trägårdh, C., Bergenståhl, B., 2009. Dynamic simulation of emulsion formation in a high pressure homogenizer. Chem. Eng. Sci. 64, 2915–2925.

17

Heath, A.R., Koh, P.T.L., 2003. Combined population balance and CFD modeling of particle aggregation by polymeric flocculent. Third International Conference on CFD in Minerals and Process Industries. CSIRO. Hermann, S., Maaß, S., Zedel, D., Walle, A., Schäfer, M., Kraume, M., 2011. Experimental and numerical investigations of drop breakage mechanism. In: Proceedings 1st International Symposium on Multiscale Multiphase Process Engineering. Kanazawa, Japan. Hossieni, S., Patel, D., Ein-Mozaffari, F., Mehrvar, M., 2010. Study of solid–liquid mixing in agitated tanks through computational fluid dynamics modeling. Ind. Eng. Chem. Res. 49, 4426–4435. Hounslow, M., Ryall, R., Marshall, V., 1988. A discretized population balance for nucleation, growth, and aggregation. AICHE J. 34, 1821–1832. Kazemzadeh, A., Ein-Mozaffari, F., Lohi, A., Pakzad, L., 2016a. Effect of the rheological properties on the mixing of Herschel–Bulkley fluids with coaxial mixers: applications of tomography, CFD, and response surface methodology. Can. J. Chem. Eng. 94, 2394–2406. Kazemzadeh, A., Ein-Mozaffari, F., Lohi, A., Pakzad, L., 2016b. Investigation of hydrodynamic performances of coaxial mixers in agitation of yield-pseudoplastic fluids: single and double central impellers in combination with the anchor. Chem. Eng. J. 294, 417–430. Kerdouss, F., Bannari, P., Proulx, P., 2006. CFD modeling of gas dispersion and bubble size in a double turbine stirred tank. Chem. Eng. Sci. 61, 3313–3322. Khajeh Naeeni, S., Pakzad, L., 2019. Experimental and numerical investigation on mixing of dilute oil in water dispersions in a stirred tank. Chem. Eng. Res. Des. 147, 493–509. Khakpay, A., Abolghasemi, H., 2010. The effects of impeller speed and holdup on mean drop size in a mixer settler with spiral-type impeller. Can. J. Chem. Eng. 88, 329–334. Kolmogorov, A., 1949. Droblenii kapelv turbulentnom potoke. Dokl. Akad. Nauk SSSR 66, 825–828. Laakkonen, M., Alopaeus, V., Aittamaa, J., 2006. Validation of bubble breakage, coalescence and mass transfer models for gas–liquid dispersion in agitated vessel. Chem. Eng. Sci. 61, 218–228. Launder, B., Spalding, D., 1974. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 3, 269–275. Lehr, F., 2001. Berechnen Von Blasengrobenverteiungen Und Stromungsfeldern in Blasensaulen. University of Hannover, Germany. Lemenand, T., Della Valle, D., Zellouf, Y., Peerhossaini, H., 2003. Droplets formation in turbulent mixing of two immiscible fluids in a new type of static mixer. Int. J. Multiphase Flow 29, 813–840. Li, D., Buffo, A., Podgorska, W., Marchisio, D.L., Gao, Z., 2017a. Investigation of droplet breakup in liquid–liquid dispersions by CFD– PBM simulations: the influence of the surfactant type. Chin. J. Chem. Eng. 25, 369–1380. Li, D., Gao, Z., Buffo, A., Podgorska, W., Marchisio, D.L., 2017b. Droplet breakage and coalescence in liquid–liquid dispersions: comparison of different kernels with EQMOM and QMOM. AIChE J. 63, 2293–2311. Liao, Y., Lucas, D., 2009. A literature review of theoretical models for drop and bubble breakup in turbulent dispersions. Chem. Eng. Sci. 64, 3389– 3406. Liu, C., Li, M., Liang, C., Wang, W., 2013. Measurement and analysis of bimodal droplet size distribution in a rotor-stator homogenizer. Chem. Eng. Sci. 102, 622–631. Liu, S., Li, D., 1999. Drop coalescence in turbulent dispersions. Chem. Eng. Sci. 54, 5667–5675. Lovick, J., Mouza, A.A., Paras, S.V., Lye, G.J., Angeli, P., 2005. Drop size distribution in highly concentrated liquid–liquid dispersions using a light back scattering method. J. Chem. Technol. Biotechnol 80, 545–552. Luo, H., 1993. Coalescence, break up and liquid circulation in bubble column reactors (Ph.D. Thesis). Norwegian Institute of Technology, Trondheim, Norway. Luo, H., Svendsen, H., 1996. Theoretical model for drop and bubble breakup in turbulent dispersions. AICHE J. 42, 1225–1233. Luo, J., Issa, R., Gosman, A., 1994. Prediction of impeller-induced flows in mixing vessels using multiple frames of reference. In: Institution of Chemical Engineers Symposium Series, 136, pp. 549–556. Maggioris, D., Goulas, A., Alexopoulos, A.H., Chatzi, E.G., Kiparissides, C., 1998. Use of CFD in prediction of particle size distribution in suspension polymer reactors. Comput. Chem. Eng. 22, 315–322. Maggioris, D., Goulas, A., Alexopoulos, A.H., Chatzi, E.G., Kiparissides, C., 20 0 0. Prediction of particle size distribution in suspension polymerization reactors: effect of turbulence nonhomogeneity. Chem. Eng. Sci. 55, 4611–4627. Marchisio, D., Picturna, J., Fox, R., Vigil, R., Barresi, A., 2003. Quadrature method of moments for population balance equations. AICHE J. 49, 1266–1276. Mirshekari, F., Pakzad, L., Fatehi, P., 2019. An investigation on the stability of the hazelnut oil-water emulsion. J. Dispers. Sci. Technol. doi:10.1080/01932691.2019. 1614.459. Nambier, D., Kumar, R., Das, T., Gandhi, K., 1992. A new model for breakage frequency of drops in turbulent stirred dispersion. Chem. Eng. Sci. 47, 2989– 3002. Nastac, L., Zhang, L., Thomas, B.G., Sabau, A., El-Kaddah, N., Powel, A.C., Combeau, H., 2012. CFD Modeling and Simulation in Material Processing. Wiley, Florida. Ohtake, T., Hano, T., Takagi, K., 1987. Effects of viscosity on drop diameter of W/O emulsion dispersed in a stirred tank. J. Chem. Eng. Jpn. 20, 443–447. Ohtake, T., Hano, T., Takagi, K., 1988. Analysis of water entrainment into dispersed W/O emulsion drops. J. Chem. Eng. Jpn. 21, 272–276. Ok, T., Ookawara, S., Yoshikawa, S., Ogawa, K., 2003. Drop size distribution in liquid–liquid mixing. J. Chem. Eng. Jpn. 36, 940–945.

18

S.K. Naeeni and L. Pakzad / International Journal of Multiphase Flow 120 (2019) 103100

Pacek, A.W., Chamsart, S., Nienow, A.W., Bakker, A., 1999. The influence of impeller type on mean drop size and drop size distribution in an agitated vessel. Chem. Eng. Sci. 5, 4211422. Pakzad, L., Ein-Mozaffari, F., Chan, P., 2008. Using electrical resistance tomography and computational fluid dynamics modeling to study the formation of cavern in the mixing of pseudoplastic fluids possessing yield stress. Chem. Eng. Sci. 63, 2508–2522. Patel, D., Ein-Mozaffari, F., Mehrvar, M., 2015. Effect of rheological parameters on non-ideal flows in the continuous-flow mixing of biopolymer solutions. Chem. Eng. Res. Des. 100, 126–134. Qi, L., Meng, X., Zhang, R., Liu, H., Liu, Z., Klusener, P., 2015. Droplet size distribution and droplet size correlation of chloroaluminate ionic liquid-heptane dispersion in a stirred vessel. Chem. Eng. J. 268, 116–124. Ramkrishna, D., 1985. The status of population balances. Rev. Chem. Eng. 3, 49–95. Ramkrishna, D., 20 0 0. Population Balances: Theory and Application to Particulate Systems in Engineering. Academic Press, San Diego. USA. Rathore, A., Sharma, C., Persad, A., 2012. Use of computational fluid dynamics as a tool for establishing process design space for mixing in a bioreactor. AICHE J. 28, 382–391. Revankar, S., 2001. Coalescence and breakup of fluid particles in multi-phase flow. International Conference on Multiphase Flow. Roudsari, S.F., Turcotte, G., Dhib, R.R., Ein-Mozaffari, F., 2012. CFD modeling of the mixing of water in oil emulsions. Comput. Chem. Eng. 45, 124–136. Rueger, P.E., Calabrese, R.V., 2013. Dispersion of water into oil in a rotor-stator mixer. Part 2: effect of phase fraction. Chem. Eng. Res. Des. 91, 2134–2141. Ruiz, M., lermanda, P., Padilla, R., 2002. Drop size distribution in batch mixer under breakage conditions. Hydrometallurgy 63 (1), 65–74. Schiller, L., Naumann, Z., 1933. Uber die grundlegenden berechnungen bei der schwekraftaubereitung. Zeitschrift Des Vereins Deutscher Ingenieure 77, 318–320. Schütz, S., Gorbach, G., Piesche, M., 2009. Modeling fluid behavior and droplet interactions during liquid–liquid separation in hydrosyclones. Chem. Eng. Sci. 64, 3935–3952. Shinnar, R., 1961. On the behaviour of liquid dispersions in mixing vessels. J. Fluid Mech. 10 (2), 259–275. Simon, M., Schmidt, A., Bart, H., 2003. The droplet population balance model- estimation of breakage and coalescence. Chem. Eng. Technol. 26, 745–750. Sjoblom, J., 2001. Encyclopedic Handbook of Emulsion Technology. Marcel Dekker, New York. USA.

Sprow, F.B., 1967. Distribution of drop sizes produced in turbulent liquid–liquid dispersion. Chem. Eng. Sci. 22, 435–442. Srilatha, C., Morab, V., Mundada, T., Patwardhan, A., 2010. Relation between hydrodynamics and drop size distributions in pump-mix mixer. Chem. Eng. Sci. 65, 3409–3426. Stamatoudis, M., Tavalarides, L., 2007. The effect of continuous phase viscosity on the unsteady state behavior of liquid–liquid agitated dispersions. Chem. Eng. J. 35, 137–143. Sutherland, K., Pakzad, L., Fatehi, P., 2019. CFD population balance modeling and dimensionless group analysis of a multiphase oscillatory baffled column (OBC) using moving overset meshes. Chemical Engineering Science 199, 552– 570. Vankova, N., Tcholakova, S., Denkov, N.D., Ivanov, I.B., Vulchev, V.D., Danner, T., 2007. Emulsification in turbulent flow: 1. Mean and maximum drop diameters in inertial and viscous regimes. J. Colloid Interface Sci. 312, 363–380. Versteeg, H.K., Malalasekera, W., 2007. An introduction to computational fluid dynamics: the finite volume method. Pearson Education. Vladisavljevic, G., Kobayashi, I., Nakajima, M., 2011. Effect of dispersed phase viscosity on maximum droplet generation frequency in micro channel emulsification using asymmetric straight-through channels. Microfluid Nanofluid J. 10, 1199–1209. Vonka, M., Soos, M., 2015. Characterization of liquid–liquid dispersions with variable viscosity by coupled computational fluid dynamics and population balances. AIChE J. 61, 2403–2414. Wang, C.Y., Calabrese, R.V., 1986. Drop breakup in turbulent stirred-tank contactors. Part II: relative influence of viscosity and interfacial tension. AIChE J. 32, 667–676. Wang, W., Liu, J., Wang, P., Duan, J., Gong, J., 2013. Evolution of dispersed drops during the mixing of mineral oil and water phases in a stirred tank. Chemical Engineering Science 91, 173–179. Wutz, J., Lapin, A., Siebler, F., Schafer, J., Wucherpfenning, T., Berger, M., Takors, R., 2016. Predictability of KL a in stirred tank reactors under multiple operating conditions using an Euler–Lagrange approach. Eng. Life Sci. 16 (7), 633– 642. Yeoh, G.H., Cheung, C.P., Tu, J., 2014. Multiphase Flow Analysis Using Population Balance Modeling: Bubbles, Drops, Particles. Elsevier Ltd, UK. Zhou, G., Kresta, S.M., 1998. Correlation of mean drop size and minimum drop size with the turbulence energy dissipation and the flow in an agitated tank. Chem. Eng. Sci. 53, 2063–2079.