Numerical simulations of vesicle and bubble dynamics in two-dimensional four-roll mill flows

Numerical simulations of vesicle and bubble dynamics in two-dimensional four-roll mill flows

Department of Mathematics, Chung-Ang University, Dongjakgu, Heukseokdong, Seoul, 156-756, Republic of Korea

Ming-Chih Lai

Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan

Yunchang Seol

National Center for Theoretical Sciences, No. 1, Sec. 4, Road. Roosevelt, National Taiwan University, Taipei 10617, Taiwan (Received 21 January 2017; revised manuscript received 29 March 2017; published 10 May 2017)

We use a computational technique based on the immersed boundary method to construct a four-roll mill device with which we can generate a broad spectrum of flow types from an extensional flow to a rotational one. We put a vesicle or a bubble in the constructed four-roll mill device to investigate their interaction with the surrounding fluid. The vesicle dynamics are determined by its bending rigidity, inextensibility, and hydrodynamical force, whereas the bubble dynamics is governed by the surface tension and the hydrodynamic interaction. Depending on the type of the flow, these suspended objects go through either a tank-treading motion or a tumbling motion. We validate our numerical method by a convergence study and discuss the transition between tank-treading and tumbling motions for the vesicles and bubbles.

DOI:10.1103/PhysRevE.95.053105

I. INTRODUCTION

We investigate the dynamics of a vesicle and a bubble in various flow types, which are generated by a two-dimensional computational model analogous to a four-roll mill device. It has been a long-time interest for many researchers in the fluid dynamics community to generate general flows such as rotational, extensional, and composite flows. In 1934, Taylor proposed a device [1] called the four-roll mill device. The device is equipped with four rotating mills at the four corners of a rectangular box and generates various flows from an extensional flow to a rotational flow by adjusting the rotational speed and direction of the four cylindrical rollers ω1and ω2; see Fig.1.

Since the four-roll mill device can generate a broad spectrum of flow types with a stagnation point at its center, it has been widely used in the rheological studies of suspended objects, such as the deformation and breakup of droplets [2–5], the coalescence of droplets [6], and the orientation and deformation of the polymer chains [7]. It attracts many researchers to develop a microfluidic analog of the four-roll mill [8,9], and several types of microfluidic devices [10–14] have been proposed to study the dynamics of single stretching DNA molecules [15,16], the flow-induced birefringence of semidilute solutions of wormlike micelles [17], and the drug-delivery systems and photonics [18–20].

Lee et al. [14] developed a microfluidic analog of the four-roll mill device that could generate all flow types from the pure extensional to pure rotational. Four pairs of inlet and outlet with the same width are arranged symmetrically around a central circular region. By adjusting the ratio of fluxes Q1/Q2, this single device can produce the entire range of flow patterns; see Fig.1. Since this microfluidic four-roll mill device enjoys a strong capability to produce various flows, it has been used

*_{[email protected]}

to study the tumbling dynamics of DNA [21], the transition among tank-treading, trembling, and tumbling motions of a vesicle [22], the complex dynamics of deforming droplets in the creeping flow [23], and the elongation and rotation of a droplet [24].

We here use computational simulations based on the immersed boundary (IB) method [25] to model a four-roll mill device with which we can generate various flows from an extensional flow to a rotational one. The IB method was developed to study flow patterns around heart valves and is a generally useful method for problems in which elastic materials interact with a viscous incompressible fluid. We put a vesicle or a bubble in the constructed four-mill device to investigate their interaction with the surrounding fluid. The dynamics of moving vesicles are determined by their membrane elasticity (bending resistance), inextensibility, and hydrodynamical force, whereas the bubble dynamics are governed by the surface tension and the hydrodynamic interaction. Thus, the IB method is a proper computational technique to model these elastic suspensions in a viscous fluid [26,27]. The novelty of the present work is that we use the unified IB framework to design a four-roll mill device and investigate the dynamics of vesicle and bubble in four-roll mill flows.

It has been known for a while that a vesicle in a shear flow has two types of motion: tank-treading (TT) and tumbling (TU). This transition is of fundamental importance since it can alter rheological properties of a vesicle solution by reducing dissipation [28,29]. The selection between these two types of vesicle motion depends on the reduced volume and the viscosity contrast between the interior and exterior fluids [30–33]. Even though we use the viscosity contrast one, we can observe both the tank-treading and tumbling motions by changing the flux ratio Q1/Q2and thus the flow type [22].

When the flux ratio increases in the negative direction, the magnitude of a vorticity tensor increases, which induces a rota-tional motion of the vesicle: either tank-treading or tumbling.

ω₁ ω₁ ω₂ ω₂ Q₁ Q₁ Q₂ Q₂ 100 μm 40 μm ω₁ ω₁ ω₂ ω₂ Q₁ Q₁ Q₂ Q₂ 100 μm 40 μm

FIG. 1. Schematic diagrams of four-roll mill device in Ref. [1] (left) and the microfluidic four-roll mill device designed here and in Ref. [14] (right). A broad spectrum of flow types can be generated by adjusting the angular velocity ratio ω1/ω2 (left) or the flux ratio

Q1/Q2(right).

Unlike the vesicle dynamics, we cannot observe a tumbling motion in the bubble dynamics. The bubble has only the surface tension, which does not resist the tangential force from the fluid, and thus the boundary of bubble follows simply the rota-tional flow without tumbling. The inextensibility of the vesicle, which resists both the tangential and normal force, plays a significant role in the tumbling motion of the vesicle. When the flux ratio increases in the positive direction, the magnitude of a deformation (or extensional) tensor increases, which induces an elongation and sometimes a breakup of a bubble [24].

In order to show that the four-roll mill device computa-tionally modeled by the IB method is a robust and efficient numerical tool to generate various flow types and to handle an elastic suspension in a viscous fluid, we first perform a convergence study, which confirms the consistent first-order convergence for the velocity field in L2 norm. Then we simulate a single vesicle or a bubble with various reduced areas in the different flows with various flux ratios. We discuss the transition between tank-treading and tumbling motions for the vesicles and the elongation and rotation of the bubbles.

II. MODEL OF VESICLE AND BUBBLE IN FOUR-ROLL MILL FLOWS

A. Equations of motion

