Chemical Engineering Journal 330 (2017) 245–261

Contents lists available at ScienceDirect

Chemical Engineering Journal

journal homepage: www.elsevier.com/locate/cej

CFD analysis of microfluidic droplet formation in non–Newtonian liquid

Somasekhara Goud Sontti, Arnab Atta ⇑
Multiscale Computational Fluid Dynamics (mCFD) Laboratory, Department of Chemical Engineering, Indian Institute of Technology Kharagpur, West Bengal 721302, India

h i g h l i g h t s g r a p h i c a l a b s t r a c t

Presented a 3D numerical
investigation of Newtonian
microfluidic droplet formation in
non–Newtonian liquids.

Identified different regimes, viz.,
squeezing, dripping and jetting in a
T-junction microchannel.

Analyzed the influence of various
controlling parameters on droplet
characteristics.

Revealed new insights on droplet
formation mechanism in non–
Newtonian liquids.

a r t i c l e i n f o a b s t r a c t

Article history: A three-dimensional, volume-of-fluid (VOF) based CFD model is presented to investigate droplet forma-
Received 22 February 2017 tion in a microfluidic T-junction. Genesis of Newtonian droplets in non–Newtonian liquid is numerically
Received in revised form 5 July 2017 studied and characterized in three different regimes, viz., squeezing, dripping and jetting. Various influ-
Accepted 15 July 2017
encing factors such as, continuous and dispersed phase flow rates, interfacial tension, and non–
Available online 18 July 2017
Newtonian rheological parameters are analyzed to understand droplet formation mechanism. Droplet
shape is reported by defining a deformation index. Near spherical droplets are realized in dripping and
Keywords:
jetting regimes. However, plug shaped droplets are observed in squeezing regime. It is found that rheo-
Non–Newtonian liquid
Droplet
logical parameters have significant effect on the droplet length, volume, and its formation regime. The
T–junction microchannel formation frequency increases with increasing effective viscosity however, the droplet volume decreases.
CFD This work effectively provides the fundamental insights into microfluidic droplet formation characteris-
Flow regime tics in non–Newtonian liquids.
Ó 2017 Elsevier B.V. All rights reserved.

1. Introduction become a significant area of research in microfluidics due to its

paramount importance in biomedical engineering. Droplets gener-
In recent years, droplet–based microfluidics offer a wide range ated in non–Newtonian medium are frequently encountered in
of applications in the fields of lab–on–a–chip, chemical, biological biochemical and drug delivery applications. However, precise
and nanomaterial synthesis [1–6]. In a microfluidic device, each control over the droplet size is required to deliver accurate dosing
droplet provides a compartment microreactor in which species of a drug or chemical reactant [16]. Therefore, monodisperse
transport or reactions can occur [7–14]. Recently, there has been droplets are highly desirable in several areas of microfluidics
a rapid development on emulsion generation in microfluidic [17]. Various microfluidic devices, such as, T–junction, flow-
devices [15]. Non–Newtonian liquid in multiphase system has also focusing, and co-flowing devices are available for generating
droplets in microchannel by shearing dispersed phase with a
continuous stream of liquid [18–22,15]. Microfluidic T–junction
⇑ Corresponding author. finds wide spread utilization compared to other devices owing to
E-mail address: arnab@che.iitkgp.ernet.in (A. Atta).

http://dx.doi.org/10.1016/j.cej.2017.07.097
1385-8947/Ó 2017 Elsevier B.V. All rights reserved.
246 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261

Nomenclature


Rew Reynolds number ¼ Dw Ulw qw Dx cell size (m)
w

 g viscosity (Pa.s)
0 KU n W 1n
Ca Modified Capillary number ¼ L c c
h contact angle (°)
U velocity (m/s)
q density (kg/m3)
v average velocity in a cell (m/s)
r surface tension (N/m)
W width of channel
c_ shear rate
H height of the droplet
jn radius of curvature
D.I droplet deformation index
geff effective viscosity ðPa:sÞ
t time (s)
s stress tensor
P pressure (Pa)
F SF volumetric surface tension force (N/m3)
K consistency coefficient ðPa:sn Þ Subscripts
Q volumetric flow rate ðlL=sÞ w water
D rate–of–deformation tensor c continuous phase
Q 0r flow rate ratio ðQ o =Q w Þ d dispersed phase
D droplet
o oil
Greek symbol
a volume fraction
Dt time step (s)

its simplicity in geometric configuration and superior control over squeezing to dripping mechanism. Fu et al. [45] also proposed scal-
droplet size. In a T–junction, cross–flowing continuous and dis- ing laws in terms of flow rate ratio and Ca for the estimation of
persed phase streams meet at the junction and the consequent bubble sizes in various regimes. Additional insights to the dis-
shear force leads to droplet formation. The droplet size and fre- persed and continuous phase pressure profiles in droplet breakup
quency are guided by continuous and dispersed phase flow rate process can be gained from the research of Sivasamy et al. [46].
ratio, as well as by adjusting viscous and interfacial tension forces. Effects of liquid viscosity and interfacial tension on droplet forma-
T–junction microchannel can be operated in three regimes namely, tion were investigated by Chen et al. [47] which showed that the
squeezing, dripping and jetting which are active functions of the period of slug formation increases with increasing interfacial ten-
primary driving forces acting on the system such as flow rate ratio, sion. Raj et al. [48] studied droplet formation in T-junction and
Capillary number (Ca), pressure gradient across the droplet, and Y–junction microchannels using volume-of-fluid (VOF) method
wetting properties of the channel surface. Thorsen et al. [23] ini- and analyzed the effect of flow rate ratio, liquid viscosity, interfa-
tially reported generation of water droplets in oil using a T–junc- cial tension, channel size, and wall adhesion properties on slug
tion microchannel. Nisisako et al. [24] experimentally length for Newtonian liquids. Fu et al. [34] investigated oil droplet
investigated the controlling parameters of droplet formation, size formation in a flow–focusing microchannel and proposed scaling
and frequency in an oil–water system. Numerous researchers also laws to predict the droplet length. Benefiting from the control over
proposed several strategies for droplet formation [25–28]. shape and size, few studies are also reported in T–junction device
Garstecki et al. [29] studied the droplet and bubble formation to produce spherical and plug shaped nanoparticles [6,49,50].
mechanism in a T-junction and proposed a power law correlation Interestingly, most of the reported research are concerned with
for predicting their sizes. Subsequently, several research works the droplet formation mechanism and flow regimes in Newtonian
were carried out to develop scaling laws in forecasting the droplet fluids, while in several applications, liquid phases are likely to
length for different inlet configuration and fluid properties [30– exhibit complex behaviors, such as non–Newtonian properties.
36]. With various surfactants and their concentrations, Xu et al. For example, it is apparent that in most cases of oil–water emul-
[37] showed that adjusting interfacial tension and the wetting sions, oil phase is usually considered as the Newtonian liquid dur-
properties could lead to ordered or disordered two-phase flow pat- ing modeling of such flow systems. However, studies on droplet
terns. Bashir et al. [38] numerically addressed similar issues and formation in non–Newtonian media is scarce. Abate et al. [51]
commented on active control of two–phase flow pattern by alter- studied monodisperse microparticle formation in non–Newtonian
ing interfacial tension and wetting properties. van der Graaf et al. polymer solutions in a flow–focusing device. Arratia et al. [52]
[39] computationally modeled droplet formation at the mesoscale reported polymeric filament thinning and breakup of Newtonian
using Lattice Boltzmann method (LBM). Wang et al. [33] and Riaud and viscoelastic liquids in a flow–focusing microchannel. The
et al. [40] investigated the droplet formation mechanism by LBM results showed different breakup mechanisms for Newtonian and
with simple modifications in a T–junction and shearing plates. polymeric liquids having same viscosity. This phenomena was
Yang et al. [41] also applied LBM to analyze the droplet formation attributed to the rheological difference between the two types of
and cell encapsulation process where three flow regimes were liquids. Qiu et al. [53] numerically investigated the droplet forma-
illustrated and the droplet shapes are reported in each regime, tion in non–Newtonian liquids in a cross–flow microchannel. It is
such as, plug shape in squeezing, and bullet shape in dripping as apparent from their findings that rheological parameters of non–
well as in jetting regimes. Wang et al. [42] discussed the Newtonian fluid significantly influence the formation mechanism
generation of monodisperse droplets using capillary embedded and size of droplets. Aytouna et al. [54] experimentally examined
T–junction device and described its dependence on Ca, viscosity the droplet pinch–off dynamics in Newtonian, yield stress, and
ratio, and dispersed phase flow rate. De Menech et al. [43] compu- shear thinning fluids. Fu et al. [55] studied the flow patterns at
tationally illustrated different flow regimes using phase-field the mesoscale in Newtonian and non–Newtonian fluids using
method. Christopher et al. [44] and Fu et al. [45] analyzed bubble T–junction millichannel. However, in-depth understanding of
formation in microfluidic T-junctions and its transition from droplet formation and flow regimes in non–Newtonian fluids is
 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261 247

