Effect of slip on circulation inside a droplet
aa r X i v : . [ phy s i c s . f l u - dyn ] N ov Effect of slip on circulation inside a droplet
Joseph J. Thalakkottor and K. Mohseni ∗ Department of Mechanical and Aerospace Engineering, Universityof Florida, Gainesville, Florida 32611, USA Department of Electrical and Computer Engineering, Universityof Florida, Gainesville, Florida 32611, USAJune 28, 2018
Abstract
Internal recirculation in a moving droplet plays an important role inseveral droplet-based microfluidic devices as it enhances mixing, chemicalreaction and heat transfer. The occurrence of fluid slip at the wall, whichbecomes prominent at high shear rates and lower length scales, resultsin a significant change in droplet circulation. Using molecular dynamics(MD) simulations, the presence of circulation in droplets is demonstratedand quantified. Circulation is shown to vary inversely with slip length,which is a measure of interface wettability. A simple circulation model isestablished that captures the effect of slip on droplet circulation. Scalingparameters for circulation and slip length are identified from the circula-tion model which leads to the collapse of data for droplets with varyingaspect ratio ( A ) and slip length. The model is validated using continuumand MD simulations and is shown to be accurate for droplets with high A . Circulation in a droplet is known to play an important role in several droplet-based applications due to the enhanced mixing, chemical reaction and heattransfer that results from it. It has an increased impact at small scales wheresurface effects start to dominate, as seen in several microfluidic applications suchas lab-on-a-chip, microreactors and digitised heat transfer. In a microchan-nel, mixing or chemical reaction occurs mostly by diffusion. Since the flowremains laminar due to viscous forces dominating over body or inertia forces,fluid circulation or vortices could enhance mixing and/or chemical reactions.This helps in improving the efficiency of miniature-sized bio/chemical analy-sis systems that use microfluidic devices, such as lab-on-a-chip or micro to-tal analysis systems ( µ TAS) (Reyes et al. , 2002; Auroux et al. , 2002). Circu-lation also increases mass transfer, which enhances diffusive penetration and ∗ Email address for correspondence: mohseni@ufl.edu et al. , 2003). In the case of digitisedheat transfer, discrete microdroplets are used to “digitally” transfer heat awayfrom the source (Mohseni et al. , 2005; Baird & Mohseni, 2008). Presence ofa fluid-fluid interface leads to the formation of a circulatory flow inside thedroplet. This results in convection normal to the wall, improving the efficiencyof the thermal management system. Some research groups have already demon-strated rapid mixing or reaction by means of shuttling a droplet in a microchan-nel (Hosokawa et al. , 1999) or transporting it through a winding microchannel(Song et al. , 2003). Thereby, the importance of circulatory flow in microfluidicdevices is observed.Although, slip at the boundary is prevalent for single phase flows, it is neg-ligibly small in most continuum and macro-scale applications. However, inmany micro- and/or nano-scale applications the first breakdown of continuumassumption often occurs at a solid boundary in the form of velocity slip. Slipat micro scale has been investigated extensively for single phase fluids, butthe research has been limited for effects of slip in two phase flows. Recently,Thalakkottor & Mohseni (2013) extended Maxwell’s slip model from steady tounsteady flows. The established unsteady model for single phase flows showeda dependence on gradient of shear rate in addition to shear rate of the flow.Also, a time-dependent universal curve for slip at the wall that encapsulatedunsteady and steady flows was found.The presence of circulatory flow inside a droplet is known for continuum scaleproblems (Paik et al. , 2003; Mohseni & Baird, 2007), but unlike some contin-uum scale hydrodynamic phenomena, circulation is observed even at nanoscales.Koplik et al. (1989) and Thompson & Robbins (1989) studied the motion of acontact line using molecular dynamics (MD) simulations and observed the pres-ence of circulatory flow. Here, these observations are replicated and also quan-tified. Changes in circulation is studied for varying degree of wetting of thewall. It is seen that the occurrence of fluid slip at the wall results in a signif-icant change in circulation. Circulation varies inversely with slip length whichcould be considered as a measure of interface wettability. In this paper a simplecirculation model is established that captures the effect of slip on circulation.Scaling circulation and slip length, the results for different droplet aspect ra-tio ( A ) collapse onto a single curve. The model is verified using results fromcontinuum and MD simulations and is shown to be accurate for droplets withhigh A . It is shown that there can be as much as reduction in circulationfor a droplet moving in a non-wetting channel. Finally, the effect of slip at thefluid-fluid interface on circulation in a droplet is briefly discussed.Details of the problem setup are specified in section 2. In section 3 thecirculation model is derived and the results are discussed in section 4. In order to study the effects of wall slip on circulation, a droplet moving incontact with a wall is simulated using Couette flow as a canonical problem.Two immiscible fluids occupying equal volume are placed in between two par-allel walls separated by a distance H . Top and bottom walls move in opposite2igure 1: Schematic of the problem. The problem simulates a Couette flow withtwo immiscible fluids.directions with a velocity U . A schematic of the problem geometry is shown infigure 1. Periodic boundary conditions are imposed along the x and z directions.The size of fluid domain used in MD simulations, for different droplet aspectratios are listed in table 1. The problem is studied primarily using MD simu-lations. Due to the high cost of MD simulations, results from a Navier Stokescontinuum solver are used to supplement the data, by performing runs for awider range of droplet aspect ratios and wall-fluid slip lengths. Two immisciblefluids with identical properties are chosen in order to simplify the problem andfocus only on how circulation is affected by slip at the wall.The molecular dynamics simulations presented in this paper are performedusing the LAMMPS package (Plimpton, 1995). The fluid’s initial state is mod-elled as a face centered cubic (fcc) structure with the x direction of the channelbeing aligned along the [11¯2] orientation of the face-centered cubic lattice. Thewall is comprised of – layers of atoms oriented along the (111) plane of fcclattice. The wall atoms are fixed to their lattice sites.The pairwise interaction of molecules separated by a distance r is modelledby the Lennard Jones (LJ) potential V LJ = 4 ǫ (cid:20)(cid:16) σr (cid:17) − (cid:16) σr (cid:17) (cid:21) , (1)where ǫ and σ are the characteristic energy and length scales. The potential iszero for r > r c , and the cutoff radius, r c is . σ .The simulations are performed for immiscible fluids having properties sameas that of liquid Argon. An interface between the two fluids is simulated bycontrolling the attractive force between the two. This is done by varying the LJparameters. Reducing the value of ǫ reduces the attractive potential betweenthe two fluids, while increasing σ allows to control the thickness of the interface.For the studies presented here the LJ parameters for fluid-fluid interaction arechosen as . ǫ and . σ .The fluid is maintained in its equilibrium state having a number density ρ = 0 . σ − and temperature T = 1 . k B /ǫ . The temperature is regulated by athermostat which simulates the transfer of heat from the system to an externalreservoir. A Langevin thermostat with a damping coefficient of Γ = 1 . τ − ,where τ = p mσ /ǫ , is used here. The damping term is only applied to the z direction to avoid biasing the flow. The equation of motion along the z H ( σ ) L ( σ ) W ( σ ) . .
152 144 . . .
152 144 . . .
222 144 . Table 1: Dimensions of the fluid domain in MD simulations, for droplets withdifferent aspect ratio. H , L and W are height, length and width of the domain,respectively.components is, m ¨ z i = X j = i ∂V ij∂z i − m Γ ˙ z i + η i . (2)Here, m is mass of the fluid molecule, P j = i denotes a sum over all interactionswith i and η i is a Gaussian distributed random force. The equations of mo-tion were integrated using the Verlet algorithm (Verlet, 1967; Allen & Tildesley,1987) with a time step ∆ t = 0 . τ .The simulation is initially run for a time of τ without any shear and thenfor an additional ∼ τ with moving walls, allowing the flow to equilibrate.After ∼ τ , the spatial averaging is performed by dividing the fluid domaininto square bins of size ∼ . σ along the x – y plane, and extending through theentire width of the channel. In addition to spatial averaging, ensemble averagingis done for a duration of τ in intervals of τ . In the case of non-wettingwall, averaging was done for an extended time of τ . A mathematical model that describes the effect of slip length on circulation, isestablished for a droplet moving in a channel where the two fluids are immiscible.The flow is driven by walls moving in opposite directions with a velocity, U , inthe frame of reference of the droplet.Velocity of fluid adjacent to the wall, u wallf , can be written from the Navierslip boundary condition u wallf = U − L s du f dy , (3)where L s is the slip length. In the case of a Couette flow, the x -component ofvelocity has a linear profile across the height of the droplet. Hence, the aboveequation can be re-written as u wallf = U (cid:18) HH + 2 L s (cid:19) . (4)Circulation inside a droplet is expressed as Γ = Z BA u f · d l + Z CB u f · d l + Z DC u f · d l + Z AD u f · d l , (5)where A, B, C and D are the four corners of the droplet, see figure 1. Assumingthe circulation contribution by the top and bottom walls and the left and right4nterfaces are the same, the equation reduces to Γ = 2 Z L u wallf · d l + 2 Z L intf u intff · d l , (6)where L and L intf are the length of droplet along the wall and interface, respec-tively and u intff is the fluid velocity along the fluid-fluid interface. Assumingthe velocity of fluid along the interface is same as that along the wall and sub-stituting u wallf onto the above equation, one obtains Γ = 2 Z L U (cid:18) HH + 2 L s ( l ) (cid:19) d l + 2 Z L intf U (cid:18) HH + 2 L s ( l ) (cid:19) d l. (7)Performing a change of variable from dl along the wall to x -direction andfrom dl along the interface to y -direction we obtain, Γ = 2 Z L U (cid:18) HH + 2 L s ( x ) (cid:19) d x + 2 Z H U (cid:18) HH + 2 L s ( y ) (cid:19) d y sin θ , (8)where θ is the respective averaged contact angle. Here, averaged contact angleis defined as the acute angle between the wall and a straight line which describesthe fluid-fluid interface.Scaling x and y with L and H respectively and doing a change of variable,circulation along a droplet can be represented as Γ2 U L = Z
11 + 2 L s ( ζ ∗ ) H d ζ ∗ (cid:18) A sin θ (cid:19) , (9)where A = L/H is the droplet aspect ratio.The variation of L s can be assumed to be symmetric about the dropletcenterline and can be approximated by a power law, Γ2 U L = 2 Z .
11 + 2 L os H ( αζ ∗ β + γ ) d ζ ∗ (cid:18) A sin θ (cid:19) . (10)Here, α, β and γ are curve fitting coefficients to L s /L os and x data, where L os is the asymptotic value of slip length as defined by Thompson & Troian (1997)for a single phase Couette flow.Scaling the droplet circulation with slip at the walls, Γ s , by the circulationwith no-slip, Γ ns , results in Γ s Γ ns = 2 Z . α ′ ζ ∗ β + γ ′ d ζ ∗ (cid:18) / A sin θ s / A sin θ ns (cid:19) , (11)where α ′ = 2 L os H α , γ ′ = 2 L os H γ + 1 and θ s , θ ns are the average contact angles forwall with slip and no-slip, respectively.For a long droplet, the decrease in fluid velocity in the vicinity of the triplecontact point has a negligible affect on the total droplet circulation. Hence, sliplength can be assumed to be a constant along the length of the droplet. In thiscase the above equation reduces to, Γ s Γ ns = 11 + 2 L os H (cid:18) / A sin θ s / A sin θ ns (cid:19) . (12)5ase ǫ wf /ǫ σ wf /σ ρ w /ρ . . . . . .
75 4 . .
75 4 . .
75 4
Table 2: Five different cases of wall-fluid properties are considered, with wallwettability decreasing from case 1 to 5.From the above expression it can be seen that circulation varies inverselywith slip length, though the extent of its effect on circulation is dependenton whether the length scale of the problem is comparable to the slip lengthscale. Also, it can be observed that for droplets with high A , the contributionto circulation from the fluid-fluid interface reduces relative to the contributionfrom the wall. Several runs are performed for five different cases of wall-fluid properties, asdescribed in table 2, with decreasing wall wettability from case 1 to case 5. Thewalls move with a speed of U = 0 . σ/τ in opposite directions and the velocitiesfrom MD presented in this paper are normalised by it. In order to computevarious values at the wall, a reference plane is defined at a distance of . σ wf away from the wall lattice site.The continuum simulations simulate a rectangular droplet with ° contactangle and a free surface boundary condition at the interface. The top andbottom walls move in opposite directions. The length of the droplet is fixed at L = 1 and the height is calculated based on the aspect ratio. Details of theNavier Stokes continuum solver is given in Zhang & Mohseni (2012). A set ofruns are made with varying slip lengths and A .The fluid-fluid interface in the droplet forces convection normal to the wall,causing the fluid to circulate inside a moving droplet. A plot of the velocityfield of a moving droplet in figure 2 shows the presence of circulation insidethe droplet. As the two fluids share the same wall they have a co-rotatingflow. The flow field is similar to what one would see in continuum scale droplets(Paik et al. , 2003; Baird & Mohseni, 2008) and was previously observed by Koplik et al. (1989) and Thompson & Robbins (1989) in MD simulations.For the Couette flow problem considered here, wall movement drives thecirculation inside a droplet, hence it is intuitive that the occurrence of slip atthe wall will affect circulation. Partial slip at the wall reduces circulation and ifthere is perfect slip the circulation reduces to zero. This reduction in circulationdirectly impacts the rate of mixing, chemical reaction, mass and heat transferessential for various droplet based microfluidic devices. Figure 3 shows how totaldroplet circulation varies inversely with slip length, which could be a measureof wall wettability. Its variation with slip length is considered for droplets withdifferent aspect ratios. 6
50 1000510152025 x( σ ) y ( σ ) Figure 2: MD result for the velocity field of a droplet moving in a channel.Each droplet exhibits a co-rotating circulatory flow. The y-axis is stretched toprovide clarity. Γ ( σ / τ ) L s ( σ ) − . − . − − − AAAAA
Figure 3: Circulation versus slip length for droplets with different A . Resultsfrom continuum simulations shows the decrease in circulation with increasingslip length which is a measure of wettability.7a) L s /H Γ s / Γ n s − − − − − . − . AAAAAA (b) L s /H Γ s / Γ n s ∼ ∼ ∼ AAA
Figure 4:
Nondimensional circulation versus slip length for, (a) continuum simulations and(b) MD simulations. Results from continuum and MD simulations for various droplet aspectratios are shown. For droplet with high A the results collapse onto a curve, confirming theproper scaling of circulation and slip length. For low A droplets the data does not collapseonto the same curve as that for high A droplets. U ∗ y ∗ − − . AA Figure 5: Velocity profile across the center of the droplet for A of 0.5 and 2.0.Here, y ∗ = y/H and U ∗ = | u/U − | , which gives the velocity in the lab frameof reference. Deviation from the linear velocity profile of a Couette flow is seenwhich also implies the formation of two circulatory flows in a droplet.The proper scaling for circulation and slip length is identified from the cir-culation model. Circulation with slip, Γ s , is scaled by circulation for a perfectlywetting (no-slip) wall, Γ ns , and slip length, L S , is scaled by the droplet height, H . For droplets with high A , the data for circulation collapse onto a singlecurve, see figure 4. From the MD results presented in figure 4(b), a decreasein circulation of as much as is seen for a non-wetting case, demonstratingthe extent of effect that slip can have on circulation. The impact of length scaleof the problem on circulation can also be seen here. For droplets where theheight of the droplet is comparable to slip length, there is a larger reduction incirculation in comparison to a taller droplet. For low A droplets, A ≤ . for continuum results and A ≤ for MD results, the data deviates from thecollapsed curve. This deviation can be attributed to the change in flow dy-namics of tall droplets. As figure 5 shows, for tall or narrow droplet a pair ofcirculatory flows form inside the droplet. This breaks down the assumption ofa linear velocity profile across the height of the droplet.The variation of slip length along the wall needs to be described in order tocalculate circulation based on the model established in the previous section. Itis found that by scaling slip length, L s , by its asymptotic value, L os , the data8 x/L L s / L o s Case 3Case 4Case 5 α x ∗ β + γ x ∗ + 1 Figure 6: Non-dimensional slip length along the wall. Plot shows the collapseof data for various cases of wall wettability. A power law of the form αx ∗ β + γ ,where x ∗ = x/L , is a good fit for this data.for different cases of wall wettability collapse onto a single curve, see figure 6.A power law of the form αx ∗ β + γ , where x ∗ = x/L , is a good fit for thisdata. The distribution of slip length can be assumed to be symmetric aboutthe center of the droplet. The coefficients of the power law are obtained bynonlinear least squares method as, α = 0 . , β = − . and γ = 1 . .They can be rounded off to α = 0 . , β = − . and γ = 1 . , to obtain L s /L os =(1 / L/x + 1 .The circulation model established in the previous section is verified by com-paring the scaled circulation ( Γ s / Γ ns ) obtained from MD and continuum sim-ulations with that calculated from the model. Figure 7 shows the model to beaccurate for droplets with high A , having an error of less than . The accu-racy of the circulation model decreases sharply for droplets with low A . Thisis mostly attributed to the breakdown of assumption that the velocity profileis linear across the height of the droplet, as discussed previously. In addition,the error also increases gradually with slip length until it reaches an asymptoticvalue. The rate of increase in error is seen to increase with decrease in aspectratio. This is a result of the assumption made that fluid velocity at the interfaceis the same as that at the wall. But, in fact the fluid moving along the fluid-fluidinterface decays due to viscous diffusion as it travels towards the center of thedroplet. Hence, decrease in fluid inertia due to decrease in wettability and/ortravelling longer distances along the fluid-fluid interface in the case of dropletswith lower A leads to reduction in accuracy of the model. From the circulationmodel given by equation 11 it is seen that the contribution of velocity along theinterface to total circulation increases with decreasing aspect ratio. Although,there is a reduction in accuracy of the model for low A droplets, it is knownthat for most practical applications, high A droplets are used for which themodel is accurate.In the current study only the case of immiscible fluids were considered butslip at a general fluid-fluid interface of two partially miscible fluid could alsohave an affect on circulation in a droplet. Immiscible fluids have a sharp, dis-tinct interface where the two fluids do not experience any transfer of tangentialmomentum at the interface. This is analogous to a fluid with perfect slip at9a) L s /H % e rr o r Γ ∗ − − − − − . − . AAAAAA (b) L s /H % e rr o r Γ ∗ ∼ ∼ ∼ AAA
Figure 7: % error of non-dimensional circulation in (a) continuum simulations and (b)MD simulations. The percentage relative error for circulation calculated using the cir-culation model is plotted for various droplet A . The error is computed as, ( | Γ ∗ model − Γ ∗ simulation | ) / Γ ∗ simulation , where Γ ∗ = Γ s / Γ ns . The error increases with decrease in droplet A and increase in slip length. For large aspect ratio droplets, the error is less than . the wall. In the case of partially miscible fluids the interaction between thetwo fluids leads to a change in momentum along the interface. As the flowin the droplets are co-rotating there could be an apparent reduction in totalcirculation in a droplet, depending on the relative molecular mass of the twofluids. But, as circulation is calculated over a material line, the occurrence ofmass transfer across the interface leads to difficulties in its calculation. UsingMD, it is currently computationally expensive to obtain sufficient statistics for acase with partial miscibility between the fluids as the interface starts stretching,eventually leading to breakup of the droplet. Circulation is an important factor that affects mixing, chemical reaction andheat transfer inside a droplet, which have many practical applications. This pa-per demonstrates the presence of circulation in droplets using MD simulationsand its inverse dependence on slip length at the wall-fluid interface. A mathe-matical model for circulation is established that shows this dependence on sliplength, Γ s Γ ns = 2 Z . L s H αζ ∗ β + γ ′ d ζ ∗ (cid:18) / A sin θ s / A sin θ ns (cid:19) . The model also shows that the extent of change in circulation caused by slipdepends on whether the length scale of the problem is of the order of slip lengthscale. On scaling circulation and slip length, the results for different droplet A collapse onto a single curve. The model is verified using results from con-tinuum and MD simulations and is shown to be accurate for droplets with high A . Also, it is seen that there can be as much as change in circulation fora droplet moving in a non-wetting channel.This research was supported by the Office of Naval Research.10 eferences Allen, M.P. & Tildesley, D.J.
Computer Simulation of Liquids . Eng-land: Clarendon Press Oxford.
Auroux, Pierre-Alain, Iossifidis, Dimitri, Reyes, Darwin R. & Manz,Andreas
Analytical Chemistry (12), 2637–2652. Baird, E. & Mohseni, K.
IEEE Transactions onComponents and packaging Technologies (1), 143–151. Burns, JR & Ramshaw, C
Lab on a Chip (1), 10–15. Dummann, Gerrit, Quittmann, Ulrich, Gröschel, Lothar, Agar,David W, Wörz, Otto & Morgenschweis, Konrad
Catalysis Today (0), 433 – 439.
Hosokawa, Kazuo, Fujii, Teruo & Endo, Isao
Analyt-ical Chemistry (20), 4781–4785. Koplik, J., Banavar, J. & Willemsen, J.F.
Physics of Fluids, A , 781–794. Mohseni, K. & Baird, E.
Nanoscale and Microscale Thermophysical Engineering (1 & 2), 109– 120. Mohseni, K., Baird, E. & Zhao, H.
Proceedings of the 2005 ASMEInternational Mechanical Engineering Congress and R & D Expo , pp. IMECE2005–79372. Orlando, Florida: ASME.
Paik, P., Pamula, V. K. & Fair, R. B.
Lab Chip (2), 253–259. Plimpton, S.
J. Comp. Physics , 1–19.
Reyes, Darwin R., Iossifidis, Dimitri, Auroux, Pierre-Alain & Manz,Andreas
Analytical Chemistry (12), 2623–2636. Song, Helen, Tice, Joshua D. & Ismagilov, Rustem F.
Angewandte ChemieInternational Edition (7), 768–772. Thalakkottor, J.J. & Mohseni, K.
Phys. Rev. E (033018).11 hompson, P. & Robbins, M. Phys. Rev. Lett. , 766–769. Thompson, P.A. & Troian, S.M.
Nature , 360–362.
Verlet, L.
Phys. Rev. (1).
Zhang, P. & Mohseni, K.10th international Conference on Nanochannel, Microchannel andMinichannels