We use the IB method [25] to simulate the dynamics of a vesicle or a bubble in various flows, which are generated by a two-dimensional device analogous to a four-roll mill [14]; see the left panel of Fig.2. Whereas the four-roll mill of Taylor [1] generates a desirable flow by controlling the angular velocities of the four rollers, our model based on the device proposed in Ref. [14] does it by varying the fluxes in the inlet and outlet regions.

Consider a two-dimensional square domain filled with an incompressible viscous fluid that contains internal fixed walls. In the IB framework, the fluid exists not only inside the internal walls but also outside the walls. In order to model the fixed walls, we utilize the “target boundary” idea. We designate target boundary points Y0(s) in the place where we want the internal walls to be. To avoid fluid leakage through the wall, the target boundary points should be spaced about half a mesh width apart (or closer). The target boundary points neither move nor interact with the fluid directly, but they are connected by a system of stiff springs to the immersed boundary points Y(s,t) that move at the local fluid velocity and apply force locally to the fluid; see Fig.3. When the wall boundary Y(s,t) moves apart from the target boundary Y0(s), a restoring force comes into play to keep them as close to each other as possible. The restoring force FK(s,t) acting on the immersed boundary

Y(s,t) is defined as

FK(s,t)= K[Y0(s)− Y(s,t)], (1) where K is a large stiffness constant. This provides a feedback mechanism for computing the boundary force needed to force the immersed boundary to stay in the stationary internal walls. This idea has been successfully used to simulate a stationary boundary in many previous works [34–36].

In order to generate the desired velocity in the domain of interest, which is enclosed by the fixed walls, we choose the

-200 -120 -20 20 120 200

x (

μm)

y (

μ

m)

x (

μm)

x (

μm)

-

Q

Ω

-

Q

Ω

Q

Ω

Q

Ω

