Motion and deformation of immiscible droplet in plane Poiseuille flow at low Reynolds number*

Ding-yi PAN （潘定一）， Yu-qing LIN （林雨青）， Ling-xin ZHANG （张凌新）， Xue-ming SHAO （邵雪明）

State Key Laboratory of Fluid Power and Mechatronic Systems， Key Laboratory of Soft Machines and Smart

Devices of Zhejiang Province， Department of Engineering Mechanics， Zhejiang University， Hangzhou 310027，China， E-mail:dpan@zju.edu.cn



Motion and deformation of immiscible droplet in plane Poiseuille flow at low Reynolds number*

Ding-yi PAN （潘定一）， Yu-qing LIN （林雨青）， Ling-xin ZHANG （张凌新）， Xue-ming SHAO （邵雪明）

State Key Laboratory of Fluid Power and Mechatronic Systems， Key Laboratory of Soft Machines and Smart

Devices of Zhejiang Province， Department of Engineering Mechanics， Zhejiang University， Hangzhou 310027，China， E-mail:dpan@zju.edu.cn

Droplet migration in plane Poiseuille flow is numerically investigated with a dissipative particle dynamics method. The single droplet deformation in the channel flow is first studied to verify the current method and the physical model. The effect of the viscosity ratio between the droplet and the solvent and the effect of the confinement are systematically investigated. The droplet is in an off-centerline equilibrium position with a specific selection of the parameters. A large viscosity ratio makes the droplet locate in a near-wall equilibrium position， and a large capillary number makes the droplet migrate to the near-centerline region of the channel. For the droplet migration at the same Capillary number， there is a critical width of the channel， which is less than twice of the droplet diameter， and the droplet will only migrate to the channel centerline if the width is less than this critical value.

droplet， migration， particle separation， microfluidics， dissipative particle dynamics

Introduction

Microfluidics is one of the promising research areas with its potential applications in biomedicine，chemical engineering， and other fields. In the biological research， the ill and infected cells have to be isolated to detect and predict diseases at an early stage，such as cancer and blood diseases［1，2］. The properties of ill cells may change during the isolation process，including the size， the deformability or the electric properties. A properly designed microfluidics device can help to isolate different cells by utilizing hydrodynamic interactions［3］. Meanwhile， in food， pharmaceutical and material industries， it is important to have approaches to filter different size droplets and manufacture monodispersed emulsion with predictable properties［4］.

From a hydrodynamic viewpoint， one of the basic mechanisms for the particle/cell separation is that the particle has an equilibrium lateral position in the channel flow due to the inertial effect. The migration of a solid particle was first studied in the 1990s and an equilibrium position between the centerline and the wall was reported. A second equilibrium position was found recently with the tube Reynolds number up to 1 000［5，6］. On the other hand， the mechanism of the cell migration is more complicated than that of the solid particle. Additional effects should be taken into account， including the deformability，the viscosity ratio， the interfacial tension as well as the viscoelasticity. The droplet suspended in another immiscible fluid can be considered as a simplified model of capsules or cells. The boundary integral method was first used to investigate the droplet motion and deformation in a parabolic flow， neglecting the inertial effect， i.e. as a Stokes flow［7，8］. The lateral migration towards the channel centerline was reported， and the size and the viscous effects were studied as well. In order to account for the inertial effect， Mortazavi and Tryggvason［9］coupled their two-dimensional Navier-Stokes solver to a front tracking algorithm to capturethe liquid-liquid interface， and the off-centerline equilibrium position was reported in their simulation. Meanwhile， the dependences of the droplet migration on the deformation， the viscosity ratio， the periodic length， and the Reynolds number were examined. Later on， this algorithm was extended to a threedimensional simulation［10］and the capsule migration with an elastic membrane model， i.e.， as neo-Hookean material［11］.