still lacking at the microscale. Therefore, analyses considering the respectively. The interface is tracked by VOF method, which solves
liquid phase as non–Newtonian fluid is imperative and desirable. a single set of conservation equations for both phases, as follows
In this study, we develop a CFD based model to understand the [56]:
droplet formation in a T–junction microchannel considering the Equation of continuity:
continuous phase as non–Newtonian liquid. We describe the role
~ ¼0
r:ðqUÞ ð1Þ
of rheological properties, interfacial tension, and flow rate ratio
on droplet characteristics, in terms of formation mechanism, dro- Equation of motion:
plet length, volume, velocity, and its shape. These understandings
~
@ðqUÞ
can significantly benefit in setting guidelines to control droplet size ~UÞ
þ r:ðqU ~ ¼ rP þ r:s þ ~
F SF ð2Þ
and shape in non–Newtonian liquids. @t
~ P; q, and s are velocity vector, pressure, volume averaged
where U;
2. Problem formulation density, and stress tensor, respectively. For incompressible Newto-
nian fluids, the shear stress is proportional to the rate–of–strain
A three–dimensional T–junction microchannel, as shown in tensor ðc_ Þ, described by:
Fig. 1a, is considered to investigate the droplet formation in a ~ þ rU
~T Þ
non–Newtonian liquid which flows through the main channel hav-
s ¼ gc_ ¼ gðrU ð3Þ
ing a cross-section of 100 lm  33 lm. Newtonian dispersed where g is the volume-averaged viscosity. The volume-averaged
phase is introduced through a perpendicular channel of properties are defined in terms of oil ðao Þ and water ðaw Þ volume
50 lm  33 lm. At the merging junction of these two liquids, con- fractions, as follows:
tinuous phase shear force acts on the dispersed phase and the con-
sequent pressure gradient results in formation of droplets that flow q ¼ ao qo þ ð1  aw Þqw ð4Þ
through downstream of the channel (Fig. 1b). Typically, the droplet
formation characteristics are governed by the complex interactions g ¼ ao go þ ð1  aw Þgw ð5Þ
between the two phases resulting from forces like viscous, interfa- Equation of volume fraction:
cial, shear and pressure gradient in the channel. These forces are The volume fraction of each liquid phase is calculated by solving
substantially influenced by the liquid properties and flow rates. the following equation:
To quantify droplet shape, a deformation index ðD:IÞ is calculated
@ aq ~
as shown in Fig. 1b, which indicates an undeformed state at þ U:raq ¼ 0 ð6Þ
D:I ¼ 0 whereas, D:I > 0 suggests deformed droplets with length @t
greater than its height. where the subscript q refers to either oil (o) or water (w) phase. In
each computational cell, the volume fractions of all phases are con-
P
2.1. Governing equations served by aq ¼ 1. For aq ¼ 0, the reference cell is assumed to be
devoid of the qth phase, and aq ¼ 1 indicates that the cell is com-
2.1.1. Equations of continuity and momentum pletely filled with qth phase. Consequently, the interface between
In this work, incompressible two–phase (oil–water) flow is con- two phases is identified by marking the cell with volume fraction
sidered, where oil and water are continuous, and dispersed phases, range 0 < aq < 1.

2.1.2. Surface tension force

The continuum surface force (CSF) model [57] is used to define
the volumetric surface tension force ðF SF Þ term in Eq. (2), as
follows:
 
~ qjN rao
F SF ¼ r ð7Þ
1
2
ðqo þ qw Þ
where jN is the radius of curvature and r is the coefficient of sur-
face tension. The interface curvature ðjN Þ is calculated in terms of
b as:
unit normal N,
1 h
N~ 
i
b ¼ N
~
jN ¼ r: Nb ¼ ~ ~ :
jNj
~  r:N
r jNj ~ ; where N
~
ð8Þ
jNj jNj
In VOF formulation, surface normal, N, is expressed as the gra-
dient of phase volume fraction at the interface which can be writ-
ten as:
~ ¼ raq
N ð9Þ
This surface tension force is implemented by the piecewise-
linear interface calculation (PLIC) scheme that provides accurate
calculation of curvatures for reconstruction of the interface front
[58,59]. Wall adhesion effect is also taken into consideration by
defining a three-phase contact angle at the channel wall ðhW Þ.
Accordingly, the surface normal at the reference cell next to the
wall is given by:

Fig. 1. (a) Considered 3D T–junction geometry, and (b) 2D schematic of droplet b ¼N

N b W sinhW
b W coshW þ M ð10Þ
formation in a T–junction microchannel.
248 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261

where N b W and Mb W are the unit vectors normal and tangential to the where K and n are the consistency and power-law indices, respec-
wall, respectively [58]. In this work, a static contact angle condition tively. The local shear rate ðc_ Þ is related to the second invariant of
is specified which is assumed to be independent of the moving con- D and is expressed as [58]:
tact line and the velocity [48]. The surface normal one cell away rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
from the wall along with the contact angle govern the local curva- 1 ~ ~ T Þ ðrU ~ þ rU ~T Þ
c_ ¼ ðrU þ rU ij ji ð13Þ
ture of the surface, which is utilized in the calculation of the surface 2
tension force (Eq. (7)) to determine the body force term in Eq. (2).
For real surfaces, the contact angle varies dynamically between an
2.2. Implementation in numerical model
advancing and a receding contact angle. If the contact angle remains
within the range of advancing and receding angles, the contact line
Aforementioned time dependent governing equations are
does not move [60]. Typically, the use of static contact angle has
solved in a CFD solver (Ansys Fluent 17.0) based on finite volume
been proved to be adequate for analyzing the flow behavior in
method. Pressure implicit with splitting of operators (PISO) algo-
microchannels [33,39,61–63]. Nonetheless, the cases of dynamic
rithm [66] is used to resolve the pressure–velocity coupling in
contact angles can be resolved by defining a level set function
momentum equation. The spatial derivatives are discretized using
[64] and coupled with the VOF model. However, such investigation
quadratic upstream interpolation for convective kinetics (QUICK)
is beyond the scope of this present work and can be found else-
scheme [67]. To avoid spurious currents as a result of mismatch
where [60,65].
between pressure and surface tension force discretization, the
pressure staggering option (PRESTO) is employed for pressure
2.1.3. Constitutive equation of continuous phase
interpolation [58,59]. The geometric reconstruction scheme is
For non–Newtonian liquids, the shear stress can be written in
adopted to solve the volume fraction equation. First-order implicit
terms of a non–Newtonian viscosity:
method is applied for the discretization of temporal derivatives.
s ¼ gðc_ Þc_ ð11Þ Subsequently, variable time step and fixed Courant number
ðCo ¼ 0:25Þ are considered for simulating the governing equations.
where g is a function of all three invariants of the rate–of–deforma
Constant velocity for both continuous and dispersed phases are
tion tensor. However, in power–law model, the non–Newtonian liq-
imposed at the inlets and outflow boundary condition is specified
uid viscosity ðgÞ is considered to be a function of only shear rate ðc_ Þ. at the outlet. It is assumed that continuous phase completely wets
gðc_ Þ ¼ K c_ n1 ð12Þ the channel wall and the solid walls are set to no–slip boundary
condition with a static contact angle of 1350 . Structured hexahe-