FIG. 2. The flow generator analogous to the four-roll mill device. Poiseuille flows are derived with the normal fluxes inward to the domain being−Q1in 1and Q2in 2. The resulting velocity vector fields in the whole domain (middle panel) and in the central region (enclosed by

FIG. 3. Schematic view of wall target boundary, wall immersed boundary, and the interface immersed boundary. The wall IB points are connected to the target points by a system of stiff springs with zero rest length, and the interface IB points represent the elastic boundary.

regions 1and 2(shaded regions in the left panel of Fig.2) that are used to prescribe the desired velocity in the following way. Let 0be the union of the four shaded regions. The way of driving a flow in 0 is to apply an external force per unit area equal to

f0(x,t)=



α[u0(x,t)− u(x,t)], x ∈ 0= 1∪ 2

0, otherwise, (2)

where u0(x,t)= [u0(x,t),v0(x,t)] is the desired velocity and αis a constant. This force will become the body force in the fluid equations. When α is large, the fluid velocity u is driven rapidly toward u0within 0.

The tangential component of the desired velocity u0(x,t) to the wall is assumed to be zero in both 1 and 2. The normal components of the velocity u0(x,t) in 1 and 2are assumed to be a Poiseuille’s flow of parabolic profiles with the normal fluxes inward to the domain being−Q1in 1and Q2 in 2, see the flow profiles in the left panel of Fig.2, which also shows the resulting velocity vector fields in the whole domain (middle panel) and in the central region (enclosed by the dashed line in the left panel) containing vesicle or bubble (right panel).

We model a vesicle and a bubble as a closed curve in the two-dimensional domain; see Figs.2 and3. We assume that the fluid inside the vesicle and bubble is the same as the outside fluid. Whereas the elasticity of the vesicle is governed by the inextensibility and the bending resistance, that of the bubble is expressed only by the surface tension. Let X(s,t) be a configuration of a vesicle or bubble at time t, then the elastic forces FE(s,t) can be written, respectively, as

FE(s,t)= cs ∂ ∂s  ∂X ∂s   − 1τ  − cb ∂4X ∂s4, (3) FE(s,t)= ∂ ∂s(γ τ), where τ(s,t) = ∂X ∂s  ∂X_∂s, (4) where cs, cb, and γ are constants [26,27]. The first term of the right-hand side of Eq. (3) expresses the extension and compression resistance of the vesicle. Thus when cs is large,

the vesicle can be almost inextensible. The second term of the right-hand side of Eq. (3) expresses the bending resistance, and Eq. (4) represents the surface tension force.

Let the fluid velocity and pressure be denoted by u(x,t) and p(x,t), respectively, where x= (x,y) is a fixed Cartesian coordinate system and t is time. Then the mathematical formulations of this problem can be summarized as follows:

ρ  ∂u ∂t + u · ∇u  = −∇p + μ∇2_{u+ f + f0, (5)} ∇ · u = 0, (6) f(x,t)= FK(s,t)δ[x− Y(s,t)]ds + FE(s,t)δ[x− X(s,t)]ds, (7) ∂Y ∂t (s,t)= u[Y(s,t),t] = u(x,t)δ[x− Y(s,t)]dx. (8) ∂X ∂t (s,t)= u[X(s,t),t] = u(x,t)δ[x− X(s,t)]dx. (9) Equations (5) and (6) are the familiar Navier-Stokes equations for a viscous incompressible fluid with the constant fluid density ρ and viscosity μ. Equations (7)–(9) involve the two-dimensional Dirac δ function δ(x)= δ(x)δ(y), which expresses the local character of the interaction between forces and velocities along the immersed boundaries Y(s,t) and X(s,t). Equation (7) expresses the relation between the Eulerian force f(x,t)dx and the Lagrangian forces FK(s,t)ds and FE(s,t)ds. Equations (8) and (9) are the equations of motion of the immersed boundaries Y(s,t) and X(s,t), which says that the boundaries move at the local fluid velocity.

We impose the open (or natural) boundary condition for the fluid equations [Eqs. (5) and (6)] for all time on the boundary of the computational domain, i.e., (μ∇u − pI) n = 0, where I is the identity matrix and n is the unit normal vector to the computational boundary [37–39]. This condition is similar to the traction-free boundary condition, which implies that the flow moves in and out of the domain without any internal stress. Thus, in our case, while the fluid flows in and out of the domain on the regions 1and 2with the prescribed velocity through Eq. (2), the fluid moves naturally with no stress through the outflow tracks, which are the four open regions that are neither 1nor 2in Fig.2.

One might ask why we derive the flows in 1and 2 by imposing the body force in Eq. (2) with the open boundary condition, instead of simply applying velocity boundary conditions on the fluid equations. We try to keep the fluid solver as simple as possible. The boundary condition used in our numerical scheme enables us to use Fourier transform methodology to solve the fluid equations; see the Sec.II Bfor the numerical implementation. The direct use of the velocity boundary conditions, which are inhomogeneous along the computational boundary, makes the Fourier transform methods inapplicable and requires that an iterative method such as multigrid be used instead.

B. Numerical implementation

For the numerical implementation of the system Eqs. (1)–(9), we employ a first-order IB method. We use a

superscript to denote the time level; thus, Xn(s) is a shorthand for X(s,n t), where t is the size of the time step, and similarly for all other variables. The computational domain is discretized by a fixed uniform lattice of meshwidth h= x= y, and we use a staggered marker-and-cell (MAC) scheme [40] in which the fluid velocity field u are defined at the edges of the cells and the pressure p is defined at the cell centers. The elastic IB variables Xn_{(s) and F}n

E(s) are defined as functions of s with meshwidth s, and similarly for the wall IB variables Yn_{(s) and F}n

K(s). Our goal is to compute un+1_{, X}n+1_{, and Y}n+1_{from given data u}n_{, X}n_{, and Y}n_. The step-by-step procedure of the numerical implementation proceeds as follows.

Step 1. Calculate the force densities FnK(s) and F n E(s) defined on the IB points Xn_{(s) and Y}n_{(s), respectively. For} the detailed description of the discrete force densities; see Refs. [26,27].

Step 2. Distribute these force densities defined on La-grangian grid points into the force at Eulerian spatial grid points to be applied in the Navier-Stokes equations. This is done by a discretization of Eq. (7) as

fn(x)= s

FnE(s) δh[x− Xn(s)] s, (10) where x= (x1,x2) is the fluid grid points, and δh(x)=

1 h2φ(

x1

h)φ( x2

h) with φ(r) of the form

φ(r)= ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 3−2|r|+√1+4|r|−4r2 8 , if|r|1 5−2|r|−√−7+12|r|−4r2 8 , if 1|r|2 0, if 2|r|. (11)

The motivation and derivation for this particular choice is discussed in Ref. [25], and we spread the force density FnK(s) in the same fashion.

Step 3. Given the computed force density fn_{(x), we solve the} discretized version of the fluid Eqs. (5) and (6) to update the velocity and pressure. To do that, we employ the first-order incremental pressure-correction projection method [37,38], which can be summarized as follows:

ρ  ˜un+1_{− u}n t + (u n_{· ∇} h)un  + ∇hpn= μ h˜un+1+ fn, (μ∇h˜un+1− pnI) n= 0 on ∂, (12) hp˜n+1= ρ t∇h· ˜u n_+1, with p˜n+1 = 0 on ∂, (13) un+1 = ˜un+1− t ρ ∇hp˜ n+1_, and pn+1= ˜pn+1+ pn, (14) where the discrete operators∇h,∇h·, and hare the standard second-order finite differences in staggered MAC grid to approximate the gradient, divergence, and Laplace operators, respectively. One can see that we need to solve one Poisson equation for the pressure-like variable ˜pn+1_{and two modified} Helmholtz equations for the velocity field ˜un+1_{. Fortunately,} these elliptic solvers can be implemented efficiently using fast Fourier transform (FFT) [41].

Step 4. Once un+1(x) is found in the fluid grid points, we can update the boundary configuration Xn+1_{(s) by the following} interpolation of un+1_(x):

Xn+1(s)= Xn(s)+ t

x

un+1(x) δh[x− Xn(s)] h2. (15) The wall boundary Yn+1_{(s) can be similarly updated, and this} completes the time step.

III. RESULTS

A. Flow generations and convergence study

Throughout this paper, we choose the computational do-main [−200,200] × [−200,200] μm2_{, which is filled with} a viscous incompressible fluid with the fluid viscosity μ= 0.01 g/(cm s) and the density ρ = 1.0 g/cm3_{. The flux} Q2 in the region 2 is fixed as 2.5× 105 μm2/s, and we vary the flux Q1 in the region 1 to generate various flow types. The mesh width is h= x = y = 400/512 μm, which is uniform and fixed all the time, and the timestep is t = 5.0 × 10−7s. The grid spacing s of the IB points for the wall and elastic boundaries is around 0.36 μm, which is less than h/2. The spring constant in Eq. (1) is set equal to K= 8 × 104_dyne/cm2_{. This choice of K makes the} max-imum deviation of the immersed boundary from the target boundary to be less than h/20; i.e., ||Y0(s)− Y(s,t)||_∞< h/20. The constant α in Eq. (2) is chosen to be 2× 10−6 dyne· s/μm3_{, which makes the error of the computed velocity} u from the desired velocity u0to be less than 4%, in the region 1∪ 2.

It is well known that the flow type in the central cavity of the four-roll mill is dependent on the flux ratio f = Q1/Q2 which goes from−1 to 1. The flow type can be assessed by the flow-type parameter ξ defined as

ξ = |E| − ||

|E| + ||, (16)

where|E| and || are the magnitudes of the deformation and vorticity tensors, respectively. In an ideal four-roll mill device, this parameter varies from−1 (rotational) to 1 (extensional).

Figure 4 shows three different types of flows; the exten-sional flow with f = 1.0 (left panels), the shear flow with f = −0.585 (middle panels), and the rotational flow with f = −1.0 (right panels). The first row shows the contours of the flow-type parameters ξ , the second and third rows draw the streamlines of the velocity in [−130,130] × [−130,130] μm2 and in the central region [−60,60] × [−60,60] μm2_, respec-tively, with some indications of velocity vectors. Figure4 is comparable to Fig. 3 in Ref. [14] in which the authors used a pseudo-3D model to investigate the dependence of the flow type on the flux ratio.

Figure5depicts the averaged flow-type parameter ¯ξ over a central disk with the radius 12.5 μm in terms of the flux ratio f . We have found that the flow-type parameter ξ deviates no more than 0.05 from the averaged value ¯ξ in this central region, which implies that the flow-type parameter ξ is almost homogeneous near the central stagnation point. We can see from Fig. 5 that the flow type varies from purely rotational (ξ = −1) to purely extensional (ξ = 1) by changing the flux

FIG. 4. The contours of the flow-type parameters (first row), the streamlines of the velocity field in [−130,130] × [−130,130] μm2_(second

row), and the streamlines with some indications of the velocity vectors in [−60,60] × [−60,60] μm2_{(third row). The columns from left to right}

show an extensional flow with the flux ratio f= 1.0, a mixed shear flow with f = −0.585, and a rotational flow with f = −1.0, respectively. ratio over the range of−1  f  1. Figure5is qualitatively

the same as Fig.2in Ref. [14]; however, our model achieves the pure shear flow (ξ = 0) at f = −0.585, while the shear flow appears at f = −0.763 in Ref. [14].

For the verification of our numerical method, we perform a convergence study of the computed solutions. We put a vesicle of a shape shown in Fig.2 in the center of the domain and generate flows with three different flux ratios f = 1.0, −0.585, and −1.0. Then we vary the mesh size of the domain as N = 128, 256, 512, and 1024 so that the mesh width becomes h= 400/N μm, correspondingly. We also choose t proportional to h with t = 6.4 × 10−7h, so that the

refinements for the fluid mesh width, the boundary mesh width, and the time-step duration are scaled by the same factor. When we refine the mesh width and time step, we increase the spring constant in Eq. (1) as K= 2.0 × (25N/64)2_dyne/cm2_.

Since we do not know the exact solution of the problem, the estimation of a convergence ratio requires three numerical so-lutions for three consecutive N ’s. Let (uN,vN) be the velocity field, and let|| · ||2be the L2norm. Figure6shows the con-vergence ratios (||uN− u2N||2₂+ ||vN− v2N||2₂)1/2/(||u2N− u_4N||2

2+ ||v2N− v4N||22)1/2 as functions of time for each of the cases N = 128 (line with “+”) and 256 (line with “◦”). The flux ratios are f = 1.0 (left), −0.585 (middle), and −1.0

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 flux ratio f=Q 1/Q2

flow type parameter

ξ

FIG. 5. Flow-type parameter ¯ξ averaged over a central disk with the radius 12.5 μm in terms of the flux ratio f . The flow type varies from purely rotational (ξ= −1) to purely extensional (ξ = 1) by changing the flux ratio over the range of −1  f  1. When the flow-type ratio f = −0.585, the flow-type parameter ξ is almost 0.

0 0.1 0.2 1.6 1.8 2 ratio time (s) 0 0.1 0.2 1.6 1.8 2 time (s) 0 0.1 0.2 1.6 1.8 2 time (s)

FIG. 6. Convergence ratios of the computed velocity field u(x,t) of the fluid. The convergence ratios (defined in the text) are plotted as functions of time for three different flux ratios: f= 1.0 (left), −0.585 (middle), and −1.0 (right). In each panel, the line with “+” is the convergence ratio obtained using the grids N= 128,256,512, and the line with “◦” is the ratio obtained with N = 256, 512, 1024. The convergence ratio is near two (first-order accuracy) for all the different cases.

(right). The convergence ratio 2 implies that the scheme has first-order accuracy, which is typical of the IB method as applied to problems with thin elastic boundaries [42,43].

B. Vesicle dynamics

We now put a vesicle in the center of the four-roll mill device as shown in Fig. 2 to investigate the vesicle motion in various flow types. The initial configuration of the vesicle is a closed curve with several different shapes. We choose the vesicles so that they have the same perimeter L= 2πR0, where R0= 10 μm, but with different internal areas A, which can define different reduced areas V = _{π R}A2

0. We use the same

computational and physical parameters as those used in the

previous section for the four-roll mill device. The stiffness and bending coefficients used in Eq. (3), which represents the elasticity of the vesicle, are cs = 32 dyne/cm and cb = 10−10 dyne· cm2_{, respectively. With these values, the relative error} of the vesicle length is found to be less than 2× 10−4.

It is known that a vesicle in a simple shear flow (ξ = 0) undergoes a tank-treading tangential motion at its equilibrium configuration [29–33]. As shown in Fig. 5, when the flux ratio f is around −0.6, the four-mill device can generate a simple shear flow with the flow-type parameter ξ = 0. If the flux ratio f is away from 1 and thus the flow is not dominantly extensional, we can also expect a tank-treading motion of the vesicle, since there is a nonzero rotational tensor of the fluid. Figure7shows the motions of the vesicle

y ( μ m) −20 0 20 −20 −10 0 10 20 −20 0 20 −20 −10 0 10 20 −20 0 20 −20 −10 0 10 20 y ( μ m) x (μm) −20 0 20 −20 −10 0 10 20 x (μm) −20 0 20 −20 −10 0 10 20 x (μm) −20 0 20 −20 −10 0 10 20

FIG. 7. The motion of the vesicle with the reduced area V = 0.7, together with the streamlines of the velocity field at t = 0.4 s. The flows are generated with different values of the flux ratio: f= 0.9 (upper-left), 0.4 (upper-middle), 0.0 (upper-right), −0.4 (lower-left), −0.6 (lower-middle), and−0.67 (lower-right). Each vesicle reaches a stable tank-treading (TT) motion with some fixed inclination angle between the long axis of the vesicle and the negative x axis.

−0.80 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 10 20 30 40 flux ratio f

inclination angle (deg)

V=0.5 0.6 0.7 0.8 0.9 −0.80 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 flux ratio f TT frequency (Hz)

FIG. 8. The inclination angles (upper) and the angular velocities (lower) of the TT motion of the vesicles in terms of the flux ratio f for various reduced areas V when the vesicles are in the steady TT motions.

with the reduced area V = 0.7, together with the streamlines of the velocity field at t= 0.4 s. The flows are generated with the flux ratio f = 0.9 (upper-left), 0.4 (upper-middle), 0.0 (upper-right),−0.4 (lower-left), −0.6 (lower-middle), and −0.67 (lower-right). Each vesicle in this figure reaches a stable TT motion with some fixed inclination angle between the long axis of the vesicle and the negative x axis. When we decrease the flux ratio further, we have observed the tumbling motion of the vesicle, which is discussed below.

When a vesicle is in a steady TT state, the motion of the vesicle can be characterized by both the inclination angle between the long axis of the vesicle and the flow direction, and the angular frequency of the TT motion. These two values have been found to be strongly dependent on the reduced area V and the flux ratio f . Figure 8 draws the inclination angle (upper) and the angular frequency (lower) in terms of the flux ratio f for various reduced areas V when the vesicles are in their steady states. As the flux ratio f increases to 1, independent of the reduced area of the vesicle, the inclination angle increases and converges around to 39◦, and the angular frequency of the TT motion decreases to 0. For a fixed flux ratio, the inclination angle and the angular frequency both increase as the reduced area V increases, which agrees well with the results in Refs. [27,44,45].

It has been known for a while that a vesicle in a shear flow has two types of motion: TT and TU. The selection between these two types of vesicle motion depends on the reduced volume (a flattened vesicle would tumble more easily than a swollen one) and the viscosity contrast (the more viscous the internal fluid in comparison to the external one

is, the easier is the tumbling) [29–33]. Even though we use the viscosity contrast 1, we can also observe a tumbling motion by decreasing the flux ratio f and thus by increasing the magnitude of the rotational tensor ||. As we can see from Fig.8, the inclination angle decreases as the flux ratio decreases, and, at some critical value fcof the flux ratio, the motion of the vesicles changes from TT to TU. Figure9shows the critical flux ratio fcfor the transition between TT and TU in terms of the reduced area V . As the reduced area gets smaller, the critical flux ratio becomes larger, i.e., the vesicle with small reduced area needs small rotational tensor for the tumbling motion, which implies that the flattened vesicle would tumble more easily than a swollen one.

When the flux ratio f gets smaller than the critical flux ratio fc, the vesicles get into a TU state. The left panels of Fig.10 depict the configurations of a vesicle with the reduced area V = 0.7 suspended in a flow with the flux ratio f = −0.8 at some chosen times, which clearly shows the tumbling motion of the vesicle. The right panel draws the time evolution of the inclination angle of the vesicle with the reduced area V = 0.7 in flows with three different flux ratios: f = −0.7 (dotted line),−0.8 (dashed line), and −0.9 (solid line). We can see that, when the flux ratio decreases, the tumbling frequency increases.

When a vesicle is in a steady tumbling state, the tumbling frequency does not only depend on the flux ratio f as shown in Fig.10, but it also depends on the reduced area of the vesicle. Figure11shows the tumbling frequency of the motion of the vesicle in terms of the flux ratio f for various reduced areas V . The tumbling frequency increases with the decreasing flux

0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 −0.8 −0.75 −0.7 −0.65 −0.6 reduced area V

critical flux ratio f

c

TU

TT

FIG. 9. The critical flux ratio fcfor the transition between TT and TU in terms of the reduced area V .

t=0.16 s t=0.32 s t=0.48 s t=0.64 s t=0.8 s t=0.96 s 0 1 2 3 4 −50 0 50 time (s)

inclination angle (deg)

FIG. 10. The configuration of a vesicle with the reduced area V= 0.7 suspended in a flow with the flux ratio f = −0.8 at some chosen times (left panels) and the time evolution of the inclination angle (right panel) of the vesicle with the reduced area V= 0.7 in flows with three different flux ratios: f = −0.7 (dotted line), −0.8 (dashed line), and −0.9 (solid line).

−10 −0.95 −0.9 −0.85 −0.8 −0.75 −0.7 −0.65 −0.6 0.2 0.4 0.6 0.8 1 flux ratios f TU frequency (Hz) V=0.5 0.6 0.7 0.8 0.9

x (μm) y ( μ m) −20 0 20 −20 0 20 x (μm) −20 0 20 −20 0 20 −10 0 10 −10 0 10 θ x (μm) 0 0.5 1 1.5 2 2.5 3 3.5 4 −200 −100 0 100 200 time (s) θ (deg)

FIG. 12. The motion of a bubble with the surface tension γ= 0.0004 dyne at t = 0.8 s (upper-left), and t = 2.4 s (upper-middle), together with the streamlines of the velocity field. The upper-right panel describes the velocity vectors on the bubble. The flow is purely rotational with the flux ratio f = −1. The lower panel shows the rotational angle, which indicates that the bubble rotates uniformly with the angular frequency ω= 0.8928 Hz.

ratio for all types of the vesicle. This is because the magnitude of rotational tensor || is inversely proportional to the flux ratio f . For a fixed flux ratio, as the reduced area gets smaller, the tumbling frequency increases. We can again conclude that the flattened vesicle would tumble more easily than a swollen one.

We here have observed the transition between TT and TU depending on the flux ratio. Recently, another stable state has been found in which the vesicle deforms significantly and the inclination angle oscillates with the inclination angle being less than π/2 [22,46,47]. This is called a trembling state (TR) or vacillating-breathing state, and the selection among the three states TT, TR, and TU depends on the viscosity contrast, the reduced area, and the bending rigidity of vesicle, the shear rate, and so on [48,49]. In our model, however, we cannot observe the trembling motion of the vesicle, even though we change the bending rigidity from cb= 10−9 dyne· cm2 to 10−15 dyne· cm2. We suspect that the trembling motion of a vesicle would occur only in 3D flow. The fluid flow in the third direction perpendicular to the 2D plane of the main flow enables the vesicle to tremble in an oscillatory manner with a large shape deformation with combination of the elasticity (bending rigidity) of the vesicle. Our 2D model does not have these effects of the 3D flow; see Ref. [50], which also argues that there is no TR of a vesicle in 2D flow.

C. Bubble dynamics

We here put a bubble in the center of the four-roll mill device as shown in Fig.2and investigate the bubble dynamics in various flow types. Since the bubble dynamics is governed by the surface tension, the motions of the bubbles with different reduced areas are almost the same, unlike the vesicles in

the previous section. Thus, we choose one type of bubble, which has the same initial configuration as the vesicle with the reduced area V = 0.7. While we use the same computational and physical parameters as those used in the previous sections, we vary the surface tension γ used in Eq. (4) from 0.0002 to 0.2 dyne to see the dependence of the bubble dynamics in various flows on the surface tension.

The motion of a bubble depends on the magnitudes of the deformation tensor|E| and the vorticity tensor ||, and the surface tension γ . When the flow is purely rotational (ξ = −1) with the flux ratio f = −1, the bubble keeps the circular shape and rotates, independent of the surface tension. Figure 12 shows the motion of a bubble with the surface tension γ = 0.0004 dyne at t = 0.8 s (upper-left) and t = 2.4 s (upper-middle), together with the streamlines of the velocity field. The upper-right panel shows that the velocity vectors on the bubble are all tangential, which indicates the rotational motion of the bubble. The lower panel shows the angle of the line connecting the bubble center and one chosen point on the bubble (circle on the bubble) from the negative x axis; see the upper-right panel of Fig.12. We can see that the bubble rotates uniformly with the angular frequency ω= 0.8928 Hz. When the surface tension varies between 0.0002 and 0.2 dyne, we have observed that the motions of the bubbles are almost the same with the same angular frequency independent of the surface tension. The only difference is that the initial transient time to reach the steady rotational motion decreases as the surface tension increases.

When the flow is purely extensional (ξ = 1) with the flux ratio f = 1, the bubble can extend with the inclination angle being fixed at around π/4. The deformation of droplets in the extensional flow has long been studied [2–4,24,51,52].

−50 0 50 −50 0 50 y ( μ m) −50 0 50 −50 0 50 −50 0 50 −50 0 50 −20 0 20 −20 −10 0 10 20 x (μm) y ( μ m) x (μm) −20 0 20 −20 −10 0 10 20 −20 0 20 −20 −10 0 10 20 x (μm)

FIG. 13. The motion of bubbles with the surface tensions γ= 0.0004 dyne (upper) and γ = 0.004 dyne (lower) at t = 0.032 s (left), 0.064 s (middle), and 0.096 s (right). The flow is purely extensional (ξ= 1) with the flux ratio f = 1. The bubble with a small surface tension cannot resist the extensional deformation of the flow and elongates in time, gets thinner, and finally breaks up (top panels). The bubble with a large surface tension approaches to a stable elliptical shape.

Figure 13 shows the motions of the bubbles with two different surface tension γ = 0.0004 dyne (upper) and γ = 0.004 (lower) at t= 0.032 s (left), 0.064 s (middle), and

0.096 s (right). The lower-middle panel of Fig.13draws the streamlines of the velocity field, and the lower-right panel indicates the velocity vectors on the bubble. When the surface

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 101 102 time (s) long axis ( μ m) 0.0004 dyne 0.0008 dyne 0.001 dyne 0.0012 dyne 0.004 dyne 0.008 dyne 0.012 dyne 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 8 9 10 11 12 13

surface tension (dyne)

long axis (

μ

m)

FIG. 14. The length of the long axis of the elliptical bubble in terms of time for various surface tensions (upper) and the constant value of the long axis in the stable state of the bubble in terms of the surface tensions (lower).

y ( μ m) −20 0 20 −20 −10 0 10 20 −20 0 20 −20 −10 0 10 20 −20 0 20 −20 −10 0 10 20 x (μm) y ( μ m) −20 0 20 −20 −10 0 10 20 x (μm) −20 0 20 −20 −10 0 10 20 −20 0 20 −20 −10 0 10 20 x (μm)

FIG. 15. The motion of bubbles with the surface tensions γ = 0.0004 dyne (upper) and γ = 0.004 dyne (lower) at t = 0.4 s (left), 0.8 s (middle), and 1.2 s (right). The flow is composite with extensional and rotational flows with the flux ratio f = −0.6.

tension is too small, the bubble cannot resist the extensional deformation of the flow and elongates in time, gets thinner, and finally breaks up (Fig.13, top panels) [2]. When the bubble has a surface tension large enough to resist the extensional deformation, it approaches to a stable elliptical shape (Fig.13, bottom panels).

Since we fix the fluxes Q1 and Q2 with the ratio Q1/Q2= 1, the magnitude of the deformation tensor is a constant, and thus the amount of the elongation of the bubble is inversely proportional to the surface tension. The upper panel of Fig. 14 shows the length of the long axis of the elliptical bubble in terms of time for various surface tensions.

0 0.5 1 1.5 2 2.5 3 3.5 4 8 10 12 14 time (s) long axis ( μ m) 0.0002 dyne 0.0004 dyne 0.0008 dyne 0.002 dyne 0.004 dyne 0.016 dyne 0 0.5 1 1.5 2 2.5 3 3.5 4 −200 −100 0 100 200 time (s) θ (deg)

FIG. 16. The length of the long axis of the elliptical bubble (upper) and the rotational angle (lower) in terms of time for various surface tensions. The flow is composite with extensional and rotational flows with the flux ratio f = −0.6.

0 0.002 0.004 0.006 0.008 0.01 0.012 8 9 10 11 12

surface tension (dyne)

long axis ( μ m) f=−1.0 f=−0.9 f=−0.8 f=−0.6 f=0.0 f=0.4 0 0.002 0.004 0.006 0.008 0.01 0.012 0.2 0.4 0.6 0.8

surface tension (dyne)

angular velocity (Hz)

FIG. 17. The length of the long axis of the elliptical bubble (upper) and the rotational angle (lower) in terms of the surface tension for various surface tensions when the bubbles are in steady rotational states.

We have found that, when the surface tension is smaller than 0.0012 dyne, the bubble elongates with the increasing long axis in time and eventually breaks up. The bubble with the surface tension from 0.0012 dyne quickly approaches to the stable state with a constant length of long axis, see the upper panel. The constant value of the long axis in the stable state of the bubble decreases as the surface tension increases, as indicated in the lower panel of Fig.14[2,24].

When the flux ratio is larger than−1 and smaller than 1, the flow is composite with extensional and rotational flows, as shown in Fig.5. Thus, the bubble in these flows goes through a steady motion with both rotation and elongation depending on the surface tension. Figure15shows the motion of bubbles with the surface tensions γ = 0.0004 dyne (upper) and γ = 0.004 dyne (lower) at t= 0.4 s (left), 0.8 s (middle), and 1.2 s (right). The streamlines of the velocity field and the velocity vectors on the bubble are also drawn in the first two columns and in the third column, respectively. We can observe that the bubble with a smaller surface tension elongates more than the bubble with a larger surface tension. The bubbles experience the rotational tank-treading motion, as we can see from the position (marked by◦) of one chosen point of the bubble and from the velocity vectors on the bubble.

Figure16depicts the length of the long axis of the elliptical bubble (upper) and the rotational angle (lower) in terms of time for various surface tensions. The flow is composite with extensional and rotational flows with the flux ratio f = −0.6. We can see from the figure that the bubbles spend some transient time to reach a steady state when they have elliptical shapes and rotate with constant angular velocities. As the

surface tension increases, the long axis gets shorter, and the angular frequency becomes larger. Note that, unlike the vesicle dynamics in the previous section, the rotational motion in the bubble dynamics is always tank-treading without tumbling motion even when the flux ratio f is−1. This is because the bubble has only the surface tension, which does not resist the tangential force from the fluid, and thus the boundary of bubble follows simply the rotational flow without tumbling.

Even though the length of the long axis is inversely proportional and the angular frequency is proportional to the surface tension of the bubble when it is small, the changes in the long axis and the angular frequency gets smaller when the surface tension increases more above some value, compare the graphs for γ = 0.002, 0.004, and 0.016 dyne in Fig.16. When we change the flow type by adjusting the flux ratio f , we have found the same behavior as shown in Fig.17, which shows the length of the long axis of the elliptical bubble (upper) and the rotational frequency (lower) in terms of the surface tension for various flow types when the bubbles are in steady rotational states. When the surface tension increases up to about γ = 0.001 dyne, the long axis gets shorter, and the angular frequency becomes larger. However, the surface tension larger than γ = 0.001 dyne makes only a small difference in the long axis and no difference in the angular frequency.

IV. CONCLUSIONS

We have introduced a computational model for a four-roll mill device and used it to investigate the dynamics of a vesicle and a bubble in four-roll mill flows. Using the immersed

boundary method, we have showed that our computational model is a robust and efficient numerical tool to generate a broad spectrum of flow types by comparing the computed flow types with those in literature and by convergence study.

Then we have simulated a single vesicle with various reduced areas, which is positioned in the central cavity of the four-roll mill device and observed the transition of the vesicle motions between TT and TU. The critical flux ratio for the transition is different for the vesicles with different reduced areas. In the TT state, the rotational frequency decreases and the inclination angle increases as the flux ratio increases for all types of vesicle. In the TU state, the tumbling frequency is also inversely proportional to the flux ratio. In both cases of TT and TU, a flattened vesicle with a small reduced area rotates faster than a swollen one with a large reduced area.

We have also simulated a single bubble in four-roll mill flows. Unlike the vesicle dynamics, the rotational motion in the bubble dynamics is always tank-treading without tumbling, since the bubble has only the surface tension, which does not resist the tangential force from the fluid. The inextensibility of the vesicle, which resists both the tangential and normal force, plays a significant role in the tumbling motion of the vesicle. When the bubble has a small surface tension, it can elongate

and break up sometimes in a flow with the flux ratio close to 1. The bubble with a large surface tension keeps its length of long axis and angular frequency of TT motion for a fixed flux ratio. Our simulation results have confirmed the existence of the transition of the vesicle motions between TT and TU depending on the flux ratio; however, there is no TR state, which was found recently in literature. We conclude that the TR motion of a vesicle would occur only in 3D flow. Thus, the extension of the present model to 3D is needed to explore more comprehensively the vesicle dynamics, which will be a future work. The transition could also depend on the difference of the fluid properties inside and outside the vesicle. The investigation of the vesicle dynamics with various viscosity contrasts and various nonuniform densities will also be the subject of future work.

ACKNOWLEDGMENTS

Y. Kim was supported by National Research Foundation of Korea Grant funded by the Korean Government (Grant No. 2015R1A2A2A01005420). The work of M.-C. Lai is supported in part by Ministry of Science and Technology of Taiwan under research Grant No. MOST-104-2115-M-009-014-MY3 and NCTS.

[1] G. I. Taylor, The formation of emulsions in definable fields of flow,Proc. R. Soc. London A 146,501(1934).

[2] B. J. Bentley and L. G. Leal, An experimental investigation of drop deformation and breakup in steady, two-dimensional linear flows,J. Fluid Mech. 167,241(1986).

[3] H. A. Stone, B. J. Bentley, and L. G. Leal, An experimental study of transient effects in the break up of viscous drops, J. Fluid Mech. 173,131(1986).

[4] H. A. Stone, Dynamics of drop deformation and breakup in viscous fluids,Annu. Rev. Fluid Mech. 26,65(1994). [5] Y. T. Hu and A. Lips, Estimating Surfactant Surface Coverage

and Decomposing its Effect on Drop Deformation,Phys. Rev. Lett. 91,044501(2003).

[6] H. Yang, C. C. Park, Y. T. Hu, and L. G. Leal, The coalescence of two equal-sized drops in a two-dimensional linear flow,Phys. Fluids 13,1087(2001).

[7] E. C. Lee and S. J. Muller, Flow light scattering studies of polymer coil conformation in solutions in extensional flow, Macromolecules 32,3295(1999).

[8] C. J. Pipe and G. H. McKinley, Microfluidicrheometry,Mech. Res. Commun. 36,110(2009).

[9] X. Hu, P. E. Boukany, O. L. Hemminger, and L. Lee, The use of microfluidics in rheology,Macromol. Mater. Eng. 296,308 (2011).

[10] H. P. Babcock, R. E. Teixeira, J. S. Hur, E. S. G. Shaqfeh, and S. Chu, Visualization of molecular fluctuations near the critical point of the coil-stretch transition in polymer elongation, Macromolecules 36,4544(2003).

[11] S. D. Hudson, F. R. Phelan, M. D. Handler, J. T. Cabral, K. B. Migler, and E. J. Amis, Micro-fluidic analog of the four-roll mill,Appl. Phys. Lett. 85,335(2004).

[12] F. R. Phelan, S. D. Hudson, and M. D. Handler, Fluid dynamics analysis of channel flow geometries for materials

characterization in microfluidic devices, Rheol. Acta 45, 59 (2005).

[13] X. Hu, S. Wang, Y. J. Juang, and L. Lee, Five-cross microflu-idic network design free of coupling between electrophoretic motion and electro-osmoticflow,Appl. Phys. Lett. 89,084101 (2006).

[14] J. Lee, R. Dylla-Spears, N. P. Teclemariam, and S. J. Mullera, Microfluidic four-roll mill for all flow types,Appl. Phys. Lett.

90,074103(2007).

[15] T. T. Perkins, D. E. Smith, and S. Chu, Single polymer dynamics in an elongational flow,Science 276,2016(1997).

[16] C. M. Schroeder, H. P. Babcock, E. S. G. Shaqfeh, and S. Chu, Observation of polymer conformation hysteresis in extensional flow,Science 301,1515(2003).

[17] J. A. Pathak and S. D. Hudson, Rheo-optics of equilibrium polymer: Solutions: Wormlike micelles in elongational flow in a microfluidic cross-slot,Macromolecules 39,8782(2006). [18] G. M. Whitesides, The origins and the future of microfluidics,

Nature 442,368(2006).

[19] S.-Y. Teh, R. Lin, L.-H. Hung, and A. P. Lee, Droplet microfluidics,Lab Chip 8,198(2008).

[20] J.-T. Wang, J. Wang, and J.-J. Han, Fabrication of advanced particles and particle-based materials assisted bydroplet-based microfluidics,Small 7,1728(2011).

[21] J. S. Lee, E. S. G. Shaqfeh, and S. J. Muller, Dynamics of DNA tumbling in shear to rotational mixed flows: Pathways and periods,Phys. Rev. E 75,040802(2007).

[22] J. Deschamps, V. Kantsler, E. Segre, and V. Steinberg, Dynamics of a vesicle in general flow,Proc. Natl. Acad. Sci. USA 106, 11444(2009).

[23] Y.-N. Young, J. Bławzdziewicz, V. Cristini, and R. H. Goodman, Hysteretic and chaotic dynamics of viscous drops in creeping flows with rotation,J. Fluid Mech. 607,209(2008).

[24] J. Wang, J. Han, and D. Yu, Numerical studies of geometry effects of a two-dimensional microfluidic four-roll mill on droplet elongation and rotation,Eng. Anal. Bound. Elem. 36, 1453(2012).

[25] C. S. Peskin, The immersed boundary method,Acta Numer. 11, 479(2002).

[26] Y. Kim, M.-C. Lai, and C. S. Peskin, Numerical simulations of two-dimensional foam by the immersed boundary method, J. Comput. Phys. 229,5194(2010).

[27] Y. Kim and M.-C. Lai, Simulating the dynamics of inextensible vesicles by the penalty immersed boundary method,J. Comput. Phys. 229,4840(2010).

[28] F. Rioual, T. Biben, and C. Misbah, Analytical analysis of a vesicle tumbling under a shear flow,Phys. Rev. E 69,061914 (2004).

[29] J. Beaucourt, F. Rioual, T. Seon, T. Biben, and C. Misbah, Steady to unsteady dynamics of a vesicle in a flow,Phys. Rev. E 69, 011906(2004).

[30] S. R. Keller and R. Skalak, Motion of a tank-treading el-lipsoldal particle in a shear flow, J. Fluid Mech. 120, 27 (1982).

[31] M.-A. Mader, V. Vitkova, M. Abkarian, A. Viallat, and T. Podgorski, Dynamics of viscous vesicles in shear flow, Eur. Phys. J. E 19,389(2006).

[32] M.-A. Mader, H. Ez-Zahraouy, C. Misbah, and T. Podgorski, On coupling between the orientation and the shape of a vesicle under a shear flow,Eur. Phys. J. E 22,275(2007).

[33] Y. Kim and M-C. Lai, Numerical study of viscosity and inertial effects on tank-treading and tumbling motions of vesicles under shear flow,Phys. Rev. E 86,066321(2012).

[34] D. Goldstein, R. Handler, and L. Sirovich, Modeling a no-slip flow with an external force field,J. Comput. Phys. 105, 354 (1993).

[35] C. Lee, Stability characteristics of the virtual boundary method in three-dimensional applications, J. Comput. Phys. 184,559 (2003).

[36] Y. Kim, J. Lee, and S. Lim, Nodal flow simulations by the immersed boundary method, SIAM J. Appl. Math. 74, 263 (2014).

[37] J. L. Guermond, P. Minev, and J. Shen, An overview of projection methods for incompressible flows,Comput. Methods Appl. Mech. Eng. 195,6011(2006).

[38] J. Liu, Open and traction boundary conditions for the incom-pressible Navier-Stokes equations,J. Comput. Phys. 228,7250 (2009).

[39] A. Poux, S. Glockner, and M. Azaiez, Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods,J. Comput. Phys. 230,4011(2011).

[40] F. H. Harlow and J. E. Welsh, Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface,Phys. Fluids 8,2181(1965).

[41] M. Frigo and S. G. Johnson, The Design and Implementation of FFTW3,Proc. IEEE 93,216(2005).

[42] B. E. Griffith and C. S. Peskin, On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems, J. Comput. Phys. 208, 75 (2005).

[43] Y. Kim and C. S. Peskin, Penalty Immersed Boundary Method with an Elastic Boundary with Mass,Phys. Fluids 19,053103 (2007).

[44] M. Kraus, W. Wintz, U. Seifert, and R. Lipowsky, Fluid Vesicles in Shear Flow,Phys. Rev. Lett. 77,3685(1996).

[45] S. K. Veerapaneni, D. Gueyffier, D. Zorin, and G. Biros, A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous fluid in 2D, J. Comput. Phys. 228,2334(2009).

[46] V. Kantsler and V. Steinberg, Transition to Tumbling and Two Regimes of Tumbling Motion of a Vesicle in Shear Flow,Phys. Rev. Lett. 96,036001(2006).

[47] C. Misbah, Vacillating Breathing and Tumbling of Vesicles under Shear Flow,Phys. Rev. Lett. 96,028104(2006). [48] V. Lebedev, K. Turitsyn, and S. Vergeles, Dynamics of Nearly

Spherical Vesicles in an External Flow, Phys. Rev. Lett. 99, 218101(2007).

[49] H. Noguchi and G. Gompper, Swinging and Tumbling of Fluid Vesicles in Shear Flow,Phys. Rev. Lett. 98,128103(2007). [50] R. Finken, A. Lamura, U. Seifert, and G. Gompper,

Two-dimensional fluctuating vesicles in linear shear flow,Eur. Phys. J. E 25,309(2008).

[51] Y. Navot, Critical behavior of drop breakup in axisymmetric viscous flow,Phys. Fluids 11,990(1999).

[52] V. Cristini, S. Guido, A. Alfani, J. Blawzdziewicz, and M. Loewenberg, Drop breakup and fragment size distribution in shear flow,J. Rheol. 47,1283(2003).