In the last decade， newly developed numerical algorithms coupled with Navier-Stokes equations were successfully applied to droplet migration studies， including the front tracking method［9］， the volume of fluid method［12］， and the lattice Boltzmann method［13，14］. The feasibilities of these algorithms for simulating the motion and deformation of a single droplet were demonstrated properly. However， certain phenomena in droplet hydrodynamics were seldom investigated with these algorithms， such as the droplet breakup and the coalescence of droplets. The topology of the droplet changes， which increases the complexity of these algorithms. The particle-based method might provide an alternative approach for handling these complex phenomena with their intrinsic features，without additional models being required to handle the topology variation. Recently， the dissipative particle dynamics （DPD） method was employed to study the droplet hydrodynamics， including the deformation and breakup［15］and the droplet-droplet interaction［16，17］in a linear shear flow. In this paper， we apply the DPD method to study the droplet motion and deformation in the plane Poiseuille flow， with consideration of the inertial effect， the viscosity ratio effect and the channel confinement effect. Several benchmark studies are presented， which may potentially provide insights for its applications in a wide range of droplet hydrodynamics.

1. Numerical algorithm and model

1.1 DPD method

The DPD method was first introduced by Hoogerbrugge and Koelman［18］， and is considered as a coarse-grained approach for the molecular dynamics simulation. The motion of each DPD particlei is governed by Newton's second law of motion

where riand viare the position and velocity vectors of particlei . Meanwhile， it is assumed that the mass of each DPD particle is unity.fiis the summation of inter-particle forces exerted by all the other particles. Espanol and Warren［19］suggested that the interparticle force fican be decomposed into three pairwise and center-to-center forces， i.e.， the conservative force， the dissipative force and the random force

Here， these forces will be active within a certain cutoff radius which is set to be unity as rC=1. These forces can be calculated as follows:

where kBTis the Boltzmann temperature of the system. In our simulation， we adopt the weight function proposed by Fan et al.［20］to have the dynamics response of the fluid with a high Schmidt number，

Furthermore， a large repulsive force is imposed between DPD particles from the solvent and the droplet， and hence the interfacial tension can be generated on the immiscible liquid-liquid interface. For more information of the interfacial tension， please refer toRef.［21］. In addition， our DPD simulation was also validated in our previous work［16］， concerning a single droplet in the shear flow and the droplet-droplet collision.

Fig.1 Sketch of moving droplet in plane Poiseuille flow

1.2 Physical model

The plane Poiseuille flow is one of the typical flows for studying the particle migration. As shown in Fig.1， in our simulation， a channel is created with two walls on top and bottom and periodic boundaries along xand y directions. The fluid flow from left to right（x -direction） driven by a pressure drop ∆pand the droplet with diameterD， is moving freely in the fluid. Frozen DPD particles are used to represent the solid wall and to satisfy the no-slip boundary condition. A parabolic velocity profile along the z-direction will be formed if there is no perturbation on the pure fluid and both the velocity profile and the maximum velocity at the centerline can be respectively calculated by

Here，LandH are the length（x）and the width（z）of the channel， andηis the viscosity of the fluid. In the DPD simulation， a gravity force along the flow direction is imposed on the DPD particles and the corresponding pressure drop can be determined as∆p=ρgL. Thus， the maximum velocity is Vmax= ρgLH2/（8η）. In addition， the gravity along z-direction is not considered in our study.

To classify the flow and the droplet deformation，the Reynolds number and the capillary number are defined， respectively， as

whereΓis the interfacial tension of the liquidliquid interface between the solvent fluid and the droplet.

2. Validation tests

2.1 Newtonian fluid in plane Poiseuille flow

The pure Newtonian fluid in a plane Poiseuille flow is first investigated in this subsection， demonstrating the capability of the current DPD simulation to capture the correct flow patterns of the Newtonian fluid. In our simulation， all variables are in the DPD unit. The gravity force/acceleration （g=0.004）acts on all fluid particles to generate the parabolic velocity profile. The velocity profile at the steady state is shown in Fig.2， and compared with the theoretical formula （Eq.（6））. Our DPD results agree well with the theoretical results.