Fig. 2. Comparison of the model predictions with the experimental results of Garstecki et al. [29](figures reprinted with permission from the publisher, Royal Society of
Chemistry) for (a) Q w ¼ 0:14 lL=s and Q o ¼ 0:124 lL=s, (b) Q w ¼ 0:14 lL=s and Q o ¼ 0:408 lL=s, (c) Q w ¼ 0:004 lL=s and Q o ¼ 0:028 lL=s, (d) Q w ¼ 0:006 lL=s and
Q o ¼ 0:028 lL=s. Comparison of droplet length with (e) experimental results of Garstecki et al. [29] and numerical predictions by Raj et al. [48] at
Q w ¼ 0:14 lL=s; r ¼ 0:0365 N=m; go ¼ 0:01 Pa:s; gw ¼ 0:001 Pa:s, and (f) experimental as well as LBM results of Graaf et al. [39] at Q w ¼ 0:055 lL=s, r = 5 mN/m,
go ¼ 6:71 mPa:s and gw ¼ 1:95 mPa:s.
 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261 249

dral meshes with Green–Gauss node based schemes [68] are used Raj et al. [48]. The developed model is further verified with the
in the computational domain to accurately calculate the gradients experimental results of van der Graaf et al. [39] for asserting its
and to overcome calculation inaccuracies resulting from spurious accuracy in droplet length prediction. In this case, oil droplet for-
currents, as recommended by Gupta et al. [69]. Moreover, the mation is simulated by changing the continuous phase (water)
effect of mesh element size is initially investigated for analyzing flow rate in a T–junction microchannel and a maximum deviation
numerical diffusion that mainly stems from poor meshing scheme. of 6% from the experimental results is observed, as shown in Fig. 2f.
Thereafter, the model is validated for the droplet length estimation This validation also establishes the efficacy of the developed model
with the reported experimental data of Garstecki et al. [29]. It has to forecast the droplet length better than the LBM simulations,
been realized that for more than 2.5 times of mesh elements (com- reported by van der Graaf et al. [39].
paring between 137,176 and 360,640 elements), the difference in
droplet length is nearly 1%. Consequently, for the considered 3. Results and discussion
geometry all simulations are performed with element size of
4 lm (137,176 elements) to optimize the computational time. Armed with fairly validated model, the droplet formation mech-
anism and its behavior in non–Newtonian liquids are systemati-
2.3. Model validation cally investigated for different operating conditions by varying
the continuous and dispersed phase velocities (Q o and Q w , respec-
Fig. 2 shows the comparison of model predictions with the tively). The influence of various rheological parameters, namely,
experimental observations of Garstecki et al. [29]. At a constant power–law index (n), consistency index (K), and, interfacial tension
water flow rate (Q w ¼ 0:14 lL=s), the droplet size variation with ðrÞ are elaborated in the following sections.
increasing oil flow rate (Q o ¼ 0:124—0:408 lL=s) is illustrated in
Fig. 2a and b. Similarly, Fig. 2c and d describe the effect of increas- 3.1. Effect of power–law index
ing dispersed phase (water) flow rate (Q w ¼ 0:004—0:006 lL=sÞ for
Q o ¼ 0:028 lL=s. Fig. 2e portrays the quantitative comparison of Fig. 3 portrays the temporal evolution of droplet formation phe-
droplet length estimation with the experimental observation by nomena in different power–law liquids, obtained by varying n from
Garstecki et al. [29] which shows the maximum error of 18% at 0.80 to 1.10. It can be understood from Fig. 3 that at a fixed oper-
the highest continuous phase flow rate. However, the model ating condition, droplet formation time decreases with increasing
results are found to be identical with the numerical prediction by power–law index. This is attributed to the interplay between vis-

Fig. 3. Droplet formation mechanism for Newtonian and power-law liquids at a fixed operating condition of Q o ¼ 0:408 lL=s; Q w ¼ 0:14 lL=s, K ¼ 0:01 Pa:sn ;
gw ¼ 0:001 Pa:s and r = 0.0365 N/m.
250 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261

cous and interfacial forces. Although the droplet formation process continuous phase and consequently, pressure gradient in continu-
is similar for n ¼ 1:0 and n ¼ 1:05 (mainly due to the small change ous phase across the droplet increases. Eventually, the formation of
in power-law index) till 0.0014 s, it shows a significant change at plug shaped droplets are formed once the gradient is sufficiently
0.003 s where three droplets are formed for n ¼ 1:05. With increas- large to overcome the pressure inside the dispersed phase droplet.
ing n, the effective viscosity ðgeff Þ of the continuous phase (oil) is Dripping regime of droplet formation is observed for Newtonian
enhanced, which in turn, imparts higher viscous resistance and and shear thickening liquids ðn  1Þ, where the droplet pinch off
helps in rapid detachment of droplets. In the cases of shear thin- occurs at the merging junction, as depicted in Fig. 4b. In the drip-
ning liquids ðn < 1Þ, interfacial force plays a dominant role and ping regime, viscous force that acts on the interface to snap off dis-
the detachment of dispersed phase at the merging junction is persed phase, typically dominates over the interfacial tension force
delayed. [15]. As a consequence of larger viscous and shear forces in Newto-
Moreover, for shear thinning liquid, squeezing regime is nian and shear thickening liquids with increasing n, Fig. 3 shows
observed, where the dispersed phase grows slowly and covers rapid formation of droplets before its enlargement up to the chan-
the entire flow area (Fig. 4a). This results in restricted flow of the nel top wall. It can be further noted that as the flow regime shifts

Fig. 4. Droplet formation mechanism (a) squeezing regime (n = 0.80) and (b) dripping regime (n = 1.10) at a fixed operating condition of Q o ¼ 0:408 lL=s; Q w ¼ 0:14 lL=s,
K ¼ 0:01 Pa:sn ; gw ¼ 0:001 Pa.s and r = 0.0365 N/m.

Fig. 5. Effect of power–law index on (a) non-dimensional droplet length (b) droplet velocity and volume, (c) deformation index ðD:IÞ, and (d) pressure profiles along the
channel centerline at fixed K ¼ 0:01 Pa:sn ; gw ¼ 0:001 Pa:s; r ¼ 0:0365 N=m; Q o ¼ 0:408 lL=s, and Q w ¼ 0:14 lL=s.
 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261 251

from squeezing to dripping, droplet shape also changes from plug To realize the effect of power-law index in terms of bulk viscos-
to spherical. ity and overall viscous force, effective viscosities of various liquids
Fig. 5a shows that non-dimensional droplet length ðLD =W c Þ are estimated from the simulations and compared with the results
decreases with increasing power –law index (n), due to increase derived from Eq. (14) [71].
in effective viscosity of continuous phase liquid. It is also apparent
 n  n1
that in shear thinning liquids, elongated droplets are formed 3n þ 1 8U L
ðLD =W c > 1Þ, whereas smaller droplets are formed for
geff ¼ K ð14Þ
4n Wc
ðLD =W c < 1Þ Newtonian and shear thickening liquids. Accordingly,
droplet volume also decreases with increasing n, as shown in where K; U L ; W c , and n are the consistency index, liquid inlet veloc-
Fig. 5b. However, the droplet velocity is found to increase with ity, width of the channel, and power–law index, respectively. Fig. 6a
increasing n. This can be ascribed to the change in film thickness shows enhanced effective viscosity with increasing power–law
and flow profile that occur in continuous phase from shear thin- index and the CFD calculations are in excellent agreement with
ning to thickening nature. Typically, the velocity profile is sharper the theoretical values. Effective viscosity distributions in the middle
for shear thickening liquids compared to a plug flow profile for of the liquid slugs are also presented in Fig. 6b. In the middle of the
shear thinning cases [70]. Based on length and height, droplet microchannel, the effective viscosity increases for shear thinning
shape is also quantified in terms of droplet deformation index liquids and decreases for shear thickening liquids which are the
ðD:IÞ, as illustrated in Fig. 5c. Plug shaped droplets are identified typical representation of non–Newtonian flow behavior. Velocity
for n < 1 and almost spherical shaped droplets are observed for profiles in the middle of the slug and the droplet are also illustrated
n = 1 when all the other parameters such as K, interfacial tension, in Fig. 6c and d, respectively. In all cases, flatter velocity profiles are
and flow rates were kept constant. For n > 1, small deformation is observed for shear thinning liquids, as expected. Due to enhanced
realized from spherical shape due to higher viscous force. The pres- film thickness near the wall and shape of the droplet, the velocity
sure evolution along the channel length for different power–law inside the droplet for shear thickening liquids is considerably higher
liquids is described in Fig. 5d. It can be seen from Fig. 5d that pres- than the shear thinning liquids.
sure drop increases with increasing power–law index and each To understand the rheological behavior of non–Newtonian liq-
peak represents the droplet position (encircled by dotted red line) uids on the droplet formation, a set of simulations are performed
where there is a pressure difference before and after droplet by adjusting the value of K which results in maintaining the iden-
detachment. tical effective viscosities for both shear thinning (n = 0.80) and

Fig. 6. (a) Effect of power–law index on bulk liquid viscosity. Distribution profiles of (b) effective viscosity, and (c) velocity in middle of the liquid slug. (d) Velocity profiles in
the middle of droplet at fixed K ¼ 0:01 Pa:sn ; gw ¼ 0:001 Pa:s; r ¼ 0:0365 N=m; Q o ¼ 0:408 lL=s, and Q w ¼ 0:14lL=s.
252 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261

Fig. 7. Distribution profiles of (a) effective viscosity, and (b) velocity in middle of liquid slug in shear thinning (n = 0.80 and K ¼ 0:01 Pa:sn ) and shear thickening liquids
(n = 1.10 and K ¼ 0:00066 Pa:sn ) at fixed effective viscosity geff ¼ 0:001669 Pa:sn ; r ¼ 0:0365 N=m; Q o ¼ 0:408 lL=s, and Q w ¼ 0:14 lL=s.

Fig. 8. Influence of K during droplet formation in (a) shear thinning liquid (n = 0.80), (b) shear thickening liquid (n = 1.10) with K 1 ¼ 0:008 Pa:sn ; K 2 ¼
0:01 Pa:sn ; K 3 ¼ 0:012 Pa:sn ; K 4 ¼ 0:014 Pa:sn ; K 5 ¼ 0:018 Pa:sn . Droplet formation in a shear thinning liquid (n = 0.80) during (c) squeezing mechanism with
K ¼ 0:018 Pa:sn , and (d) dripping mechanism with K ¼ 0:16 Pa:sn at a fixed operating condition of Q o ¼ 0:408 lL=s; Q w ¼ 0:14 lL=s; gw ¼ 0:001 Pa:s and r ¼ 0:0365 N=m.
 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261 253

shear thickening (n = 1.10) liquids. The effective viscosity distribu- phase is pulled downstream of the main channel before the dro-
tion profiles in the middle of the liquid slug is shown in Fig. 7a. plets breakup due to dominant interfacial tension forces, as shown
Similar to the earlier observation, typical shear dependent viscosity in Fig. 8a at K ¼ 0:008Pa:sn . For n < 1, interfacial and shear forces
profiles are also apparent for both the liquids in this case. Dimen- together provide squeezing action due to lower viscous forces. A
sionless droplet length for shear thinning liquid ðLD =W C ¼ 1:860Þ is long and thick layer of dispersed phase leads to the generation of
found to be higher than shear thickening liquid ðLD =W C ¼ 1:575Þ elongated droplets for n = 0.80 at K = 0.008–0.014 Pa:sn . However,
and the corresponding phase contours are shown in the inset of further increase in consistency index value for n = 0.80 and 0.90
Fig. 7b. The velocity profiles illustrated in Fig. 7b also depict the results in flow regime transition from squeezing to dripping. Smal-
combined influence of n and K even when the effective viscosities ler droplets are observed for n > 1 and flow regime shifts from
are same for both the liquids. squeezing to dripping due to the increase in effective viscosity,
as discussed earlier. In a power law liquid, the droplet formation
3.2. Effect of consistency index through squeezing as well as dripping mechanism can be realized
by varying the consistency index value, as shown in Fig. 8c and d
In this section, effect of consistency index (K) on droplet forma- for n ¼ 0:80.
tion mechanism, length, and velocity has been methodically The droplet length is found to decrease with increasing K as
explored by altering K = 0.008–0.018 Pa:sn of the power–law liq- shown in Fig. 9a, due to increase in effective viscosity of the contin-
uids. Fig. 8 shows the droplet formation mechanism in squeezing uous phase liquid. For higher K values, droplet length hardly
and dripping regimes. For the considered range of K, squeezing changes for shear thickening liquids. From Fig. 8b, it can be noticed
regime is experienced for n < 1, and dripping mechanism is real- that under such scenario, droplet shape gradually changes toward
ized for nP1, as depicted in Fig. 8a and b, respectively. Typically, nearly spherical (see for K = 0.01 Pa:sn ) with reduced core diameter.
the forces acting on the dispersed phase at the merging junction Beyond this value of K, droplets are observed to fuse, attributing to
are viscous, pressure difference, and interfacial tension forces. Vis- extremely high viscous stress. Droplet velocity is also estimated to
cous force is caused by the viscous stress acting on a liquid–liquid understand the effect of K on droplet dynamics.
interface and is proportional to the dispersed phase area along Fig. 9b depicts the gradual increase in droplet velocity with
with the velocity gradient. For lower values of K, the dispersed increasing K which is more pronounced in shear thickening liquids.

Fig. 9. Effect of consistency index on droplet (a) length, (b) velocity (c) volume, and (d) deformation index at Q o ¼ 0:408 lL=s; Q w ¼ 0:14 lL=s; gw ¼ 0:001 Pa:s, and
r ¼ 0:0365 N=m.
254 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261

This is mainly due to increase in liquid film thickness around the 3.3. Effect of interfacial tension
droplets and reduction in droplet height (H in Fig. 1). Conse-
quently, the droplet velocity marginally changes in shear thinning To understand the effect of interfacial tension on the droplet
cases. However, a noticeable change in droplet velocity is observed length and velocity, a set of simulation is carried out keeping other
for shear thickening and Newtonian liquids. Droplet volume is properties constant. All the results are analyzed based on modified


found to decrease with increase in K, as shown in Fig. 9c. For higher KU n w1n
Capillary number Ca0 ¼ L r c . It can be seen from Fig. 11a that
K values, droplet volume remains unchanged for shear thickening
for shear thinning liquid, flow regime shifts from dripping to jet-
liquids. Deformation in shape is further quantified in Fig. 9d which
ting with increasing interfacial tension. At higher interfacial ten-
shows that, at lower consistency index (K = 0.008–0.010 Pa:sn ),
near spherical shaped droplets are observed for Newtonian liquids. sion ðr ¼ 0:0565 N=m and ðCa0 ¼ 0:005ÞÞ, a thick layer of
For shear thinning liquid, plug shaped droplets are observed and dispersed phase covers the channel length and droplet formation
the deformation decreases with increase in K. In the range of this is not observed. Interestingly, for shear thickening (n = 1.10) liquid,
study, it is observed that the droplet deformation index (D.I) hardly the flow regime shifts from small beads that are linked by a fila-
varies for Newtonian and shear thickening liquids. However, there ment of the dispersed phase to dripping with increase in interfacial
is a dependence of D.I in cases of shear thinning liquids which tension and is illustrated in Fig. 11b. Similar observation was
becomes stronger with lower effective viscosity. This can be reported by Arratia et al. [52] in their experimental findings. This
attributed to the lower shear force for droplet detachment from can be attributed to the fact that at lower interfacial tension, the
the neck of the junction in cases of shear thinning liquids. The growth of dispersed phase at the merging junction is hindered by
non-dimensional size of the formed droplets ðLD =W c Þ are scaled the higher shear force, which in turn results into smaller droplet
length for all cases.
with the modified Capillary number ðCa0 Þ as a power–law
Fig. 12a indicates the change in non-dimensional droplet length
relationship and are illustrated in Fig. 10. The proposed scaling
laws for Newtonian and non–Newtonian liquids in the range of with varying Ca0 . It can be seen that with decreasing interfacial ten-
0:008 6 K 6 0:018 Pa:sn are summarized in Table 1 which sion, the droplet length and height also decrease in all cases. Sim-
have a maximum deviation of 5%. It can be observed that the ple power law relations for various n are proposed for the droplet
proposed scaling laws to predict droplet size are similar to those length as a function of Ca0 which are shown in Fig. 12a. The appli-
suggested for the droplet/bubble formation in Newtonian and cability range of these relations are listed in Table 2 which show a
non–Newtonian liquids by several researchers [30,35,36,72] in maximum error of 8%. Consequently, considerable change in dro-
various microfluidic devices, but with different prefactors and plet velocity is observed, as shown in Fig. 12b. As droplet blocks
exponents. the entire channel cross-section with increasing interfacial tension,
its velocity decreases and this phenomena is pronounced in shear
thickening cases unlike shear thinning liquid where the droplet
height (H) remains constant. In line with earlier discussions, dro-
plet volume is found to decrease with decreasing interfacial ten-
sion (Fig. 12c). Pressure drop across the droplet increases with
increasing interfacial tension owing to change in curvature of the
dispersed phase, as shown in Fig. 12d. However, the overall pres-
sure drop across the channel length is constant.

3.4. Effect of flow rate

3.4.1. Continuous phase

At a fixed operating condition of K = 0.01 Pa:sn ; r ¼ 0:0365 N/m,
and Q w ¼ 0:14 lL=s, the effect of continuous phase flow rate on
droplet formation mechanism, length and deformation index is
investigated. For a range of continuous flow rates from
Q o ¼ 0:297 lL=s to Q o ¼ 0:693 lL=s, only squeezing regime is
observed for shear thinning liquid (Fig. 13a and b). However, tran-
sition from squeezing to dripping regime is observed for Newto-
nian liquid, as described in Fig. 13c. From Fig. 13d and e, it is
apparent that in shear thickening liquid, dripping mechanism is
dominant with increase in continuous phase flow rate. At lower
Fig. 10. The scaling relation for the non-dimensional droplet size with the modified flow rates of shear thinning continuous phase, the dispersed phase
Capillary number ðCa0 Þ for various power-law liquids at Q o ¼ 0:408 lL=s; Q w ¼ pushes easily into the main channel due to lower resistance in con-
0:14 lL=s; gw ¼ 0:001 Pa:s, and r ¼ 0:0365 N=m. tinuous phase and completely blocks the cross section leading to

Table 1
Scaling of the droplet size with consistency index for various power-law liquids.

Power-law index (n) Consistency index K (Pa:sn ) Flow rates and fluid properties Scaling Laws Modified Capillary number (Ca0 )
0.80 0 0:36 0.006–0.0146
LD =W C ¼ 0:311Ca
0.90 Q o ¼ 0:408 lL=s, LD =W C ¼ 0:250Ca0
0:38 0.0132–0.0297
1.0 0.008–0.018 Q w ¼ 0:14 lL=s and LD =W C ¼ 0:462Ca 0 0:18 0.0269–0.0606
1.05 r ¼ 0:0365 N=m LD =W C ¼ 0:270Ca0
0:35 0.03847–0.0865
1.10 0:36 0.0549–0.1235
LD =W C ¼ 0:266Ca0
 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261 255

Fig. 11. Influence of interfacial tension on droplet formation for (a) shear thinning liquid ðn ¼ 0:80Þ having Ca0 ¼ 0:017 ðr ¼ 0:0165 N=mÞ, Ca0 ¼ 0:011 ðr ¼ 0:0265 N=mÞ,
Ca0 ¼ 0:008 ðr ¼ 0:0365 N=mÞ, Ca0 ¼ 0:006 ðr ¼ 0:0465 N=mÞ, and Ca0 ¼ 0:005 ðr ¼ 0:0565 N=mÞ, and (b) shear thickening liquid ðn ¼ 1:10Þ having Ca0 ¼ 0:151
ðr ¼ 0:0165 N=mÞ, Ca0 ¼ 0:095 ðr ¼ 0:0265 N=mÞ, Ca0 ¼ 0:068 ðr ¼ 0:0365 N=mÞ, Ca0 ¼ 0:054 ðr ¼ 0:0465 N=mÞ and Ca0 ¼ 0:044 ðr ¼ 0:0565 N=mÞ at K ¼ 0:01 Pa:sn ;
gw ¼ 0:001 Pa:s; Q o ¼ 0:408 lL=s, and Q w ¼ 0:14 lL=s.

Fig. 12. Effect of interfacial tension on droplet (a) length, (b) velocity, (c) volume and (d) pressure profiles along the channel length at
Q o ¼ 0:408 lL=s; Q w ¼ 0:14 lL=s; K ¼ 0:01 Pa:sn , and gw ¼ 0:001 Pa:s in various power-law liquids.
256 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261

Table 2
Scaling of the droplet size with interfacial tension for various power-law liquids.

Power-law index (n) Interfacial tension r (N/m) Flow rates and fluid properties Scaling Laws Modified Capillary number ðCa0 Þ
0.80 0 0:34 0.0052–0.01796
LD =W C ¼ 0:354Ca
0.90 Q o ¼ 0:408 lL=s, LD =W C ¼ 0:236Ca0
0:37 0.0106–0.0365
1.0 0.0165–0.0565 Q w ¼ 0:14 lL=s and LD =W C ¼ 0:431Ca 0 0:21 0.0217–0.0745
1.05 K ¼ 0:01 Pa:sn LD =W C ¼ 0:434Ca0
0:20 0.0310–0.1063
1.10 0:28 0.0443–0.1518
LD =W C ¼ 0:33Ca0

Fig. 13. Effect of continuous phase to dispersed phase flow rate ratio ðQ 0r Þ on droplet formation in (a) n = 0.80, (b) n = 0.90, (c) n = 1.0, (d) n = 1.10, and (e) n = 1.10 for different
flow rate ratios of Q 0r1 ¼ 0:212 lL=s; Q 0r2 ¼ 2:915 lL=s; Q 0r3 ¼ 3:536 lL=s; Q 0r4 ¼ 4:243 lL=s and Q 0r5 ¼ 4:95 lL=s at a fixed operating condition of
K ¼ 0:01 Pa:sn ; gw ¼ 0:001 Pa:s, r ¼ 0:0365 N=m and Q w ¼ 0:14 lL=s.

squeezing flow regime. The dispersed phase thread gradually However, with increasing n and continuous phase flow rate,
decreases with increasing shear stress and the detachment occurs increased viscous and inertia forces result in rapid detachment of
on attainment of critical thickness at the neck. the droplet at the neck, as evident from some cases of Newtonian
 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261 257

and all cases of shear thickening liquid. For shear thickening liquid, flow rate. The proposed relations are listed in Table 3 for
greater resistance is imparted by the continuous phase liquid caus- 0.297 6 Q o 6 0:693 lL=s and predicts with a maximum deviation
ing the dispersed phase to grow slowly into the main channel until of 3%.
it balances all forces at the junction and eventually leads to
reduction in droplet length (Fig. 14a). Similar to the previous 3.4.2. Dispersed phase
observation, near spherical and plug shaped droplets are observed Keeping the continuous phase flow rate constant at
in dripping and squeezing regimes, respectively (Fig. 14b). Fig. 14c Q o ¼ 0:408 lL=s, dispersed phase flow rate ðQ w Þ is varied from
shows power law relationship between the non-dimensional dro- 0.14 lL/s to 0.539 lL/s to understand its effect on droplet forma-
plet size and Ca0 resulting from the variation of continuous phase tion characteristics.

Fig. 14. Effect of continuous phase flow rate on (a) non-dimensional droplet length, and (b) droplet deformation index at a fixed operating conditions of
K ¼ 0:01 Pa:sn ; r ¼ 0:0365 N=m, gw ¼ 0:001 Pa:s and Q w ¼ 0:14 lL=s. (c) The scaling of non-dimensional droplet length with Ca0 for different continuous phase flow rates
in power-law liquids.

Table 3
Scaling of the droplet size with continuous phase flow rate for various power-law liquids.

Power-law index (n) Flow rates Fluid properties Scaling Laws Modified Capillary number ðCa0 Þ
0.80 0:17 0.001–0.0124
LD =W C ¼ 0:692Ca0
0.90 0:16
LD =W C ¼ 0:537Ca0 0.0015–0.0267
1.0 Q o ¼ 0:297  0:693 lL=s r ¼ 0:0365 N=m LD =W C ¼ 0:561Ca0
0:10 0.0024–0.0575
1.05 Q w ¼ 0:14 lL=s K ¼ 0:01Pa:sn LD =W C ¼ 0:502Ca0
0:11 0.003–0.0843
1.10 0 0:15 0.0038–0.1236
LD =W C ¼ 0:388Ca
258 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261

Fig. 15. Effect of dispersed phase flow rate ðQ w Þ on droplet formation in (a) n = 0.80 (b) n = 0.80 (c) n = 1.05 and (d) n = 1.10 for different dispersed phase Reynolds numbers:
Rew1 =4.24, Rew2 =6.30, Rew3 =8.45, Rew4 =10.3 and Rew5 =16.36 at K ¼ 0:01 Pa:sn ; gw ¼ 0:001 Pa:s, r ¼ 0:0365 N=m and Q o ¼ 0:408 lL=s.

The results are explained in terms of dispersed phase Reynolds phase Reynolds number ðRew Þ as shown in Fig. 16c. The range of


number Rew ¼ Dw Ug w qw . From Fig. 15, it can be observed that for applicability is provided in Table 4.
w

shear thinning cases, squeezing regime prevails with increasing

Q w (Figs. 15a and 15b). However, for Newtonian liquid, the regime 4. Conclusion
shifts from dripping to squeezing with increasing Q w due to
increased amount of dispersed phase being pushed into the chan- A comprehensive computational study of Newtonian droplet
nel at a fixed operating condition (Fig. 15c). In case of shear thick- formation in non–Newtonian power-law liquids is carried out in
ening liquid, dripping regime is observed which moves toward a T-junction microchannel using VOF method. New insights are
jetting regime with increased Q w (Figs. 15d and 15e) due to the obtained regarding the droplet formation process in non–Newto-
combined effect of higher inertia force (from dispersed phase) nian liquids. With increasing power-law and consistency index,
and viscous stress (in continuous phase). Consequently, the droplet droplet length is found to decrease as effective viscosity increases.
length increases significantly in shear thinning liquid with increas- The droplet length also decreases with increasing continuous
ing dispersed phase velocity (Fig. 16a) as the inertia force is signif- phase flow rate. Squeezing, dripping, and jetting regimes are found
icant enough to resist the opposing continuous phase shear stress. to strongly depend on the flow rate ratio, interfacial tension, and
However, in Newtonian and shear thickening liquid, that inertia rheological properties. A parameter defined as droplet deformation
effect is suppressed by continuous phase shear and viscous stresses index to identify the shape of droplets, shows that near spherical
resulting in small change in droplet length. Due to increase in dro- droplets are typically formed in all cases of dripping and jetting
plet volume with increasing Q w , the droplet shape also changes regimes. However, plug shaped droplets are obtained in squeezing
from near spherical to plug type, as indicated by droplet deforma- regime. Like Newtonian medium, it is observed that interfacial ten-
tion index in Fig. 16b. For various n, linear scaling relationships are sion has significant influence on the droplet formation pattern and
proposed between non-dimensional droplet size and the dispersed size in non–Newtonian liquids. With increasing interfacial tension,
 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261 259

Fig. 16. Effect of dispersed phase Reynolds number on droplet (a) length and (b) droplet formation index ðD:IÞ at K ¼ 0:01 Pa:sn ; gw ¼ 0:001 Pa:s, r ¼ 0:0365 N=m and
Q o ¼ 0:408 lL=s. (c) The scaling of non-dimensional droplet length with Rew for different dispersed phase flow rates in power-law liquids.

Table 4
Scaling of the droplet size with dispersed phase flow rate for various power-law liquids.

Power-law Flow rates Fluid properties Scaling Laws Reynolds number ðRew Þ
index (n)
0.80 LD =W C ¼ 0:152Rew þ 1:11
0.90 LD =W C ¼ 0:174Rew þ 0:202
1.0 Q w ¼ 0:14  0:539 lL=s r ¼ 0:0365 N/m LD =W C ¼ 0:085Rew þ 0:47 4.24–16.36
1.05 Q o ¼ 0:408 lL=s K ¼ 0:01Pa:sn LD =W C ¼ 0:065Rew þ 0:54
1.10 LD =W C ¼ 0:053Rew þ 0:46

droplet size increases in all cases, however, the regime shifts from Acknowledgments
filament linked small beads to dripping in shear thickening case.
The development of microfluidic methods to generate and This work is supported by Sponsored Research & Industrial Con-
manipulate monodisperse droplets have led to an increasing sultancy (SRIC), IIT Kharagpur under the scheme for ISIRD (Code:
number of potentially interesting applications. Although the com- FCF).
putational study can not completely diminish the necessity of
exhaustive and expensive, at times, experimental investigation; a
References
well-validated CFD model can certainly complement with various
aspects of physical phenomena that are not attainable by experi- [1] S. Li, J. Xu, Y. Wang, G. Luo, Controllable preparation of nanoparticles by drops
ments. Our findings are expected to provide better understanding and plugs flow in a microchannel device, Langmuir 24 (2008) 4194–4199.
and experimental guidelines in terms of controlling parameters [2] H. Breisig, J. Hoppe, T. Melin, M. Wessling, On the droplet formation in hollow-
fiber emulsification, J. Membr. Sci. 467 (2014) 109–115.
on the formation of desired Newtonian droplet of different shape [3] A. Knauer, A. Thete, S. Li, H. Romanus, A. Csaki, W. Fritzsche, J. Köhler, Au/ag/au
and size in a non–Newtonian medium. double shell nanoparticles with narrow size distribution obtained by
260 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261

continuous micro segmented flow synthesis, Chem. Eng. J. 166 (2011) 1164– [36] T. Fu, Y. Ma, H.Z. Li, Breakup dynamics of slender droplet formation in shear-
1169. thinning fluids in flow-focusing devices, Chem. Eng. Sci. 144 (2016) 75–86.
[4] J. Xu, J. Tan, S. Li, G. Luo, Enhancement of mass transfer performance of liquid– [37] J. Xu, S. Li, J. Tan, Y. Wang, G. Luo, Controllable preparation of monodisperse o/
liquid system by droplet flow in microchannels, Chem. Eng. J. 141 (2008) 242– w and w/o emulsions in the same microfluidic device, Langmuir 22 (2006)
249. 7943–7946.
[5] E.Y. Basova, F. Foret, Droplet microfluidics in (bio) chemical analysis, Analyst [38] S. Bashir, J.M. Rees, W.B. Zimmerman, Simulations of microfluidic droplet
140 (2015) 22–38. formation using the two-phase level set method, Chem. Eng. Sci. 66 (2011)
[6] L. Xu, J. Peng, M. Yan, D. Zhang, A.Q. Shen, Droplet synthesis of silver 4733–4741.
nanoparticles by a microfluidic device, Chem. Eng. Process. Process Intensif. [39] S. van der Graaf, T. Nisisako, C. Schroen, R. Van Der Sman, R. Boom, Lattice
102 (2016) 186–193. boltzmann simulations of droplet formation in a t-shaped microchannel,
[7] T. Nisisako, T. Torii, T. Higuchi, Novel microreactors for functional polymer Langmuir 22 (2006) 4144–4152.
beads, Chem. Eng. J. 101 (2004) 23–29. [40] A. Riaud, K. Wang, G. Luo, A combined lattice-boltzmann method for the
[8] P. Plouffe, D.M. Roberge, J. Sieber, M. Bittel, A. Macchi, Liquid–liquid mass simulation of two-phase flows in microchannel, Chem. Eng. Sci. 99 (2013)
transfer in a serpentine micro-reactor using various solvents, Chem. Eng. J. 285 238–249.
(2016) 605–615. [41] H. Yang, Q. Zhou, L.-S. Fan, Three-dimensional numerical study on droplet
[9] D. Tsaoulidis, P. Angeli, Effect of channel size on mass transfer during formation and cell encapsulation process in a micro t-junction, Chem. Eng. Sci.
liquid–liquid plug flow in small scale extractors, Chem. Eng. J. 262 (2015) 87 (2013) 100–110.
785–793. [42] K. Wang, Y. Lu, J. Xu, J. Tan, G. Luo, Generation of micromonodispersed droplets
[10] J. Wan, L. Shi, B. Benson, M.J. Bruzek, J.E. Anthony, P.J. Sinko, R.K. Prudhomme, and bubbles in the capillary embedded t-junction microfluidic devices, AIChE J.
H.A. Stone, Microfluidic generation of droplets with a high loading of 57 (2011) 299–306.
nanoparticles, Langmuir 28 (2012) 13143–13148. [43] M. De Menech, P. Garstecki, F. Jousse, H. Stone, Transition from squeezing to
[11] A.M. Nightingale, J.C. deMello, Segmented flow reactors for nanocrystal dripping in a microfluidic t-shaped junction, J. Fluid Mech. 595 (2008) 141–
synthesis, Adv. Mater. 25 (2013) 1813–1821. 161.
[12] L. Zhang, G. Niu, N. Lu, J. Wang, L. Tong, L. Wang, M.J. Kim, Y. Xia, Continuous [44] G.F. Christopher, N.N. Noharuddin, J.A. Taylor, S.L. Anna, Experimental
and scalable production of well-controlled noble-metal nanocrystals in observations of the squeezing-to-dripping transition in t-shaped microfluidic
milliliter-sized droplet reactors, Nano Lett. 14 (2014) 6626–6631. junctions, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 78 (2008)
[13] G. Niu, A. Ruditskiy, M. Vara, Y. Xia, Toward continuous and scalable 036317.
production of colloidal nanocrystals by switching from batch to droplet [45] T. Fu, Y. Ma, D. Funfschilling, C. Zhu, H.Z. Li, Squeezing-to-dripping transition
reactors, Chem. Soc. Rev. 44 (2015) 5806–5820. for bubble formation in a microfluidic t-junction, Chem. Eng. Sci. 65 (2010)
[14] G. Niu, M. Zhou, X. Yang, J. Park, N. Lu, J. Wang, M.J. Kim, L. Wang, Y. Xia, 3739–3748.
Synthesis of pt- ni octahedra in continuous-flow droplet reactors for the [46] J. Sivasamy, T.-N. Wong, N.-T. Nguyen, L.T.-H. Kao, An investigation on the
scalable production of highly active catalysts toward oxygen reduction, Nano mechanism of droplet formation in a microfluidic t-junction, Microfluid.
Lett. (2016). Nanofluid. 11 (2011) 1–10.
[15] P. Zhu, L. Wang, Passive and active droplet generation with microfluidics: a [47] S. Chen, K. Liu, C. Liu, D. Wang, D. Ba, Y. Xie, G. Du, Y. Ba, Q. Lin, Effects of
review, Lab. Chip 17 (2017) 34–75. surface tension and viscosity on the forming and transferring process of
[16] Y. Park, T.A. Pham, C. Beigie, M. Cabodi, R.O. Cleveland, J.O. Nagy, J.Y. Wong, microscale droplets, Appl. Surf, Sci, 2016.
Monodisperse micro-oil droplets stabilized by polymerizable phospholipid [48] R. Raj, N. Mathur, V.V. Buwa, Numerical simulations of liquid- liquid flows in
coatings as potential drug carriers, Langmuir 31 (2015) 9762–9770. microchannels, Ind. Eng. Chem. Res. 49 (2010) 10606–10614.
[17] Y. Li, K. Wang, J. Xu, G. Luo, A capillary-assembled micro-device for [49] L. Zhang, Y. Wang, L. Tong, Y. Xia, Synthesis of colloidal metal nanocrystals in
monodispersed small bubble and droplet generation, Chem. Eng. J. 293 droplet reactors: The pros and cons of interfacial adsorption, Nano Lett. 14
(2016) 182–188. (2014) 4189–4194.
[18] X. Wang, K. Wang, A. Riaud, X. Wang, G. Luo, Experimental study of liquid/ [50] Y. Li, D.G. Yamane, S. Li, S. Biswas, R.K. Reddy, J.S. Goettert, K. Nandakumar, C.S.
liquid second-dispersion process in constrictive microchannels, Chem. Eng. J. Kumar, Geometric optimization of liquid–liquid slug flow in a flow-focusing
254 (2014) 443–451. millifluidic device for synthesis of nanomaterials, Chem. Eng. J. 217 (2013)
[19] G.F. Christopher, S.L. Anna, Microfluidic methods for generating continuous 447–459.
droplet streams, J. Phys. D: Appl. Phys. 40 (2007) R319. [51] A.R. Abate, M. Kutsovsky, S. Seiffert, M. Windbergs, L.F. Pinto, A. Rotem, A.S.
[20] C.N. Baroud, F. Gallaire, R. Dangla, Dynamics of microfluidic droplets, Lab. Chip Utada, D.A. Weitz, Synthesis of monodisperse microparticles from non–
10 (2010) 2032–2045. newtonian polymer solutions with microfluidic devices, Adv. Mater. 23
[21] H. Gu, M.H. Duits, F. Mugele, Droplets formation and merging in two-phase (2011) 1757–1760.
flow microfluidics, Int. J. Mol. Sci. 12 (2011) 2572–2597. [52] P.E. Arratia, J.P. Gollub, D.J. Durian, Polymeric filament thinning and breakup in
[22] A. Huerre, V. Miralles, M.-C. Jullien, Bubbles and foams in microfluidics, Soft microchannels, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 77 (2008)
Matter 10 (2014) 6888–6902. 036309.
[23] T. Thorsen, R.W. Roberts, F.H. Arnold, S.R. Quake, Dynamic pattern formation in [53] D. Qiu, L. Silva, A.L. Tonkovich, R. Arora, Micro-droplet formation in non–
a vesicle-generating microfluidic device, Phys. Rev. Lett. 86 (2001) 4163. newtonian fluid in a microchannel, Microfluid. Nanofluid. 8 (2010)
[24] T. Nisisako, T. Torii, T. Higuchi, Droplet formation in a microchannel network, 531–548.
Lab. Chip 2 (2002) 24–26. [54] M. Aytouna, J. Paredes, N. Shahidzadeh-Bonn, S. Moulinet, C. Wagner, Y.
[25] J.D. Tice, A.D. Lyon, R.F. Ismagilov, Effects of viscosity on droplet formation and Amarouchene, J. Eggers, D. Bonn, Drop formation in non–newtonian fluids,
mixing in microfluidic channels, Anal. Chim. Acta 507 (2004) 73–77. Phys. Rev. Lett. 110 (2013) 034501.
[26] Y.-C. Tan, J.S. Fisher, A.I. Lee, V. Cristini, A.P. Lee, Design of microfluidic channel [55] T. Fu, L. Wei, C. Zhu, Y. Ma, Flow patterns of liquid–liquid two-phase flow in
geometries for the control of droplet volume, chemical concentration, and non–newtonian fluids in rectangular microchannels, Chem. Eng. Process.
sorting, Lab. Chip 4 (2004) 292–298. Process Intensif. 91 (2015) 114–120.
[27] L.S. Roach, H. Song, R.F. Ismagilov, Controlling nonspecific protein adsorption [56] V. Ranade, Computational Flow Modeling for Chemical Reactor Engineering,
in a plug-based microfluidic system by controlling interfacial chemistry using Process systems engineering, Academic Press, 2002.
fluorous-phase surfactants, Anal. Chem. 77 (2005) 785–796. [57] J. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface
[28] D. Dendukuri, K. Tsoi, T.A. Hatton, P.S. Doyle, Controlled synthesis of tension, J. Comput. Phys. 100 (1992) 335–354.
nonspherical microparticles using microfluidics, Langmuir 21 (2005) 2113– [58] A. Fluent, Ansys fluent 17.0 user’s guide, ANSYS FLUENT Inc., 2017.
2116. [59] Z. Guo, D.F. Fletcher, B.S. Haynes, Implementation of a height function method
[29] P. Garstecki, M.J. Fuerstman, H.A. Stone, G.M. Whitesides, Formation of to alleviate spurious currents in cfd modelling of annular flow in
droplets and bubbles in a microfluidic t-junctionscaling and mechanism of microchannels, Appl. Math. Modell. 39 (2015) 4665–4686.
break-up, Lab. Chip 6 (2006) 437–446. [60] J. Choi, G. Son, Numerical study of droplet motion in a microchannel with
[30] J. Xu, S. Li, J. Tan, G. Luo, Correlations of droplet formation in t-junction different contact angles, J. Mech. Sci. Technol. 22 (2008) 2590.
microfluidic devices: from squeezing to dripping, Microfluid. Nanofluid. 5 [61] D. Hoang, L. Portela, C. Kleijn, M. Kreutzer, V. Van Steijn, Dynamics of droplet
(2008) 711–717. breakup in a t-junction, J. Fluid Mech. 717 (2013).
[31] A. Gupta, S.S. Murshed, R. Kumar, Droplet formation and stability of flows in a [62] Z. Zhang, J. Xu, B. Hong, X. Chen, The effects of 3d channel geometry on ctc
microfluidic t-junction, Appl. Phys. Lett. 94 (2009) 164107. passing pressure–towards deformability-based cancer cell separation, Lab.
[32] A. Gupta, R. Kumar, Effect of geometry on droplet formation in the squeezing Chip 14 (2014) 2576–2584.
regime in a microfluidic t-junction, Microfluid. Nanofluid. 8 (2010) [63] Y. Kagawa, T. Ishigami, K. Hayashi, H. Fuse, Y. Mino, H. Matsuyama, Permeation
799–812. of concentrated oil-in-water emulsions through a membrane pore: numerical
[33] W. Wang, Z. Liu, Y. Jin, Y. Cheng, Lbm simulation of droplet formation in micro- simulation using a coupled level set and the volume-of-fluid method, Soft
channels, Chem. Eng. J. 173 (2011) 828–836. Matter 10 (2014) 7985–7992.
[34] T. Fu, Y. Wu, Y. Ma, H.Z. Li, Droplet formation and breakup dynamics in [64] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions
microfluidic flow-focusing devices: from dripping to jetting, Chem. Eng. Sci. 84 to incompressible two-phase flow, J. Comput. Phys. 114 (1994) 146–159.
(2012) 207–217. [65] W. Lee, G. Son, Bubble dynamics and heat transfer during nucleate boiling in a
[35] T. Fu, Y. Ma, D. Funfschilling, C. Zhu, H.Z. Li, Breakup dynamics of slender microchannel, Numer. Heat Tr. A - Appl. 53 (2008) 1074–1090.
bubbles in non–newtonian fluids in microfluidic flow-focusing devices, AIChE [66] R.I. Issa, Solution of the implicitly discretised fluid flow equations by operator-
J. 58 (2012) 3560–3567. splitting, J. Comput. Phys. 62 (1986) 40–65.
 S.G. Sontti, A. Atta / Chemical Engineering Journal 330 (2017) 245–261 261

[67] B.P. Leonard, A stable and accurate convective modelling procedure based on [70] M. Yoshino, Y. Hotta, T. Hirozane, M. Endo, A numerical method for
quadratic upstream interpolation, Comput. Methods in Appl. Mech. Eng. 19 incompressible non–newtonian fluid flows based on the lattice boltzmann
(1979) 59–98. method, J. Non–Newtonian Fluid Mech. 147 (2007) 69–78.
[68] M. Seifollahi, E. Shirani, N. Ashgriz, An improved method for calculation of [71] M.J. Rhodes, Introduction to Particle Technology, John Wiley & Sons, 2008.
interface pressure force in plic-vof methods, Eur. J. Mech. B-Fluid 27 (2008) 1– [72] Y. Lu, T. Fu, C. Zhu, Y. Ma, H.Z. Li, Scaling of the bubble formation in a flow-
23. focusing device: role of the liquid viscosity, Chem. Eng. Sci. 105 (2014) 213–
[69] R. Gupta, D.F. Fletcher, B.S. Haynes, On the cfd modelling of taylor flow in 219.
microchannels, Chem. Eng. Sci. 64 (2009) 2941–2950.