UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Simulation of selected interfacial dynamic problems using Cahn-Hilliard diffuse-interface method Mehrabian, Hadi 2014

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2014_spring_mehrabian_hadi.pdf [ 2.63MB ]
JSON: 24-1.0166873.json
JSON-LD: 24-1.0166873-ld.json
RDF/XML (Pretty): 24-1.0166873-rdf.xml
RDF/JSON: 24-1.0166873-rdf.json
Turtle: 24-1.0166873-turtle.txt
N-Triples: 24-1.0166873-rdf-ntriples.txt
Original Record: 24-1.0166873-source.json
Full Text

Full Text

Simulation of Selected InterfacialDynamic Problems usingCahn-Hilliard Diffuse-InterfaceMethodbyHadi MehrabianB.Sc., Amirkabir University of Technology, 2003M.Sc., Sharif University of Technology, 2005A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Chemical and Biological Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)March 2014c? Hadi Mehrabian 2014AbstractUsing the Cahn-Hilliard diffuse-interface model, I have studied three inter-facial dynamic problems for incompressible immiscible two-phase flows. Asthe first problem, capillary instability of a liquid torus is computed. Themain differences between the torus and a straight thread are the presenceof an axial curvature and an external flow field caused by the retraction ofthe torus. We show that the capillary wave initially grows linearly as on astraight thread. The axial curvature decreases the growth rate of the cap-illary waves while the external flow enhances it. Breakup depends on thecompetition of two time scales: one for torus retraction and the other forneck pinch-off. The outcome is determined by the initial amplitude of thedisturbance, the thickness of the torus relative to its circumference, and theviscosity ratio.The second problem concerns interfacial dynamics and three-phase con-tact line motion of wicking through micropores of two types of geometries:axisymmetric tubes with contractions and expansions of the cross section,and two-dimensional planar channels with a Y-shaped bifurcation. Resultsshow that the liquid meniscus undergoes complex deformation during itspassage through contraction and expansion. Pinning of the interface at pro-truding corners limits the angle of expansion into which wicking is allowed.Capillary competition between branches downstream of a Y-shaped bifur-cation may result in arrest of wicking in the wider branch.As the third problem, auto-ejection of drops from capillary tubes is stud-ied. This study focuses on two related issues: the critical condition for auto-ejection, and the role of geometric parameters in the hydrodynamics. Fromanalyzing the dynamics of the meniscus in the straight tube and the nozzle,we develop a criterion for the onset of auto-ejection based on a Weber num-ber defined at the exit of the nozzle and an effective length that encompassesthe geometric features of the tube-nozzle combination. In particular, thiscriterion shows that ejection is not possible in straight tubes. With steepercontraction in the nozzle, we predict two additional regimes of interfacialrupture: rapid ejection of multiple droplets and air bubble entrapment.iiPrefaceThis PhD thesis entitled ?Simulation of Selected Interfacial Dynamic Prob-lems using the Cahn-Hilliard Diffuse-Interface? presents the main featuresof the research that I led and carried out during my PhD study under su-pervision of Professor J. J. Feng. In this preface, the contributions andcollaborations to the papers published or submitted for publication fromcurrent thesis are briefly explained.? A version of chapter 3 has been published. H. Mehrabian, and J.J. Feng (2013), Capillary breakup of liquid torus, Journal of FluidMechanics 717: 281-292. Under supervision of J. J. Feng, I studiedthe capillary instability of a liquid torus and drafted the paper. J. J.Feng helped me to prepare the final version of the paper. In addition,I acknowledge discussions with Giovanni Ghigliotti, Ekapop Pairamand Qiming Wang.? A version of Chapter 4 has been published. H. Mehrabian, P. Gao,and J. J. Feng (2011), Wicking flow through microchannels. Phys. ofFluids 23: 122108. Under supervision of J. J. Feng, I studied wickingflow inside two type of geometries: axisymmetric area change and Y-shape 2D branching and drafted the paper. J. J. Feng helped me toprepare the final version of the paper.? A version of Chapter 5 has been prepared for submission. I conductedthis study under supervision of J.J. Feng. The final version of themanuscript was prepared with help of J. J. Feng. I would thank An-drew Wollman for generously providing data and video of his experi-ments to me.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiiList of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . xviiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Methodology of Research . . . . . . . . . . . . . . . . . . . . . 72.1 Governing equations . . . . . . . . . . . . . . . . . . . . . . . 72.1.1 Simulation of moving contact lines . . . . . . . . . . . 92.2 Parameters of the Cahn-Hilliard model . . . . . . . . . . . . 102.3 Discretization of governing equations and computational do-main . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.4 Limitations on resolving small lenght scales . . . . . . . . . . 122.4.1 Bubble production in meniscus sessile droplet collision 132.4.2 Dependence of bubble dynamics on meniscus shape . 142.4.3 Dependence of bubble volume on cut-off length . . . 153 Capillary Breakup of a Liquid Torus . . . . . . . . . . . . . . 163.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 163.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . 183.3 Results: linear growth of capillary waves . . . . . . . . . . . 19ivTable of Contents3.3.1 Quasi-static retraction: effect of axial curvature . . . 193.3.2 Faster retraction: effect of external flow . . . . . . . . 213.4 Results: nonlinear growth and breakup . . . . . . . . . . . . 223.4.1 Fastest mode . . . . . . . . . . . . . . . . . . . . . . . 223.4.2 Pinch-off time versus retraction time . . . . . . . . . 233.5 Comparison with experiment . . . . . . . . . . . . . . . . . . 263.6 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 284 Wicking Flow through Microchannels . . . . . . . . . . . . . 304.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . 324.2.1 Choice of Cahn-Hilliard parameters . . . . . . . . . . 334.3 Wicking in a tube with contraction or expansion . . . . . . . 364.3.1 Contraction . . . . . . . . . . . . . . . . . . . . . . . 364.3.2 Regularization of corner singularity . . . . . . . . . . 384.3.3 Expansion . . . . . . . . . . . . . . . . . . . . . . . . 394.3.4 Penetration time . . . . . . . . . . . . . . . . . . . . . 414.4 Wicking in Y-shaped branches: capillary competition . . . . 454.4.1 Flow in both branches . . . . . . . . . . . . . . . . . 464.4.2 Flow in one branch . . . . . . . . . . . . . . . . . . . 474.4.3 Flow reversal due to spatially inhomogeneous hydrophilic-ity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 504.5 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 525 Auto-ejection of Liquid Drops from Capillary Tubes . . . 545.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 545.2 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . 565.3 Physical model and numerical algorithm . . . . . . . . . . . 585.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 605.4.1 Meniscus dynamics . . . . . . . . . . . . . . . . . . . 605.4.2 Ejection criterion . . . . . . . . . . . . . . . . . . . . 655.4.3 Rapid ejection and air entrapment . . . . . . . . . . . 695.4.4 Contact line depinning at nozzle lip . . . . . . . . . . 715.5 Summary and conclusions . . . . . . . . . . . . . . . . . . . . 746 Conclusions and Recommendations . . . . . . . . . . . . . . 766.1 Summary of key findings . . . . . . . . . . . . . . . . . . . . 766.1.1 Capillary breakup of a liquid torus . . . . . . . . . . 766.1.2 Wicking flow through microchannels . . . . . . . . . . 77vTable of Contents6.1.3 Auto-ejection of liquid jets and drops from capillarytubes . . . . . . . . . . . . . . . . . . . . . . . . . . . 776.2 Significance and limitations . . . . . . . . . . . . . . . . . . . 786.2.1 Recommendations for future work . . . . . . . . . . . 80Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81viList of Figures1.1 (a) Schematic of a liquid torus suspended in a bath of surroundingliquid (b) Experimental results of [57] on the breakup of silicone-oil tori suspended in glycerine. The thin torus breaks down intomultiple droplets while the fat torus shrinks into one droplet. . . . 21.2 Schematic of an axisymmetric contraction and expansion in tubearea. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.3 A sequence of snapshots showing spontaneous capillary rise andauto-ejection of droplets in the experiment of [88] under micro-gravity. The inner diameter of the glass tube is 9.2 mm in thestraight section, and the liquid is PDMS of viscosity 0.65 cs. Thedrop volume is roughly 20 ?l. Ohnesorge number Oh = 0.0015,static contact angle ? = 0o , contraction angle of nozzle ? = 17o,contraction ratio of nozzle C = 0.42, tube length L = 8.0 (see ? 5.2for detailed definition of the terms ?,C,Oh, and L). The photosare taken 0.1 s apart. Adapted from [88] with permission, c?Springer. 52.1 (a) A meniscus with a velocity U and contact angle of ? will impacta sessile droplet with a equivalent radius R (b) Three dimensionalschematic of the problem setup. . . . . . . . . . . . . . . . . . . 142.2 Effect of diffusion length scale on the coalescence type. . . . . . . 153.1 (a) A quarter of the top half of a liquid torus for simulating thecapillary growth of an even mode, i.e. with an even number ofwavelengths around the torus. For odd modes, a half of the tophalf must be included. (b) The interface on the symmetric mid-plane with and without a sinusoidal disturbance. . . . . . . . . . . 17viiList of Figures3.2 (a) Dispersion relation on a shrinking torus compared to that fora straight filament. The latter is computed by our diffuse-interfacemethod and agrees with the Tomotika formula within 4%. Thewavelength l and the growth rate ? are made dimensionless by theinstantaneous 2?a and tc, respectively. (b) The linear growth ratedecreases with the axial curvature for a prescribed dimensionlesswave length l0= 2. The point at 1/? = 0 corresponds to a straightfilament. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.3 (a) Ratio of growth rates on a torus as a function of the viscosityratio for a capillary wave of dimensionless wavelength l = 2. (b) Ra-tio of growth rates on a straight filament under uniform extensionalflow, calculated from the theoretical result of [52]. . . . . . . . . . 213.4 Nonlinear evolution of three modes of instability, with wave numberk = 2, 3 and 4, for ? = 6.7, m = 0.033 and initial amplitude?0= 0.02. ? is the instantaneous amplitude of the capillary waves.The curves for k = 2 and k = 3 end in breakup, with the onset ofsecondary necking also marked on the latter. The k = 4 mode endsin complete retraction. . . . . . . . . . . . . . . . . . . . . . . . 233.5 Snapshots of the evolving interface on the mid-plane of the torusfor ? = 6.7, m = 0.033 and ?0= 0.02. The interface is given by thelevel set of ? = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . 243.6 (a) The pinch-off time decreases with increasing initial amplitudeof disturbance. ? = 5.3, m = 0.033, k = 2. The solid curve isthe best fitting by Eq. (3.2). (b) The critical initial amplitude ?cdecreases with the initial aspect ratio ?. The solid curve is the bestfitting by Eq. (3.3). . . . . . . . . . . . . . . . . . . . . . . . . . 243.7 Effect of the viscosity ratio on the growth of disturbance. ? = 5.3,?0= 0.02 and k = 2. . . . . . . . . . . . . . . . . . . . . . . . . 253.8 (a) Determining the initial amplitude of perturbation ?0from thevariation of the thickest radius at versus the thinnest radius an onthe torus. Both radii are normalized by the initial value a0. (b)Determining the interfacial tension ? from the temporal variationof an. ? = 5.3, k = 2 and m = 0.033. . . . . . . . . . . . . . . . . 273.9 Comparison between the predicted and observed number of primarydrops after breakup, for tori with five initial aspect ratios. m =0.033 and ?0= 0.01. N = 0 and 1 refer to, respectively, completeretraction with no breakup and breakup at a single primary neck. . 284.1 Schematic of the flow geometry for wicking into a capillary tubewith contraction. . . . . . . . . . . . . . . . . . . . . . . . . . . 32viiiList of Figures4.2 Sharp-interface limit for computing capillary rise. (a) Contact linemotion indicated by the rise of Hw in time for two Cahn numbersCn = 0.01 and 0.02. (b) Variation of the effective meniscus curva-ture ? with time for the same two Cn values. S = 0.04, ? = 60?,the tube radius R = 1, total length Ht = 20 and the initial columnheight H0= 15. . . . . . . . . . . . . . . . . . . . . . . . . . . 344.3 Comparison between the diffuse-interface simulation and the ana-lytical Lucas-Washburn formula at different S values. (a) Imbibi-tion with ? = 60?, Cn=0.01, m = 0.02 and H0= 15. (b) Drainagewith ? = 120?, Cn=0.01, and H0= 19. Now the less viscous com-ponent is wetting, and the non-wetting-to-wetting viscosity ratiom = 50. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.4 Meniscus movement through a contraction with ? = 45? repre-sented by (a) the variation of the base point with the wall point,and (b) the effective curvature defined in Eq. (4.1). ? = 60?,Cn = 0.005 and Hu = 10. In (a) the insets correspond to thefour points marked by squares on the curve. In (b) the dashed lineindicates the curvature expected of a quasi-static spherical meniscus. 374.5 Gray-scale contours of ? depicting the interface traversing a concavecorner through Cahn-Hilliard diffusion. The light line indicates thecontour of ? = 0. The contraction angle ? = 75? is greater thanthe wetting angle ? = 60?. Hu = 10 and initially the meniscus isat H0= 9.8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.6 Schematic of an expansion illustrating the pinning criterion. ?b =? + ? is the breakthrough angle, and ?m = 90? is the maximumangle that the interface may reach at the corner. . . . . . . . . . . 404.7 Wicking through an expansion with ? = 25?, ? = 60?, Cn = 0.005and Hu = 5. The insets correspond to the four points a?d on thecurve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.8 Permanent pinning of the interface at the entrance to an expansionwith ? = 30?. ? = 60?, Cn = 0.005, Hu = 5 and H0 = 4.9. . . . . . 414.9 (a) Comparison of wicking speed through 2:1 contractions at twocontraction angles ? = 15? and 25?. The inset illustrates the flowgeometry. Hu = 10, H0 = 9, ? = 60?. (b) Similar comparison for1:2 expansions at ? = 15? and 25?. Hu = 5, H0 = 4, ? = 60?. . . . 424.10 Wicking through multiple contraction-expansion cycles. ForN = 1,2 and 3, ? = 26.6?, 45? and 56.3?. The wetting angle ? = 30?,Cn = 0.02. The total length Ht = 21, and the meniscus starts atH0= 10 at the beginning. . . . . . . . . . . . . . . . . . . . . . 44ixList of Figures4.11 Schematic of a planar microchannel with a Y-shaped bifurcation,showing three stages of wicking: (a) the meniscus reaches the ex-pansion; (b) the meniscus breaks into two at the bifurcation; (c)wicking continues in each branch under suitable conditions. . . . . 454.12 Wicking in both branches with a relatively wide root tube: D0=1.6. The origin of time is when the meniscus first touches the tip atthe junction (cf. Fig. 4.11b). The inset shows that wicking is fasterin the narrow branch initially (t < 11) but the wide branch winnsfor longer times. . . . . . . . . . . . . . . . . . . . . . . . . . . 474.13 Capillary competition between two branches with a relatively nar-row root tube, D0= 0.6. Wicking proceeds in the narrow branchbut is suppressed in the wide branch until tc = 194, marked by adot on both curves. The origin of time is when the meniscus firsttouches the tip at the junction (cf. Fig. 4.11b). . . . . . . . . . . . 484.14 Evolution of the interfacial morphology for the simulation depictedin Fig. 4.13. (a) The meniscus touches the salient corner at t = 0.(b) The meniscus relaxes toward the equilibrium curvature insideeach branch. (c) After a brief retraction, the meniscus is immo-bilized in the wide branch. (d) After the restarting of flow in thewide branch, the menisci advance in both branches. . . . . . . . . 494.15 Flow reversal in the narrow branch when the meniscus in the widebranch moves onto a more hydrophilic portion with ? = 20? atLw = 1.25. Elsewhere ? = 60?. D0 = 0.6, L0 = 4, D2 = 0.9. . . . . 515.1 Schematic of meridian plane of the axisymmetric computationaldomain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 575.2 With a wall relaxation parameter ? = 0.4, the simulation approx-imates experimental results closely in terms of (a) the position ofthe center of the meniscus, and (b) the centerline velocity of themeniscus. The arrows indicate the moment when the contact linereaches the start of the nozzle, and the curves end when a droppinches off, indicated by a filled square. The geometric and physi-cal parameters match the experiment of [86]: Oh = 0.011, L = 5.98,C = 0.493, ? = 23.8? and ? = 0?. In addition, S = 8 ? 10?3 andCn = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 595.3 Comparison of the dynamic contact angle ?d in a straight tube be-tween our numerical simulation and two experimental correlationsdue to [41] and [11]. The model parameters are the same as inFig. 5.2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60xList of Figures5.4 (a) Evolution of the contact line velocity Vw , meniscus center veloc-ity Vc, and average velocity at nozzle exit Vn in time. (b) Evolutionof the dynamic contact angle ?d. (c) Snapshots showing the posi-tion and shape of the meniscus at significant moments marked inthe velocity plot. Oh = 0.01, ? = 0?, L = 5, C = 0.5 and ? = 30?. . 625.5 Temporal variation of the instantaneous velocity Vn(t) at the nozzleexit, starting from point f at tf = 4.85. L = 5, Oh = 0.01, ? =30?, C = 0.5, ? = 0, Le = 1.68, We = 7.0. . . . . . . . . . . . . . . 645.6 (a) Number of drops produced as a function of We. The wideoverlaps between different outcomes indicate that We does not pro-vide an adequate criterion for auto-ejection of droplets. These datacover most of the parameter ranges studied: 0.005 ? Oh ? 0.02,0.25 ? C ? 1, 1 ? L ? 10 and 0 ? ? ? 40?. (b) A short tube(L = 1.5) fails to produce drop ejections under identical conditionsto Fig. 5.4, where a longer tube (L = 5) does produce ejection. Thejet reaches maximum length at t? = 1.07 and then retracts. . . . . 665.7 (a) Criterion for self-ejection: number of droplets plotted as a func-tion of We and the effective length Le of equation (5.5). Thethree outcomes are demarcated by We = 3.4(1+0.8/Le and We =5.5(1+0.8/Le), shown as the solid and dashed curves, respectively.(b) The gray band, representing 5 ? Lj/Rn ? 7 for the jet lengthof equation (5.9), indicates a rough threshold for auto-ejection. . . 675.8 (a) Pressure field inside nozzles with ? = 45? and ? = 30? when jetstarts protruding from the nozzle exit. There is a two-dimensionalpressure filed with a high pressure region around the centreline forthe nozzle with ? = 45?. (b) Comparison of the axial velocityprofile at the nozzle exit between ? = 30? and ? = 45?. . . . . . . 685.9 The regime of rapid ejection at contraction angle ? = 45?, otherconditions being identical to those in Fig. 5.4. (a) t? = 0; (b)ejection of the first drop at t? = 0.15; (c) ejection of the seconddrop at t? = 0.2; (d) ejection of the third drop at t? = 2.37. (e)retraction of the filament. . . . . . . . . . . . . . . . . . . . . . 695.10 Auto-ejection under gravity for large contraction angle. Bo = 0.4,Oh = 0.01, ? = 50?, C = 0.25, L = 2. After ejecting a singledroplet at t? = 0.34, the jet pinches off at its base (t? = 0.47), andlater breaks up into two more drops (t? = 0.5). . . . . . . . . . . 71xiList of Figures5.11 Air entrapment at contraction angle ? = 55?, other conditionsbeing identical to those in Fig. 5.4. (a) t? = 0; (b) formation of airfinger at t? = 0.069; (c) pinch-off of the air bubble at t? = 0.079;(d) jet continues out of the nozzle (e) droplet pinch-off at t? = 1.684after bubble entrapment. Air bubble disappears due to diffusion ofnumerical (f ) filament retraction. . . . . . . . . . . . . . . . . . 725.12 (a) Effect of contact line depinning on the growth of the jet anddrop ejection. The ordinate is the length of the jet measured fromthe exit of the nozzle, and the abscissa is time starting from themoment of the contact line reaching the inner corner of the lip.Drop pinch-off is indicated by a round dot. The width of the uppersurface is fixed at Wl = 0.25Rn, and ?l is the contact angle on theupper surface. The other parameters of the simulation are C = 0.5,? = 30?, L = 4, ? = 0 (on the inner surface) and Oh = 0.01. (b)Snapshots of the interface for ?l = 45? at points marked on the curve. 73xiiNomenclatureSymbol Units (SI) DescriptionA m2 Solid surfacean ?? Radius of thickest part of the torus cross sectionat ?? Radius of thinnest part of the torus cross sectiona0m Initial radius of torus cross sectiona?10m?1 Azimuthal curvature of the torusa(t) m Instantaneous radius of torus cross sectionBo ?? Bond numberC ?? Contraction ratio of the nozzleCa ?? Capillary numberCn ?? Cahn numberd m distance between the droplet and the meniscusDe m Effective widthD1m Diameter of wide branchD2m Diameter of narrow branchD0m Diameter of root tubeE N.m Energy inside the tube-nozzle-reservoir combinationF N.m Total free energyfmix N/m2 Mixing energy of fluid-fluid systemfw N/m Wall energyG N/m2 Chemical potentialGp ?? Grading parameterHb m Position of the center of the meniscusHc m Position of the center of the meniscusHt m Total lengthHu m Length of upstream sectionH(t) m Instantaneous height of the liquid columnHw m Position of the meniscus wall pointH0m Initial position of the meniscush1?? Mesh size at the interfaceh2, h3?? Mesh sizes in two bulksk ?? Wavenumber of the disturbancexiiiNomenclaturel?(t) m Instantaneous wavelength of the disturbancel?0m Initial wavelength of the disturbancel0?? Dimensionless initial wavelength of the disturbancel ?? Dimensionless wavelength of the disturbancelc m Characteristic length scaleld m Diffusion length scaleLe ?? Effective lengthLj ?? Length of the jetLcn ?? Critical length in narrow branchLn ?? Length of liquid column in narrow branchLw ?? Position of the contact line (Chapter 5)Lw ?? Length of liquid column in wide branch (Chapter 4)L0m Initial position of the interfacem ?? Viscosity ration ?? Unit normal vectorN ?? Number of produced droplets (Chapter 3)N ?? Number of contraction-expansion cycles (Chapter 4)Oh ?? Ohnesorge numberp N/m2 Pressurepa N/m2 Ambient pressurepj N/m2 Pressure at junctionpn N/m2 Pressure at the narrower branch in Chapter 4pn N/m2 Pressure at the nozzle exit in Chapter 5pw N/m2 Pressure at the wider branchR m Equivalent drop radius (Chapter 2)Rb m Radius of entrapped bubbleR m Radius of the tubeRn m Radius of the nozzleR(t) m Instantaneous axial radius of the torusR0m Initial axial radius of the torusR?10m?1 Initial axial curvature of the torusR1m Radius of narrower sectionR2m Radius of the wider sectionr, z, ? Cylindrical coordinatesS ?? Dimensionless diffusion lengtht ?? timetc ?? Critical time for flow in wide branchtc sec Capillary time scaletci sec Capillary-inertial timexivNomenclaturetp sec Pinch-off time of the torus (Chapter 3)tp ?? Passage time (Chapter 4)ts sec Shrinkage time of the torust? ?? Time with respect to contact line pinning momentU m/sec Impact velocityv m/sec Velocity vectorvci m/sec Capillary-inertial velocityVc ?? Velocity of the center of the meniscusVn ?? Average velocity at the nozzle exitVw ?? Velocity of the contact linevWL m/sec Lucas-Washburn velocityV0m/sec Average velocity in root tubeV1m/sec Average velocity in wider branchV2m/sec Average velocity in narrower branchWe ?? Weber numberWl m Width of the lipGreek symbolsSymbol Units (SI) Description? ?? Torus aspect ratio (Chapter 1 and 3)? ?? Expansion angle of the Y-branching (Chapter 4)? Pa.sec Viscosity?a Pa.sec Viscosity of the air?m Pa.sec Viscosity of the surrounding medium?t Pa.sec Viscosity of the torus? kg.m?3 Density of the fluid?a kg.m?3 Density of the air?s m Radius of spherical surface? N/m Interfacial tension?w1, ?w2 N/m Solid-liquid Interfacial tensions? m3 fluid domain? ?? Growth rate of the disturbance (Chapter 3)? ?? Contraction or expansion angle (Chapters 4 and 5)?m sec?1 Growth rate of fastest growing wavelength?M ?? Ratio of growth rate terms in Eq. 49 of Ref. [53]?r ?? Ratio between growth rate on a torus and straightfilament? ?? Instantaneous amplitude of the disturbance?c ?? Critical amplitude of the disturbancexvNomenclature?m m Height of the meniscus?0?? Initial amplitude of the disturbance? 1/m Effective curvature? ?? Phase field parameter? N Mixing energy density? ?? Static contact angle?b ?? Breakthrough angle?d ?? Dynamic contact angle?l ?? Static contact angle on the lip?m ?? Maximum angle that interface can reach in the corner m Interfacial thickness? ?? Wall relaxation? m4/N.sec Mobility parameter? m/N.sec Rate of wall relaxationxviList of AbbreviationsGRUMMP Generation and Refinement of Unstructured Mixed-Element Meshes in ParallelAMPHI Adaptive Meshing using a Phase Field (?)GDM Gas Diffusion MediumxviiAcknowledgementsI would like to express my special appreciation and thanks to my advisorprofessor James J. Feng, who has been a tremendous mentor for me. Hisexpertise, devotion to research, sense of responsibility and kindness helpedme to grow as a research scientist.A special thanks to my family. Words cannot express how grateful I amto my mother and father for all of the sacrifices that they have made on mybehalf and for all the prayers that sustained me through my life; and to mymother-in law and father-in-law, whose support and encouragement was abig asset for me during these years. At the end, I would like express mysincerest appreciation to my beloved wife and best friend Faezeh whose loveand patience enabled me to finish this thesis.xviiiChapter 1Introduction1.1 MotivationThe category of problems in which the interaction between forces createdby the surface and bulk flow is important is called interfacial hydrodynam-ics. Beside being scientifically interesting, in recent years due to increasinginterest in small scale applications (like micro-engineering), interfacial flowhas assumed greater significance. Using the Cahn-Hilliard diffuse interfacemethod, we are going to numerically study the following problems:? Capillary breakup of a liquid torus? Wicking flow through microchannels? Auto-ejection of liquid drops from capillary tubesThese problems have been selected based on several considerations. First, inthe framework of immiscible mixtures of incompressible Newtonian fluids,we are interested in fundamental scientific questions about the dynamicsof the interface. Second, these questions, while scientifically significant,must be computationally amenable to our numerical algorithms and codes.Typically, the thickness of the interface in diffuse-interface method shouldbe two orders of magnitude smaller than the bulk to produce physicallymeaningful results. This gives rise to high computational cost and limitsthe amount of inetrfacial area/length of the problem especially for threedimensional configurations. Third, for the problems with moving contactlines, there should be an experimental or a theoretical datum to calibratethe parameters of the numerical model (see section ?2.2 for more details).Here, a brief introduction for each problem and numerical method ispresented. More detailed description of the studied problems and methodis presented in the corresponding sections.As the first problem, we will study the capillary instability of a Newto-nian liquid torus suspended in a surrounding Newtonian liquid. This studyis motivated by the recent experiment of Pairam and Fernandez-Nieves [57]11.1. Motivation(a) (b)Figure 1.1: (a) Schematic of a liquid torus suspended in a bath of surroundingliquid (b) Experimental results of [57] on the breakup of silicone-oil tori suspendedin glycerine. The thin torus breaks down into multiple droplets while the fat torusshrinks into one droplet.on the breakup of silicone-oil tori suspended in glycerine. As shown inFig.1.1(b), they observed that a torus can either shrink to a single dropletor break down into multiple droplets depending on the ratio of its axial ra-dius R0to the radius of its cross section a0(Fig.1.1(a)). There is a clearconnection to the classical problem of Rayleigh-Tomokita instability [65, 81]of a straight, infinitely long filament for which the growth rate depends onthe wavenumber, defined relative to the radius of the filament and also itsviscosity ratio. With the torus, several complications arise. First, due to pe-riodicity of the torus geometry only a few discrete wavelengths are possibleon a torus of given aspect ratio ? = R0/a0. Second, the torus has an axialcurvature which may affect the growth of the capillary wave. In addition,axial curvature produces shrinkage which induces flow in the surroundingfluid which may modify the capillary instability as well [52, 81]. Finallythe presence of the shrinkage puts a limitation on the available time for thedisturbance growth. Here we are going to explore the effect of these differ-ences on the initial growth of small disturbances as well as final breakupand number of produced droplets.The second problem concerns wicking flow through microchannels. Wick-ing is the suction of a liquid by the negative capillary pressure due to themeniscus curvature. It is a key mechanism for flow in porous media [21] like21.1. Motivationwater transport in the gas diffusion medium of proton-exchange-membranefuel cells [54], and in microfluidics for chemical analysis and biological as-say [10, 17]. Porous media and microfluidics devices can have a complexmicrochannel structure through which flow happens in a range of differentgeometrical configurations where the interfacial morphology and motion de-termine to a large degree the efficacy and efficiency of the devices [10, 17].Two types of common geometrical features of microchannels are contractionor expansion in the tube area and interconnectivity of the channels. Thefocus of this study is to investigate the interfacial dynamics in a simplifiedform of these two geometries.Effect of area change on the meniscus dynamics has been studied before[47, 60, 66, 68, 73]. These studies are generalizations of the Lucas-Washburnsolution to tubes and channels of gradually varying cross sections. However,a common feature of all these studies is that they ignore dynamics at themeniscus. The capillary pressure is quasi-statically equilibrated along themeniscus, thus giving it a spherical shape whose curvature is used to computethe capillary pressure via the Young-Laplace equation. Such assumptionis not valid in relatively sharp area changes and the transient dynamicsof the meniscus needs to be considered. As is shown in Fig.1.2, here wewill consider an axisymmetric contraction, expansion and combination ofthem and study the dynamics of an interface through these geometries.Competition among interconnected pores is considered the key mechanismin developing a tree-like morphology of water transport in the gas diffusionmedium (GDM) of fuel cells [54, 59]. To the best of our knowledge suchdynamics has not been studied for wicking flow. Therefore we are going toconsider a two-dimensional Y-shaped branching configuration and study theinterfacial dynamics on and after the junction.As the third problem, we have studied the auto-ejection of drops fromcapillary tubes. The term auto-ejection is used by Wollman and coworkersfor a new mode of drop formation that relies on wicking in a capillary tube[86, 87, 88]. As shown in Fig. 1.3, a glass tube with a tapered end is putinto contact with a reservoir of silicone oil, which wets the glass perfectly.The liquid meniscus rises with sufficient momentum such that a jet is ejectedfrom the nozzle, and later disintegrates into droplets. The sequence of photosshown here were captured under microgravity in a drop tower [88]. Similarexperiments have been done in the International Space Station and undernormal gravity on earth [86, 87]. The process is interesting in that it involvesno external forces or flux, and is entirely autonomous.The critical condition for ejection is the most important question aboutthe auto-ejection process. Wollman et al. [88] have suggested an instan-31.1. MotivationFigure 1.2: Schematic of an axisymmetric contraction and expansion in tube area.taneous Weber number at the end of the nozzle as the parameter whichcontrols the ejection. Experimental data of Wollman and Weislogel [88]shows that such Weber number is not sufficient to quantify the number ofejected droplets. There is a large scatter of data in terms of the instanta-neous Weber number.Another important question is to elucidate the role of the geometric pa-rameters on the auto-ejection. Dynamics of auto-ejection happens throughtwo main stages. The first stage is meniscus rise inside the tube and noz-zle. It is governed by interplay between contact line dynamics, inertia andviscous forces which are themselves a function of the geometric parametersand liquid-solid properties. The second stage of liquid ejection from thetube is governed by the interplay between incoming inertia and capillaryforces. Unlike pulsed ejection in which there is an imposed source of forceor velocity, in capillary ejection velocity decreases as liquid moves out of thetube. Variation of velocity as well as its profile at the nozzle exit depends onthe geometric parameters and fluid-solid properties. The aim of this studyis to develop a criterion for auto-ejection and to study the role of geometricparameters as typically used in experiments.There are theoretical and numerical difficulties in computing interfacialdynamics. These include the lack of a good model for the moving contactline, the need to capture dynamically a moving and deforming interface,morphological singularities in coalescence and rupture of interfaces, andthe complex flow geometries in practically interesting problems like flow41.1. MotivationFigure 1.3: A sequence of snapshots showing spontaneous capillary rise and auto-ejection of droplets in the experiment of [88] under microgravity. The inner diameterof the glass tube is 9.2 mm in the straight section, and the liquid is PDMS ofviscosity 0.65 cs. The drop volume is roughly 20 ?l. Ohnesorge number Oh =0.0015, static contact angle ? = 0o , contraction angle of nozzle ? = 17o, contractionratio of nozzle C = 0.42, tube length L = 8.0 (see ? 5.2 for detailed definition ofthe terms ?,C,Oh, and L). The photos are taken 0.1 s apart. Adapted from [88]with permission, c?Springer.in porous medium. The Cahn-Hilliard diffuse-interface method is a power-ful method in treating such difficulties. As real interfaces are actually thinmixing layers, in this method an interface has a finite thickness and stores amixing energy. The two phases are distinguished by a phase field parameterthat varies smoothly through the interface. When the interface thicknessgoes to zero the model approaches a sharp-interface level set formulation.An advantage of using the diffuse-interface method is that it regularizesthe singular events on the interfaces like breakup, coalescence and movingcontact lines. Thus, this formulation allows us to capture the moving in-terface and its morphological changes accurately and naturally, includingpinning at sharp corners and otherwise singular interfacial breakup at bifur-cations without manual interventions. In the aforementioned problems, thephysical origin of the method can also shed additional light on the underlyingphysics.51.2. Thesis outline1.2 Thesis outlineIn Chapter 2, the Cahn-Hilliard diffuse-interface method is introduced andthe way that the moving contact line is captured in this model is explained.The key features of the finite element solver are explained and finally, twocommon difficulties in using the Cahn-Hilliard diffuse-interface method areexplained in the context of an example problem. Chapters 3 to 5, deal withthe three interfacial dynamics problems to be studied. In Chapter 3, wereport simulation results for the dynamics of a Newtonian torus suspendedin a surrounding Newtonian liquid in three dimensions (3D). This study hasthree man parts. First we will study the linear growth of a sinusoidal distur-bance on the torus and investigate the effect of the retraction and the axialcurvature on the growth rate. Then we will examine the nonlinear insta-bility and the final breakup into droplets in terms of competition betweenretraction and pinch-off mechanism. Finally the numerical results will becompared with the experiment.In Chapter 4, we report numerical simulations of wicking through micro-pores of two types of geometries, axisymmetric tubes with contractions andexpansions of the cross section, and two-dimensional planar channels with aY-shaped bifurcation. At a contraction and an expansion, the dynamics ofthe meniscus at concave and convex corners is illustrated and its effect onpassage time is discussed. For branching geometries, dependence of the flowtrajectory on geometrical and wetting properties of conduits is explained.In Chapter 5, the auto-ejection of drops and jets from capillary tubes isinvestigated. A criterion for ejection is developed in terms of the instanta-neous Weber number at the nozzle exit and the effective tube length. Effectsof a large contraction angle on the meniscus dynamics inside the nozzle andon the auto-ejection are studied. Finally, the effect of the thickness of thetube wall on ejection is explained.Finally, Chapter 6 summarizes the thesis, outlines the significance andlimitations of the current work, and makes recommendations for future work.6Chapter 2Methodology of Research2.1 Governing equationsThe problems that we have chosen to study all fall in the category of im-miscible two-phase flows with both fluids being Newtonian. The flow hy-drodynamics is governed by the continuity equation (2.1) and momentumequation (2.2) for incompressible flows:? ? v = 0, (2.1)?DvDt = ??p+? ? (??v) +G??, (2.2)where v is the velocity, p is the pressure, ? is the density, and ? is viscosity.The term G?? in equation (2.2) represents the role of interfacial tensionin the momentum equation [39, 97] and is obtained by means of variationalcalculus. This is the way that a diffuse-interface model handles the interface,as will be explained below. For the first two problems the Bond number isvery small and the third problem, on auto-ejection from a capillary tube,concerns mostly microgravity experiments. Thus, gravity is ignored in mostof the simulations to be presented. In the first two problems consideredhere the resulting Reynolds number is typically much below unity and theflow dynamics is governed by capillary and viscous forces. Therefore, wewill neglect the inertia term (left hand side of equation 2.2) for the first twoproblems and solve the modified form of the Stokes equation. For the thirdproblem, viscous forces are negligible and the flow dynamics is determinedby the balance between capillary and inertial forces.To capture the interface, the Cahn-Hilliard diffuse-interface method isused. In the diffuse-interface model the two fluid components are viewedas mixing to a limited extent in a narrow interfacial layer. A scalar phasefield ? is introduced to distinguish the components such that ? = 1 in oneliquid, ? = ?1 in the other liquid or gas, and ? = 0 gives the position ofthe interface. In this framework all flow parameters are continuous throughthe interface and instead of solving equations for two parts of the domainseparated by the interface, the same equations are solved for the whole72.1. Governing equationsdomain including the interface itself. In the interfacial region, v may beviewed as a volume-averaged velocity and ? as an average viscosity betweenthose of the two components:? = 1 + ?2?1+1? ?2?2. (2.3)The Cahn-Hilliard model is an energy-based approach to the diffuse-interface method. The total free energy of the system can be written in thefollowing formF =??fmix(?,??)d? +?Afw(?)dA, (2.4)where ? is the fluid domain and A is the solid surface, and fmix is the mixingenergy of the fluid-fluid systemfmix(?,??) =? |??|22+?42 (?2? 1)2, (2.5)where  is the interfacial thickness and ? is the mixing energy density. Inthe limit of thin interfaces, the classical concept of interfacial tension ? canbe recovered from the mixing energy:? = 2?23? . (2.6)fw in 2.4 is the wall energy [40]:fw(?) = ?? cos ??(3? ?2)4+?w1 + ?w22, (2.7)At ? = ?1, i.e. away from the contact line, fw should give the fluid-wallinterfacial tension for the two fluids, ?w1 and ?w2. This requirement leadsto Young?s equation that prescribes the static contact angle ?:cos ? = ?w2 ? ?w1? . (2.8)Using variational calculus, we can calculate the chemical potential as thevariation of the mixing energy with respect to ?:G = ???2?+ ?2?(?2? 1). (2.9)82.1. Governing equationsUsing the gradient of the chemical potential to drive the diffusive fluxesacross the interface, the Cahn-Hilliard equation for the evolution of ? isderived [13]:???t + v ? ?? = ? ? (??G), (2.10)where ? is the mobility parameter and assumed to be constant. Flow hy-drodynamics and interface dynamics is obtained by solving equation (2.10)along with equation (2.1) and (2.2) as a boundary value problem, withboundary conditions specified below.2.1.1 Simulation of moving contact linesProblems studied in Chapter 4 and Chapter 5 involve contact line dynamics.Interaction of molecular and large scale dynamics at the three phase con-tact line makes it a stress singularity point for Navier-Stokes formulation.Different methods use different approaches to handle it. Traditional sharp-interface models typically impose a slip velocity on the wall [102]. Therefore,it is important to explain how the contact line dynamics is captured in theCahn-Hilliard diffuse-interface method.In the diffuse-interface method diffusion across the interface driven bythe gradient of chemical potential (G) allows the contact line movementwithout imposing a slip velocity on the wall. This has two implications, firstit is possible to use no-slip boundary condition on the wall (equation 2.13).Second, the mobility parameter in equation (2.10) has an important role incontact line movement which will be discussed in the next section.Physically it is known that the interface at the wall is almost at equilib-rium. In the Cahn-Hilliard model, fluid-fluid and solid-fluid intermolecularforces are reflected in mixing energy and wall energy respectively. In mostsituations, the fluid layer next to the wall is assumed to be always in equilib-rium with the wall. Mathematically this is reflected by a natural boundarycondition [40, 96]n ? ?? = ?f?w(?)? , (2.11)where the normal vector n points into the solid wall. It is also possibleto allow the wall energy to relax at a finite rate. This leads to anotherboundary condition (equation 2.12) in which ? is the rate constant of wallrelaxationn ? ?? = ?f?w(?)? ?1??(???t + v ? ??). (2.12)92.2. Parameters of the Cahn-Hilliard modelFor ? ? ?, it can be shown that to the leading order, equation (2.12)constrains the dynamic contact angle to be equal to the static one [40, 93].Either of the two equations above, along with the no-slip condition and zeromass flux into the wall,v = 0, (2.13)n ? ?G = 0, (2.14)gives the complete set of boundary conditions for the flow problem.2.2 Parameters of the Cahn-Hilliard modelThe Cahn-Hilliard model formulated above has three parameters that haveno counterparts in conventional Navier-Stokes problems. The interfacialthickness , mobility parameter ?, and rate of wall relaxation ?. Theygive rise to three additional dimensionless parameters: the Cahn numberCn = /lc, diffusion length scale S = ld/lc, and wall relaxation ? = 1/??lc.Cn is the ratio between the interfacial thickness and the macroscopic lengthscale (lc), S is the ratio of the diffusion length ld = ?1/2(?1?2)1/4 to themacroscopic length [96], and ? shows how fast flow close to the wall equi-librates with the wall. The choice of these three parameters (Cn, S, and?) is informed by their physical meanings and the requirement of achievingthe sharp-interface limit. Real interfaces, a few nanometers in thickness, aretypically not resolvable in macroscopic flow simulation. Thus the diffuse-interface method uses an artificial  that may be much larger than the realvalue. This is allowable if the sharp-interface limit is achieved: when and Cn are sufficiently small such that the results are not affected by theunrealistic thickness of the interface [12, 96]. For interfacial flows with-out contact lines such as the capillary instability of a suspended torus, thesharp-interface limit is typically approached at Cn ? 0.01 [97]. In such flowsthe Cahn-Hilliard diffusion across the interface, represented by ld or S, isimmaterial as long as the sharp-interface limit is achieved.With moving contact lines, Yue et al. [96] have shown that the achieve-ment of sharp-interface limit is also dependent on the diffusion length. ForCouette and Poiseuille flows with a transverse interface, they found that theshape of the meniscus converges to a unique solution after Cn falls belowa threshold Cn ? 4S and suggested this as the criterion for achieving thesharp-interface limit. On the other hand, the motion of the contact lines isaffected by interfacial diffusion and wall relaxation. Thus, S and ? mustbe chosen judiciously, e.g., to coincide with an experimental measurement102.3. Discretization of governing equations and computational domain[96]. We use the following procedure proposed by Yue and Feng [93] toset values of S and ?. Choose as small an Cn value as computationallyaffordable, and then pick a value for S to ensure the sharp-interface limit isachieved. Finally, determine the wall relaxation parameter ? by fitting anexperimental datum. For studying the wicking flow, contact line movementhappens very slowly and meniscus keeps its equilibrium shape. Thereforewe have set ? equal to zero for the problem studied in Chapter 4. For theauto-ejection problem, due to fast movement of the meniscus, there is aconsiderable deviation from the static shape and ? should have a non-zerovalue.2.3 Discretization of governing equations andcomputational domainThe governing equations are solved using a finite-element package calledAMPHI in planar 2D, axisymmetric and 3D geometries. [97] and [100]have described the numerical algorithm in detail, presented numerical ex-periments on grid and time-step refinements and validated the methodologyagainst numerical benchmarks. The numerical package has been applied toa number of interfacial-dynamics simulations [1, 27, 28, 96, 101]. Here, wewill only mention a few important features.The discretization of the governing equations follows Galerkin formalism.The fourth-order Cahn-Hilliard equation is decomposed into two second-order equations:???t + v ? ?? =??2 ?2(? + s?), ? = ?2?2?+ (?2 ? 1? s)?, (2.15)where s is a positive number to enhance the stability of the numerical methodand set to 0.5 [69]. Piecewise quadratic (P2) elements are used for velocity,? and ? and piecewise linear (P1) elements are used for pressure. Time-stepping is done in a second-order implicit scheme. The nonlinear algebraicsystem generated from the weak form of the equations is solved using New-ton?s method. The equations are discretized on an unstructured grid andthe finite thickness of the interface is resolved using adaptive meshing, withthe grid being dynamically refined and coarsened, respectively, upstreamand downstream of the moving interface. This is done by a general-purposemesh generator called GRUMMP that produces a triangular mesh in 2Dand a tetrahedral mesh in 3D [26] by Delaunay triangulation. It controlsthe grid size by using a scalar field, which can easily be computed from the?? field of the diffuse-interface method.112.4. Limitations on resolving small lenght scalesFour parameters are used to specify the domain discretization: meshsize at the interface (h1), mesh sizes at the two bulks (h2and h3), and agrading parameter (Gp) which controls the transition between the bulk andinterface mesh size. For all of the results to be presented, we have usedsufficiently refined spatial and temporal discretization to ensure numericalconvergence. This is confirmed by grid and time-step refinements to verifythat the results change little with further refinement. A few typical spatialand time steps are indicated below as examples. For capillary breakupof liquid torus (Chapter 3), mesh sizes inside the torus, at the interface,and in the surrounding liquid are chosen to be equal to 1/5 of the torusthickness, interfacial thickness (), and 1/2 of torus thickness, respectively.Grading parameter Gp is chosen to be 3. For wicking flow in microchannel(Chapter 4) and auto-ejection of drops from capillary tubes (Chapter 5), thebulk mesh sizes inside the tube and nozzle are chosen to be 1/8 of the tubeand nozzle exit radii, respectively. Mesh size at the interface is chosen tobe /2 and grading parameter for these problems is equal to 5. The typicaltime-step of the simulations for problems in Chapters 3, Chapter 4, andChapter 5 is 10?2, 5? 10?3, and 5? 10?4, respectively.2.4 Limitations on resolving small lenght scalesIn this section two types of difficulties in using the Cahn-Hilliard diffuse-interface model will be discussed. These difficulties happen due to limi-tations on resolving the small length scales which are not specific to theCahn-Hilliard model. The real diffusion length of the contact line is aroundsix orders of magnitude smaller than the bulk length scale which is not com-putationally resolvable. Therefore, there are two parameters in the model?, ? whose values are not known and they are important to the contactline dynamics. To capture the interfacial dynamics using a large diffusionlength, it is required to calibrate the diffuse-interface parameters using anexperimental or a theoretical data point. Such data points are not availablefor all physical problems.The second difficulty, is the high computational cost of the diffuse-interface method. Numerically capturing an interfacial region which isaround six orders of magnitude smaller than bulk region is almost impos-sible with current computational resources. Therefore, the sharp-interfaceconcept is used to show that with much thicker interfaces, it is possible toget physically meaningful and numerically converged results. Despite beinga remedy for many interfacial problems, such an approach fails when there122.4. Limitations on resolving small lenght scalesare disparate length scales in the problem which are important to the finaloutput of the simulation. An example is to capture the physical quantitieswhich are very sensitive to coalescence/pinch-off phenomenon. Physicallysuch a process depends on the intermolecular forces. Sharp-interface meth-ods face a singularity in dealing with coalescence/pinch-off which is usuallyhandled by imposing a cut-off length. In the diffuse-interface method, thereis no need to impose an extra cut-off condition and such length is determinedby the diffusive nature of the method and finite thickness of the interface.For some problems such natural treatment of the coalescence/pinch-off is de-sirable like the filament breakup for capillary instability of a torus (Chapter3), and droplet pinch-off for the auto-ejection (Chapter 5). But sometimesthe desired physical quantity is strongly dependent on the inherent cut-offlength of the diffuse-interface method. One such quantity may be, for ex-ample, the duration of a drop pinch-off process. As the characteristic lengthscale eventually goes to zero, no numerical scheme can follow the process?exactly? till the end. In reality, short-range molecular effects come in atsome point to dominate the rest of the process. Then it is a subtle pointas to how a diffuse-interface model might represent such a process. In thefollowing, we illustrate both difficulties by a concrete example.2.4.1 Bubble production in meniscus sessile dropletcollisionLet?s apply the Cahn-Hilliard diffuse-interface method to study air-entrapmentwhen a nearly flat liquid meniscus impacts a sessile droplet. Such an impactmay happen in dip-coating [6, 19], immersion lithography [75] and resultsin a production of small air bubbles which are not desirable. Keij et al. [43]have studied this problem experimentally and showed that when an inter-face impacts a sessile droplet, there are two possible scenarios: coalescenceof the droplet and meniscus on or close to the contact line and coalescencefarther above the substrate. The first case will result in either no air bubbleformation or appearance of a floating air bubble. The second case producesa bubble that sticks to the substrate. It is shown that the size of the floatingbubble increases with increasing the capillary number and the size of stick-ing bubbles is independent from the impact velocity and is about an orderof magnitude smaller than the first.We aimed to carry out numerical simulations in both two (Fig. 2.1(a))and three dimensional setup (Fig. 2.1(b)) in which a meniscus with a velocityU and contact angle of ? will impact a sessile droplet with equivalent radiusR. As is shown in Fig. 2.1(a) a simple shear flow is imposed at the inlet132.4. Limitations on resolving small lenght scales(a) (b)Figure 2.1: (a) A meniscus with a velocity U and contact angle of ? will impacta sessile droplet with a equivalent radius R (b) Three dimensional schematic of theproblem setup.and it is assumed that the meniscus has reached its steady shape beforethe impact. The nondimensioanl parameters of the problem are capillarynumber, viscosity ratio, contact angle, and gap size. The main desiredoutput is how the coalescence, bubble shape and its volume vary with thegoverning parameters. We also assume rapid relaxation of the wall energyso ? is equal to zero. We were unsuccessful in simulating both floating andsticking bubbles, reasons for which are explained in the following sections.2.4.2 Dependence of bubble dynamics on meniscus shapeThe main result that we want to capture in this study is the bubble sizeand its position, which depend on whether the coalescence between themeniscus and the drop happens on the wall or at some point above the wall[43]. Position of the coalescence is dependent on the shape of the meniscuswhich itself is dependent on the value of the diffusion length scale. It isshown by Yue et al. [96] that in shear flows the meniscus shape or interfaceinclination is highly sensitive to S. In Fig. 2.2, meniscus shapes for twodifferent values of the diffusion length scale are shown. It can be easilyseen that at S=0.0025, collision on the substrate never happens while forS=0.04, drop and meniscus will collide on some point above the meniscus.In addition, the amount of the air trapped in the bubble depends on theshape of the meniscus. This means that a slight change in the value ofthe diffusion length can change the shape of the meniscus and hence thetype of the coalescence and the size of the bubble considerably. Thereforein order to get physically meaningful results, it is very important to havea good estimate for the value of S based on an experimental datum. Theexperimental data point is available for a case in which the bubble size is at142.4. Limitations on resolving small lenght scalesxy4 5 6 7 8 9 100123S=0.0025S=0.04Figure 2.2: Effect of diffusion length scale on the coalescence type.least ten times smaller than the drop radius. This brings us to the seconddifficulty which is discussed in the next section.2.4.3 Dependence of bubble volume on cut-off lengthAs mentioned before, the important quantity that we want to capture is thevolume for the bubble, which equals the amount of air entrapped betweenthe meniscus and the droplet. Therefore, the distance between droplet andmeniscus (d) at the moment of pinch-off determines the volume of the bub-ble. If we assume that air entrapped between the droplet and interface is anair sheet with the thickness d and area equal to the front face of the droplet,then the radius of the bubble can be roughly related to the distance betweenthe meniscus and droplet asRbR ? (dR )13 (2.16)where Rb is the radius of the bubble. The largest reported bubble sizein the experiment is around one tenth of the droplet size. Therefore, tocapture such bubble size d/R ratio should be at least 1 ? 10?3. Accordingto numerical results, two interfaces start affecting each other when theyare approximately 5 apart ( is the interface thickness) for S = 0.01. Thismeans that the Cn number should be around 2?10?4. The smallest possibleCn achievable in 3D simulation is 0.015. Therefore, capturing such smallbubble size is impossible with current computational resources.Problems to be studied in Chapter 3 to Chapter 5, are not subjectedto the mentioned difficulties where obtained results are not sensitive to thepinch-off/coalescence time/length scale. Problems in Chapters 4 and 5 in-volve contact line movement for which there are theoretical and experimentaldata points respectively.15Chapter 3Capillary Breakup of aLiquid Torus3.1 IntroductionAs discussed in Chapter 1, we want to study the capillary breakup of aliquid ring suspended in a surrounding liquid which has a main additionalgeometric characteristic compared to a straight filament: the presence ofaxial curvature around the center of the torus. The axial curvature modifiesthe curvature driven capillary instability and produces shrinkage which im-poses an external flow around the torus, changes the instantaneous geometryof torus and puts a time limitation on disturbance growth.The initial growth rate of the disturbances can be related to linear growthrate of the liquid filament which is a classical problem in fluid mechanics[24, 71]. A long cylindrical liquid thread becomes linearly unstable to distur-bances with a wavelength longer than the circumference of the thread 2?a,a being the radius of the filament. The most unstable wavelength is 9.02afor an inviscid filament [65], and longer and dependent on the viscosity ratiofor a viscous thread in a viscous surrounding fluid [80]. The capillary wavesgrow into the nonlinear regime and ultimately lead to breakup, and satellitedrops may appear depending on the viscosity ratio [79]. Thus, capillarybreakup of long straight filaments is well understood.In comparison, we have a rather limited knowledge of the stability ofcurved filaments. Experimentally, Pairam and Fernandez-Nieves [57] stud-ied the retraction and breakup of Newtonian tori in a Newtonian surround-ing liquid. McGraw et al. [50] and Wu et al. [89] further considered thebreakup of nano-scale polymer and liquid metal rings on solid substrates.Several theoretical and numerical studies have appeared in the literature,and most of these have dealt with the more complicated situation of a liq-uid ring or torus in contact with solid substrates. For instance, Wu [90]computed the Rayleigh modes on a liquid ring spreading on a solid afterimpingement. Bostwick and Steen [9] considered the static stability of the163.1. Introduction(a) (b)Figure 3.1: (a) A quarter of the top half of a liquid torus for simulating thecapillary growth of an even mode, i.e. with an even number of wavelengths aroundthe torus. For odd modes, a half of the top half must be included. (b) The interfaceon the symmetric mid-plane with and without a sinusoidal disturbance.so-called torus lift, a liquid ring constrained by a solid ribbon in contactwith part of the liquid surface. Nguyen et al. [55] carried out molecular-dynamics and long-wave continuum simulations of the capillary breakup ofa nano-scale liquid metal ring on a solid surface. Gomes [32] computedthe stability of a rotating toroidal gas bubble constrained between two con-centric cylinders. The baseline situation, of a freely suspended torus in aquiescent medium, seems to have been studied only in Yao and Bowick [92];they solved the Stokes flow during the contraction of the torus but did notinvestigate its capillary instability.In this study we simulate the dynamics of a Newtonian torus suspendedin a surrounding Newtonian liquid in three dimensions (3D). First we willstudy the linear growth of a sinusoidal disturbance on the torus and investi-gate the effect of the retraction and the axial curvature on the growth rate.Then we will examine the nonlinear instability and the final breakup intodroplets. Finally the numerical results will be compared with the experi-ment.173.2. Problem setup3.2 Problem setupConsider a Newtonian liquid torus of viscosity ?t suspended in an immiscibleNewtonian medium of viscosity ?m. Initially the cross section of the torusis a circle of radius a0, and the axis through the centre of the cross section isa circle of radius R0. Hereafter, we refer to the curvatures due to R?10anda?10as the axial curvature and azimuthal curvature, respectively. Althoughnon-varicose modes of instability are possible under external forcing, theexperiments showed only varicose necking and breakup. Thus, we assumesymmetry about the mid-plane of the torus, and only need to consider itstop half. Furthermore, we can compute a half or a quarter of the top halffor the growth of odd and even sinusoidal modes (Fig. 3.1). A sinusoidalperturbation of wavelength l?0is imposed on the torus at the start:(r ?R0)2 + z2 = a20[1 + ?0cos(2?R0l?0)?]2, (3.1)where r, z, and ? show the surface of the torus in cylindrical coordinates,k = 2?R0/l?0is the number of waves along the circumference 2?R0, and?0is the initial dimensionless amplitude. In presenting results, k will becalled the wave number, though it differs from the usual sense of the word(2?a0/l?0). We use the subscript 0 to indicate the initial condition. Withcontraction of the torus and growth of the disturbance, a(t), R(t), l?(t) and?(t) all change in time.The subsequent fluid flow is governed by the Stokes equation; inertia andbuoyancy are negligible in the experiment and will be neglected in the com-putations. For boundary conditions, we assume symmetry on the bottomand planar side walls of the domain of Fig. 3.1(a). The top wall is 11a0abovethe top of the torus, on which we impose zero stresses. The outer cylindricalwall is at least 10a0from the torus, and is solid with vanishing velocity. Theouter boundaries are sufficiently removed from the torus that they do notaffect the retraction and capillary instability on the latter. Toward the endof the chapter, when trying to match the experimental geometry of [57], wewill bring the side wall closer to the torus.Two dimensionless numbers quantify the physical problem: the torus-to-medium viscosity ratio m = ?t/?m and the initial aspect ratio of thetorus ? = R0/a0. The Cahn-Hilliard model introduces two more parame-ters, the Cahn number Cn = /a0and a diffusion length scale S = ld/a0.We have used S = 0.02 and Cn = 0.05 throughout this chapter; this en-sures the attainment of the sharp interface limit during torus retraction.The final breakup involves length scales shrinking to zero, and the finite183.3. Results: linear growth of capillary wavesthickness of the interface and the diffusion within will eventually manifestthemselves. With Cn = 0.05, numerical experiments show that the pinch-off time increases by less than 5% when S decreases from 0.02 to 0.004. Inpresenting results, we use a0as the characteristic length and the capillarytime tc = a0?t/? as the characteristic time, ? being the interfacial tension.The wavelength l, however, will be scaled by the instantaneous circumfer-ence of the cross section of the torus 2?a to facilitate comparison with thestraight-filament results. Note that tc characterizes the capillary waves onthe torus. Its retraction in the presence of a viscous external fluid is on thetime scale (R0? a0)?m/? = tc(? ? 1)/m.3.3 Results: linear growth of capillary wavesCompared with the Rayleigh-Tomotika instability on a straight filament,several complications arise on the torus. First, due to the finite circumfer-ence of the torus, only a number of discrete wavelengths are possible for agiven aspect ratio ?. Second, the torus has an axial curvature (R?1) whichmay affect the growth of the capillary wave. Finally, the contracting torusinduces a flow in the surrounding fluid which may modify the capillary in-stability as well [52, 81]. Under the constraint of quantized wavelengths, thelast two effects will be explored separately.3.3.1 Quasi-static retraction: effect of axial curvatureBy choosing a large initial aspect ratio ? and a small viscosity ratio m, wecan separate the time scales for the growth of the capillary wave and theretraction of the torus. In physical terms, this corresponds to a thin torusretracting slowly in a highly viscous bath. The speed of retraction dR/dtdecreases in time. As an indication of its magnitude, dR/dt = ?0.0036at R = 4 for m = 0.033. For larger m, the retraction speed increases inproportion as expected. Such a quasi-static process is convenient in that wecan probe the effect of the axial curvature on the linear instability of thetorus while excluding the dynamic effect of the retraction-induced externalflow. Furthermore, if we use a small enough initial perturbation and carryout the simulations on the time scale of torus retraction tc(??1)/m, we canrecord the linear growth rate at different axial curvatures and wavelengths.Thus a dispersion relation can in principle be generated in one simulation.For one such torus with initial aspect ratio ? = 5.3 and viscosity ratiom = 0.033, we impose two wave forms on it (k = 2). Different initialamplitudes (?0= 0.005 and 0.01) are tested, and ln(?/?0) initially grows193.3. Results: linear growth of capillary waves(a) l?0.5 1 1.5 2 2.5 300.0010.0020.0030.0040.0050.006Torus (k=2)Straight filament(b) 1/??0 0.05 0.1 0.15 0.2 0.25 0.30.0040.00450.0050.00550.006k=2k=4k=6Figure 3.2: (a) Dispersion relation on a shrinking torus compared to that for astraight filament. The latter is computed by our diffuse-interface method and agreeswith the Tomotika formula within 4%. The wavelength l and the growth rate ? aremade dimensionless by the instantaneous 2?a and tc, respectively. (b) The lineargrowth rate decreases with the axial curvature for a prescribed dimensionless wavelength l0= 2. The point at 1/? = 0 corresponds to a straight filament.linearly in time with a slope ? that is independent of ?0. This confirms thatwe are in the linear regime, with ? being the growth rate. Over longer times(on the order of tc(?? 1)/m ? 100tc), the growth rate remains independentof ?0but starts to change in time. This is an effect of the torus retractioneven though the instability is still in the linear regime. Since the wavenumber k = 2 is fixed, the wavelength shrinks with the retraction, not onlyin dimensional terms, but also relative to the thickening filament radius a.Thus, recording the growth rate as a function of the changing wavelengthproduces the dispersion relation in Fig. 3.2(a). The growth rate on the torusis some 15% below that on the straight filament, although the difference isexpected to diminish for larger ?. For ? ? 10, the difference narrows downto within 5%. In the limit of R0 a0, of course, one recovers the growthrate on a straight circular cylinder. Therefore, the axial curvature on thetorus tends to hinder the growth of the capillary waves. Note also that boththe minimum wavelength for instability and the fastest growing wavelengthhave shifted slightly to longer waves from those for the straight filament.The simulation above is not ideal in quantifying the effect of the ax-ial curvature R?1 on the growth rate ? since the former cannot be pre-scribed but continues to increase in time. For this purpose, we have con-ducted a series of simulations with tori of the same initial a0, but different203.3. Results: linear growth of capillary waves(a) m?r0 0.05 0.1 0.15 0.2 0.25 m?M0 0.1 0.2 0.311. 3.3: (a) Ratio of growth rates on a torus as a function of the viscosity ratiofor a capillary wave of dimensionless wavelength l = 2. (b) Ratio of growth rates ona straight filament under uniform extensional flow, calculated from the theoreticalresult of [52].initial aspect ratio ? in proportion to the wave number k. Thus, thesecapillary waves have the same initial wavelength [in dimensionless forml0= (2?R0/k)/(2?a0) = ?/k = 2], and differ only in the axial curvatureR?10. Fig. 3.2(b) plots the initial linear growth rate ? as a function of1/? = a0R?10, the non-dimensionalized axial curvature. It shows unequivo-cally that the instantaneous growth rate decreases with the axial curvature.3.3.2 Faster retraction: effect of external flowTo examine the effect of the external flow field on capillary instability ofthe torus, we have gradually decreased the viscosity of the suspending fluidto produce faster retraction of the torus. Even on a straight filament, inthe absence of the flow effect being examined, the ambient viscosity wouldhave affected the growth rate. To remove this effect and isolate that ofthe retraction-induced external flow, we compute the ratio ?r between thegrowth rate on a retracting torus and that on a straight filament, the latterbeing calculated from the Tomotika formula using the same viscosities andthe instantaneous filament diameter and wavelength of the torus. This ra-tio, as a function of m, demonstrates how the flow affects the growth of theinstability. Note that the torus viscosity ?t remains unchanged in this pro-cess; it gives a fixed time scale tc against which the growth rate is measured.The faster retraction is then indicated by an increasing viscosity ratio m.213.4. Results: nonlinear growth and breakupFig. 3.3(a) plots the ratio of growth rates ?r against the viscosity ratiom for a dimensionless wavelength l = 2. With increasing m and hence in-creasing retraction speed, the growth rate ratio increases. This implies thatthe external flow induced by the torus retraction has the effect of enhancingthe growth of instability. That ?r is below unity reflects the quasi-staticeffect of the axial curvature discussed in the preceding subsection.It is interesting to compare this flow effect with that on a straight fila-ment. [52] computed the effect of a uniform extensional flow on the capillaryinstability on a straight filament. The growth rate is written as the sum oftwo terms (see their Eq. 59). The first, due to the thinning of the filamentand advective lengthening of the wavelength, had previously been computedby [81]. This effect is quasi-static in nature, and its counterpart on the torushas been included in the analysis of the last subsection. The second term,proportional to the strain rate G, explicitly accounts for the flow effect.From our torus retraction simulation, we extract a negative G from the rateof filament thickening, and then compute the two terms for the same wave-length l = 2. We take the ratio between the total growth rate and the firstterm, and plot it as a ratio of growth rates ?M in Fig. 3.3(b). This is notthe same ratio as that in Fig. 3.3(a) since there is no axial curvature. Nev-ertheless, the qualitative trend is clear and confirms our observations on theretracting torus: the compression of a straight filament enhances the growthof capillary instability.3.4 Results: nonlinear growth and breakupThe nonlinear instability and breakup of the torus must take place beforethe torus contracts onto itself. In this process, the quantized wavelengthavailable and the initial amplitude of the perturbation are both importantfactors. Besides, the initial aspect ratio of the torus and the viscosity ratioare key parameters.3.4.1 Fastest modeOn a retracting torus, with the wavelength and filament thickness changingcontinually, the initially dominant mode does not necessarily persist tillbreakup. In fact, the torus retraction should favor initially longer waves andthis is illustrated in Fig. 3.4, with ? = 6.7, m = 0.033 and ?0= 0.02. Basedon the dispersion relation for the torus, the linearly dominant wavelength isl = 2.03 and corresponds to a wave number k = 3.3. Thus, k = 3 or k = 4should initially produce the fastest growth. Indeed, the two modes grow223.4. Results: nonlinear growth and breakupt?0 200 400 600 800 1000 120000. 3.4: Nonlinear evolution of three modes of instability, with wave numberk = 2, 3 and 4, for ? = 6.7, m = 0.033 and initial amplitude ?0= 0.02. ? is theinstantaneous amplitude of the capillary waves. The curves for k = 2 and k = 3end in breakup, with the onset of secondary necking also marked on the latter. Thek = 4 mode ends in complete retraction.at comparable rates at the beginning. But as the torus shrinks, the k = 3mode maintains a high growth rate while the growth rate for k = 4 declines,leading eventually to retraction, not breakup. This can be rationalized bynoting that for a retracting torus with a fixed wave number k, the wavelengthgets shorter in time, in dimensional terms and especially relative to thegrowing thickness a. Thus the initially longer wave (k = 3) is favored overthe shorter one (k = 4). The k = 2 mode grows more slowly but doeslead to breakup. The breakup of the torus into droplets is depicted bysnapshots in Fig. 3.5 for k = 3, starting from an initial perturbation ofamplitude ?0= 0.02. Primary necking proceeds at three points around thecircumference of the torus until t = 678, when two secondary necks emergearound each primary neck. At t = 748 the torus breaks down into threeprimary drops and three satellite droplets. In time these all relax towardthe spherical shape.3.4.2 Pinch-off time versus retraction timeFrom the preceding discussion, it is clear that the breakup of the torusdepends on the competition of two time scales: tp needed for the neck to233.4. Results: nonlinear growth and breakup(a)xy2 4 6 8-8-6-4-202468t=0(b) xy2 4 6 8-8-6-4-202468t=678(c) xy2 4 6 8-8-6-4-202468t=748(d)xy2 4 6 8-8-6-4-202468t=784Figure 3.5: Snapshots of the evolving interface on the mid-plane of the torus for? = 6.7, m = 0.033 and ?0= 0.02. The interface is given by the level set of ? = 0.(a) ?0t p0 0.04 0.08 0.12 0.16300400500600700800900no pinchoff(b) ?? c3 3.5 4 4.5 5 5.5 6 6.5 700. 3.6: (a) The pinch-off time decreases with increasing initial amplitude ofdisturbance. ? = 5.3, m = 0.033, k = 2. The solid curve is the best fitting byEq. (3.2). (b) The critical initial amplitude ?c decreases with the initial aspectratio ?. The solid curve is the best fitting by Eq. (3.3).pinch off, and ts needed for the torus to shrink onto itself. This competitioncan be affected by multiple factors. For example, the k = 4 mode of Fig. 3.4can survive till breakup if the initial perturbation has a sufficiently largeamplitude; ?0defines tp. Besides, the breakup depends on the initial aspectratio ? and the viscosity ratio m, each having a role in ts. These threefactors will be examined in turn.Fig. 3.6(a) demonstrates the dependence of the pinch-off time tp on the243.4. Results: nonlinear growth and breakupt?0 200 400 600 800 100000. 3.7: Effect of the viscosity ratio on the growth of disturbance. ? = 5.3,?0= 0.02 and k = 2.initial amplitude ?0for ? = 5.3, m = 0.033 and k = 2, which is the initiallydominant mode. If ?0is below a critical value ?c ? 0.02, no breakup occurs.For ?0> ?c, the torus breaks up into two principal drops and two satellitedroplets, and tp decreases with increasing ?0 as expected. Besides, the fasterthe breakup, the larger the satellite droplets. The critical amplitude ?cdecreases with increasing initial aspect ratio ?, as shown in Fig. 3.6(b). Thethinner, longer torus offers a longer ts within which breakup can take place.In the ? range shown, k = 2 persists till breakup from all ?0> ?c; no othermodes emerge from noise to overtake the imposed k = 2 mode.The viscosity ratio m = ?t/?m is another parameter that modulatesthe competition between pinch-off and retraction. Our results show thatthe torus retraction is more influenced by the matrix viscosity ?m while thenecking and pinch-off more by the torus viscosity ?t. As m increases from0.033 to 0.05 and 0.1, the critical amplitude ?c increases from 0.02 to 0.03and 0.07. For m = 0.5 even ?0= 0.18 is unable to break down the relativelyviscous torus before it contracts into a single drop, often entrapping a dropletof the ambient fluid in the centre [97].Fig. 3.7 illustrates the effect of m on the growth of an initial distur-bance with k = 2, which is the initially dominant mode for all the m valuesconsidered here. Since time is scaled by tc = a0?t/?, using the torus vis-cosity, increasing m can be conveniently thought of as due to a decreasing?m. As ?m decreases, the initial growth rate of the capillary wave increases.However, the retraction of the torus becomes faster as well. Numerical ex-253.5. Comparison with experimentperiments show that the latter has the upper hand. Thus, for lower ?m, ?reaches a maximum quickly and then declines, due to the thickening of thetorus and the effective shortening of the wavelength. It is for the largestmatrix viscosity, at m = 0.033, that the slow retraction offers the capillarydisturbance sufficient time to grow till breakup, despite the slower lineargrowth rate.The competition between time scales can be represented by scaling argu-ments. As noted earlier, the shrinkage time ts ? tc(?? 1)/m. The pinch-offtime can be taken as that required for the disturbance to grow from the non-dimensionalized initial amplitude ?0to 1: tp = ? ln ?0/?m, where the fastestgrowth rate ?m can be estimated from the Tomotika solution: ?m ??m/tc[18]. Therefore, we can writetp = tcc1?m ln(1?0), (3.2)and c1= 40.6 gives a reasonably good fitting to the numerical data inFig. 3.6(a). Furthermore, equating this tp with the shrinkage time ts givesus the critical initial amplitude for breakup:?c = exp(?c2? ? 1?m), (3.3)which fits the data in Fig. 3.6(b) well with c2= 0.16. Given that much of thenecking and pinch-off is nonlinear, these linearly based scaling relationshipswork remarkably well.3.5 Comparison with experimentAs far as we know, the only prior experiment on the breakup of a freely sus-pended torus is that of Pairam and Fernandez-Nieves [57]. With Newtonianglycerol tori in a Newtonian oil bath, these authors reported that thick torishrink to one droplet while thin ones break down into a number of dropletsthrough Rayleigh-Tomotika instability. We match the liquid viscosities andflow geometry in the experiment, where the torus is confined in a cylindricaldrum, with the top and side walls being some 6a away from the outer edgeof the torus. Our numerical experimentation shows that this confinement isessential for slowing down the torus retraction and allowing breakup. Stilltwo uncertainties complicate a direct comparison. The first is the initialamplitude of perturbation ?0. In the experiment, the torus is generated byreleasing a glycerine jet into silicone oil while the drum rotates. There is263.5. Comparison with experiment(a) ana t0.1 0.3 0.5 0.7 0.9 ta n0 5 10 15 20 25 3000. 3.8: (a) Determining the initial amplitude of perturbation ?0from thevariation of the thickest radius at versus the thinnest radius an on the torus. Bothradii are normalized by the initial value a0. (b) Determining the interfacial tension? from the temporal variation of an. ? = 5.3, k = 2 and m = 0.033.a complex flow history, and it is not obvious how to gauge the magnitudeof the initial perturbation. The second is the interfacial tension ? in theexperiment. It was not reported and cannot be made available to us. Wedetermine ?0and ? first by fitting the experimental data.First note that the capillary time tc is the only time-scale of the problem,and the only role of ? is to lengthen or compress tc. Thus, in Fig. 3.8(a) weplot the radius at of the thickest part of the torus against the thinnest radiusan at the neck. Such a curve should be independent of tc. Among numericalresults starting from different ?0values, ?0= 0.01 agrees very closely withthe experiment. So we take ?0= 0.01 to be the initial amplitude for thiscase. Now plotting the temporal variation of the neck radius in Fig. 3.8(b)gives us a fitting of ? = 31.8 mN/m, which is within the 5 percent of thehandbook values [63].With the ?0and ? values determined, we compare the number of pri-mary drops N between the simulation and the experiment for a range oftorus aspect ratio ? (Fig. 3.9). All the simulations have started with thefastest linear mode for the ? value. The results agree with the experimentexcept for ? = 4, where the simulation predicts complete retraction, whilethe experiment reported N = 1, breakup at a single primary neck for thek = 1 mode. We cannot explain this at present; possibly this experimenthad a different ?0from that fitted in Fig. 3.8(a) for ? = 5.3. Numericalexperimentation indicates that ?0= 0.02 would lead to breakup at a single273.6. Summary and conclusions?N0 2 4 6 8 10012345Numerical resultsExperimental resultsFigure 3.9: Comparison between the predicted and observed number of primarydrops after breakup, for tori with five initial aspect ratios. m = 0.033 and ?0= 0.01.N = 0 and 1 refer to, respectively, complete retraction with no breakup and breakupat a single primary neck.neck. In all the cases leading to breakup, N corresponds to the fastest linearmode. Even though the wavelength and filament thickness both change dur-ing the retraction, we have never seen the linearly dominant mode yieldingto a nascent mode in the nonlinear stage. This reflects the fact that thereis a limited time window for growth and it is too short for another mode toemerge spontaneously from random noise.3.6 Summary and conclusionsCapillary instability of a Newtonian liquid torus is studied by imposing asinusoidal disturbance on an initially stationary torus in quiescent surround-ing Newtonian liquid. The main findings of this study are(a) If the initial disturbance has sufficiently small amplitude, the capillaryinstability initially grows linearly, at a growth rate independent of theamplitude of the disturbance.(b) The geometry of the torus can affect the capillary instability throughits axial curvature and the external flow field.283.6. Summary and conclusions(c) Effect of axial curvature on the disturbance growth is studied by choos-ing a very viscous surrounding liquid and hence producing very slowretraction. A dispersion relation is constructed for the torus by look-ing at the growth rate of the disturbance during the shrinkage. It isshown that axial curvature decreases the growth rate of the distur-bance compared to straight filament.(d) Effect of the external flow field on capillary instability of the torusis examined by gradually decreasing the viscosity of the surroundingfluid to produce faster retraction of the torus. By finding the ratio ofthe growth rate of a certain wavelength on a torus to its counterpart ona straight filament, it is shown that shrinkage enhances the capillaryinstability. This is further confirmed by comparing the flow inducedterm and quasi-static term in the equation proposed by [52] for theflow effects in a straight filament.(e) The final shape of torus is a result of competition between two time-scales called shrinkage and pinch-off timescale. Shrinkage can be scaledas ts ? tc(? ? 1)/m which shows that it is inversely proportional toviscosity ratio and is proportional to aspect ratio of the torus as wellas capillary timescale. The pinch-off time depends on the wavelength.If we choose the fastest one and assume that the growth rate is almostequal to the initial growth rate during the retraction. The follow-ing scaling can be proposed for the pinch-off time: tp ? tc?m ln(1?0).Competition of these two time-scales determines the final shape of thetorus.29Chapter 4Wicking Flow throughMicrochannels4.1 IntroductionAs mentioned in Chapter 1, capillary force is a key mechanism to moveflow inside narrow channels of porous media and microfluidic devices thattypically consist of complex geometric features. Dynamics of interface isdependent on the geometry of the conduit, and it is usually ignored in pre-vious studies of wicking flow, which is typically one-dimensional. There aretheoretical and numerical difficulties in computing wicking flows throughcomplex geometries. These include the lack of a good model for the mov-ing contact line, the need to capture dynamically a moving and deforminginterface, morphological singularities in coalescence and rupture of inter-faces, and the complex flow geometries in practically interesting problems.The first three are generic to simulation of interfacial flows. The last, ongeometry, is especially pertinent to flow in porous medium. The geomet-ric features of a pore include changes in the cross-sectional area betweenwide pore chambers and narrow pore throats, branching and intersection ofpores, and the appearance of sharp edges and corners on which a gas-liquidinterface can be pinned.The classic work on wicking flows is that of Lucas [49] and Washburn[85], who computed capillary rise in straight tubes. This solution is notablefor its simplicity. Dynamics of the meniscus is completely ignored. In itsplace, a static interfacial shape is assumed such that the interface merelysupplies a constant suction pressure. In addition, the flow in the tube istaken to be fully developed Poiseuille flow. Now the capillary rise can becomputed by balancing the capillary pressure against the viscous friction.Later work has sought to include inertia, dynamic contact angle and entryeffects [8, 23, 42, 70, 98]. More recently, tubes of non-circular cross sectionshave also been considered [35].Of more relevance to our work are generalizations of the Lucas-Washburn304.1. Introductionsolution to tubes and channels of gradually varying cross sections. Using thelubrication approximation, one-dimensional (1D) solutions have been ob-tained for sinusoidal tubes [68, 73], sinusoidal tubes with tortuosity [60] andtubes and channels with convergent, divergent and power-law cross sections[66]. Liou et al. [47] extended the previous solutions to 2D axisymmetricflows by using approximate velocity profiles. This allowed them to includeinertia as well as viscous stresses that vary with the cross-sectional area.The only numerical study of wicking flow inside tubes with convergent-divergent cross section is done by Erickson et al. [25] They have computedthe 2D axisymmetric flow inside tubes by finite elements, but excluded thehydrodynamics of the interface as in the analytical studies mentioned above.By assuming a spherical shape of the meniscus, they update its positionfrom the liquid volume flow rate. The capillary pressure is introduced bythe Young-Laplace equation along with a dynamic contact angle. Perhapsto minimize the disturbance to the meniscus and to make the sphericalshape a more accurate approximation, they have used very long tubes withexceedingly mild contractions and expansions, with a contraction/expansionangle of 0.5?. A surprising prediction is that if the total lengths of the widerand narrower portions of the tube are fixed, the time for the meniscus to passthrough the tube is independent of the number of the contraction/expansioncycles along the tube?s length.To summarize this brief review of the literature, previous studies haveignored the hydrodynamics of the meniscus. The motion of the contactline is unaccounted for, and the meniscus shape is always prescribed to bespherical. Little can be found in the literature that deals with the detailedmorphological changes of the meniscus during its motion through complexgeometries. As such, their applicability to two-phase transport in porousmedia is quite limited. For one, pinning of the interfaces on sharp corners ofthe pore is responsible for pore blockage and, if external pressure is applied,eventual capillary breakthrough [22, 48]. Competition among interconnectedpores will depend on the dynamics of the interfaces in complex geometries,including rupture and coalescence as the meniscus negotiates bifurcationsand junctions. Finally, the wettability of porous medium is often manipu-lated to enhance two-phase transport [58]. The underlying mechanism hasto be sought from the hydrodynamics of the interface. In this context, afundamental understanding of interfacial dynamics during wicking throughcomplex geometry is essential. In view of the various difficulties mentionedabove, however, a rigorous study of the dynamics of the meniscus based onhydrodynamic principles has yet to be done.This chapter presents an initial effort toward addressing these issues. We314.2. Problem setupFigure 4.1: Schematic of the flow geometry for wicking into a capillary tube withcontraction.examine the two quintessential geometric features of a porous medium: theareal changes between pore throats and chambers and the branching of flowconduits. Specifically, we simulate the wicking flow in axisymmetric tubeswith non-uniform cross sections and 2D planar channels that bifurcate intotwo branches. Using a Cahn-Hillard description of the contact line dynam-ics allows us to capture the moving interface and its morphological changesaccurately and naturally, including pinning at sharp corners and otherwisesingular interfacial breakup at bifurcations. We show that the meniscus un-dergoes complex deformations through contractions and expansions, withcontact line pinning at protruding corners and turning of the interface atconcave corners. Capillary competition between bifurcating channels maysuppress wicking in the wider branch in favor of the narrower one. Manip-ulating the wettability in the branches can even produce flow reversal.4.2 Problem setupWicking is significant in small capillary tubes and pores, and the resultingReynolds and Bond numbers are typically much below unity. Therefore,we will neglect inertia and gravity throughout this work, and highlight theroles of capillarity and viscosity. The hydrodynamics is governed by thecontinuity equation and a modified Stokes equation.The computational domain is illustrated in Fig. 4.1. Natural boundaryconditions are employed at the inlet and the outlet. In addition, the pressureis set to be equal between the inlet and the outlet such that the motion ofthe liquid column is driven entirely by wicking, i.e. by the capillary pressure324.2. Problem setupgenerated by the meniscus.The wicking flow consists in a column of hydrophilic fluid of viscosity?1displacing a hydrophobic one of viscosity ?2(Fig. 4.1). The geometricalparameters include the contraction or expansion angle ?, the total lengthHt, length of the upstream section Hu, and the larger and smaller tube radiiR1and R2. In Fig. 4.1, Hb denotes the position of the center of the meniscus,hereafter called its base point. Hw marks the position of the contact lineon the wall, hereafter called the wall point of the meniscus. Initially, thereis no flow and the liquid column is at Hw = H0. Later, both Hw and Hbvary in time as the wicking proceeds. In this dynamic process, the interfaceshape is determined by the viscous and capillary forces, and is in general notspherical. However, it will prove convenient to use an effective curvature ?,defined for a spherical surface, in discussing the evolution of the interface.From the height of the meniscus ?m = Hw ? Hb and the local tube radiusR, we can calculate the radius of the spherical surface that passes throughthe wall and base points of the meniscus: ?s = R2+?2m2?m, from which we candefine:? = 2?mR2 + ?2m. (4.1)Note that ? is an overall indication of the meniscus curvature, and does notreflect the local deformation of the interface. It varies along the axis in anexpansion or contraction as R does. In presenting results in dimensionlessform, we scale length by R1, curvature by R?11, velocity by ?/?1and time by?1R1/?. Throughout this study we have set R2/R1= 0.5 and Ht/R1 = 20except in Fig. 4.10 where Ht/R1 = 21.The physical parameters of the problem can be combined into two di-mensionless parameters: the static contact angle ? and the viscosity ratiom = ?2/?1. Unless noted otherwise, m is set to 0.02 to represent the vis-cosity ratio between air and water at room temperature.4.2.1 Choice of Cahn-Hilliard parametersAs discussed in Chapter 2, the contact line velocity is very small for theconsidered wicking flow so we assume fast wall relaxation and set ? equalto zero. Therefore, the diffuse-interface model introduces two additionallengths [96]: the interfacial thickness  and a diffusion length ld = ?1/2(?1?2)1/4.They produce two dimensionless parameters: the Cahn number Cn = /R1and S = ld/R1.  or the Cahn number Cn should be sufficiently small sothat the numerical results no longer depend on it; this is known as the sharpinterface limit. To ensure that the sharp-interface limit for moving contact334.2. Problem setup(a) tHw0 200 400 600151617181920S=0.04, Cn=0.02S=0.04, Cn=0.01(b) t?0 200 400 6000.450.460.470.480.490.5S=0.04, Cn=0.02S=0.04, Cn=0.01Figure 4.2: Sharp-interface limit for computing capillary rise. (a) Contact linemotion indicated by the rise of Hw in time for two Cahn numbers Cn = 0.01 and0.02. (b) Variation of the effective meniscus curvature ? with time for the sametwo Cn values. S = 0.04, ? = 60?, the tube radius R = 1, total length Ht = 20 andthe initial column height H0= 15.lines is achieved, we use the criterion of Yue et al. [96]: (Cn ? 4S). In ourwicking problem, the criterion turns out to be more stringent than that ofYue et al. [96]. We tested a range of Cn values for S = 0.04, and found theresults to be essentially independent of Cn once it is below 0.02. Fig. 4.2shows that the contact line motion and the meniscus shape agree closely be-tween Cn = 0.02 and Cn = 0.01. Thus, the sharp-interface limit is achievedby using Cn = 0.02 in this case.As discussed in Chapter 2, the diffusion length ld or the parameter S,should be chosen to match a single experimental measurement. In thissection, we examine these issues in the simple geometry of imbibition anddrainage in a straight capillary tube. Our aim is twofold: to select a suitablevalue for S, and to validate the numerical results for a straight circular tubeagainst the Lucas-Washburn formula [49, 85].In the capillary rise problem, the Lucas-Washburn formula [49, 85]H(t) =?H20+?R cos ?(t? t0)2? (4.2)is widely accepted as an accurate representation of the interfacial movement,as long as the liquid column is long enough such that the ?end effect? isnegligible and the flow can be approximated by the Poiseuille flow. Here,R is the tube radius, ? is the viscosity of the liquid and that of the gas is344.2. Problem setup(a) (b)Figure 4.3: Comparison between the diffuse-interface simulation and the analyt-ical Lucas-Washburn formula at different S values. (a) Imbibition with ? = 60?,Cn=0.01, m = 0.02 and H0= 15. (b) Drainage with ? = 120?, Cn=0.01, andH0= 19. Now the less viscous component is wetting, and the non-wetting-to-wetting viscosity ratio m = 50.neglected.Fig. 4.3 compares our diffuse-interface calculation of capillary imbibitionand drainage at several S values with the Lucas-Washburn formula. Afew observations can be made. First, in both imbibition and drainage, thenumerical result approaches the analytical formula as S increases. While theanalytical solution neglects the dynamics at the meniscus and the contactline completely, the Cahn-Hilliard model includes a friction at the contactline in terms of an additional dissipation [97]. As the diffusion length or Sincreases, the effective slippage at the contact line increases, thus reducingthe influence of this friction. For S = 0.04, the effective curvature of theinterface is 5% lower than that expected of a spherical surface. This isdue to the flow effects at the meniscus that the Lucas-Washburn formuladisregards.Second, the contact line speed is insensitive to S. With a tenfold changein S, the contact line speed changes by some 6%. This forms an interestingcontrast to the situation studied by Yue et al. [96], where in shear flowsthe meniscus shape or interface inclination is highly sensitive to S. Thiscan be rationalized by the fact that the contact line speed is determined byequating the viscous dissipation to the surface energy gained by wetting ordewetting. Thus, insofar as most of the dissipation occurs in the bulk of thecolumn, the effect of S is mild. In shear flows, in contrast, the contact line354.3. Wicking in a tube with contraction or expansionspeed is prescribed, and the amount of Cahn-Hilliard diffusion affects theshape of the interface greatly.Third, the same S produces larger deviation from the Lucas-Washburnformula for imbibition than for drainage. This reflects the fact that theviscosity of the displacing and the displaced components contributes to thecontact line motion differently. But this asymmetry is not reflected by thedefinition of S used here. Finally, there is an upper limit to reasonable Svalues. Using too large a diffusion length ld exaggerates the area that isdirectly affected by the contact line. Our numerical experiment shows thatfor S = 0.15, for example, the overall features of the flow are distortedby the interfacial diffusion and the solution becomes very inaccurate. Intypical flow situations, the slip length is orders of magnitude smaller thanthe macroscopic length scale [20, 94].To sum up this section, we have demonstrated how the sharp-interfacelimit can be achieved by using a small enough Cn and how S can be selectedby comparing with the Lucas-Washburn formula. Most of the results to bepresented are for S = 0.04 and Cn = 0.01. For the wicking through multiplecontraction-expansion combinations (Fig. 4.10), we have used Cn = 0.02.To better resolve contact line pinning and turning at corners (e.g., Figs. 4.4and 4.7), we have used a smaller Cn = Wicking in a tube with contraction orexpansionWe consider the wicking flow of a liquid column into an axisymmetric tube,with a contraction as shown in Fig. 4.1 or with an expansion. The goal isto elucidate the detailed hydrodynamics of the moving interface, in partic-ular how the contact line negotiates concave and convex corners. Also ofinterest is the passage time as a function of the flow geometry, with a singlecontraction-expansion combination or multiple cycles of it.4.3.1 ContractionFig. 4.4 illustrates the wicking of a liquid through a 2:1 contraction at con-traction angle ? = 45?. The wall is hydrophilic to the liquid, with a wettingangle ? = 60?. The evolution of the interface is punctuated by severalcritical points marked on the Hb ? Hw and ? ? Hw curves as well as bythe insets. In the first stage of the process (Hw < Hu), the meniscus moveswith a constant shape within the wide tube before it reaches the contraction.364.3. Wicking in a tube with contraction or expansion(a) HwHb9.8 10 10.2 10.4 10.6 Hw?9.8 10 10.2 10.4 10.6 10.800.511.52abcdFigure 4.4: Meniscus movement through a contraction with ? = 45? representedby (a) the variation of the base point with the wall point, and (b) the effectivecurvature defined in Eq. (4.1). ? = 60?, Cn = 0.005 and Hu = 10. In (a) the insetscorrespond to the four points marked by squares on the curve. In (b) the dashedline indicates the curvature expected of a quasi-static spherical meniscus.The base point and wall point advance at equal speed and the trajectoryin Fig. 4.4(a) is a straight line with slope 1. The meniscus is not spherical,however. Viscous forces distort it so that the constant effective curvature? is 6% below that expected of a spherical meniscus at equilibrium in theupstream portion of the tube.As the contact line reaches the corner at the beginning of the contraction,marked by point a in the plot, a new behavior sets in. First, the contact linequickly moves past the concave corner. Once it is on the inclined wall of thecontraction, the interface must rotate by ? to maintain the same contactangle ? with the wall. This rotation first occurs locally at the contact line,elevating the local curvature. Then the interfacial distortion propagatestoward the center by interfacial tension, causing the central portion of themeniscus to pull back upstream. This process is reflected in Fig. 4.4(a) bythe downturn of the trajectory and the sharp upturn of ? in Fig. 4.4(b). Theinterfacial adjustment is completed by point b, when the base point of theinterface is at a minimum. During this highly dynamic transition, ? falls farbelow what one would expect by assuming a quasi-static spherical meniscus.After point b, the base point moves forward again, and at a higher speedthan the wall point because R is shrinking and the interface is continuouslybecoming more curved. Hence the continued increase in ?. The next mile-post is when the contact line reaches the convex corner marking the end of374.3. Wicking in a tube with contraction or expansionthe contraction (point c). The contact line is pinned at the corner [37] whilethe base point continues to move forward. Thus the interface rotates towardthe downstream as if hinged at the corner, and in the mean time straightenswith a steep decline in ?. The pinning ends when the angle between theinterface and an extension of the downstream wall reaches the contact angle?, in accordance with Gibbs? pinning criterion, and the meniscus as a wholemoves into the narrow channel. This moment is marked by d in Fig. 4.4.The effective curvature ? settles into a steady value roughly 14% below thatfor a perfectly spherical meniscus at equilibrium.4.3.2 Regularization of corner singularityThe turning of the interface at the concave corner (point a) deserves acloser examination. The case illustrated in Fig. 4.4 has a relatively mildcontraction with ? < ?. Thus, as the contact line advances from the corneronto the ramp, the interface rotates locally by ? to form a tight curve, whichis subsequently smoothed out over the rest of the interface. Now imagine astronger contraction with ? > ?. If the interface rotates by ? at the corner,it would have to penetrate the wall upstream. Hence, a concave cornerwith a contraction angle larger than the wetting angle appears to present asingularity to the contact line.This singularity is not real, of course. It arises because the foregoingargument is made in the classical sharp-interface framework, with the in-terface being viewed as a mathematical surface of zero thickness. This is agood representation of real interfaces as long as the length scale of interestis much larger than the interfacial thickness. At a concave corner of suffi-ciently large ?, the turning of the interface entails intersection with the solidwall, which in reality would invoke physics on the molecular length scale.Little surprise that an apparent singularity should appear. In fact, a strictimplementation of the sharp-interface model would encounter difficulty evenfor a contact line moving on a flat substrate [61, 62].This is where the diffuse-interface model presents a distinct advantage.By preserving the reality that interfaces are diffuse mixing layers ratherthan discontinuities, the model circumvents the traps of singularity on flatsubstrate as well as at corners. On a flat substrate, Cahn-Hilliard diffusionallows a contact line to move and predicts a dynamic contact angle [94, 96].Inside a concave corner, diffusion allows the interface to turn a large ? in anatural manner as demonstrated in Fig. 4.5.We have to point out that the diffuse-interface model introduces a lo-cal length scale , the interfacial thickness. If the physical process being384.3. Wicking in a tube with contraction or expansionFigure 4.5: Gray-scale contours of ? depicting the interface traversing a concavecorner through Cahn-Hilliard diffusion. The light line indicates the contour of? = 0. The contraction angle ? = 75? is greater than the wetting angle ? = 60?.Hu = 10 and initially the meniscus is at H0 = 9.8.studied involves a length scale that shrinks indefinitely, as occurs here inFig. 4.5 and during interfacial pinch-off or rupture [29, 95], the finite- effectmanifests itself eventually, and is intrinsic to the diffuse-interface formalism.Therefore, the negotiation of the corner in Fig. 4.5 occurs more slowly withdecreasing . The question of choosing suitable Cahn-Hilliard parametershas been discussed elsewhere [94]. For the current problem, the diffuse-interface model regularizes the singularity at the corner and captures thequalitative features of the process, but cannot foretell what  value wouldpredict reality quantitatively.4.3.3 ExpansionWicking through an expansion, schematically depicted in Fig. 4.6, differsfrom wicking through a contraction in that the contact line first encounters aconvex corner, and then a concave one. The process is illustrated in Fig. 4.7.When the meniscus is entirely inside the narrower channel upstream, Hw andHb advance with the same speed. As the contact line reaches the corner atthe start of the expansion (point a), it is pinned temporarily according toGibbs? pinning criterion [37]. Meanwhile the base point moves forward veryquickly until point b, when the interface reaches an angle of ? + ? = 85?with respect to the upstream wall. It depins from the corner, and the entiremeniscus advances through the expansion. This corresponds to the segmentbetween points b and c. At point c, the meniscus reaches the end of theexpansion with the contact line at the concave corner. As the wall rotatescounterclockwise by ? at this corner, so must the interface before it couldmarch downstream onto the straight portion of the tube. This causes a largelocal curvature of the interface, which propagates toward the center, causing394.3. Wicking in a tube with contraction or expansionFigure 4.6: Schematic of an expansion illustrating the pinning criterion. ?b = ?+?is the breakthrough angle, and ?m = 90? is the maximum angle that the interfacemay reach at the corner.HwHb4 4.5 5 5.5 6 6.5 744.555.566.5abcdFigure 4.7: Wicking through an expansion with ? = 25?, ? = 60?, Cn = 0.005and Hu = 5. The insets correspond to the four points a?d on the curve.the base point of the meniscus to retreat, as illustrated by the decline of Hbbeyond point c. Once this interfacial adjustment is completed at point d,the entire meniscus moves down the wider straight channel, again with thebase point and wall point advancing at the same speed.Naturally one contrasts the above process with wicking through a con-traction (Fig. 4.4). The behavior at the concave corner at the end of theexpansion, between points c and d in Fig. 4.7, is essentially the same as404.3. Wicking in a tube with contraction or expansionFigure 4.8: Permanent pinning of the interface at the entrance to an expansionwith ? = 30?. ? = 60?, Cn = 0.005, Hu = 5 and H0 = 4.9.appears at the start of the contraction, from a to b in Fig. 4.4. If the ex-pansion is too abrupt, with ? > ?, the corner would present a singularity toa sharp-interface model but not to our Cahn-Hilliard model. On the otherhand, the convex corner at the start of the expansion, point a in Fig. 4.7,differs fundamentally from its counterpart, point c in Fig. 4.4; here it hasthe potential of permanently pinning the interface. This would happen ifthe expansion is abrupt enough such that the required breakthrough angle?b = ? + ? is beyond the maximum achievable ?m = 90?:?+ ? ? 90?. (4.3)Such a situation is illustrated by the snapshots in Fig. 4.8 for ? = 60? and? = 30?. After the contact line gets pinned at t = 6.6, the interfacial tensionacts to move the rest of the interface forward so as to minimize the interfacialarea. This continues till t = 101.6, when the meniscus becomes a flat surfaceand can move no further. At this point, the angle between the interface andthe upstream wall is ?m = 90?, barely equal to the breakthrough angle?b. The contact line cannot depin and the flow is arrested permanently.Thus, the geometric constraint of Eq. (4.3), based on the Gibbs pinningcondition [37], specifies a degree of expansion beyond which a hydrophilicfluid cannot enter. This may be contrasted with the convex corner at the endof the contraction (point d of Fig. 4.4). As long as the fluid is hydrophilic(? < 90?), the contact line always depins before the meniscus becomes flatat 90? angle with the downstream wall.4.3.4 Penetration timeThe speed of wicking and the time required to penetrate a given depth areof practical significance in various applications [10, 54, 60], and have been414.3. Wicking in a tube with contraction or expansion(a)tHw0 100 200 300 400 500911131517?=15?=25acb(b)Figure 4.9: (a) Comparison of wicking speed through 2:1 contractions at twocontraction angles ? = 15? and 25?. The inset illustrates the flow geometry. Hu =10, H0= 9, ? = 60?. (b) Similar comparison for 1:2 expansions at ? = 15? and25?. Hu = 5, H0 = 4, ? = 60?.studied by a few groups [25, 73]. In this subsection, we will examine how thespeed of wicking through contractions and expansions is affected by the flowgeometry. Consider a tube of total length Ht. The two radii R1 and R2 areprescribed, as is the upstream length Hu. The rest of the length consists of acontraction or expansion and possibly a straight downstream segment. Thequestion is what contraction or expansion angle gives the fastest wickingthrough the total length. We have tested a range of contraction/expansionangles and wetting angles. A clear trend emerges, and is illustrated inFig. 4.9 by comparing ? = 15? and 25? for ? = 60?.In Fig. 4.9(a), the two trajectories coincide prior to reaching the start ofthe contraction, point a. Afterwards the wicking accelerates faster throughthe sharper contraction at ? = 25?, evidently because of the faster increasein curvature and capillary pressure. But the sharper contraction is shorter,and the acceleration ends at point b, after which the meniscus enters thenarrow downstream segment and decelerates. In comparison, the mildercontraction sees a more gradual acceleration that lasts longer, till point c.Downstream of point c, wicking decelerates as well, but at a gentler rate thanin the sharper contraction. This is because the sharper contraction incursmore viscous friction. As a result, the meniscus eventually overtakes thatin the sharper contraction, at Hw ? 14.8. Therefore, the question of whichgeometry gives faster wicking depends on the length of the downstreamsegment. If it is long enough, a gentler contraction wins. If Ht and Hu are424.3. Wicking in a tube with contraction or expansionprescribed, then there is an optimal ? that gives the shortest penetrationtime through the entire length. For example, for Ht = 13, Hu = 10, and? = 60?, we have tested ? values from 15? to 65?, and ? = 25? gives theshortest transit time.For expansion, the story is simpler (Fig. 4.9b). Wicking is slower inthe sharper expansion because the driving force, the capillary pressure, de-creases more steeply with the expanding tube radius. This effect is so strongthat the sharper expansion (from a to b) takes longer time to traverse thanits milder counterpart (from a to c) despite its shorter length. Note thesudden surge of Hw at b and c when the contact line rapidly traverses theconcave corner. Upon entering the downstream segment, wicking acceler-ates to more or less the same speed in both geometries. This speed willgradually decline in the downstream tube. Overall, the sharper expansionalways causes a longer penetration time. Besides, comparing Fig. 4.9(a) and(b), the expansion takes much longer time than the contraction of the samelength and same ?, by 15-fold for ? = 15?, and 53-fold for 25?. This impliesthat in a contraction-expansion combination, the latter takes up most of thepenetration time.The penetration or passage time tp, as it turns out, sheds unique lighton the validity of the quasi-static assumption widely used in the literature.By using the formula of Liou et al. [47], based on a quasi-static sphericalinterface, we have calculated tp through expansions at different angles. Theformula overpredicts tp by 2.5% for ? = 5?, and by 0.81% for ? = 25?. Withincreasing ?, wicking becomes slower, giving the interface more time to relaxtoward equilibrium. A different picture emerges for contractions. As ? in-creases from 5? to 25?, the underestimation of tp by the quasi-static methodincreases from 16% to 32%. Evidently, a faster moving interface deviatesmore from the spherical shape, and renders the quasi-static assumption lessaccurate. At large contraction angles, however, another factor comes intoplay. The strong radial flow tends to restore the meniscus toward spherical(cf. Fig. 4.4b).To better reflect the flow geometry in porous media, Erickson et al. [25]studied wicking through multiple contraction-expansion cycles. They cameto the surprising conclusion that as long as the total lengths of the wide seg-ments and narrow segments are each kept constant, the penetration time tpremains the same regardless of the number of contraction-expansion cycles.This implies that adding additional contraction and expansion pairs costsno delay in the wicking, something inconsistent with our observations inFig. 4.9. To probe this further, we compare wicking through three channelswith N = 1, 2 and 3 contraction-expansion cycles (Fig. 4.10). The total434.3. Wicking in a tube with contraction or expansiontHw0 1000 2000 3000101214161820N=1N=2N=3Figure 4.10: Wicking through multiple contraction-expansion cycles. For N = 1,2 and 3, ? = 26.6?, 45? and 56.3?. The wetting angle ? = 30?, Cn = 0.02. Thetotal length Ht = 21, and the meniscus starts at H0 = 10 at the beginning.lengths of the straight segments are the same among the three, 16 for thewider part and 3 for the narrower part. The sloping segments also add tothe same length of 2, and the contraction/expansion angle then increaseswith N . This geometric setup is modeled after Erickson et al. [25]According to Fig. 4.10, the total penetration time tp is not the sameamong the three; it increases by 17% from N = 1 to 2 and by another70% to N = 3. As expected, additional contraction-expansion pairs do costpenetration time, more so for larger N as ? increases. The discrepancyis mainly because Erickson et al. [25] used a much smaller ? (? 0.5?) andstraight sections much longer than the contractions and expansions. Thustraversing the contraction and expansion takes up only a small fraction ofthe total tp. Moreover, they ignored the local fluid dynamics at the meniscusand replaced it by a quasi-static spherical surface. To probe smaller ? inour model, we have computed gentler slopes with ? increasing from 1? atN = 1 to 15.6? at N = 16, with ? = 30?, R2/R1= 0.75, Hu = 10 andHt = 41. Compared with N = 1, tp increases by a mere 2.1% for N = 8and 11% for N = 16. Since Erickson et al. [25] only investigated N up to3, they would not have noticed much change in tp even if they had not usedthe quasi-static assumption.444.4. Wicking in Y-shaped branches: capillary competition(a) (b) (c)Figure 4.11: Schematic of a planar microchannel with a Y-shaped bifurcation,showing three stages of wicking: (a) the meniscus reaches the expansion; (b) themeniscus breaks into two at the bifurcation; (c) wicking continues in each branchunder suitable conditions.4.4 Wicking in Y-shaped branches: capillarycompetitionConnectivity between pores is an important attribute of porous media thathas not been considered in the above. When the meniscus reaches the bifur-cation where one pore branches into two, will it split into two and go throughboth branches, or will one branch dominate the other? What parametersdetermine the interfacial dynamics at and after the bifurcation? These arethe questions that we turn to in this section.Consider the wicking flow in the 2D planar geometry of Fig. 4.11. Thesame ambient pressure pa exists at the far-upstream of the root channeland the far downstream of both branches. When the interface reaches thebranching point, it breaks up into two smaller menisci (Fig. 4.11b), eachthen quickly adjusting to the size of the branches. A bifurcation into twoidentical branches is a trivial case; wicking proceeds in each branch withequal velocity. If the two branches differ in size, then there is a potentialfor capillary competition governed by the following three factors. (i) Thepressure behind each meniscus depends on its curvature and hence the sizeof the branch. Although the interface is generally non-spherical, we canroughly speak of the capillary pressure in the wide branch pw being higherthan that in the narrow one pn: pw > pn. The narrow channel engenders alower capillary pressure and is thus more conducive to wicking flow. (ii) At454.4. Wicking in Y-shaped branches: capillary competitionthe junction, we can roughly think of a pressure pj that is shared by bothbranches. The pressure drops pj ? pw and pj ? pn drive the flow in eachbranch (Fig. 4.11c). (iii) pj is determined by the viscous friction in the roottube, and continuously rises in time. This is because as wicking proceeds,the flow in one or both branches slows down and so does the flow in the roottube.Depending on whether pj is greater than pw and pn, we can differentiatetwo situations: wicking in both branches and wicking in one branch only.The former happens if pw and pn differ little, or if the pressure drop expendedin the root tube is small such that pj is initially high. The latter happensif the two branches are disparate in size, or if there is a long and thin roottube to yield a weak pj. In discussing these two regimes in the following twosubsections, we have found it convenient to fix D2= 0.5D1and L0= 4D1,and vary the width of the root tube D0relative to D1. In addition, thecontact angle is set at ? = 60?. Length will be scaled by D1and time by?1D1/?.4.4.1 Flow in both branchesWith a wide root tube, the viscous dissipation in it is small and it is likeconnecting the branches directly to a reservoir. In this simple situation,wicking occurs through both branches, though at different speeds dependingon their size (Fig. 4.12).When the meniscus reaches the end of the root tube (Fig. 4.11a), it facesan expansion at angle ?, and the discussion of pinning in ?4.3.3 applies. Inparticular, we require ? < 90??? such that the meniscus can depin from thecorner and proceed beyond this point. Throughout this section we have used? = 20?. When the interface reaches the point of bifurcation (Fig. 4.11b),it breaks into two smaller menisci, whose curvature, at this point, reflectsthe larger dimension of the junction. Thus they are not at equilibrium withthe smaller size of each branch. A short period of equilibration ensues, withthe wall points on the outside walls pulling back and those on the walls inthe middle surging ahead. This is why in Fig. 4.12 the curves appear tostart from a small positive L value at t = 0. In the inset, points a and bmark when the equilibration is completed in the narrow and wide branches,respectively. Note that the bifurcation angle ? determines the size of themeniscus in Fig. 4.11(b) and the equilibration process. But it has littleeffect on the subsequent wicking in each branch. Once the equilibration iscompleted, each meniscus is orientated symmetrically with respect to theaxis of its branch. The geometric setup is such that the pressure pj is higher464.4. Wicking in Y-shaped branches: capillary competitiontL0 50 100 150 200 250 3000246LwLnFigure 4.12: Wicking in both branches with a relatively wide root tube: D0=1.6. The origin of time is when the meniscus first touches the tip at the junction(cf. Fig. 4.11b). The inset shows that wicking is faster in the narrow branch initially(t < 11) but the wide branch winns for longer times.than both pn and pw, and wicking proceeds in both branches.Initially, the narrow tube enjoys faster wicking because the pressure droppj ?pn driving the flow is larger than its counterpart in the wide tube. Thislasts till t ? 11, marked by point c in the inset of Fig. 4.12. As the liquidcontinues to invade both branches, the viscous dissipation increases withthe column height and the flow speed declines. This effect is stronger forthe narrow branch since, as the flow approaches the Poiseuille flow, theviscous wall stress scales with the meniscus velocity divided by the channelwidth. Thus, the wider tube has faster wicking for later times, similar tothe prediction of the Lucas-Washburn equation.Once the flow starts in either branch, it does not stop in finite time.This is because pj must exceed the pressure pn or pw for there to be flowin the branch. In time pj and hence the pressure drop only increase as theflow and pressure drop in the root tube declines. Thus, the flow continuesin both branches, and gradually slows down toward zero in time.4.4.2 Flow in one branchWith thinner root tubes, pj may initially fall below the capillary pressurepw in the wide branch such that wicking occurs only in the narrow branch.474.4. Wicking in Y-shaped branches: capillary competitiontL0 100 200 300 400 500 600-1012345t cLncLnLwFigure 4.13: Capillary competition between two branches with a relatively narrowroot tube, D0= 0.6. Wicking proceeds in the narrow branch but is suppressed inthe wide branch until tc = 194, marked by a dot on both curves. The origin of timeis when the meniscus first touches the tip at the junction (cf. Fig. 4.11b).This behavior is demonstrated forD0= 0.6 by the trajectories of the menisciin Fig. 4.13 and by the snapshots of the interface in Fig. 4.14. After theinterface splits into two at the bifurcation, they reorient with respect tothe axes of the branches and adjust their curvature to the local tube size(Fig. 4.14b). Afterwards wicking starts in the narrow branch, and the flowin the root tube entails a pressure drop. The junction pressure pj thusproduced turns out to be lower than the capillary pressure pw in the widebranch, and no wicking occurs there. In fact, the negative pressure pj ? pwcauses the interface to retreat until the contact line becomes pinned at theinner corner at t = 42. Fig. 4.14(c) depicts a moment soon afterwards withthe meniscus immobilized in the wide branch. But the arrest of flow in thewide branch is necessarily temporary. As the flow slows down in the narrowtube, pj rises continually, eventually surpassing pw to produce wicking flowin the wide branch as well. This is marked in Fig. 4.13 by tc = 194 when theliquid column in the narrow tube is at Lcn = 2.56. After that the situationbecomes qualitatively the same as in Fig. 4.12, and Fig. 4.14(d) shows asnapshot in this stage. Eventually wicking slows down toward zero in bothbranches.The onset of wicking in the wide branch, indicated by tc or Lcn, is ofpractical interest. For instance, in a porous medium of finite thickness, thecritical value Lcn for the small pores will determine whether the bigger pores484.4. Wicking in Y-shaped branches: capillary competition(a) (b)(c) (d)Figure 4.14: Evolution of the interfacial morphology for the simulation depictedin Fig. 4.13. (a) The meniscus touches the salient corner at t = 0. (b) Themeniscus relaxes toward the equilibrium curvature inside each branch. (c) Aftera brief retraction, the meniscus is immobilized in the wide branch. (d) After therestarting of flow in the wide branch, the menisci advance in both branches.will contribute to liquid transport at all. If the liquid traverses the entirelength of the smaller pores before wicking even starts in the bigger ones, thelatter will be dead ends, which have been observed in experiments [33] andconsidered a major hinderance to water transport through the GDM of fuelcells [48]. Note that all the ideas and qualitative arguments discussed so farin this section apply as well to 3D flows in real porous media.In the spirit of the Lucas-Washburn analysis (Eq. 4.2), we can esti-mate the onset of wicking in the wide branch by neglecting dynamics at themenisci and assuming fully developed Poiseuille flow in the root and narrowtubes. Let us denote the instantaneous average velocity in the narrow tube494.4. Wicking in Y-shaped branches: capillary competitionby V2and that in the root tube by V0. Then the junction pressure pj canbe estimated either from the force balance on the liquid in the root tube orthat in the narrow tube:pj = pa ?12?L0V0D20= pa ?2? cos ?D2?12?LnV2D22. (4.4)The critical condition for wicking in the wide tube is pj being equal to thecapillary pressure behind the meniscus in the wide tube:pj = pw = pa ?2? cos ?D1. (4.5)In addition, volume conservation requires V0D0= V2D2. Eliminating V0and V2from the above gives the following critical condition on the liquidcolumn Lcn:Lcn = L0(D1D2? 1)(D2D0)3. (4.6)Recall our previous argument that wicking in the wide tube dependson the viscous friction in the root tube and the dissimilarity between thetwo branches. It is no surprise that Lcn turns out to depend on the lengthand diameter of the root tube as well as the size difference between thetwo branches. For the conditions in Fig. 4.13, Eq. (4.6) predicts Lcn = 2.3,reasonably close to the numerical result of 2.56. Numerical experimentationwith narrower D0values has confirmed further delays in the wide branch inagreement with Eq. (4.6). Finally, we note that the above calculation canbe easily generalized to 3D circular tubes, and the formula has the exponenton (D2/D0) changed from 3 to Flow reversal due to spatially inhomogeneoushydrophilicityInsofar as the Young-Laplace equation gives a capillary pressure in the formof ? cos ?/D, varying the contact angle ? in a branch is in a way tantamountto varying the tube size D. Thus, capillary competition between branchescan be controlled by varying ? as well as D. Suppose that in Fig. 4.11, wemake the downstream portion of the wide branch more hydrophilic, witha smaller contact angle. Then a flow reversal may occur in the narrowchannel, as illustrated in Fig. 4.15.In this geometry, ? = 60? throughout the Y-branch except for the down-stream portion of the wide branch starting from Lw = 1.25 that features504.4. Wicking in Y-shaped branches: capillary competitiontL0 400 800 120000.511.522.533.54LnLwFigure 4.15: Flow reversal in the narrow branch when the meniscus in the widebranch moves onto a more hydrophilic portion with ? = 20? at Lw = 1.25. Else-where ? = 60?. D0= 0.6, L0= 4, D2= 0.9.a smaller ? = 20?. The geometric and physical parameters of the setupare such that wicking occurs initially only in the narrow branch, and startslater in the wide channel around t = 110. When the meniscus encountersthe more hydrophilic portion in the wide branch (t = 942), the wickingsuddenly accelerates, causing a flow reversal in the narrow tube. This isbecause the elevated flow rate in the root tube depresses the pressure at thejunction so much that it falls below the capillary suction pressure pn in thenarrow tube.Depending on the physical and geometric parameters, the liquid columnmay retreat entirely from the narrow tube, with the interface pinned at thecorner of the bifurcation, or reverse its course again before that. Thereafter,the situation becomes similar to Fig. 4.13 or Fig. 4.12. Based on the Young-Laplace equation, one may view the wicking in the more hydrophilic portionof the wide branch as occurring in a tube with ? = 60? but a smaller effectivewidth De = D1 cos 60?/ cos 20? = 0.53D1, which is narrower than D2. (Theviscous friction will be different, of course.) Thus, the wicking continueswith dwindling speed in the wide branch until the junction pressure hasagain risen above the capillary pressure in the narrow tube to restart wickingthere.Such flow reversal has been observed experimentally. Litster et al. [48]reported that in a model GDM for fuel cells, a sudden acceleration in oneflow path, due to breakthrough from the GDM into open space, causes the514.5. Summary and conclusionsliquid to retreat in a neighboring connected path. The underlying principle issimple and robust, and suggests how surface properties can be manipulatedto control the flow pattern in porous media. Indeed, the GDM of fuel cellsis often surface-treated in a spatially inhomogeneous way to enhance watertransport [46]. In addition, a more hydrophilic micro-porous layer with finerpores is often attached to the GDM to create a jump in wettability along theflow direction [78]. Another potential application for capillary competitionand flow reversal is as a precise switching mechanism in microfluidic devices[7, 16, 44]. By careful choice of the root and branch sizes it is possible todesign a flow loop in which different branches are impregnated by liquidat precise moments. The mechanism of capillary competition works formultiple branches as well, and one may design microfluidic manifolds usingthe same principle.4.5 Summary and conclusionsThis work aims for a detailed and rational understanding of two-phase trans-port through micropores in porous media. Using finite-element computa-tions, we capture the evolving morphology of the interfaces in geometriesthat retain the salient features of real pores, including expansion, contrac-tion and branching. From a fundamental viewpoint, the most importantfindings are the following:(a) The meniscus undergoes complex deformation during transit throughmicropores, governed by the dynamic balance among fluid-solid andgas-liquid interfacial tensions and viscous friction. Such flow effectstend to distort the meniscus away from a spherical shape.(b) The dynamics of the contact line plays a central role. It pins at pro-truding corners, potentially barring wicking into expansions with toosteep a slope. The contact line negotiates inner corners thanks to thediffuseness of the interface.(c) Capillary competition between connected branches depends on thecapillary pressure due to meniscus curvature inside each, and in turnon the size of the branches and surface wettability. Under suitableconditions, wicking can be arrested in wider branches in favor of anarrower one, and the flow may even reverse course when wickingaccelerates in a neighboring path.524.5. Summary and conclusionsWe have hinted at the relevance of these insights to technological ap-plications, e.g., in proton-membrane exchange fuel cells. Against this back-ground, however, the work reported here must be seen as a preliminary step.Real 3D flow through porous media includes many complicating factors thathave not been accounted for, including 3D connectivity, pore size distribu-tion and tortuosity of the flow path. Nevertheless, this serves as a startingpoint for an approach to two-phase flow in porous media that is more ratio-nal and accurate than the traditional one centered on an empirical relativepermeability.53Chapter 5Auto-ejection of LiquidDrops from Capillary Tubes5.1 IntroductionDroplet production is a fluid dynamical process of considerable importancein engineering applications. The rapid development of microfluidic tech-nology has given new impetus to the study of controlled drop productionin miniaturized devices [31]. A common method for drop production is topump liquid through a tube such that a jet issues from the end, and breaksup due to capillary instability. In microfluidics, this is typically realizedby flow focusing [4, 77], and two regimes, jetting and dripping, have beenidentified [2, 84, 99]. Jet breakup can be actively promoted and controlledby a pressure pulse, as in drop-on-demand devices [91]. In these schemes ofdrop formation, the jet is always fed by an externally controlled flow rate.As mentioned in Chapter 1, Wollman and coworkers have demonstrateda novel method of drop formation that relies on wicking in a capillary tube[86, 87, 88]. Two interesting questions can be asked about this process:what the critical condition is for ejecting one or more drops, and how geo-metric parameters of the problem affect the ejection. The ejection process isgoverned by inertia as well as capillarity, much like for Worthington jets [30]and cavity jets [5]. Regarding the first question, it seems reasonable to arguethat auto-ejection occurs when the upward momentum of the liquid columnovercomes capillary restriction of the liquid surface. As will be shown later,viscous friction is negligible under typical experimental conditions. How-ever, it is difficult to quantify this idea in terms of a Weber number. This isbecause both the liquid momentum and the capillary restriction vary in timeas complex functions of several factors, including the dynamic contact angle,the shape of the nozzle and contact line pinning. In particular, auto-ejectionhas never been recorded at the end of a straight tube; the converging nozzleseems to be necessary [70, 88].To analyze this intricate process, it seems appropriate to divide it into545.1. Introductiontwo stages, the acceleration of the meniscus inside the tube, including thenozzle at the end, and the protrusion and possible breakup of the jet outsidethe nozzle. In the following, we will briefly summarize the current state ofknowledge on each phenomenon.Capillary rise inside straight tubes has been extensively studied before;see for example [72]. In the absence of gravity, the dynamics is mainlygoverned by the interplay among capillary, viscous and inertial forces. Atthe initial stage of the rise, viscous forces are negligible and the balancebetween capillary and inertial forces yields a constant rise velocity [64]:vci =(2? cos ?d?R)12, (5.1)where R is the tube radius, ? is the liquid density, ? is the surface tensionand ?d is the dynamic contact angle. This is known as the capillary-inertialvelocity. As the imbibition proceeds, the liquid column increases in lengthand mass. Viscous friction becomes important and the meniscus velocitystarts to decline. Eventually inertia becomes unimportant and the dynamicsenters the Lucas-Washburn regime where capillary pressure balances theviscous friction [49, 85]. Denoting the liquid viscosity by ?, we can writethe velocity of rise asvLW =R? cos ?d4?H , (5.2)which decreases with the length of the liquid column H. The above are twolimiting behaviors for short and long times. In the auto-ejection process,however, it is not clear a priori if the meniscus velocity follows either equa-tion. What is more, these simple models disregard the contact line dynamics.At high velocities, the dynamic contact angle ?d may deviate considerablyfrom the static one ? [11, 38]. Thus, the capillary force driving the meniscuschanges with its velocity, adding another subtlety to the problem.As the nozzle is essential for auto-ejection, the meniscus accelerationinside the nozzle is a key aspect of the process. For inertialess flows, [51]have investigated the meniscus dynamics inside contractions, including thetransient turning of the interface, its evolving curvature as well as the overallacceleration of the liquid column. Auto-ejection requires a high incomingmomentum with a large inertia, and the meniscus dynamics inside the nozzleremains to be studied.In the second stage of auto-ejection, a jet emanates from the nozzle,and one or more droplets form through a capillary mechanism known as555.2. Problem setupend-pinching [74]. Essentially, capillary retraction at the tip produces a bul-bous end, whose neck then becomes susceptible to capillary pinch-off. End-pinching has been studied by linear instability analysis [45], one-dimensionallubrication model [3], experiments [15, 83] and numerical simulations [36,56, 67, 82]. These studies have assumed either zero incoming flow at the baseof the jet or a constant flow rate. Work of Gordillo and Gekle [34] appearsto be the only one that allows a transient incoming flow; a linearly decreas-ing incoming velocity is assumed for Worthington jets. The auto-ejectionproblem differs in that the jet is being fed by a time-dependent flow ratethat is governed by the morphology of the jet and the physical conditionsinside the tube and nozzle. Thus, spatial and temporal variations of theliquid velocity determine the fate of the jet and the number and size of anydroplets that may form. Prior studies have indicated additional geometriccomplications related to the shape and wettability of the lip of the nozzle[2]. How the interface may depin from the inner edge of the lip and movealong its width turns out to have a strong influence on drop pinch-off.The review of prior work suggests the criterion for auto-ejection to bethe most prominent question. To answer this question, one must study themeniscus dynamics in the tube and the nozzle, as well as the jet behav-ior outside. In particular, the criterion should predict how auto-ejectiondepends on geometric factors: tube length, contraction angle, and even thewidth of the lip at the exit of the nozzle. We undertake such an investigationusing numerical simulations that captures detailed features of the contactline dynamics.5.2 Problem setupThe axisymmetric computational domain consists of a capillary tube con-nected to a liquid reservoir at the bottom and ambient air at the top(Fig. 5.1). In most of the simulations the tube has a contracting nozzleat its upper end. The contraction angle is ? and the radius shrinks from thetube radius R to Rn at the end of the nozzle. The total length of the tube,including the nozzle, is L. Thus, the flow geometry is completely specifiedby three dimensionless quantities: the contraction angle ?, the contractionratio C = Rn/R and the aspect ratio L/R. Initially the air-liquid interfaceis assumed flat at a small distance L0inside the tube. For the most part, L0represents the capillary climb under normal gravity before the drop-towerexperiment commences [88]. There is also a numerical incentive for placingthe interface inside the tube to avoid complications at the corner.565.2. Problem setup  Reservoirinitial position   of interfaceAmbient AirL0RRn?rLz4R3R4RnWLFigure 5.1: Schematic of meridian plane of the axisymmetric computational do-main.The liquid and air reservoirs are sufficiently large that their boundarieshave no effect on dynamics of the meniscus, liquid jet and drops. Basedon numerical experiments, we have chosen the liquid reservoir to be 3R inradius and 4R in height. On its bottom and side walls, we impose zero nor-mal stress and zero tangential velocity as boundary conditions. Its top wallis taken to have zero shear stress and zero normal velocity. This boundarycondition avoids the computational cost of tracking the slight deformationof the liquid-air interface outside the tube. [72] have shown that this sim-plification has little effect on the meniscus motion, and we have reached thesame conclusion by benchmarking our simulation of capillary rise againstexperiments. The air reservoir on top is 4Rn in radius, and its height rangesfrom 12Rn to 30Rn depending on the length of the ligament in differentsimulations. Zero stress boundary conditions are used on the top, bottomand side of the air reservoir. On the sloping walls of the nozzle, no-slipconditions are imposed. The upper surface of the nozzle (or the ?lip?) isa horizontal ring of width Wl. For most of the simulations, this surface isassigned a contact angle ?l = 180? to ensure that the contact line remainspinned at the inner corner of the lip. Smaller ?l values are used in ?5.4.4 to575.3. Physical model and numerical algorithmexplore depinning of the interface from the sharp corner.In addition to the geometric ratios, the problem is characterized by fourdimensionless groups based on material properties: the liquid-air densityratio ?/?a and viscosity ratio ?/?a, the Ohnesorge number Oh = ?/??R?,and the static contact angle ? inside the tube and nozzle. On the innersurface of the tube, we impose the no-slip condition, and model the motionof the three-phase contact line by Cahn-Hilliard diffusion to be discussedbelow. Gravity is neglected in all presented results except in Fig. 5.10.This is because most of the experimental data have been collected undermicrogravity, and gravity tends to inhibit auto-ejection. We will fix theseparameters: ? = 0? (perfect wetting), ?/?a = 200 and ?/?a = 100. Incomparison with the silicone oils used in the experiments [88], the densityratio is too low but the viscosity ratio is within the range of experimentalvalues. In view of the numerical difficulties in computing larger densityratios, we are satisfied that the air has little influence on the liquid jet anddrops [28]. We will vary the three geometric ratios C, ? and L/R along withthe Ohnesorge number Oh. We will use R as the characteristic length, thecapillary-inertial time tci =??R3/? as the characteristic time, and R/tcias the characteristic velocity, and present the results in dimensionless form.5.3 Physical model and numerical algorithmFrom a computational viewpoint, the auto-ejection process is difficult tosimulate as the interface moves, deforms and eventually breaks up, and theprocess features a prominent role for the moving contact line.As discussed in Chapter 2, our diffuse-interface model leaves us withthree new model parameters, say , ? and ?. We follow the procedurerecommended by [93, 94] to choose their values. For smallest value of which is computationally affordable, ? is chosen to ensure that the sharp-interface limit is achieved. Then value for the wall relaxation parameter ?is determined by fitting an experimental datum.To implement this procedure, we make the parameters dimensionless us-ing a characteristic length lc: Cn = /lc, S = ld/lc, and ? = 1/(??lc). Cnis commonly known as the Cahn number [100]. One needs to be careful inchoosing lc. Accurate simulation using the diffuse-interface model requiresthat the interfacial thickness  and the diffusion length ld both be muchsmaller than the global length scale [93]. As the meniscus advances throughthe nozzle, the effective global length scale is shrinking. We find it neces-sary to reduce  and ld accordingly to maintain accuracy of the simulation.585.3. Physical model and numerical algorithm(a) tHc0 2 4 6 80246810 SimulationExperimentpinch-off(b) HcVc2 4 6 8 10024681012SimulationExperimentpinch-offFigure 5.2: With a wall relaxation parameter ? = 0.4, the simulation approximatesexperimental results closely in terms of (a) the position of the center of the meniscus,and (b) the centerline velocity of the meniscus. The arrows indicate the momentwhen the contact line reaches the start of the nozzle, and the curves end when a droppinches off, indicated by a filled square. The geometric and physical parametersmatch the experiment of [86]: Oh = 0.011, L = 5.98, C = 0.493, ? = 23.8? and? = 0?. In addition, S = 8? 10?3 and Cn = 0.01.Therefore, when the contact line is in the straight portion of the tube, wetake lc = R. When it is in the nozzle, we take lc to be the local radius ofthe nozzle at the contact line. After the contact line reaches the lip, we fixlc = Rn. Thus, with fixed values of Cn, S and ?, the microscopic lengths and ld shrink inside the nozzle as required.We determine the model parameters by the experiment of [87] on cap-illary ejection. The geometric and material parameters are matched suchthat Oh = 0.011, L = 5.98, C = 0.493, ? = 23.8?. The static contact angleis 0? inside the tube and 40? degrees outside. The initial height of the liquidcolumn L0= 0.08L matches the experiment condition at the start. Follow-ing [94], we choose a small Cahn number Cn = 10?2 that is comfortablycomputable, and a corresponding S = 8?10?3. Then we found that ? = 0.4gives the closest agreement with the experimental results. This is illustratedin Fig. 5.2 in terms of the position and velocity of the center of the meniscus.A notable feature of this simulation is the evolution of the dynamic contactangle ?d. The wall energy relaxation in equation (2.12) allows ?d to deviatefrom ? [94]. Fig. 5.3 compares our computed ?d for capillary rise in a straighttube with two experimental correlations. In our computation, the meniscusrises with an essentially constant speed V , with which we define a capillarynumber Ca = ?V/?. The correlation of Ref. [11] is for solid strips drawn595.4. ResultsCa? d0 0.005 0.01 0.0150102030405060708090Jiang et al.Bracke et al.SimulationFigure 5.3: Comparison of the dynamic contact angle ?d in a straight tube betweenour numerical simulation and two experimental correlations due to [41] and [11].The model parameters are the same as in Fig. 5.2.into a pool of liquid, while that of Ref. [41] is based on the experimentsof Hoffman [38] on pushing non-polar liquids through glass capillary tubes.The numerical and experimental results all indicate an increase of ?d withCa, but the former exhibits a somewhat steeper slope than the experiments.One reason for the difference is that the Cahn-Hilliard model is phenomeno-logical, and the mechanism of wall energy relaxation cannot be expected tocapture quantitatively the dynamic contact angle. Moreover, in our simula-tions ?d is measured from the slope of the interface where it intersects thewall. In the experiments, it is estimated from fitting a circular arc to thecentral portion of the meniscus. This introduces some discrepancy as well.5.4 Results5.4.1 Meniscus dynamicsWe begin with an overview of the dynamics of the meniscus as it advancesthrough the straight portion of the tube and the contracting nozzle, andforms a jet outside the nozzle. For this purpose we select a typical set ofphysical and geometric parameters: Oh = 0.01, L = 5, C = 0.5, ? = 30?and ? = 0?. To describe the motion and deformation of the meniscus, wetrack the contact line velocity along the wall Vw and the velocity at the605.4. Resultscenter of the meniscus Vc in time.Fig. 5.4(a) plots Vw and Vc, as well as the average velocity across thenozzle exit Vn, as functions of time. Before the liquid meniscus arrives at thenozzle exit, Vn is computed from the velocity profile of air. From an initiallyflat shape (Fig. 5.1), the meniscus experiences an acceleration and adjust-ment phase at the start of the imbibition. The contact line immediatelymoves upward at a roughly constant speed, while the center of the meniscusoscillates several times before settling into a steady shape and speed of rise(point a in Fig. 5.4). This marks the start of the capillary-inertial regime.The meniscus velocity Vw = Vc = 1.09 agrees closely with the theoreticalresult vci = 1.07 (cf. equation (5.1)). This steady rise persists till point b,when the contact line arrives at the start of the nozzle. It is remarkable thatthe meniscus velocity stays roughly constant so far, showing little decreasedue to viscous dissipation. This can be rationalized by an estimation of theviscous effect in a straight tube. By balancing the capillary, viscous and in-ertial forces, Bosanquet [8] derived an analytical solution for the rise of themeniscus. For short times (Oh ? t  1), this solution predicts the followingvariation of the meniscus velocity along the tube:1VcdVcdHc= 2.4 Oh, (5.3)where Vc and Hc are dimensionless. Our numerical simulation verifies theproportionality to Oh, but with a milder slope of 1.8. For Oh = 0.01 andan axial distance of about 2.5 for the capillary-inertial regime in Fig. 5.4,therefore, viscous reduction of the meniscus velocity Vc is only about 5%. Infact, viscosity never plays an appreciable role throughout the entire process,and will be disregarded for the rest of the chapter. In the experiments, Ohis typically on the order of 0.01 [88], and viscosity is generally immaterial.Once the contact line reaches the nozzle, the interface immediately ro-tates at the contact line so as to adjust its orientation relative to the taperingwall of the nozzle. This rotation pushes the central portion of the meniscusbackward by capillarity, thus reducing the centerline velocity Vc to a min-imum at point c in Fig. 5.4. Afterwards, the meniscus accelerates rapidlyupward, mainly because of the upward momentum of the liquid column be-ing channeled through a narrowing conduit. Capillarity also contributes tothe acceleration since the meniscus is trailing the spherical shape at themoment, having been delayed by the rotation of the interface from b to c.This is illustrated in the snapshots of the interfaces in Fig. 5.4. This stagecontinues until point d, when the upward acceleration has moved the cen-tral portion of the meniscus ahead of the spherical surface dictated by the615.4. Results(a) tVelocity0 1 2 3 4 5 6 7024681012VcVnVwabcdfghe(b) t? d0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 55060708090100bcdae(c)gb fc d e h-1 0 10123456789aFigure 5.4: (a) Evolution of the contact line velocity Vw, meniscus center velocityVc, and average velocity at nozzle exit Vn in time. (b) Evolution of the dynamiccontact angle ?d. (c) Snapshots showing the position and shape of the meniscusat significant moments marked in the velocity plot. Oh = 0.01, ? = 0?, L = 5,C = 0.5 and ? = 30?.local dynamic contact angle. Thus capillary forces now pull the meniscusbackward, causing the reduction in Vc until point e, when the contact linereaches the lip of the exit. Note that with the contact line inside the nozzle,the contact line speed Vw is still measured by the axial position of the con-tact line, not by the distance traveled along the wall. Fig. 5.4(b) shows thatfrom point a till point e, the dynamic contact angle ?d closely tracks theevolution of the contact line speed Vw, in accordance with the observations625.4. Resultsin Fig. 5.3.At point e, the contact line becomes pinned at the sharp inner corner ofthe lip according to Gibb?s criterion [27, e.g.]. Constraining the upward flownear the nozzle wall results in a high pressure at the nozzle exit that thruststhe central portion of the meniscus out in the form of a jet (point f). As thejet is ejected and lengthens against surface tension, its tip velocity declinestoward point g, when capillary necking commences on the jet, eventuallyleading to a droplet pinching off at the tip (point h).It is interesting to contrast the behavior of the average velocity at thenozzle, Vn, with that of the meniscus velocity Vc. Note that because ofincompressibility, Vn gives the average liquid velocity in the tube (subject toa factor C2 due to area contraction) even before the meniscus reaches the lipof the nozzle. During the initial acceleration of the meniscus, prior to pointa, Vn monotonically increases to a constant level that corresponds to thecapillary-inertial regime. It starts to decline near point d, when capillaritystarts to oppose the upward motion of the liquid. The decline continuesmonotonically even as the jet rapidly issues from the nozzle. This can berationalized from how the interfacial tension, acting on the pinned contactline, continually depletes the upward momentum of the liquid column.The decline of Vn(t) in time after the jet formation (roughly from pointf onward) can be quantified from an energy argument. Consider a controlvolume that encloses the inside of the tube and the nozzle, as well as theliquid reservoir. The kinetic energy of the liquid inside, E, decreases becauseof the energy eflux at the nozzle as well as the pressure work there:dEdt = ?12?R2n?V 3n ? pn?R2nVn, (5.4)where pn is the liquid pressure at the exit of the nozzle. Note that thepressure in the reservoir equal that in the ambient, and has been put to zero,and that the incoming energy flux has been neglected as the velocity at theboundary the reservoir is much smaller than that inside the capillary tube. Eis the sum of the kinetic energy in the tube, the nozzle and the reservoir. Toestimate the fluid velocity inside the nozzle, we make a one-dimensional plugflow approximation by assuming that the liquid velocity is axial, and changesfrom Vn at the nozzle exit to C2Vn inside the capillary tube. Similarly, theflow in the reservoir is assumed to be radial and uniform on spherical surfacescentered at the entry of the tube, with a velocity that can be related to Vnthrough mass conservation [76]. Thus, E can be expressed in terms of Vn:635.4. ResultstVn5 5.5 6 6.5 7 7.500.511.522.533.54Eqn 5.6SimulationFigure 5.5: Temporal variation of the instantaneous velocity Vn(t) at the nozzleexit, starting from point f at tf = 4.85. L = 5, Oh = 0.01, ? = 30?, C = 0.5, ? = 0,Le = 1.68, We = 7.0.E ? ?2?R2nLeV 2n , where the effective lengthLe = RC(1? C)2tan? +(L+ 76R)C2 (5.5)is a purely geometric parameter. To estimate the exit pressure pn, wenote that the capillary pressure decreases from 2?/Rn to ?/Rn as the liquidinterface inflates from a semi-spherical shape to a cylinder with radius Rn.Taking pn = 2?/Rn, plugging E into equation (5.4) and integrating in time,we obtainVn(t) = u tan[?u(t? tf )2Le+ tan?1(Vfu)]= Vf1? uVftan[u(t?tf)2Le]1 + Vfu tan[u(t?tf)2Le] , (5.6)where u = ( 2??Rn)1/2, tf is the starting time for the integration, at point f ,when the average velocity across the nozzle exit is Vf = Vn(tf ).Fig. 5.5 compares Vn(t) predicted from the simple one-dimensional modeland the numerical solution. At the start, the model slightly underestimatesthe rate of deceleration. Toward the end, however, it overestimates it as thecapillary pressure at the exit falls below 2?/Rn and approaches ?/Rn. Insimulations the velocity profile is not a perfect plug flow. Such deviationsfrom plug flow, have compensated for the overestimation of the capillarypressure.645.4. ResultsBut over the course of the deceleration of the jet, the simple modelcaptures reasonably well the dynamics of the jet velocity. In particular,note how Le dictates the rate of decrease of Vn in time; a larger Le impliesa longer liquid column moving with a greater kinetic energy. Thus, thedeceleration will be slower, and a longer jet will likely be produced, in favorof auto-ejection. This point will be revisited in the next subsection. Finally,we have also confirmed that for the small Oh tested, viscous dissipationmakes a very small contribution to the energy balance of equation (5.4),consistent with previous arguments on the unimportance of viscosity in theprocess.5.4.2 Ejection criterionNaturally, we think of a Weber number to represent the idea that the upwardmomentum must overcome the capillary restriction. However, there are twodifficulties in constructing such a Weber number. First, there is no obviouscharacteristic velocity. The meniscus velocity is itself determined by thewicking inside the tube, and in turn by the contact angle and geometry(especially length) of the tube and nozzle. It also changes in time andin space. Wollman et al. [88] suggested an instantaneous Weber numberdefined using the liquid velocity at the exit of the nozzle when the meniscusfirst reaches that point. This corresponds to our point e in Fig. 5.4. Let ustake this point as the nominal start of the jet-formation process t? = 0, witht? = t ? te measuring the time from this point onward. Using the velocityVe = Vn(te) at this point, we can define an instantaneous Weber number:We = ?V2e Rn? . (5.7)Second, the instantaneous velocity Ve or We does not completely deter-mine the fate of the jet and breakup. Fig. 5.6(a) shows that We does notdelineate sharply the boundaries separating no-pinchoff and pinchoff, noramong different number of droplets produced. This inadequacy is not hardto appreciate. Roughly speaking, the eventual length of the jet is deter-mined by converting the kinetic energy of the entire liquid column at t? = 0into surface energy. Thus, the length of the liquid column should matteras well. Fig. 5.6(b) shows that under conditions that are otherwise iden-tical to Fig. 5.4, a shorter capillary tube (L = 1.5) produces a short jetand no breakup, whereas the longer tube (L = 5) of Fig. 5.4 does lead toauto-ejection. The Weber number We = 7 in both cases. Since we havepreviously introduced an effective tube length Le (equation 5.5), it seems655.4. Results(a) WeN0 5 10 15 20012(b) rz-1 -0.5 0 0.5 10.511.522.533.5t*=1.07t*=0.0t*=0.28t*=1.76Figure 5.6: (a) Number of drops produced as a function of We. The wide overlapsbetween different outcomes indicate that We does not provide an adequate crite-rion for auto-ejection of droplets. These data cover most of the parameter rangesstudied: 0.005 ? Oh ? 0.02, 0.25 ? C ? 1, 1 ? L ? 10 and 0 ? ? ? 40?. (b)A short tube (L = 1.5) fails to produce drop ejections under identical conditionsto Fig. 5.4, where a longer tube (L = 5) does produce ejection. The jet reachesmaximum length at t? = 1.07 and then retracts.natural to use it to account for the total amount of kinetic energy prior to jetformation. Fig. 5.7 plots the outcome of jet breakup against two parameters,We and Le, where Le has been made dimensionless by R. The overlaps inthe We plot (Fig. 5.6a) have now been sorted out by Le. This plot suggeststhe following criterion for predicting drop formation in auto-ejection:N =???0 if We < 3.4f(Le)1 if 3f(Le) < We < 5.5f(Le)2 if We > 5.5f(Le)(5.8)where f(Le) = 1+0.8/Le. For the range of parameters tested here, ejectionof 3 and more droplets has been observed mainly for large contraction angles,which produce a different flow regime to be considered in ?5.4.3. Thus wedo not include these cases here.A few remarks about this criterion seem in order. First, the criterion isgeneral as it encompasses almost the entire parameter ranges explored in oursimulations. The material and geometric parameters of the problem havebeen included through We and Le. The only exception is large contractionangles ? that induce new flow patterns. This is to be dealt with separately665.4. Results(a) WeL e0 5 10 15 2002468no dropone droptwo dropsWe=3.4 (1+0.8/Le)We=5.5 (1+0.8/Le)(b) WeL j /Rn0 5 10 15 200246810121416no dropone droptwo dropsFigure 5.7: (a) Criterion for self-ejection: number of droplets plotted as a func-tion of We and the effective length Le of equation (5.5). The three outcomes aredemarcated by We = 3.4(1+ 0.8/Le and We = 5.5(1+ 0.8/Le), shown as the solidand dashed curves, respectively. (b) The gray band, representing 5 ? Lj/Rn ? 7for the jet length of equation (5.9), indicates a rough threshold for auto-ejection.in ?5.4.3. Second, the auto-ejection criterion is in terms of We and Le, anddoes no explicitly account for the jet dynamics outside the nozzle, includ-ing the process of end-pinching. This is because the later dynamics are inprinciple dictated by these two control parameters. More specifically, Weindicates the instantaneous upward momentum of the liquid column beforea jet is produced, and Le governs how that momentum decays in time (cf.eqution 5.6). Taken together, they determine the ultimate length of the jetthat can be produced, which in turn determines whether end-pinching oc-curs and how many drops result. Equating the kinetic energy and meniscussurface energy at t? = 0 to the surface energy of a cylindrical jet of radiusRn, ?2R2n?LeV 2e + 2?R2n? = 2?RnLj?, we estimate the eventual length ofthe jet Lj once the kinetic energy has been completely converted to surfaceenergy:Lj =We4Le +Rn. (5.9)Now the numerical results of Fig. 5.7(a) can be reinterpreted in terms of Ljin Fig. 5.7(b). Roughly speaking, the transition from non-ejection to ejectionoccurs over the range of 5 ? Lj/Rn ? 7. This coincides with the critical jetlength that is determined for end-pinching on a initially stationary filament,Lj/Rn = 6 ? 1 [14]. Thus Lj provides a connection between auto-ejection,in which the mass flux at the nozzle exit varies in time, and end-pinching675.4. Results(a)?=45? ?=30?(b) rV0 0.1 0.2 0.3 0.4 0.505101520?=30??=45?Figure 5.8: (a) Pressure field inside nozzles with ? = 45? and ? = 30? when jetstarts protruding from the nozzle exit. There is a two-dimensional pressure filedwith a high pressure region around the centreline for the nozzle with ? = 45?. (b)Comparison of the axial velocity profile at the nozzle exit between ? = 30? and? = 45?.on a stationary filament where that flux is nil. The correspondence is notperfect, of course, since our jet shape can differ considerably from a perfectcylinder. At small Weber numbers, the strong capillary force makes theshape of the jet more spherical and hence increases the critical value of Lj .At high Weber numbers, the decaying velocity field at the exit producesa conical jet shape with a tapering tip. This amounts to an effectivelythinner jet diameter, and consequently a smaller critical aspect ratio Lj forbreakup. Still, the criterion of equation (5.8) is somewhat unsatisfactoryin that it is expressed in terms of We based on the instantaneous velocityVe, which is not one of the material or geometric parameters but a complexfunction of them. We have found no straightforward way to model Ve. Thisis because the acceleration of the meniscus in the nozzle depends on thedynamic contact angle ?d, which depends on the meniscus velocity in turn(cf. Fig. 5.3). Thus, we have to content ourselves for the moment with anejection criterion in terms of an instantaneous Weber number.The criterion appears consistent with the experimental data of [88].These data were presented in terms of a Weber number at the exit, sim-ilar to our We except that the local velocity was estimated using scalingarguments. Similar to our Fig. 5.6, different outcomes overlap considerablyin terms of We values. Non-ejection was observed for We from around 2up to nearly 20. The ejection of 1 or 2 droplets occurred for 6 < We < 20,while three or more drops were seen for We above 10. Since the geometricparameters were not reported for the individual data points, we are unableto compute Le and use it to untangle the data as we have done in Fig. 5.7.Thus, we can only observe that the experimental data suggest threshold We685.4. Resultst*=0.2t*=0.0 t*=0.15 t*=2.89t*=2.37(a) (b) (c) (d) (e)Figure 5.9: The regime of rapid ejection at contraction angle ? = 45?, otherconditions being identical to those in Fig. 5.4. (a) t? = 0; (b) ejection of the firstdrop at t? = 0.15; (c) ejection of the second drop at t? = 0.2; (d) ejection of thethird drop at t? = 2.37. (e) retraction of the filament.values that are consistent with our results in Fig. 5.7.Finally, the auto-ejection criterion makes an interesting prediction aboutthe possibility of auto-ejection in a straight capillary tube. The maximummeniscus velocity in a straight capillary tube is the capillary-inertial velocityvci (equation 5.1), which yields a Weber number We = 2cos ? ? 2. Thisis smaller than the minimum We for auto-ejection We = 3.4f(Le) > 3.4.Thus, auto-ejection cannot occur in straight tubes, as has been suggestedby empirical observations [70, 88].5.4.3 Rapid ejection and air entrapmentThis subsection deals with two new flow regimes encountered at large valuesof the contraction angle ?. In constructing the pinch-off criterion, we haveencoded all geometric effects into Le. The contraction induces an inwardradial flow, one consequence of which is to accelerate the average velocityof the liquid and increase the total kinetic energy. Using a one-dimensionalplug flow assumption, we have represented the acceleration effect in Le. Forlarger contraction angles, however, the two-dimensional nature of the flowbecomes important, and the radial flow tends to modify the meniscus shapeand the dynamic contact angle, thus producing new regimes of interfacialbreakup. As a baseline, we take the simulation depicted in Fig. 5.4 at695.4. Resultscontraction angle ? = 30?. Note that when the contact line reaches the exit(point e), the meniscus as a whole arrives at the exit as well, with a more orless flat interface and uniform velocity profile. Subsequently, a more or lesscylindrical jet is formed (point g), which grows to a maximum length around5Rn before end-pinching produces a single large droplet. Considering thisbaseline scenario as ?regular ejection?, we encounter two new regimes athigher ?, termed rapid ejection and air entrapment.Rapid ejection is illustrated in Fig. 5.9 for ? = 45?. The stronger con-traction leads to faster acceleration of the contact line speed, as well as alarger and increasing contact angle. Since capillarity cannot keep up withthe rapid contact line movement, the meniscus deviates markedly from aspherical shape, and a deep depression forms in the center (Fig. 5.9(a),t? = 0). Afterwards, the strong radial flow converges toward the center,while surface tension lifts the meniscus rapidly. The pressure field inside thenozzle slightly after jet?s protrusion from nozzle exit is shown in Fig. 5.8(a).These two effects conspire to produce a highly non-uniform velocity profilewhen the meniscus as a whole reaches the exit. As shown in Fig. 5.8(b), thecenterline velocity is much higher than the average velocity Vn for ? = 45?,as compared with the baseline case of ? = 30?. As a result, the first drop isejected quickly at t? = 0.15, with drop radius r = 0.07 and velocity v = 14.8.This is followed by a second small drop (r = 0.046, v = 6.9) at t? = 0.20,and a much larger third one (r = 0.54, v = 0.06) after a much longer intervalat t? = 2.37. In contrast, the baseline case has its first and only ejection att? = 1.66, producing a much larger and slower drop (r = 0.61, v = 1.02).After the third drop, the jet retracts. In view of the rapid ejection of high-speed droplets, higher ? may help induce auto-ejection under normal-gravityconditions. Indeed, the ancillary video of [87] depicts auto-ejection undernormal gravity using a large contraction angle ? = 50?. We have carriedout a limited exploration of such scenarios, and an example is depicted inFig. 5.10 for Bond number Bo = ?R2g/? = 0.4 at ? = 50?. After the ejec-tion of one droplet, the jet grows a bulb at the tip while forming a neck atthe base (t? = 0.42). Shortly afterwards, the neck pinches in and the bulbdetaches, producing two drops of disparate size (t? = 0.5). Under the sameconditions, contraction angles below 40? do not produce auto-ejection at all.It is thanks to the stronger radial flow that a thin jet forms and breaks upinto droplets.Air entrapment occurs at an even larger contraction angle of ? = 55?(Fig. 5.11). At t? = 0, the interface forms a depression as in the rapid-ejection regime. Subsequently, however, the radial flow is so strong as tocause the depression to narrow and deepen, producing an air finger. At705.4. Resultst*=0.34t*=0.0 t*=0.47 t*=0.5t*=0.42Figure 5.10: Auto-ejection under gravity for large contraction angle. Bo = 0.4,Oh = 0.01, ? = 50?, C = 0.25, L = 2. After ejecting a single droplet at t? = 0.34,the jet pinches off at its base (t? = 0.47), and later breaks up into two more drops(t? = 0.5).t? = 0.079, the neck of the air finger pinches off, entrapping a bubble inthe liquid. Given the relatively short length of the air finger, the pinch-off is dynamically driven by the inward liquid flow rather than interfacialtension as in Rayleigh-Plateau instability. After this, the strong momentumof the liquid continues to propel the jet forward, much like the later stageof Fig. 5.9. This leads to the ejection of a large drop (r = 0.49, v = 0.71) att? = 1.68. Eventually the jet retracts. Experimentally, [87] demonstratedthe possibility of the air entrapment regime under normal gravity at ? ? 50?and Oh = 0.005. This provides direct evidence for this unusual flow regime.To conclude the investigation of the contraction angle ?, we note that auto-ejection favors an intermediate range of ? values. Too gentle a contractiondoes not provide sufficient flow focusing to produce a long jet. Too abrupta contraction stifles the upward momentum of the liquid column, againsuppressing drop ejection.5.4.4 Contact line depinning at nozzle lipSo far, we have assumed the nozzle exit to be a horizontal surface of widthWlthat is completely non-wettable by the liquid (?l = 180?). Thus, the contactline is pinned at the inner corner of the lip. Under certain experimentalconditions, the contact line has been observed to depin and move outward[87]. This effectively broadens the base of the jet and changes the outcomeof drop ejection [2]. This has motivated us to relax the pinning conditionby imposing a smaller ?l so that the effect of contact line depinning can be715.4. Resultst*=0.0 t*=1.684t*=0.079 t*=2.037t*=0.069 t*=0.182(a) (b) (c) (d) (e) (f )Figure 5.11: Air entrapment at contraction angle ? = 55?, other conditions beingidentical to those in Fig. 5.4. (a) t? = 0; (b) formation of air finger at t? = 0.069;(c) pinch-off of the air bubble at t? = 0.079; (d) jet continues out of the nozzle (e)droplet pinch-off at t? = 1.684 after bubble entrapment. Air bubble disappears dueto diffusion of numerical (f ) filament retraction.investigated.Fig. 5.12 depicts the effect of contact line depinning by tracking theposition of the centreline of interface in time for several values of ?l. Asthe jet emanates from the nozzle, the interfacial slope never exceeds 90?relative to the upper surface of the lip. Thus, for ?l ? 90?, the contact lineremains pinned at the inner corner of the lip and ?l has no effect. These casesare represented by the ?l = 90? curve in Fig. 5.12(a). The geometric andphysical conditions for these runs correspond to We = 6.9 and Le = 1.46,and thus auto-ejection of a single drop occurs according to Fig. 5.7. As ?lreduces to 80? and 70?, the contact line depins and moves outward. Thishampers the lengthening of the jet and delays the pinch-off. The dropproduced is also somewhat larger. At the point of pinch-off, the contactline is somewhere on the flat part of the upper surface, not having reachedthe outer corner. For ?l ? 60?, the length of the jet is further stunted anddrop ejection is completely suppressed. For these cases, the contact linereaches the outer edge of the lip and stays pinned there, at least until thejet retracts.Fig. 5.12(b) analyzes the suppression of drop ejection for ?l = 45?. De-pinning of the contact line takes place at point a when the interface makes anangle of 45? with respect to the upper surface of the exit. After de-pinning,the contact line moves radially outward, broadening the base of the jet. This725.4. Results(a) t*Hc-L0 0.5 1 1.5 2 2.500.511.522.53?l=90?l=80?l=70?l=60?l=45baecdpinchoff(b)45 ?t*=0acdebFigure 5.12: (a) Effect of contact line depinning on the growth of the jet and dropejection. The ordinate is the length of the jet measured from the exit of the nozzle,and the abscissa is time starting from the moment of the contact line reaching theinner corner of the lip. Drop pinch-off is indicated by a round dot. The width ofthe upper surface is fixed at Wl = 0.25Rn, and ?l is the contact angle on the uppersurface. The other parameters of the simulation are C = 0.5, ? = 30?, L = 4, ? = 0(on the inner surface) and Oh = 0.01. (b) Snapshots of the interface for ?l = 45?at points marked on the curve.reduces the upward liquid velocity through mass conservation. Moreover,the curvature of the meniscus is moderated (point b), resulting in a lowercapillary pressure at the base of the jet. Both effects conspire to restrain thelengthening of the jet. The contact line reaches the outer corner of the lip atpoint c, and the jet length peaks at point d some time later. This maximumjet length, at 2R = 4Rn (Fig. 5.12a), is about 25% shorter than the casewithout contact line depinning (?l ? 90?). It is too short for drop ejection(cf. Fig. 5.7b). Thus, the jet retracts and flattens afterwards, driving thecontact line past the outer corner, producing the nearly spherical interfaceof point e.Insofar as the contact line becomes pinned at the outer corner of thelip during the growth phase of the jet, the width Wl of the lip should alsoaffect the jet behavior. For a fixed ?l = 45?, we have examined the effectof increasing Wl from 0.05Rn to 2Rn. As expected, a wider lip broadensthe base of the jet, inhibits the lengthening of the jet, and suppresses thepotential for drop ejection.735.5. Summary and conclusions5.5 Summary and conclusionsAs far as we know, this study represents the first numerical computation ofthe process of auto-ejection. In interpreting the numerical results, we havealso developed simple models to describe various aspects of the process.The parameter range captures most of the experimental conditions, and wereproduce all the salient features of the experimental observations. Themain results of the study can be summarized as follows.(a) The meniscus quickly attains the capillary-inertial regime in the straighttube, and advances with a mostly constant velocity until it enters thecontraction in the nozzle, where it accelerates. The dynamic contactangle increases with the meniscus speed. Viscosity has a negligiblerole in the entire process.(b) With the contact line pinned at the inner corner of the exit, a jet issuesinto the ambient air. The lengthening of the jet is accompanied bydeceleration of the liquid column inside the tube, with kinetic energybeing converted into surface energy. An energy balance model capturesthe temporal decay of the liquid velocity at the nozzle quite accurately.This rate of decay is dictated by an effective length that embodies thegeometric features of of the tube-nozzle combination.(c) A two-parameter criterion for auto-ejection of droplets is developedusing the instantaneous Weber number when the contact first arrivesat the nozzle exit and the effective length. Together they determinethe length of the jet that may be produced when the available kineticenergy is converted into surface energy. This critical length agrees withprior studies of end-pinching on an initially stationary filament, thusdemonstrating our criterion as being rooted in essentially the samehydrodynamics.(d) With increasing contraction angle, we predict new regimes of rapidejection of multiple drops and air bubble entrapment. When the con-traction is too mild, auto-ejection is suppressed. In particular, auto-ejection is impossible in a straight tube.(e) To the extent that comparisons can be made, the numerical resultsagree with experimental observations.One limitation of the study is that the criterion for auto-ejection is givenin terms of an instantaneous Weber number, rather than in terms of the ma-terial and geometric parameters. We attempted to model the instantaneous745.5. Summary and conclusionsvelocity at the nozzle in terms of these parameters, with little success. Ascompared with other microfluidic drop-forming procedures, auto-ejectionis unique in that it involves no external force or flux, and is entirely au-tonomous. From this standpoint, it will be desirable to devote future workto refining the current criterion into one expressed in the geometric andmaterial parameters.75Chapter 6Conclusions andRecommendationsIn this thesis, the Cahn-Hilliard diffuse-interface model is used to numeri-cally study three interfacial dynamic problems. The diffusive interface re-moves the contact line singularity and singular topological events duringthe pinch-off and coalescence. In addition, the finite thickness of the inter-face with the energy-based formulation of the Cahn-Hilliard model enablesus to capture the contact line dynamics and interfacial tension more natu-rally. In the following, a summary of key findings for each studied problemis presented. Then the significance and limitations are discussed. Finallyrecommendation for future works are made.6.1 Summary of key findings6.1.1 Capillary breakup of a liquid torusThe capillary breakup of a Newtonian liquid torus suspended in a surround-ing Newtonian liquid is studied. Starting from an externally imposed sinu-soidal disturbance, the initial stage of the growth is linear and the additionalcurvature of the torus around its axis inhibits the growth of the imposedwavelength compared to its counterpart on the straight filament. The con-traction of the torus toward its center amplifies the disturbance growth byproducing an external flow field.The final shape of the torus breakup and the number of produced dropletsare an outcome of the competition between the contraction of the torus andcapillary instability. This competition is controlled by three parameterswhich are torus-to-medium viscosity ratio, torus aspect ratio, and initialamplitude of the disturbance. A large aspect ratio for the torus lengthensthe shrinkage time. In addition, a large aspect ratio increases the distur-bance growth by reducing the axial curvature. Therefore a large aspectratio favors the capillary pinch-off mechanism. The torus-to-medium vis-cosity ratio is important for both processes. Higher viscosity ratios make766.1. Summary of key findingsthe shrinkage faster and increase the growth rate of the disturbance. Thedependence of shrinkage on viscosity ratio is stronger. Therefore, increasingthe viscosity ratio favors the shrinkage mechanism. Increasing the initialamplitude of the disturbance shortens the time for the capillary pinch-off.6.1.2 Wicking flow through microchannelsWicking flow is a key mechanism for flow movement inside the complexmicropore geometries of porous media. The dynamics of a meniscus in apore is dependent on its shape, connectivity with other pores, liquid-solid,liquid-air, and solid-air interfacial tensions, and flow properties. For inertia-less flows, an axisymmetric contraction, expansion and their combinationare used to study the effect of pore shape on the meniscus dynamics. Itis shown that the interface moves through a contraction or an expansionthrough three main steps: shape adjustment at inward corners, movementof spherical meniscus, and pinning at outward corners. At inward corners,first the contact line moves in and rotates. The tendency of the meniscus tokeep a spherical profile pushes the center of the meniscus downward. Thediffusive nature of the interface enables the meniscus to negotiate the largecontraction angles at inward corners, which presents a singularity for a sharpinterface formulation. Then the meniscus moves in with nearly a sphericalshape until it reaches the end of the contraction. At outward corners, thecontact line gets pinned and the meniscus center moves up to minimize itssurface energy until it makes enough angle with the next section to depin andmove forward. For large contraction angles, dynamics of meniscus at inwardcorner can change the passage time considerably. In expansion geometries,due to a usually slow meniscus movement, neglecting the meniscus dynamicat the corner does not affect the passage time considerably.A simple 2D Y-branching geometry is used to study the connectivity ofpores. The flow trajectory depends on the capillary forces in the branchesas well as viscous dissipation inside the root channel. Counter-intuitively,flow moves into the narrower channels if there is large dissipation inside theroot conduit.6.1.3 Auto-ejection of liquid jets and drops from capillarytubesAuto-ejection is studied for low Ohnesorge number liquids which perfectlywet the solid in a zero-gravity environment with negligible surrounding aireffects.776.2. Significance and limitationsThe dynamics of the meniscus inside the tube and the nozzle are studied.It is shown that the effect of viscosity is negligible and the meniscus dynamicscan be understood qualitatively by assuming plug flow. Inside the nozzle,the interaction of accelerating flow field and contact line dynamics producescomplex meniscus dynamics which is carefully analyzed.A two-parameter ejection criterion is developed. It is shown that theejection criterion depends on the momentum of the ejecting liquid and itsdecay rate. The momentum of the ejecting liquid is quantified in terms ofthe Weber number at the nozzle exit when the meniscus first gets pinnedthere. Its decay rate is dependent on the kinetic energy. Using a one dimen-sional flow inside the tube and nozzle, and sink flow inside the reservoir, itis possible to relate the velocity at each point inside tube-nozzle reservoircombination to the velocity at the nozzle exit. Then total kinetic energy canbe expressed in terms of an effective length and velocity at the nozzle exit. Itis shown that such an effective length is important in categorizing differentejection regimes. The importance of the effective length is further shownby relating the auto-ejection data to the critical aspect ratio for breakupof stationary filaments during the retraction. It is shown that the effectivelength is related to the length of the jet that can be produced by convertingthe liquid kinetic energy into surface energy.To obtain the effective length one dimensional plug flow assumption isused inside the tube and nozzle. Such an assumption is not valid for highcontraction angles where the strong radial flow produces two-dimensionalflow inside the nozzle. It is shown that at large contraction angles, a strongradial flow produces a highly curved meniscus inside the nozzle. This leadsto two new regimes, rapid ejection and air entrapment, with increasing con-traction angles.6.2 Significance and limitationsNowadays due to miniaturization and micro-engineering, interface dynamicand its interaction with bulk flow take on increasing scientific and practicalimportance. In auto-ejection under normal gravity, for example, the de-tailed information about the shape adjustment stage of meniscus movementthrough contractions can be used to promote droplet ejection.Information on contact line dynamics is difficult to gain through exper-imental and numerical studies. Experiments always face unwanted factorssuch as surface roughness and heterogeneity and visualization problems.Numerically, there is a singularity in macroscopic equations for the contact786.2. Significance and limitationsline movement for which there is at present no satisfactory model. TheCahn-Hilliard diffuse-interface method provides a physically motivated for-mulation to investigate contact line dynamics. By tuning its parameters fora certain fluid-solid combination, it is possible to study the dynamics whichare usually hard to capture. An example is the variation of the dynamiccontact angle inside the nozzle for the auto-ejection problem.Understanding the pore scale interfacial dynamics is the building blockfor developing better models for porous media flow. Such understanding canbe used to design more efficient gas diffusion layers for fuel cells. Two mainfeatures of pores are their non-uniform cross section and connectivity. Theperformed research helps to understand the interface dynamics in these twogeometries.Capillary instability plays a main role in most of the droplet productionmechanisms. Knowledge of capillary instability is further extended by takinginto account the effect of filament curvature and compressive flow field. Inaddition, liquid rings are unstable configurations and will eventually contractonto themselves or breakup into droplets. Therefore, the simplified twotimes-scale model will give an insight into the fate of the torus.We have compared our numerical simulations with experimental andtheoretical results. There is a good agreement between physical observationsand numerical results which further demonstrates the ability of the Cahn-Hilliard model to produce physically meaningful results.There are also limitation for this research, which we summarize below:? Computational limitation. Realization of the sharp interface limit iscomputationally very expensive. This becomes a more severe issuefor three dimensional interfaces and when there are disparate lengthscales in the problems. Besides, our AMPHI algorithm uses a fullyimplicit time-updating scheme. Although this improves stability, thecomputational cost in inverting the matrix is high. More efficient splitalgorithms might help alleviate the problem.? Limitation of the Cahn-Hilliard contact line model. The phenomeno-logical nature of the contact line model implies that model parametersneed to be determined from experimental data, which are not alwaysavailable for the geometry and materials required. In addition, re-solving the diffusion length which is around six orders of magnitudesmaller than the bulk length scale is numerically challenging.796.2. Significance and limitations6.2.1 Recommendations for future workThe current work can be further extended in certain aspects. In the analysisof the torus instability, we have not considered the effect of non-symmetricdisturbances. In reality the geometry of the torus is not a prefect ring andit can have defects on its surface. How such defects modify the breakupof liquid torus can be studied. In addition, we showed that the initiallyfastest mode may not proceed to the end. How a combination of differentmodes grow on the torus and how they morph into each other needs to bestudied. It is also interesting to study the capillary instability of the ringin the presence of gravity where thick and thin parts of the disturbed toruswill experience different buoyancy forces and the growth of the disturbancewill be affected.For wicking flow one can experimentally study the flow branching andalso extend it to more realistic geometries for pore-connectivity, including3D branches and networks. Such studies will provide a benchmark for thepopular pore-network models. In addition, it can be extended to includethe effect of gravity on meniscus dynamics in Y-branches and also throughcontraction or expansion.For auto-ejection, it will be desirable to devote future work to refiningthe current criterion into one expressed in the geometric and material pa-rameters. In addition, one can explore the dependency of droplet size ongeometric and solid-fluid properties. We have done a limited test for a 1-gcondition, which can be explored in more detail. Finally, the modeling ofthe dynamic contact angle in the Cahn-Hilliard model needs to be studiedmore carefully, especially in regards to the energy dissipation at the contactline.80Bibliography[1] M. Ahmadlouydarab, Z.-S. Liu, and J. J. Feng. Interfacial flows incorrugated microchannels: flow regimes, transitions and hysteresis.Int. J. Multiphase Flow, 37:1266?1276, 2011.[2] B. Ambravaneswaran, H. J. Subramani, S. D. Phillips, and O. A.Basaran. Dripping-jetting transitions in a dripping faucet. Phys. Rev.Lett., 93:34501, 2004.[3] B. Ambravaneswaran, E. D. Wilkes, and O. A. Basaran. Drop forma-tion from a capillary tube: Comparison of one-dimensional and two-dimensional analyses and occurrence of satellite drops. Phys. Fluids,14:2606?2621, 2002.[4] S. L. Anna, N. Bontoux, and H. A. Stone. Formation of dispersionsusing ?flow focusing? in microchannels. Appl. Phys. Lett., 82:364?366,2003.[5] A. Antkowiak, N. Bremond, J. Duplat, S. Le Dize`s, and E. Villermaux.Cavity jets. Phys. Fluids, 19(9):091112?091112, 2007.[6] H. Benkreira and M. I. Khan. Air entrainment in dip coating underreduced air pressures. Chemical Engineering Science, 63(2):448 ? 459,2008.[7] G. Blankenstein and U. D. Larsen. Modular concept of a laboratoryon a chip for chemical and biochemical analysis. Biosens. Bioelectron.,13:427?438, 1998.[8] C. H. Bosanquet. On the flow of liquids into capillary tubes. Philos.Mag. Ser. 6, 45:525?531, 1923.[9] J. B. Bostwick and P. H. Steen. Stability of constrained cylindri-cal interfaces and the torus lift of Plateau-Rayleigh. J. Fluid Mech.,647:201?219, 2010.81Bibliography[10] S. Bouaidat, O. Hansen, H. Bruus, C. Berendsen, N. K. Bau-Madsen,P. Thomsen, A. Wolff, and J. Jonsmann. Surface-directed capillarysystem; theory, experiments and applications. Lab Chip, 5:827?836,2005.[11] M. Bracke, F. De Voeght, and P. Joos. The kinetics of wetting: thedynamic contact angle. In Trends in Colloid and Interface Science III,pages 142?149. Springer, 1989.[12] G. Caginalp and X. Chen. Convergence of the phase field model to itssharp interface limits. Euro. J. Appl. Math., 9:417?445, 1998.[13] J. W. Cahn and J. E. Hilliard. Free energy of a non-uniform system.Part I: interfacial free energy. J. Chem. Phys., 28:258?267, 1958.[14] A. Castrejo?n-Pita, J. Castrejo?n-Pita, and I. Hutchings. Breakup ofliquid filaments. Phys. Rev. Lett., 108(7), 2012.[15] J. R. Castrejo?n-Pita, A. A. Castrejo?n-Pita, E. J. Hinch, J. R. Lister,and I. M. Hutchings. Self-similar breakup of near-inviscid liquids.Phys. Rev. E, 86(1):015301, 2012.[16] R. Chein and S. H. Tsai. Microfluidic flow switching design usingvolume of fluid model. Biomed. Microdevices, 6:81?90, 2004.[17] J. M. Chen, P. C. Huang, and M. G. Lin. Analysis and experimentof capillary valves for microfluidics on a rotating disk. Microfluid.Nanofluid., 4:427?437, 2008.[18] I. Cohen, M. P. Brenner, J. Eggers, and S. R. Nagel. Two fluid dropsnap-off problem: Experiments and theory. Phys. Rev. Lett., 83:1147?1150, 1999.[19] O. Cohu and H. Benkreira. Air entrainment in angled dip coating.Chemical Engineering Science, 53(3):533 ? 540, 1998.[20] R. G. Cox. The dynamics of the spreading of liquids on a solid surface.Part 1: Viscous flow. J. Fluid Mech., 168:169?194, 1986.[21] H. Czachor. Modeling the effect of pore structure and wetting angleson capillary rise in soils having different wettabilities. J. Hydrol.,328:604?613, 2006.82Bibliography[22] N. Djilali. Computational modelling of polymer electrolyte membrane(PEM) fuel cells: Challenges and opportunities. Energy, 32:269?280,2007.[23] M. Dreyer, A. Delgado, and H. J. Rath. Capillary rise of liquid betweenparallel plates under microgravity. J. Colloid Interface Sci., 163:158?168, 1994.[24] J. Eggers and E. Villermaux. Physics of liquid jets. Rep. Prog. Phys.,71:036601, 2008.[25] D. Erickson, D. Li, and C. B. Park. Numerical simulation of thecapillary driven flows in non-uniform cross sectional capillaries. J.Colloid Interface Sci., 250:422?430, 2002.[26] L. A. Freitag and C. Ollivier-Gooch. Tetrahedral mesh improvementusing swapping and smoothing. Int. J. Numer. Meth. Engng., 40:3979?4002, Jul 1997.[27] P. Gao and J. J. Feng. Enhanced slip on a patterned substrate due todepinning of contact line. Phys. Fluids, 21:102102, 2009.[28] P. Gao and J. J. Feng. A numerical investigation of the propulsion ofwater walkers. J. Fluid Mech., 668:363?383, 2011.[29] P. Gao and J. J. Feng. Spreading and breakup of a compound dropon a partially wetting substrate. J. Fluid Mech., 682:415?433, 2011.[30] S. Gekle and J. M. Gordillo. Generation and breakup of worthing-ton jets after cavity collapse. part 1. jet formation. J. Fluid Mech.,663:293?330, 2010.[31] T. Goldmann and J. Gonzalez. DNA-printing: utilization of a stan-dard inkjet printer for the transfer of nucleic acids to solid supports.Journal of biochemical and biophysical methods, 42(3):105?110, 2000.[32] D. A. Gomes. Stability of rotating liquid films. Q. J. Mech. Appl.Math., 55:327?343, 2002.[33] R. C. Goodknight, W. A. Klikoff, and I. Fatt. Non-steady-state fluidflow and diffusion in porous media containing dead-end pore volume.J. Phys. Chem., 64(9):1162?1168, 1960.83Bibliography[34] J. M. Gordillo and S. Gekle. Generation and breakup of worthingtonjets after cavity collapse. part 2. tip breakup of stretched jets. J. FluidMech., 663:331?346, 2010.[35] G. J. Gutierrez, A. Medina, and F. J. Higuera. Equilibrium height of aliquid in a gap between corrugated walls under spontaneous capillarypenetration. J. Colloid Interface Sci., 338:519?522, 2009.[36] J. Ha and L. G. Leal. An experimental study of drop deformation andbreakup in extensional flow at high capillary number. Phys. Fluids,13:1568, 2001.[37] S. Herminghaus, M. Brinkmann, and R. Seemann. Wetting and dewet-ting of complex surface geometries. Annu. Rev. Mater. Res., 38:101?121, 2008.[38] R. L. Hoffman. A study of the advancing interface. I. Interface shapein liquid-gas systems. J. Colloid Interface Sci., 50:228?241, 1975.[39] D. Jacqmin. Calculation of two-phase Navier-Stokes flows using phase-field modeling. J. Comput. Phys., 155:96?127, 1999.[40] D. Jacqumin. Contact-line dynamics of a diffuse fluid interface. J.Fluid Mech., 402:57?67, 2000.[41] T.-S. Jiang, O.-H. Soo-Gun, and J. C. Slattery. Correlation for dy-namic contact angle. J. Colloid Interface Sci., 69, 1979.[42] P. Joos, P. Van Remoortere, and M. Bracke. The kinetics of wettingin a capillary. J. Colloid Interface Sci., 136:189?197, 1990.[43] D. L. Keij, K. G. Winkels, H. Castelijns, M. Riepen, and J. H. Snoeijer.Bubble formation during the collision of a sessile drop with a meniscus.Phys. Fluids, 25(8):082005, 2013.[44] G. B. Lee, C. I. Hung, B. J. Ke, G. R. Huang, and B. H. Hwei.Micromachined pre-focused 1?N flow switches for continuous sampleinjection. J. Micromech. Microeng., 11:567, 2001.[45] S. Leib and M. Goldstein. Convective and absolute instability of aviscous liquid jet. Phys. Fluids, 29, 1986.[46] G. Lin and T. V. Nguyen. Effect of thickness and hydrophobic poly-mer content of the gas diffusion layer on electrode flooding level in aPEMFC. J. Electrochem. Soc., 152:A1942?A1948, 2005.84Bibliography[47] W. W. Liou, Y. Peng, and P. E. Parker. Analytical modeling of cap-illary flow in tubes of nonuniform cross section. J. Colloid InterfaceSci., 333:389?399, 2009.[48] S. Litster, D. Sinton, and N. Djilali. Ex situ visualization of liquidwater transport in pem fuel cell gas diffusion layers. J. Power Sources,154:95?105, 2006.[49] R. Lucas. Ueber das zeitgesetz des kapillaren aufstiegs vonflu?ssigkeiten. Kolloid Z., 23:15?22, 1918.[50] J. D. McGraw, J. Li, D. L. Tran, A.-C. Shi, and K. Dalnoki-Veress.Plateau-Rayleigh instability in a torus: formation and breakup of apolymer ring. Soft Matter, 6:1258?1262, 2010.[51] H. Mehrabian and J. J. Feng. Wicking flow through microchannels.Phys. Fluids, 23:122108, 2011.[52] T. Mikami, R. G. Cox, and S. G. Mason. Breakup of extending liquidthreads. Int. J. Multiphase Flow, 2:113?138, 1975.[53] T. Mikami and S. G. Mason. The capillary break-up of a binary liquidcolumn inside a tube. Can. J. Chem. Eng., 53:372?377, 1975.[54] J. H. Nam and M. Kaviany. Effective diffusivity and water-saturationdistribution in single- and two-layer pemfc diffusion medium. Int. J.Heat Mass Transfer, 46:4595?4611, 2003.[55] T. D. Nguyen, M. Fuentes-Cabrera, J. D. Fowlkes, J. A. Diez, A. G.Gonza?lez, L. Kondic, and P. D. Rack. Competition between collapseand breakup in nanometer-sized thin rings using molecular dynamicsand continuum modeling. Langmuir, 28:13960?13967, 2012.[56] P. K. Notz and O. A. Basaran. Dynamics and breakup of a contractingliquid filament. J. Fluid Mech., 512(1):223?256, 2004.[57] E. Pairam and A. Ferna?ndez-Nieves. Generation and stability oftoroidal droplets in a viscous liquid. Phys. Rev. Lett., 102:234501,Jun 2009.[58] G.-G. Park, Y.-J. Sohn, T.-H. Yang, Y.-G. Yoon, W.-Y. Lee, and C.-S. Kim. Effect of PTFE contents in the gas diffusion media on theperformance of PEMFC. J. Power Sources, 131:182?187, 2004.85Bibliography[59] U. Pasaogullari and C. Y. Wang. Liquid water transport in gas dif-fusion layer of polymer electrolyte fuel cells. J. Electrochem. Soc.,151:A399?A406, 2004.[60] D. Patro, S. B., and V. Jayaram. Flow kinetics in porous ceramics:Understanding with non-uniform capillary models. J. Am. Ceram.Soc., 90:3040?3046, 2007.[61] L. M. Pismen. Nonlocal diffuse interface theory of thin films and themoving contact line. Phys. Rev. E, 64:021603, Jul 2001.[62] L. M. Pismen. Mesoscopic hydrodynamics of contact line motion.Colloids Surf., A, 206:11?30, 2002.[63] A. Pizzi and K. L. Mittal. Handbook of Adhesive Technology. MarcelDekker, New York, 2003.[64] D. Que?re?, E?. Raphae?l, and J. Ollitrault. Rebounds in a capillary tube.Langmuir, 15(10):3679?3682, 1999.[65] L. Rayleigh. On the instability of jets. Proc. Lond. Math. Soc., 10:4?13, 1878.[66] M. Reyssat, L. Courbin, E. Reyssat, and H. A. Stone. Imbibition ingeometries with axial variations. J. Fluid Mech., 615:335?344, 2008.[67] R. M. S. M. Schulkes. The contraction of liquid filaments. J. FluidMech., 309:277?300, 1996.[68] R. Sharma and D. S. Ross. Kinetics of liquid penetration into period-ically constricted capillaries. J. Chem. Soc., Faraday Trans., 87:619?624, 1991.[69] J. Shen. Efficient spectral-Galerkin method. II. direct solvers of secondand fourth order equations by using Chebyshev polynomials. SIAMJ. Sci. Comput., 16:74?87, 1995.[70] R. Siegel. Transient capillary rise in reduced and zero-gravity fields.J. Appl. Mech., 83:165?170, 1961.[71] W. A. Sirignano and C. Mehring. Review of theory of distortion anddisintegration of liquid streams. Prog. Energy Combust. Sci., 26:609?655, 2000.86Bibliography[72] M. Stange, M. E. Dreyer, and H. J. Rath. Capillary driven flow incircular cylindrical tubes. Phys. Fluids, 15:2587, 2003.[73] T. L. Staples and D. G. Shaffer. Wicking flow in irregular capillaries.Colloids Surf., A, 204:339?350, 2002.[74] H. A. Stone and L. G. Leal. Relaxation and breakup of an initiallyextended drop in an otherwise quiescent fluid. J. Fluid Mech., 198:399?427, 1989.[75] M. Switkes. Bubbles in immersion lithography. Journal of vacuumscience and technology. B Microelectronics processing and phenomena,23:217 ? 236, 2005.[76] J. Szekely, A. W. Neumann, and Y. K. Chuang. The rate of capillaryoenetration and the applicability of the Washburn equation. J. ColloidInterface Sci., 35:273?278, 1971.[77] S. Takeuchi, P. Garstecki, D. B. Weibel, and G. M. Whitesides. Anaxisymmetric flow-focusing microfluidic device. Adv. Mater., 17:1067?1072, 2005.[78] T. Tanuma. Innovative hydrophilic microporous layers for cathode gasdiffusion media. J. Electrochem. Soc., 157(12):B1809?B1813, 2010.[79] M. Tjahjadi, H. A. Stone, and J. M. Ottino. Satellite and subsatelliteformation in capillary breakup. J. Fluid Mech., 243:297?317, 1992.[80] S. Tomotika. On the instability of a cylindrical thread of a viscousliquid surrounded by another viscous fluid. Proc. R. Soc. London,Ser. A, 150:322?337, 1935.[81] S. Tomotika. Breaking up of a drop of viscous liquid immersed inanother viscous fluid which is extending at a uniform rate. Proc. R.Soc. London, Ser. A, 153:302?318, 1936.[82] A. Y. Tong and Z. Wang. Relaxation dynamics of a free elongatedliquid ligament. Phys. Fluids, 19:092101, 2007.[83] A. Umemura. Self-destabilizing mechanism of a laminar inviscid liquidjet issuing from a circular nozzle. Phys. Rev. E, 83(4):046307, 2011.[84] A. S. Utada, E. Lorenceau, D. R. Link, P. D. Kaplan, H. A. Stone,and D. A. Weitz. Monodisperse double emulsions generated from amicrocapillary device. Science, 308:537?541, 2005.87Bibliography[85] E. W. Washburn. The dynamics of capillary flow. Phys. Rev., 17:273?283, 1921.[86] A. Wollman. Capillarity-driven droplet ejection. Master?s thesis, Port-land State University, USA, 2012.[87] A. Wollman, T. Snyder, D. Pettit, and M. Weislogel. Spontaneouscapillarity-driven droplet ejection. arXiv:1209.3999, 2012.[88] A. Wollman and M. Weislogel. New investigations in capillary fluidicsusing a drop tower. Exp. in Fluids, 54(4):1?13, 2013.[89] Y. Wu, J. D. Fowlkes, P. D. Rack, J. A. Diez, and L. Kondic. Onthe breakup of patterned nanoscale copper rings into droplets viapulsed-laser-induced dewetting: Competing liquid-phase instabilityand transport mechanisms. Langmuir, 26:11972?11979, 2010.[90] Z.-N. Wu. Approximate critical Weber number for the breakup of anexpanding torus. Acta Mech., 166:231?239, 2003.[91] Q. Xu and O. A. Basaran. Computational analysis of drop-on-demanddrop formation. Phys. Fluids, 19:102111, 2007.[92] Z. Yao and M. Bowick. The shrinking instability of toroidal liquiddroplets in the Stokes flow regime. Euro. Phys. J. E, 34:1?6, 2011.[93] P. Yue and J. J. Feng. Can diffuse-interface models quantitativelydescribe moving contact lines? Eur. Phys. J - Spec. Top., 197:37?46,2011.[94] P. Yue and J. J. Feng. Wall energy relaxation in the Cahn-Hilliardmodel for moving contact lines. Phys. Fluids, 23:012106, 2011.[95] P. Yue, C. Zhou, and J. J. Feng. A computational study of the coales-cence between a drop and an interface in Newtonian and viscoelasticfluids. Phys. Fluids, 18:102102, 2006.[96] P. Yue, C. Zhou, and J. J. Feng. Sharp interface limit of the cahn-hilliard model for moving contact lines. J. Fluid Mech., 645:279?294,2010.[97] P. Yue, C. Zhou, J. J. Feng, C. F. Ollivier-Gooch, and H. H. Hu. Phase-field simulations of interfacial dynamics in viscoelastic fluids usingfinite elements with adaptive meshing. J. Comput. Phys., 219:47?67,2006.88Bibliography[98] B. V. Zhmud, F. Tiberg, and K. Hallstensson. Dynamics of capillaryrise. J. Colloid Interface Sci., 228:263?269, 2000.[99] C. Zhou, P. Yue, and J. J. Feng. Formation of simple and compounddrops in microfluidic devices. Phys. Fluids, 18:092105, 2006.[100] C. Zhou, P. Yue, J. J. Feng, C. F. Ollivier-Gooch, and H. H. Hu.3D phase-field simulations of interfacial dynamics in Newtonian andviscoelastic fluids. J. Comput. Phys., 229:498?511, 2010.[101] D. Zhou, P. Yue, and J. J. Feng. Viscoelastic effects on drop deforma-tion in a converging pipe flow. J. Rheol., 52:469?487, 2008.[102] M.-Y. Zhou and P. Sheng. Dynamics of immiscible-fluid displacementin a capillary tube. Phys. Rev. Lett., 64:882?885, 1990.89


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items