Fig.2 Stream-wise velocity profile in plane Poiseuille flow: A comparison between DPD results and theoretical estimations

Fig.3 Comparison between theoretical estimation （solid）， numerical simulation with VOF （dashed） and our DPD simulation （spherical particles） of droplet deformation，droplet is moving along the centerline of plane Poiseuille flow

2.2 Single droplet deformation in plane Poiseuille flow A Newtonian droplet deforming and moving along the centerline of a long channel between two parallel plates is first studied. A comparison between our DPD results and the analytical solutions as well as other numerical results is shown in Fig.3. The analytical formula given by Nadim and Stone［22］is as follows

wherer is the radial coordinate，ϕis the angular coordinate，λis the ratio of the droplet viscosity to the solvent viscosity and the shape correction function is

Lan and Khismatullin［12］numerically simulated this problem by solving the Navier-Stokes equation and tracking the liquid-liquid interface with the VOF algorithm. Two cases with different Capillary numbers，Ca =0.28and Ca=1.40， were presented in their paper. Therefore， two sets of parameter shown in Table 1 are used in our DPD simulation to make sure that the Reynolds number and the Capillary number are close to their simulation cases. The viscosity ratio λis kept to 1.0.

Table 1 Parameters for droplet shape verification

For the low Cacase shown in Fig.3（a）， our DPD result agrees well with both the analytical and numerical results that the droplet is almost spherical. WhenCais high， the DPD result is closer to the numerical result than the analytical estimation. A bullet-like droplet is formed， with a flatter back than the analytical one.

Fig.4 Droplet trajectories in plane Poiseuille flow with different initial lateral positions （Z0）， droplet and solvent have the same viscosity and density，Re≈10，Ca =1.37，g=0.006

3. Equilibrium position and viscosity ratio effect

In this section， the equilibrium position for the droplet migration in the Poiseuille flow is investigated by the DPD simulation. To account for the inertial effect， theRe number is set around 10 in our simulation. Meanwhile， to eliminate the droplet-droplet interaction along the stream direction， the length of the channel in x-direction should be as large as possible，and we set the length to 80 within the computation limit. The width of the channel is set to 40， to reduce the confinement effect. The cases of the droplet and the solvent with the same viscosity and the same density are first considered. The droplets are released at different lateral positions （z-direction） initially.

Figure 4 shows the trajectories of the droplets moving in the Poiseuille flow released from three initial lateral positions. It is observed that the droplet trajectories collapse to the same equilibrium position eventually. The droplet released from the near-wall position（Z/ H=-0.3）migrates to the equilibrium position faster than the droplet released from the centerline. In addition， there are certain fluctuations when the droplet comes to its equilibrium position， and this is caused by the random nature of the force in the DPD simulation.

The snapshots of the droplet in the process of migration are shown in Fig.5. The droplets are initially spherical， then deformed according to their local flow patterns. The droplet from a near-wall position first deforms extremely and then recovers back to an equilibrium shape. The droplet from the centerline deforms gradually to its equilibrium shape and moves to the equilibrium position. Due to the nature of the Poiseuille flow， the local shear rate is larger at the near-wall place than that near the centerline. Therefore， the near wall droplet first deforms largely， and then recovers as it moves to the equilibrium position.

Fig.5 Droplet deformation in plane Poiseuille flow with different initial lateral positions，Re≈10，Ca =1.37，g= 0.006

The effects of the viscosity ratio between the droplet and the solvent， which is defined as λ=ηD/η，are also studied. The dissipative coefficientγis varied to produce different viscosities. One issue that should be taken into account is the dissipative interaction across the interface， i.e.， the dissipative coefficient between particles from the droplet and the solvent should be selected in a reasonable way. Visser et al.［23］made series of tests on the periodic Poiseuille flow［24］， and suggested that the dissipative coefficientbetween the droplet particle and the solvent particle could be set as

