Deformation and breakup of a liquid droplet past a solid circular cylinder: A lattice Boltzmann study Qiuxiang Li,1 Zhenhua Chai,2 Baochang Shi,2,* and Hong Liang1 1

National Laboratory of Coal Combustion, Huazhong University of Science and Technology, Wuhan 430074, China 2 School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China (Received 9 December 2013; published 21 October 2014) In this paper, we present a numerical study on the deformation and breakup behavior of liquid droplet past a solid circular cylinder by using an improved interparticle-potential lattice Boltzmann method. The effects of the eccentric ratio β, viscosity ratio λ between the droplet and the surrounding fluid, surface wettability, and Bond number (Bo) on the dynamic behavior of the liquid droplet are considered. The parameter β represents the degree that the solid cylinder deviates from the center line, and Bo is the ratio between the inertial force and capillary force. Numerical results show that there are two typical patterns, i.e., breakup and no breakup, which are greatly influenced by the aforementioned parameters. When β increases to a critical value βc , the droplet can pass the circular cylinder without a breakup, otherwise, the breakup phenomenon occurs. The critical eccentric ratio βc increases significantly with increasing Bo for case with λ > 1, while for the case with λ < 1, the viscosity effects on the βc is not obvious when Bo is large. For the breakup case, the amount of deposited liquid on the tip of the circular cylinder is almost unaffected by β. In addition, the results also show that the viscosity ratio and wettability affect the deformation and breakup process of the droplet. For case with λ < 1, the viscosity ratio plays a minor role in the thickness variations of the deposited liquid, which decreases to a nonzero constant eventually; while for λ > 1, the increase of the viscosity ratio significantly accelerates the decrease of the deposited liquid, and finally no fluid deposits on the cylinder. In term of the wettability, there occurs continuous gas phase trapped by the wetting droplet, but this does not happen for nonwetting droplet. Besides, for λ < 1, the time required to pass the cylinder (tp ) decreases monotonically with decreasing contact angle, while a nonmonotonic decrease appears for λ > 1. It is also found that tp decreases monotonically with increasing Bo and is less sensitive to λ at a large Bo. DOI: 10.1103/PhysRevE.90.043015

PACS number(s): 47.61.−k, 02.60.−x, 68.08.−p, 68.05.−n

I. INTRODUCTION

The droplet dynamics past solid obstacles is frequently encountered in a broad range of industry applications, such as the immiscible displacement during oil recovery, inkjet printing, droplet-based microfluidics, and so on. In such systems, the dynamics of drop deformation and breakup are usually very complicated for the complex interactions between components [1]. On the other hand, the study of the droplet past a solid target has recently been an active field of research for its great importance in many fields. Pasandideh-Fard [2] presented experimental results of the droplet with varied diameters falling on a tube and disintegrating into smaller drops. It is found that the impact region on the tube surface significantly influences the number of daughter droplets. To obtain high-throughput droplets with precise size, Link et al. [3] introduced two possible ways, an obstacle-meditated geometry and a T junction, to break large drops into minor ones. They showed the size ratio of daughter drops can be controlled by changing the transverse position of the obstacle, or the length ratio of branch channels of the T junction. Chung et al. [4] studied the dynamics of both Newtonian and viscoelastic droplets passing through a cylinder obstruction in a confined microchannel, and they found that the breakup of the Newtonian droplet past a solid obstruction depends on the droplet size, while for viscoelastic droplet, the normal stress difference is a key factor in inducing breakup and extension. Then Chung et al. [5] extended their research

*

[email protected]

1539-3755/2014/90(4)/043015(10)