where γFDis the dissipative coefficient between particles from the droplet and the solvent，γDDis the dissipative coefficient among droplets and γFFis the dissipative coefficient among solvents. Their suggestion is followed in our simulation.

Fig.6 Droplet equilibrium position in plane Poiseuille flow versus Capillary number， compared with previous simulation results with viscosity ratio λ=1.0and λ=5.0

The equilibrium position results for the same viscosity cases are first compared with the results of Lan and Khismatullin［12］obtained with the VOF algorithm. The Capillary number varies from 0.1 to 3.0，since the droplet will break up if the Capillary number is higher and the compressibility effect should also be taken into account［25］. The comparison results are shown in Fig.6， where our DPD results are in good agreement with previous simulations. As the Capillary number increases， the equilibrium position approaches to the centerline. Two additional cases with the viscosity ratio up to 5 are also shown in Fig.6， the trend of the equilibrium position versus the Capillary number is the same with that of the case of the same viscosity. However， the variation range becomes smaller. For the same Capillary number cases， the droplet with a large viscosity has a near-wall position compared to the case of a small viscosity droplet. We does not consider cases with a higher viscosity ratio， e.g.，λ= 10， since the current model （Eq.（10）） can not give stable and correct results of the droplet equilibrium position， and a further improvement is required in the future.

4. Confinement effect

The effect of the ratio of the channel width to the droplet diameter， i.e.， the confinement effect， is one of the important aspects of the droplet hydrodynamic in the channel flow， since the interaction between the droplet and the solid wall may play a crucial role. The confinement factor can be defined as ξ=H/ D. Sibillo et al.［26］first studied the confinement effect on the droplet deformation and breakup in a linear shear flow. Later on， Guido and Preziosi［27］reviewed studies of the droplet deformation and breakup in a confined Poiseuille flow， however， the effect on the droplet equilibrium position has not been investigated comprehensively. In this section， we present our DPD simulation results of the droplet equilibrium position in a confined plane Poiseuille flow.

Fig.7 Dimensional equilibrium position of droplet at different Capillary numbers， the ratio of channel width to droplet diameter ranges from 1.5 to 5.0

Fig.8 Non-dimensional equilibrium position of droplet with different ratios of channel width to droplet diameter， the Capillary number ranges from 0.5 to 2.3

Figure 7 shows the equilibrium position of the droplet as a function of the Capillary number and each curve corresponds to one specific channel width. Here，ξvaries from 1.5 to 5.0， and obviously the trend of each curve is similar to the curve presented in Fig.6，i.e.， the equilibrium position approaches to the centerline as the Capillary number increases. A comparison between different curves shows that， with the same Capillary number， the smaller the width of the channel，the closer the equilibrium position to the centerline. For the narrowest case，ξ=1.75， the equilibriumpositions at different Capillary numbers are all located around the centerline very closely， which means the restrictions from both bound walls （top and bottom）are quite strong that lock the droplet at the centerline. Accordingly， we also plot the non-dimensional equilibrium position as functions of the confinement factor ξin Fig.8. It is shown that each curve in Fig.8 has a critical point， and for cases with the confinement factor less than this point， the droplet equilibrium position is always at the centerline. However， this critical point also depends on the Capillary number. It is shown that，the large the Capillary number， the larger the critical confinement factor.

Fig.9 Shapes of deformed droplet at its equilibrium position with different confinement factors

Three typical shapes of the deformed droplet are shown in Fig.9. In a wide channel， the droplet is elongated to take an ellipsoid-like shape. In a medium width channel， the droplet is less elongated and is balanced with a blunt head near the centerline and a sharp back on the other side. Whereas in a narrow channel， since the droplet migrates to the centerline，its shape turns to be symmetric with a small concave back. The deformation of the droplet is closely related to its local flow patterns， e.g.， the local shear rate， and meanwhile it is correlated to its internal flow structure，and thus a close inspection of this structure in the future is worthwhile.

5. Conclusions

The particle-based DPD method is used to study the droplet deformation and motion in the plane Poiseuille flow. Compared to the traditional simulation tools based on the Navier-Stokes equations， the current DPD method has its intrinsic advantages to handle the complex topology variation. In this paper，the single droplet deformation in the channel flow is first investigated for verification and validation. The results agree well with previous theoretical and numerical results， that confirms the reliability of the current DPD method.

The droplet migration phenomenon is of interest in biomedical and chemical engineering studies. The droplet might have an equilibrium lateral position due to the inertial effect in the solvent fluid， and this could be helpful for particle separation， infected cell detection， etc.. The viscosity ratio effect on the equilibrium position is discussed. The simulation results show that a large viscosity ratio leads to a near-wall equilibrium position， and meanwhile， a large Capillary number makes the droplet migrate to the channel centerline. The effect of the ratio of the channel width to the droplet diameter is investigated as well. For the droplet migration at the same Capillary number， there is a critical width of the channel， and the droplet will only locate at the centerline in the channel if the channel width is less than this critical value.

References

［1］ SOLLIER E.， GO D. E. and CHE J. et al. Size-selective collection of circulating tumor cells using Vortex technology［J］. Lab on a Chip， 2014， 14（1）: 63-77.

［2］ HUR S. C.， HENDERSON-MACLENNAN N. K. and MCCABE E. R. B. et al. Deformability-based cell classification and enrichment using inertial microfluidics［J］. Lab on a Chip， 2011， 11（5）: 912-920.

［3］ XUE Chun-dong， CHEN Xiao-dong and LIU Chao et al. Lateral migration of dual droplet trains in a double spiral microchannel［J］. Science China-Physics， Mechanics and Astronomy， 2016， 59（7）: 674711.

［4］ TAN Y. C.， LEE A. P. Microfluidic separation of satellite droplets as the basis of a monodispersed micron and submicron emulsification system［J］. Lab on a Chip， 2005，5（10）: 1178-1183.

［5］ MATAS J. P.， MORRIS J. F. and GUAZZELLI É. Inertial migration of rigid spherical particles in Poiseuille flow［J］. Journal of Fluid Mechanics， 2004， 515: 171-195.

［6］ SHAO X.， YU Z. and SUN B. Inertial migration of spherical particles in circular Poiseuille flow at moderately high Reynolds numbers［J］. Physics of Fluids， 2008，20（10）: 103307.

［7］ COULLIETTE C.， POZRIKIDIS C. Motion of an array of drops through a cylindrical tube［J］. Journal of Fluid Mechanics， 1998， 358: 1-28.

［8］ GRIGGS A. J.， ZINCHENKO A. Z. and DAVIS R. H. Low-Reynolds-number motion of a deformable drop between two parallel plane walls［J］. International Journal of Multiphase Flow， 2007， 33（2）: 182-206.

［9］ MORTAZAVI S.， TRYGGVASON G. A numerical study of the motion of drops in Poiseuille flow. Part 1. Lateralmigration of one drop［J］. Journal of Fluid Mechanics，2000， 411: 325-350.

［10］ NOURBAKHSH A.， MORTAZAVI S. A three-dimensional study of the motion of a drop in plane Poiseuille flow at finite Reynolds numbers［J］. Iranian Journal of Science and Technology， 2010， 34（2）: 179-196.

［11］ DODDI S. K.， BAGCHI P. Lateral migration of a capsule in a plane Poiseuille flow in a channel［J］. International Journal of Multiphase Flow， 2008， 34（10）: 966-986.