to a more complicated case where more obstructions are included, considering the effects of droplet size, obstruction shape, and capillary number. By using experimental and theoretical approaches, Bakshi et al. [6] investigated the droplet impacts onto a small spherical target and identified three regimes: the initial drop deformation phase, the inertiadominated phase, and the viscosity-dominated phase. Proti`ere et al. [7] studied the dynamics of a droplet past a circular cylinder and focused on the capillary effects. In this work, the authors developed a one-dimensional model to explain the experimental phenomena and found that the breakup is caused by the surface-tension-driven instability as well as variation of the axial gaps. Salkin et al. [8] investigated the similar problem except that the circular cylinder is replaced with a linear obstacle, and with a focus on the relation between drop size and critical capillary number for drop breakup. They found whether the relation is monotonic or not depends on the viscosity contrast. Recently, a better understanding of the impact behavior between a drop and a target with various shapes was obtained by Gac and Grado´n [9] and Juarez et al. [10]. The former work identified two kinds of splashing regimes which depend on the number of the target vertices, while the latter work found that the behavior of a head-to-head collision is almost independent of the target shape. To provide guidance to high-quality prints, Quang et al. [11] numerically studied an oil droplet hitting a paper-like structure that is composed of a random arrangement of solid cylinders. The results show that the behavior of the oil droplet is greatly influenced by the wettability of the solid cylinder: the droplet may deposit on the deeper layer of solid cylinder after passing through adjacent cylinders and breaking by the air, or bounce

043015-1

©2014 American Physical Society

LI, CHAI, SHI, AND LIANG

PHYSICAL REVIEW E 90, 043015 (2014)

off and retain the complete shape. From the perspective of promoting mixing in microfluidics, Yin and Luo [12] studied the flow patterns of a drop past protrusions located on one side of the wall. Different from the above mentioned studies, no breakup is observed in this work. In addition, the results also show that the flow pattern inside the drop largely depends on the viscosity ratio and geometry of the obstructions, and the low viscosity or small capillary number tends to enhance the mixing and vice versa. Although there are some works on the dynamic behavior of a droplet past one or more obstacles, the deformation and breakup of a droplet past a solid circular cylinder has not been investigated comprehensively, and also, the effects of the impact region, surface wettability, and viscosity ratio on the droplet behaviors are still not well understood. Therefore, we will perform a detailed study on the deformation and breakup of a droplet past a solid circular cylinder by using the lattice Boltzmann method. The lattice Boltzmann method (LBM), as a promising numerical approach, has received more and more attention in studying complex fluid flows for its simplicity and computational efficiency [13,14]. Originating from the lattice gas automaton (LGA) method, LBM can also be viewed a spatial, time, and velocities discretization of the continuous Boltzmann equation. By simulating the collision and streaming procedure of the fluid particles over all the grids of the flow regions, and evaluating hydrodynamic moments of the distribution function, macroscopic quantities then can be obtained. Based on its kinetic characteristics, LBM is very attractive in modeling interfacial phenomena involved in multiphase and multicomponent flows [15,16]. Among those existing multiphase LBM models, the interparticle-potential multiphase model proposed by Shan and Chen (S-C) [17,18] has been widely used for its simplicity and remarkable versatility, such as easy to incorporate complex interactions and implement different wettability conditions [19]. However, the model also receives criticisms, in which the most important one is the numerical instability for low viscosity. To overcome this shortcomings, Porter et al. [20] proposed an improved interaction-potential model for multicomponent and/or multiphase flows, which is

σ (eq) fi

developed by incorporating the external force into the discrete lattice Boltzmann equation (LBE). The newly explicit forcing (EF) model coupled with a multiple-relaxation time (MRT) collision scheme [21–23] can not only improve numerical stability for fluid with a low viscosity, but also maintain grid-independent fluid properties. Therefore in the present work, the EF interparticle-potential model is used to study the deformation and breakup of a liquid droplet past a solid circular cylinder, aiming to provide a deeper understanding on this problem. The rest of the paper is organized as follows. In Sec. II the EF interparticle-potential model for multiphase flows is briefly reviewed, followed by model validation. In Sec. III, numerical results and discussion are presented, and finally, a summary of the main conclusions is given in Sec. IV.

II. MATHEMATICAL FORMULATION A. The EF interparticle-potential lattice Boltzmann Method (LBM) model

In the EF interparticle-potential model, the force term is directly added to the discrete lattice Boltzmann equation, and the evolution of this model with a multiple-relaxation-time (MRT) collision operator can be written as [20,24] f σ (x + cδt,t + δt) − f σ (x,t) = −T−1 S¯ σ T[f σ (x,t) − f σ (eq) (x,t)] S¯ σ −1 I− TFσ (x,t), + δtT 2

(1)

where f σ (x,t) is the density distribution function of the σ th component with velocity c at position x and time t, Fσ represents the external force term, and the boldface symbols f σ and Fσ denote nine-dimensional column vectors, i.e., f σ (x,t) = [f0σ (x,t),f1σ (x,t), . . . , f8σ (x,t)]T , Fσ (x,t) = σ (eq) [F0σ (x,t),F1σ (x,t), . . . , F8σ (x,t)]T . fi (x,t) is the corresponding equilibrium distribution function, and given by

eq 2 eq 2 eq uσ ci · uσ ci · uσ , = ωi ρσ 1 + + − 2 4 cs 2cs 2cs2

(2)

where ωi are the weighting coefficients and cs is the speed of sound. For the D2Q9 lattice model (one special model with nine velocities in the two-dimensional space) considered in the present work, ωi are set as ω0 = 4/9, ω1−4 = 1/9, ω5−8 = 1/36, and the discrete velocities ci are given as ⎧ ⎪ ⎨(0,0)c, ci = ( cos[(i − 1)π/2], sin[(i − 1)π/2])c, ⎪ ⎩ 2( cos[(i − 5)π/2 + π/4], sin[(i − 5)π/2 + π/4])c,

i=0 i = 1,2,3,4,

(3)

i = 5,6,7,8

where c = δx/δt (set to be 1 in this work) with δx and δt representing the lattice spacing and time step, respectively. Note that c √ is related to the speed of sound based on the relation cs = c/ 3. 043015-2

DEFORMATION AND BREAKUP OF A LIQUID DROPLET . . .

T is a transformation matrix and defined as [21] ⎡ 1 1 ⎢−4 −1 ⎢ ⎢ 4 −2 ⎢ 1 ⎢ 0 ⎢ T = ⎢ 0 −2 ⎢ 0 0 ⎢ ⎢ 0 0 ⎢ ⎣ 0 1 0 0

1 −1 −2 0 0 1 −2 −1 0

1 −1 −2 −1 2 0 0 1 0

which can be used to project f σ , f σ (eq) , and Fσ onto the moment space according to mσ = Tf σ , mσ (eq) = Tf σ (eq) , and F¯ σ = TFσ [25]. The matrix S¯ σ in Eq. (1) is a diagonal relaxation matrix and given by S¯ σ = diag scσ ,seσ ,sεσ ,scσ ,sqσ ,scσ ,sqσ ,sυσ ,sυσ , (4) where siσ is the inverse of the relaxation time for fiσ relaxing to its equilibrium state. The parameter scσ corresponds to conserved moments and set to be one in our simulations, whereas seσ , sεσ , sqσ , and svσ corresponding to the nonconserved moments that can be adjusted independently to improve the accuracy and stability of the MRT model [21–23]. The rest of parameters are taken as follows: seσ = 1.64, sεσ = 1.54, and sqσ = svσ = 1.0. It is worth noting that svσ is related to the viscosity of the fluid based on the following relation: 1 (5) v σ = cs2 σ − 0.5 δt. sv The second term Fiσ in the right-hand side of Eq. (1) accounts for the force term imposed on the σ th component and defined as [20,26] Fiσ =

Gσ · (e − ueq ) σ (eq) fi , ρσ cs2

PHYSICAL REVIEW E 90, 043015 (2014)

1 −1 −2 0 0 −1 2 −1 0

1 2 1 −1 −1 −1 −1 0 1

⎤ 1 2⎥ ⎥ 1⎥ ⎥ 1⎥ ⎥ 1⎥ , −1⎥ ⎥ −1⎥ ⎥ 0⎦ −1

where gσ w is the interaction strength between a fluid phase and a solid phase, a positive value of gσ w is used to describe the fluid is nonwetting to the surface and negative for wetting to the surface. S is a constant at the solid phase and zero elsewhere [27]. The wettability is described by the static contact angle θ , which can be changed ranging from 0 to 180 degrees through adjusting the parameter gσ w conveniently, and it is found that the static contact angles are approximately a linear function of gσ w . Usually the contact angle is not the unique static contact angle when the droplet moves, but in the following simulations, the contact angle only means the static contact angle. Different from the original S-C model, the ueq in Eq. (6) is calculated as eq σ u = sc ρσ uσ scσ ρσ , (10) σ

ρσ =

fiσ , ρσ uσ =

i

ci fiσ +

δt σ G . 2

(11)

The total density and velocity of the fluid mixture are σ ρσ uσ σ ρ= fi , u = . (12) σ ρσ σ i

(7)

where σ¯ denotes the other fluid components. The quantity ψσ (x) generally called the “effective density” is a function of the local density, and different forms of ψσ (x) will give different equations of state; here it is taken as ρσ (x) [19]. Gσ σ¯ is the interaction strength between fluids, and it should be noted that the value of Gσ σ¯ (σ = σ¯ ) usually varies in a special region such that the relatively nondiffusing interface can be preserved.

σ

where ρσ and uσ are the density and velocity of the σ th component, and calculated by

i

σ¯

(9)

x

(6)

where Gσf , Gσs , and Gσe represent the fluid-fluid force, fluidsolid force, and other possible external force such as the gravity force. Based on the original S-C model [17,18], the fluid-fluid interaction force acting on the σ th component at site x can be given by [19] ωi Gσ σ¯ ψσ¯ (x )(x − x), (8) Gσf (x) = −ψσ (x) x

1 2 1 −1 −1 1 1 0 −1

Similarly, the fluid-solid force is defined as ωi S(x )(x − x ), Gσs (x) = −ψσ (x)gσ w

where Gσ is the external force. In the EF interparticle-potential model, the total force Gσ in Eq. (6) is composed of three parts, Gσ = Gσf + Gσs + Gσe ,

1 2 1 1 1 1 1 0 1

Generally, by performing a Taylor expansion on Eq. (1) in time and space, and employing the Chapman-Enskog analysis, the macroscopic Navier-Stokes equations can be recovered from LBE. We refer readers to Refs. [15,28] for more details. Note that the MRT model will reduce to the singlerelaxation-time (SRT) model when all the nine relaxation factors are equal to each other. Similar to the SRT model, the evolutionary process of the MRT model can be divided into the following two steps: (1) The collision in the momentum space

043015-3

f σ + (x + cδt,t + δt) = f σ (x,t) − (T−1 S¯ σ )[mσ (x,t) − mσ (eq) (x,t)] S¯ σ ¯ σ −1 I− F (x,t), (13a) + δtT 2

LI, CHAI, SHI, AND LIANG

PHYSICAL REVIEW E 90, 043015 (2014)

−3

2

−3

x 10

3

x 10

2.5 1.5

1

u

u

2 1.5 1 0.5

Analytical solution LBM

0 −50 −40 −30 −20 −10

0 y

10

20

0.5

30

40

Analytical solution LBM

0 −50 −40 −30 −20 −10

50

(a)

0 y

10

20

30

40

50

(b)

FIG. 1. Steady state velocity profiles as a function of the channel width for (a) M = 1/150, (b) M = 150.

(2) The propagation in the velocity space f σ (x + cδt,t + δt) = f σ + (x,t).

(13b)

away from the interface. In the following, we will adopt the EF interparticle-potential model to investigate the complex dynamics of the droplet past a solid circular cylinder.

B. Model validation

III. NUMERICAL RESULTS AND DISCUSSION

The two-phase Poiseuille flow driven by a body force is studied in a 100 × 200 lattice system and used to validate the EF interparticle-potential model. Initially bulk phase 1 is in the center region of 0 < |y| < a, and bulk phase 2 fills the region of a < |y| < b, the dynamic viscosity ratio is M = ρρ12 υυ12 . In our simulations, the periodic boundary condition is applied to the inlet and outlet, and no-slip boundary conditions are applied to the up and bottom walls. For this special problem, there is an analytic solution of the velocity [20]:

A. Physical problem

Figure 2 shows the schematic of the physical problem with Lx × Ly = 150 × 900 (all the variables are used in LB unit). Initially a droplet with radius R = 50 is placed in the axis

2R

u1 = A1 y 2 + C1 , 0 < |y| < a,

g

d

u2 = A2 y 2 + B2 y 2 + C2 , a < |y| < b,

s

where

s F1 F2 , A2 = − , A1 = − 2ρ1 υ1 2ρ2 υ2

2r

B2 = −A2 a + 2MA1 a,

d

Ly C1 = (A2 − A1 )a − B2 (b − a) − A2 b , 2

2

C2 = −A2 b2 − B2 b. Figure 1 shows the steady state velocity profiles at the viscosity ratio M = 1/150 and M = 150, respectively. It can be seen that the numerical results agree with the corresponding analytical solutions, except for a little discrepancy near the interface, which may be caused by the force imbalance at the interface [29] and the diffused interface feature of the LBM. This is because the interface in the diffused interface method is assumed to have a finite width, whereas the condition of a sharp interface is used to derive the analytical solution, so the error near the interface seems inevitable. However, the discrepancy has no significant effects on the velocity profiles

Lx FIG. 2. (Color online) Schematic illustration of the flow geometry.

043015-4

DEFORMATION AND BREAKUP OF A LIQUID DROPLET . . .

(a)

(b)

PHYSICAL REVIEW E 90, 043015 (2014)

(c)

(a)

(b)

(c)

FIG. 3. (Color online) Effects of the grid size on the deformation and breakup of liquid droplet past a solid circular cylinder for (left) λ = 0.02, (right) λ = 50, and grid size with Lx×Ly (a) 300×1800, (b) 150×900, and (c) 74×450.

of the channel, followed by a solid circular cylinder with radius r = 40. The movement of the droplet is driven by the gravity g. The vertical distance (s) between the droplet and circular cylinder is fixed, whereas the transverse position of the circular cylinder is varied. To quantify the deviation of the circular cylinder from the axis, an eccentric ratio β = 1 − d/( Lx − r) is used, where d is the minimum distance from the 2 right wall. Note that it represents a head-to-head impact when d is set as Lx − r, otherwise it represents the eccentric impact. 2 In the following simulations, periodic boundary condition is applied at the inlet and outlet of the channel, and the half-way bounce-back scheme is applied to all solid points, including the solid cylinder and the two channel walls. The bounce-back boundary condition means a boundary that lies halfway between the boundary nodes and the neighboring fluid nodes, the particle streams to a wall node and returns back to the node it came from. In this work both fluids are taken to be the Newtonian fluid, and the densities are set to be 1, the ratio of the droplet viscosity μd to the surrounding fluid viscosity μs is λ = μμds , and the ratio of gravitational force to capillary force is 2

the Bond number, defined as Bo = ρd gR , where σ is the surface σ tension force. The evolution time t0 is nondimensionalized with t = t0R∗U , where U is the characteristic velocity defined 2

by U = gρμd dR [19]. The custom C code is written to complete all the following simulations, and the code for each case is run on an Intel Core i5-2400 PC. Before performing the numerical simulations, two important things need to be done. The first is the gridindependence convergence test, and the second is to provide a comparison between present LBM results and some previous experiments, which can be further used to validate present LBM for such complex flows. For the first issue, three different grid sizes Lx × Ly = 300 × 1800, 150 × 900, 74 × 450 have been used perform a grid-independence test, and the results in

Fig. 3 show that the droplet dynamics obtained with 150 × 900 is similar to the case of 300 × 1800, but different from the case with 74 × 450. In the following, the grid size 150 × 900 will be used through considering the computational cost and numerical accuracy. For the second issue, we carried out a comparison between present LBM results and some previous experimental results [5], and presented the results in Fig. 4 where a droplet past a square and a circular cylinder are considered. As seen from this figure, the droplet gradually deforms during the collision with the solid cylinder, and the liquid film hanging on the cylinder becomes thin, and finally the breakup phenomena can be observed. However, it should be noted that the breakup positions of the droplets past the square and circular cylinder are quite different: for the former case, the droplet breaks at the center of the circular cylinder, and thus only two equal-sized daughter droplets are formed; for the latter case, the droplet breaks at the two corners of the square. These results qualitatively agree with experimental results reported in Ref. [5] and also support that the present LBM can be used to study the breakup dynamics of the droplet passing a solid cylinder. Here we also would like to point out that there are still some quantitative discrepancies of droplet dynamics between two-dimensional simulations and three-dimensional experiments, which may be caused by the interface curvature and viscous drag in three-dimensional problems [5,30], but the two-dimensional simulations are helpful to understand the droplet breakup behavior. For this reason, we first study a 2D droplet past a circular cylinder, and the 3D problem will be studied in the future. In addition, we also note that, although some related work have been conducted in Ref. [5], the effects of the eccentric ratio, surface wettability, and viscosity are not considered in their work. To fill the gap, a comprehensive study will be carried out to gain a clear and deep understanding of the droplet dynamics past a solid cylinder.

043015-5

LI, CHAI, SHI, AND LIANG

PHYSICAL REVIEW E 90, 043015 (2014)

t = 0.0

t = 2.09

t = 4.18

t = 5.43

t = 6.26

t = 0.0

t = 1.25

t = 1.67

t = 2.51

t = 2.92

(a)

(b)

t=1.67

t=2.92

t=3.33

t=4.17 (a)

t=5.0

t=5.42

t=5.83

t=1.67

t=2.5

t=3.33

t=4.17 (b)

t=5.0

t=5.83

t=6.25

FIG. 4. (Color online) Comparisons between the present simulations and previous experiments [5] on the droplet past an obstruction Ca = μdσU = 0.01, λ = 0.2; (a) the droplet past a square cylinder (left: numerical, right: experimental), (b) the droplet past a circular cylinder (left: numerical, right: experimental).

B. Effect of the eccentric ratio β

We first performed simulations to investigate the effect of the eccentric ratio β by changing the transverse position of the circular cylinder. Figure 5 presents the time evolution of the droplet past a circular cylinder with fixed Bo = 0.307, θ = 170◦ , λ = 0.02, and varying β. For the first case with β = 0 [see Fig. 5(a)], the droplet initially deforms and spreads radially, shaped like a meniscus hanging on the circular cylinder at t = 2.92. Later, the meniscus is pulled down by the gravity and fills the two gaps between the cylinder and solid walls. The layer of liquid film between the droplet rear and tip of the circular cylinder becomes uniformly thinner and thinner. In the meanwhile, the surface tension force enhances its contraction to two sides [31]. The combination of different forces eventually results in two break points of the stretched droplet at the left and right fronts of the circular cylinder. The formed daughter droplets continue to fall and a small blob of liquid deposits on the tip of the circular cylinder. The results are consistent with those reported in Ref. [5]. When increasing β [see Fig. 5(b), β = 1/7], the droplet breakups into two asymmetric daughter droplets in that, as the gap between the circular cylinder and the right wall decreases, a larger pressure is built up at the right side, and thus more liquid is inclined to intrude into the left gap. By further increasing β up to 2/7 [see Fig. 5(c)], the leading meniscus in the left gap moves faster and the minor meniscus in right gap is pulled back. As a result, no breakup phenomenon is observed as the droplet passes the circular cylinder. Besides that, it is also

t=1.67

t=2.5

t=3.33

t=4.17

t=5.0

t=5.83

(c)

FIG. 5. (Color online) Time evolution of the droplet past a circular cylinder at θ = 170◦ , λ = 0.02 (a) head-to-head impact: β = 0, (b) eccentric impact: β = 1/7, (c) eccentric impact: β = 2/7.

found that the nonbreak droplet falls faster than the other two cases. We further studied the effect of β on the size ratio of the right daughter droplet to the left one and show the results in Fig. 6. As seen from this figure, the size ratio decreases with the increase of β, whereas the amount of residual liquid on the tip of the circular cylinder for different β is almost same since the droplet breakup happens at a similar position. In addition, it should be noted that, when β reaches a critical value, which is about 0.257, the whole droplet can pass the circular cylinder without breakup. C. Effect of the viscosity ratio

To investigate the effects of viscosity ratio, several simulations of head-to-head impact, the same as that in Fig. 5(a) were carried out, in which the Bond number and surface wettability are set to be constants, with the sole exception of

043015-6

DEFORMATION AND BREAKUP OF A LIQUID DROPLET . . .

PHYSICAL REVIEW E 90, 043015 (2014) 1

λ=0.02 λ=0.2 λ=5 λ=50

0.9

1

0.8 0.7 0.6

0.6

h 2R

size ratio

0.8

0.5 0.4

0.4 0.3 0.2

0.2

0.1

0 0

0.05

0.1

0.15 β

0.2

0.25

0 0

0.3

FIG. 6. (Color online) Effect of the eccentric ratio β on the size ratio of daughter droplets.

a different viscosity ratio (in the following, the head-to-head impact is considered unless stated otherwise). Figure 7 shows the dynamic behavior of the droplet with λ = 50; similar as reported previously, the droplet first deforms symmetrically as it approaches the circular cylinder. However, a different behavior appears at t = 2.83, and both ends of the meniscus elongate under the gravity, forming two long and thin liquid fingers in the left and right gaps. The formed daughter droplets are also inclined to gather to the axes of the channel. The thin film on the circular cylinder stretches and breaks up at at the center tip finally, and there is no liquid deposited on the circular cylinder. The breakup dynamics is different from the case with a small λ that has two breakup points and a liquid film leaving on the circular cylinder [see Fig. 5(a)]. Actually, the difference of the flow behavior with a small [see Fig. 5(a)] and a large viscosity ratios (see Fig. 7) can be explained by the forces acting on the droplet. It is known that fluid with higher viscosity always produces a higher shear force when acted on by external forces [32]. The shear force as well as gravity in the system of this work is the main deformation force, but on the contrary, the surface tension force tries to prevent such deformation [33]. For a fixed surface tension force, the high-viscosity

0.5

1

1.5

2

2.5

3

3.5

t

4

4.5

5

5.5

6

6.5

7

FIG. 8. Time evolution of the normalized thickness of the deposited film with different λ.

droplet always suffers a higher deformation force relative to a low-viscosity droplet. When the viscosity ratio is large, the high shear force and gravity force take over the surface tension force and lead to an elongation of the daughter droplets during the impact process. Hence, eventually results in two filaments. In contrast, the low-viscosity droplet is easier to maintain a circular shape since the force responsible for deformation is relatively weak. The above results indicate that the viscosity ratio has an influence on the variations of the thickness evolution of the deposited film. To quantitatively show the results, the film thickness h, defined as the distance from the top interface of the droplet to the bottom interface on the cylinder is measured. Without loss of generality, four cases with viscosity ratio λ ranging from 0.02 to 50 are considered. The evolution of the h normalized thickness 2R of the deposited film at different λ are shown in Fig. 8. It is found that the viscosity ratio plays a minor role before the droplet impacting with the circular cylinder. Afterwards, the effect of the viscosity ratio begins to appear around t = 1.5. For a droplet with λ > 1, the thickness of a high-viscosity droplet is smaller than that of the low-viscosity droplet at the same moment, so it requires less time for a high-viscosity droplet to break up. For a droplet with λ < 1, however, the thickness is less sensitive to the variations of the viscosity ratio. In addition, the thickness does not decreases monotonically to a minimum, but oscillates in the process and finally reaches a constant value. This is because the liquid film is thin at the early stage of breakup and then contracts to a tiny droplet under the effect of the surface tension force. However, no deposition is observed when λ > 1, so h falls to zero finally. D. Effect of the surface wettability

t=1.67

t=1.92

t=2.17

t=2.83

t=3.25

t=3.50

t=3.83

FIG. 7. (Color online) Time evolution of the droplet past a circular cylinder at λ = 50, θ = 170◦ .

In both cases of Figs. 5 and 7, the droplet is nonwetting to the solid circular cylinder, and for the purpose of comparison, we also studied the dynamics of the wetting droplet past a circular cylinder. Figure 9 shows the evolutionary process of a wetting droplet with the static contact angle θ = 70◦ and low-viscosity ratio λ = 0.02. We found that the dynamics of

043015-7

LI, CHAI, SHI, AND LIANG

PHYSICAL REVIEW E 90, 043015 (2014)

FIG. 9. (Color online) Time evolution of the droplet past a circular cylinder at θ = 70◦ , λ = 0.02.

t=2.5

t=3.33

t=4.17

t=5.0

t=5.83

t=6.67

the wetting droplet is more complicated than the nonwetting one. First, two bubbles between the drop and the cylinder form during the descending process; this may be because the surface tension force tends to contract the droplet to round [33]. The continuous phase is entrapped as both ends of the meniscus are absorbed by the circular cylinder surface at t = 3.33. Then a third blob of the continuous phase is entrapped during the coalescence process, which is not observed in the nonwetting case. Second, the coalescence of daughter droplets occurs earlier than that of the nonwetting case, and third, more liquid deposits on the circular cylinder, and the deposition place is near the bottom surface of circular cylinder, which are different from those of the nonwetting case. We further conducted more simulations over a wide range of the contact angle and found that the deposition place is changed from the upward to downward along the cylinder surface with the decrease of the contact angle. This may be because a wetting droplet is inclined to attach to the surface tightly, and the thin layer of liquid deposition can easily be dragged away by the gravity force. For a higher viscosity droplet with λ = 50, the dynamic process of the wetting droplet past the cylinder is shown in Fig. 10. Compared to the nonwetting droplet with same viscosity ratio in Fig. 7, the wetting droplet undergoes more complex deformation, breakup, and coalescence. After the breakup, two ends of the daughter droplets first coalesce,

t=1.67

t=2.5

t=3.75

t=4.42

t=7.5 stretch to their maximum length at t = 5.0, and gradually recoil due to the surface tension force. It can also seen that some continuous phase is entrapped by the droplet during the coalescence process, and finally, the droplet breaks into three small droplets at t = 8.33, which is later than the nonwetting case. To further investigate the effect of surface wettability, we measured the passing time tp , defined as the time for the bulk of the droplet passing the circular cylinder. The numerical results reveal that the evolutions of tp , as a function of the contact angle, are different for cases of λ > 1 and λ < 1. To see these results more clearly, the variation of tp corresponding to the two cases are separately plotted in Fig. 11. As seen from this figure, the passing time tp decreases with the increase of the contact angle for all cases, but the decreasing extent of cases with λ < 1 is less than that of cases with λ > 1. It is easy to understand that the increasing hydrophobicity of the solid cylinder makes the movement of drops faster. For the droplet with λ < 1, tp of the droplet with a small viscosity is less than that of the droplet with a large viscosity, and a steep decrease of tp is observed as θ is varied from 30◦ to 40◦ . This is because when the contact angle is small (θ = 30◦ ), the daughter droplets can merge into a large droplet quickly, and with no continuous phase entrapped, the neck of the merged droplet is so thick that it requires more time to breakup. For a large contact angle (θ = 40◦ ), on the other hand, there is some

t=5.0

t=5.83

t=6.67

t=7.5

t=8.33

FIG. 10. (Color online) Time evolution of the droplet past a circular cylinder at θ = 70◦ , λ = 50. 043015-8

DEFORMATION AND BREAKUP OF A LIQUID DROPLET . . .

PHYSICAL REVIEW E 90, 043015 (2014)

9

9 λ=0.2 λ=0.02

8

8

7

tp

tp

7

6

6

5

5

4

4

3 20

λ=5 λ=50

40

60

80

100

θ

120

140

160

3 20

180

40

60

80

100

θ

(a)

120

140

160

180

(b)

FIG. 11. (Color online) Effect of the surface wettability on the passing time tp : (a) λ < 1, (b) λ > 1.

continuous phase entrapped that intensifies a break in the neck of the stretched droplet. For λ > 1, the passing time tp of the droplet with a small viscosity is not always larger than that of the droplet with a large viscosity. Besides, tp of the droplet with a large viscosity is more sensitive to the surface wettability variations. E. Effect of the Bond number

The Bond number (Bo) represents the relative importance of the gravitational force to the surface tension force. In this section we first discuss how Bo affects the passing time tp , a range of Bo changing from 0.3 to 1.54 is considered through varying the gravitational factor g. Figure 12 shows variations of tp in Bo. It is observed that, with the increase of Bo, tp decreases monotonically at different contact angles and viscosity ratios. Besides, we also found that the effects of surface wettability and viscosity ratios on the dynamics of

the droplet with a large Bo is smaller than that of the droplet with a small Bo. This is because, at a large Bo, the inertial force is more significant than the surface adhesion force and viscous force [33]. However, it should be noted that there is a great difference when Bo is equal to 0.301. This is because the inertia force is relatively small, and that the droplet dynamics, including the droplet shape and the way the droplet adhering to the solid, are significantly affected by the effects of the viscous force and adhesion force. Take the case of λ = 5 and θ = 70◦ as an example: for a fixed Bo, tp is less than the case of λ = 0.2 and θ = 70◦ , it is because some dispersed fluid are entrapped and the mixture becomes easy to break at the neck point. Furthermore, Bo also has a significant effect on the critical eccentric ratio βc , above which the droplet can pass the circular cylinder without breakup, and below βc , the droplet will be split into more daughter droplets. For the nonwetting droplet, we calculated βc with Bo under different viscosity ratios and Bond numbers and showed the results in Fig. 13. It can be

8

0.9

λ = 0.2, θ = 170 λ = 0.2, θ = 70

7

0.8

λ = 5, θ = 170 λ = 5, θ = 70

0.7 0.6

5

βc

tp

6

λ=0.2 λ=0.02 λ=5 λ=50

0.5

4

0.4

3 0.3

2 0.3

0.5

0.7

0.9 Bo

1.1

1.3

1.5

FIG. 12. Effect of Bo on the passing time tp .

0.2 0.2

0.4

0.6

0.8

Bo

1

1.2

1.4

FIG. 13. Effect of the Bo on the critical eccentric ratio βc . 043015-9

1.6

LI, CHAI, SHI, AND LIANG

PHYSICAL REVIEW E 90, 043015 (2014)

seen from this figure that βc increases with the increase of Bo. Also, for a fixed Bo, βc increases with viscosity ratio, which is more obvious when Bo is small. IV. CONCLUSIONS

We have adopted an improved interparticle-potential lattice Boltzmann multiphase model to investigate the complex deformation and breakup features of liquid droplet past a circular cylinder. The effects of eccentric ratio, viscosity ratio, surface wettability, and Bond number on the dynamics of the droplet are considered. Based on the present results, the main conclusions can be summarized as follows: (1) Two typical flow patterns of the droplet are obtained: breakup and no breakup, which is determined by the eccentric ratio β. In the breakup region, the relative size of the daughter droplets can be controlled by varying β, whereas the amount of deposited liquid on the tip of the circular cylinder is almost unaffected by β. (2) The viscosity ratio strongly affects the shape of the daughter droplets, breakup positions, and thickness of the deposited liquid film. Additionally, for case of a large viscosity ratio, the thickness of the deposited liquid decreases faster than that of a small viscosity ratio. (3) The wettability has an important effect on the dynamic behavior of the droplet, deposition place of liquid film, and the

[1] M. Marengo, C. Antonini, I. V. Roisman, and C. Tropea, Curr. Opin. Colloid Interface Sci. 16, 292 (2011). [2] F. M. Passandideh, M. Bussmann, and S. Chandra, Atomization Sprays 11, 397 (2001). [3] D. R. Link, S. L. Anna, D. A. Weitz, and H. A. Stone, Phys. Rev. Lett. 92, 054503 (2004). [4] C. Chung, K. H. Ahn, and S. J. Lee, J. Non-Newton. Fluid Mech. 162, 38 (2009). [5] C. Chung, M. Lee, K. Char, K. H. Ahn, and S. J. Lee, Microfluid Nanofluid 9, 1151 (2010). [6] S. Bakshi, I. V. Roisman, and C. Tropea, Phys. Fluids 19, 032102 (2007). [7] S. Proti`ere, M. Bazant, D. Weitz, and H. A. Stone, Europhys. Lett. 92, 54002 (2010). [8] L. Salkin, L. Courbin, and P. Panizza, Phys. Rev. E 86, 036317 (2012). [9] J. M. Gac and L. Grado´n, Colloids Surf. A: Physicochem. Eng. Aspects 441, 831 (2012). [10] G. Juarez, T. Gastopoulos, Y. Zhang, M. L. Siegel, and P. E. Arratia, Phys. Rev. E 85, 026319 (2012). [11] M. Do-Quang, A. Carlson, and G. Amberg, Fluid Dyn. Mater. Process 7, 389 (2011). [12] B. Yin and H. Luo, Comput. Fluids 82, 14 (2013). [13] S. Chen and G. D. Doolen, Annu. Rev. Fluid Mech. 30, 329 (1998).