［12］ LAN H.， KHISMATULLIN D. B. A numerical study of the lateral migration and deformation of drops and leukocytes in a rectangular microchannel［J］. International Journal of Multiphase Flow， 2012， 47: 73-84.

［13］ KATAOKA Y.， INAMURO T. Numerical simulations of the behaviour of a drop in a square pipe flow using the two-phase lattice Boltzmann method［J］. Philosophical Transactions of the Royal Society of London A: Mathematical， Physical and Engineering Sciences，2011， 369（1945）: 2528-2536.

［14］ SHAN Ming-lei， ZHU Chang-ping and ZHOU Xi et al. Investigation of cavitation bubble collapse near rigid boundary by lattice Boltzmann method［J］. Journal of Hydrodynamics， 2016， 28（3）: 442-450.

［15］ CHEN S.， PHAN-THIEN N. and FAN X. J. et al. Dissipative particle dynamics simulation of polymer drops in a periodic shear flow［J］. Journal of Non-Newtonian Fluid Mechanics， 2004， 118（1）: 65-81.

［16］ PAN D.， PHAN-THIEN N. and KHOO B. C. Dissipative particle dynamics simulation of droplet suspension in shear flow at low Capillary number［J］. Journal of Non-Newtonian Fluid Mechanics， 2014， 212: 63-72.

［17］ ZHANG J. Movement of dispersed droplets of W/O emulsion in a uniform DC electrostatic field: Simulation on droplet coalescence［J］. Chinese Journal of Chemical Engineering， 2015， 23（9）: 1453-1459.

［18］ HOOGERBRUGGE P. J.， KOELMAN J. Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics［J］. EPL （Europhysics Letters）， 1992，19（3）: 155-160.

［19］ ESPANOL P.， WARREN P. Statistical mechanics of dissipative particle dynamics［J］. EPL （Europhysics Letters），1995， 30（4）: 191-196.

［20］ FAN X.， PHAN-THIEN N. and CHEN S. et al. Simulating flow of DNA suspension using dissipative particle dynamics［J］. Physics of Fluids， 2006， 18（6）: 063102.

［21］ PAN D.， PHAN-THIEN N. and KHOO B. C. Studies on liquid-liquid interfacial tension with standard dissipative particle dynamics method［J］. Molecular Simulation，2015， 41（14）: 1166-1176.

［22］ NADIM A.， STONE H. A. The motion of small particles and droplets in quadratic flows［J］. Studies in Applied Mathematics， 1991， 85（1）: 53-73.

［23］ VISSER D. C.， HOEFSLOOT H. C. J. and IEDEMA P. D. Modelling multi-viscosity systems with dissipative particle dynamics［J］. Journal of Computational Physics，2006， 214（2）: 491-504.

［24］ BACKER J. A.， LOWE C. P. and HOEFSLOOT H. C. J. et al. Poiseuille flow to measure the viscosity of particle model fluids［J］. The Journal of Chemical Physics， 2005，122（15）: 154503.

［25］ PAN D.， PHAN-THIEN N. and MAI-DUY N. et al. Numerical investigations on the compressibility of a DPD fluid［J］. Journal of Computational Physics， 2013，242（3）: 196-210.

［26］ SIBILLO V.， PASQUARIELLO G. and SIMEONE M. et al. Drop deformation in microconfined shear flow［J］. Physical Review Letters， 2006， 97（5）: 054502.

［27］ GUIDO S.， PREZIOSI V. Droplet deformation under confined Poiseuille flow［J］. Advances in Colloid and Interface Science， 2010， 161（1）: 89-101.

10.1016/S1001-6058（16）60673-X

April 10， 2016， Revised June 7， 2016）

* Project supported by the National Natural Science Foundations of China （Grant Nos. 11402230， 11332009 and 11272284）.

Biography: Ding-yi PAN （1984-）， Male， Ph. D.

Xue-ming SHAO，

E-mail:mecsxm@zju.edu.cn

2016，28（4）:702-708