passing time tp of the droplet passing cylinder. It is found the deposition place is changed from the upward to the downward surface of the circular cylinder with increasing wettability. Additionally, the evolutions of tp , as a function of the contact angle, are different for λ > 1 and λ < 1. When λ < 1, tp of the droplet with a small viscosity ratio is usually less than that of the droplet with a large viscosity ratio, whereas a different phenomenon is observed for the droplet with λ > 1. (4) The Bond number influences tp and the critical eccentric ratio βc . tp decreases monotonically with increasing Bo and is less sensitive to the viscosity ratio at a large Bo. And generally, βc increases in Bo, but the increment is not obvious under a large Bo and λ > 1.

ACKNOWLEDGMENTS

The authors would like to thank Dr. Qin Lou and Changsheng Huang for many helpful suggestions and discussion during the research. This study is supported by the National Natural Science Foundation of China (Grant Nos. 11272132, 51306133, 51006040), the National Funds for Distinguished Young Scientists of China (Grant No. 51125024), and the China Postdoctoral Science Foundation (Grant No. 2012M521424).

[14] C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech. 42, 439 (2010). [15] Z. Chai and T. Zhao, Acta Mech. Sin. 28, 983 (2012). [16] Q. Li, K. H. Luo, and X. J. Li, Phys. Rev. E 87, 053301 (2013). [17] X. Shan and H. Chen, Phys. Rev. E 47, 1815 (1993). [18] X. Shan and H. Chen, Phys. Rev. E 49, 2941 (1994). [19] Q. Kang, D. Zhang, and S. Chen, Phys. Fluids 14, 3203 (2002). [20] M. L. Porter, E. T. Coon, Q. Kang, J. D. Moulton, and J. W. Carey, Phys. Rev. E 86, 036701 (2012). [21] P. Lallemand and L. S. Luo, Phys. Rev. E 61, 6546 (2000). [22] K. N. Premnath and J. Abraham, Phys. Fluids 17, 122105 (2005). [23] M. E. McCracken and J. Abraham, Phys. Rev. E 71, 036701 (2005). [24] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E 65, 046308 (2002). [25] Z. Chai and T. S. Zhao, Phys. Rev. E 86, 016705 (2012). [26] X. He, X. Shan, and G. D. Doolen, Phys. Rev. E 57, R13 (1998). [27] N. S. Martys and H. Chen, Phys. Rev. E 53, 743 (1996). [28] Q. Li, K. H. Luo, and X. J. Li, Phys. Rev. E 86, 016709 (2012). [29] Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E 83, 036707 (2011). [30] D. Harvie, J. Cooper-White, and M. Davidson, J. NonNewtonian Fluid Mech. 155, 67 (2008). [31] A. Carlson, M. Do-Quang, and G. Amberg, Int. J. Multiphase Flow 36, 397 (2010). [32] M. Martinez and K. Udell, J. Fluid Mech. 210, 565 (1990). [33] P. Dimitrakopoulos, Phys. Fluids 19, 122105 (2007).

043015-10