Crossover from Shear-Driven to Thermally Activated Drainage of Liquid-Infused Microscale Capillaries
Carlos E. Colosqui, Jason S. Wexler, Ying Liu, Howard A. Stone
CCrossover from Shear-Driven to Thermally Activated Drainageof Liquid-Infused Microscale Capillaries
Carlos E. Colosqui, ∗ Jason S. Wexler,
2, 3
Ying Liu, and Howard A. Stone † Department of Mechanical Engineering, Stony Brook University, Stony Brook, NY 11794, USA. Otherlab, San Francisco, CA 94110, USA. Department of Mechanical and Aerospace Engineering,Princeton University, Princeton, New Jersey 08544, USA.
The shear-driven drainage of capillary grooves filled with viscous liquid is a dynamic wettingphenomenon relevant to numerous industrial processes and novel lubricant-infused surfaces. Priorwork has reported that a finite length L ∞ of the capillary groove can remain indefinitely filledwith liquid even when large shear stresses are applied. The mechanism preventing full drainage isattributed to a balance between the shear-driven flow and a counterflow driven by capillary pressurescaused by deformation of the free surface. In this work, we examine closely the approach to the finalequilibrium length L ∞ and report a crossover to a slow drainage regime that cannot be describedby conventional dynamic models considering solely hydrodynamic and capillary forces. The slowdrainage regime observed in experiments can be instead modeled by a kinetic equation describinga sequence of random thermally activated transitions between multiple metastable states caused bysurface defects with nanoscale dimensions. Our findings provide new insights on the critical rolethat natural or engineered surface roughness with nanoscale dimensions can play in the imbibitionand drainage of capillaries and other dynamic wetting processes in microscale systems. I. INTRODUCTION
Dynamic wetting processes such as spreading, imbibition, and drainage are ubiquitous in natural, agricul-tural, and industrial processes that are crucial to modern technology. Engineering applications ranging from oilrecovery and water treatment to microfluidics and bioanalytical systems have been enabled by a fundamental un-derstanding of wetting that is embodied in mathematical descriptions such as the Young-Dupre, Young-Laplace,and Lucas-Washburn equations [1, 2]. These classical wetting models are derived in the framework of continuumthermodynamics under the assumption of perfectly smooth and homogeneous surfaces and predict dynamic be-haviors that are governed by deterministic forces due to capillary action and hydrodynamic effects. Althoughthese assumptions can reasonably describe wetting phenomena in macroscale systems, random thermal fluctua-tions and the microscopic details of the surface must be properly considered to understand interfacial transportprocesses at micro- and nanoscales. With the advent of micro- and nanofabrication techniques a comprehen-sive understanding of dynamic wetting has become essential to improve traditional industrial processes suchas surface coating and spraying and to fully exploit the potential of modern fabrication techniques such asmicro/nanolitography and additive manufacturing (or 3D printing).As the system dimensions shrink to micrometer scales and below, roughness and chemical heterogeneitiesinherent to natural and artificial surfaces pose a major challenge in modeling wetting processes [1, 3]. Giventhe multiscale nature of the microscopic structure of solid surfaces it is not always feasible to define a singlecharacteristic dimension. Nevertheless, surface roughness and heterogeneities are usually characterized by a“defect” size s d , determined by some relevant dimension given by the root-mean-square (rms) roughness, heightautocorrelation length, or other topographic parameters. For “macroscopic” defect sizes s d >
100 nm, thermalfluctuations can be neglected and for low Capillary numbers the dominant forces are due to elastic deformation ofthe interface and pinning at localized defects [1, 4–6]. These elastic and pinning forces are merely the consequenceof changes in interfacial energies as the contact line moves over random surface heterogeneities of physical and/orchemical nature. When multiple “macroscopic” defects collectively distort and pin the contact line, the energybarriers preventing net displacement give rise to contact angle hysteresis [7–11]. The conventional approachto consider the effects of random surface defects with macroscopic ( s d >
100 nm) or mesoscopic ( s d (cid:39) θ Y , which is determined by minimization of energy on a perfectly smooth surface. Despiteavailable predictive models based on the Wenzel [12] and Cassie-Baxter [13] equations, no analytical approach ∗ [email protected] † [email protected] a r X i v : . [ c ond - m a t . s o f t ] M a y has been established to quantitatively predict the degree of contact angle hysteresis from topographic parameterscharacterizing the surface [3, 14–16]. As result, receding and advancing contact angles for static and dynamicconditions for different surfaces and liquid pairs must often be determined empirically.It is necessary to model the effects of random thermal motion when surface defects have dimensions smallerthan 100 nm and become comparable to the nanoscale thermal fluctuations of the liquid interface. The interplaybetween thermal motion and nanoscale surface features can lead to nontrivial wetting processes that are inducedby thermal fluctuations of the contact line [17–26]. A few different approaches have been proposed to modelthe effect thermal motion and nanoscale surface defects s d ≤ E between adsorption sites separated by a distance λ (cid:39) s d . Thevirtual frictional force proposed in MKT scales linearly with viscosity and its magnitude is often comparable tohydrodynamic forces, which can make it difficult to distinguish between damping due to pinning at nanoscaledefects or hydrodynamic effects [2, 16, 29, 31]. Energy barriers ∆ E = W a in MKT are determined by the“work of adhesion” W a (cid:39) γA d (1 + cos θ Y ) at localized sites, here γ is the liquid-vapor surface tension γ and A d ∼ s d the area of the adsorption site. Predictions from MKT show agreement with experimentally observeddisplacement rates for different liquid pairs by assuming nanoscale defect sizes s d = 0.2–1 nm (e.g., see Ref. [16]).For consistency with the model assumptions of MKT the defect size must be smaller than 1 nm ( A d ∼ ),which yields energy barriers ∆ E (cid:46) k B T (here k B is the Boltzmann constant and T the system temperature).Notably, a series of recent experimental studies on diverse systems indicate that even larger defect sizes ofthe order of 10 nm can induce wetting processes that are thermally activated. For example, experimentalobservations report that single colloidal particles at water-oil interfaces exhibit surprisingly slow adsorption rateswith time scales to reach equilibrium conditions on the order of several hours or even days [32, 33]. According toconventional wetting models for perfectly spherical particles [34, 35], the adsorption dynamics of single particles isa fast monotonic decay to stable equilibrium conditions where the system energy is a global minimum. The slowadsorption rates observed for diverse microparticles were attributed to thermally activated processes induced bysurface defects with sizes ranging from 1 to 5 nm [32]. Studies of the spreading dynamics of low viscosity liquidson surfaces with defect sizes of 10 nm report that the contact line displacement is governed by thermally activatedprocesses [19, 20, 36–38]. These studies [37, 38] indicate that energy barriers prescribing the displacement rate ofthe contact line are significantly smaller than the work of adhesion, and thus energy barriers ∆ E (cid:28) γs d inducedby mesoscopic defects are smaller than predicted from the defect size.The “kinetics” of contact line displacement on surfaces with mesoscopic defects s d = 1–100 nm can be describedby wetting models based on Kramers theory of thermally activated transitions [24, 39, 40]. In this approach, theenergy barrier ∆ E and separation distance λ between long-lived metastable states can have a nontrivial relationwith the defect size s d since these quantities are determined by projecting the multidimensional energy landscapeparametrized by molecular positions and velocities onto a one-dimensional energy profile along the “reaction”coordinate describing the contact line displacement [24, 26, 40]. Theoretical models recently proposed by Colosqui et al. [24] support the idea that kinetic rates determined via Kramers theory [41, 42] can predict the displacementrates of contact lines in the presence of mesoscopic defects ( s d = 1–10 nm). According to these models [24] itis possible to observe both a fast dynamic regime, governed by capillary forces and hydrodynamic friction, or amuch slower kinetic regime governed by thermally activated processes. The distance from equilibrium at whichthe regime crossover takes place is determined by the energy barrier magnitude and defect size, as well as thelength of the contact line perimeter [24].Previous studies by Wexler et al. [43, 44] have reported the shear-driven drainage of oil-infused microgroovesand identified conditions where a finite volume of oil is retained for indefinitely long time. The observed steadystates were analytically predicted by establishing a balance between capillary forces and the applied shear stress[43]. The drainage dynamics far from equilibrium was approximately described by a Lucas-Washburn-typeequation where thermal motion is neglected and the microgroove surfaces are assumed to be macroscopicallysmooth but having a receding contact angle significantly different from the Young contact angle. Given thatthe drainage of the microgrooves involves the displacement of a contact line perimeter of microscale dimensions,similar phenomena observed in the adsorption of microparticles at water-oil interfaces [32, 33] is expected to affectthe drainage dynamics. Indeed, experimental observations by Wexler et al. show that the drainage dynamics closeto steady-state conditions presents deviations from analytical predictions from the proposed Lucas-Washburn-type equation [43].In the present work we extend the Lucas-Washburn-type equation for shear-driven drainage in order to con-sider thermal motion and the presence of nanoscale surface roughness, by following the approach proposed byColosqui et al. for microparticle adsorption [24]. Atomic force microscopy (AFM) is employed to characterize thesurface roughness and thus determine the defect dimensions used in the proposed wetting model for thermallyactivated wetting. While the rms roughness seems to determine the magnitude of the energy barriers ∆ E , theheight autocorrelation length appears to determine the separation distance λ between metastable states. Theproposed model employing mesoscopic defect sizes (3–30 nm) determined via AFM describes the drainage dynam-ics observed close to equilibrium conditions for different oil viscosities and applied shear rates. The agreementbetween the observed contact line displacements and analytical predictions indicate that the drainage close toequilibrium is dominated by thermally activated transitions between metastable states. Moreover, we propose acriterion for estimating the crossover point where the drainage transitions from dynamics governed by capillaryand hydrodynamic forces to a kinetic regime dominated by thermally activated processes. II. SYSTEM DESCRIPTION
The experimental system consists of a rectangular microfluidic cell fabricated from Norland epoxy and sealedwith a transparent glass lid for visualization purposes (see Fig. 1(a)). The microfluidic flow cell has width W cell = 7 mm, height H cell = 0 .
18 mm, and length L cell = 45 mm and is filled with a 1:1 weight mixture ofglycerol and water (i.e., the outer aqueous phase) with viscosity µ aq = 5 . ρ aq = 1150 kg/m .There is one additional port that is 10 mm downstream of the outlet slot; this port is used for filling the oil atthe beginning of the experiment, and is closed when the experiment is performed. A syringe pump maintainsconstant volumetric flow rates ( Q =1–2 mL/min) in the aqueous phase via injection of fluid through an inletport upstream of the microgrooves.As illustrated in Figs. 1(a)–(b), on one wall of the microfluidic cell there is a parallel array of 50 rectangularmicrogrooves of width w = 9 µ m, height h = 10 µ m, and length (cid:96) = 36 mm, which are infused with a siliconeoil that is immiscible with the aqueous phase. Two different silicone oils are used to infuse the microgrooves: 1)1,1,5,5-Tetraphenyl-1,3,3,5-tetramethyltrisiloxane (Gelest PDM-7040), with viscosity µ o = 42 . ρ = 1061 kg/m , and interfacial tension (with the aqueous solution) γ = 29 mN/m; and 2) 1,1,3,5,5-Pentaphenyl-1,3,5-trimethyltrisiloxane (Gelest PDM-7050) with viscosity µ o = 201 mPa-s, density ρ = 1092 kg/m , andinterfacial tension (with the aqueous solution) γ = 28 . T (cid:39) ± ◦ C. (d) (a) (c) (nm)1.5 - (𝜇m) (e) 𝐿 ∞ = 5.6 mm
18 min 𝑡 = 33 min ti m e 𝐿(𝑡) 𝑄 𝐻 𝑐𝑒𝑙𝑙 𝐿 𝑐𝑒𝑙𝑙 silicone oilwater/glycerol 𝑥𝑦 𝑤ℎ 𝐿(𝑡)ℓ 𝑥 𝑦 z (𝜇m) (b) FIG. 1. Experimental configuration. (a) Schematic of the microfluidic flow cell (not to scale). An array of 50 microgrooves(bottom wall) is infused with silicone oil (green) and connected to an oil reservoir at the flow cell terminus. (b) Schematic ofthe geometry of a single groove. (c) Image sequence (3 min between images) of a sample shear-driven drainage experiment( Q = 2 mL/min, µ o = 42 . µ m × µ m) of the groove surface. After the syringe pump starts to inject the water/glycerol mixture, a finite time t S must elapse before reachingsteady flow conditions with the prescribed volumetric rate Q . A time t S = ρ aq l /µ aq (cid:39)
150 s can be estimated byconsidering solely diffusive effects; this time is in good agreement with experimental observations for all the flowrates studied in this work. As shown in the image sequence in Fig. 1(c), the outer flow drives the gradual dewettingof the oil infused in the microgrooves until reaching a final finite length L ∞ , after which the microgrooves remainpartially filled indefinitely; the time to reach the final length L ∞ is on the order of thousands of seconds underthe studied conditions. Assuming plane Poiseuille flow and a large viscosity ratio µ o /µ aq (cid:29)
1, and given that themicrogrooves are aligned with the outer flow, the shear stress applied at the oil-water interface is estimated as τ xy = 6 µ aq Q/W cell H cell . The predicted stress τ xy is employed to describe experimental observations except forthe case of low viscosity oil and high flow rate where the shear stress employed is 15% smaller than analyticallyestimated; this deviation is attributed to the finite viscosity ratio ( µ o /µ aq = 7 .
9) for the latter case. The Reynoldsnumber in the aqueous phase is Re = (3 / ρ aq Q/W cell µ aq (cid:39) τ xy can be attributed to deviations from plane Poiseuille flow and end effects. Since the Reynoldsnumber in the oil phase is O (10 − ) and the Bond number is O (10 − ), inertial and gravitational effects can beneglected inside the microgrooves.The microfluidic device is molded from Norland Optical Adhesive (NOA 81) using the “sticker” technique[43, 45]. The array of microgrooves is molded from PDMS that is in turn molded from an etched silicon waferwith the nominal cross-section profile shown in Fig. 1(d). The cross-section profile of the microgroove arraypresents micron-scale deviations from the nominal geometry that are below 5% and can be observed by opticalmicroscopy. This small “error of form” is expected to cause small deviations from the flow conditions predictedfor the nominal microgroove geometry (see Fig. 1(d)). Analysis of the microgroove surfaces is performed witha scanning probe microscope (Bruker Dimension Icon) operating in AFM tapping mode (PeakForce Tapping (cid:114) )with a height resolution of 0.1 nm and lateral spatial resolution of 2 nm. Topographic imaging via AFM (seeFig. 1(e)) reveals a complex random topography with nanoscale physical features resembling peaks and valleyswith maximum heights and depths on the order of 3 nm and lateral dimensions reaching up to 50 nm. Asdiscussed in detail in the next section, the presence of nanoscale roughness is expected to cause pinning of thecontact line and thermally activated processes that lead to significant deviations from the dewetting dynamicspredicted for a perfectly smooth surface. III. THEORETICAL MODELING
As in previous work by Wexler et al. [43], we begin by assuming unidirectional creeping flow in the oil insidethe microgrooves so that the streamwise fluid velocity u ( y, z, t ) satisfies the governing equations ∂u/∂x = 0 and µ o ∇ u − dp/dx = 0 for mass and linear momentum balances; here, µ o is the dynamic viscosity of the oil and p ( x, t )is the pressure in the oil phase. For the studied experimental configuration and given that the oil is much moreviscous than the aqueous solution we will assume a constant pressure p o in the external aqueous phase. Under theassumed incompressible flow conditions the pressure inside the microgroove must vary linearly ( dp/dx = const.)and so must the curvature of the top free surface κ = 1 /r ( x ) since a pressure drop ∆ p = − γ/r ( x ) (for r (cid:28) L ) isinduced by capillary effects. Hence, the pressure inside the oil is p ( x, t ) = p + ( γ/r min )( x/L ) where [43] r min = w/ (2 cos θ ) for wh ≤ θ + tan θ ) h (cid:0) w/ h ) (cid:1) for wh > θ + tan θ ) , (1)is the minimum radius of curvature at the downstream end ( x = (cid:96) − L ( t )) determined by the receding contactangle θ (see Fig. 1(b)). A receding contact angle θ = 56 ± ◦ has been previously determined from experimentalmeasurements [43] and since w/h = 0 . r min = w/ (2 cos θ ) according to Eq. (1). For the assumedcurvature profile of the oil-water interface the oil volume inside the microgroove is V ( t ) = c d whL ( t ) where [43] c d = 1 − r min h (cid:32) − (cid:115) − w r min (cid:33) + r min wh arcsin (cid:18) w r min (cid:19) . (2)Conservation of mass determines that the rate of change of oil volume c d wh dLdt = − ( q s + q p ) (3)inside the grooves is determined by the volumetric flow rates q s driven by the applied shear force F s = τ xy wL ,and q p induced by the force F p = − ( γ/r min ) wh due to capillary pressure. Assuming creeping flow conditionsand a rectangular cross-section for the liquid-filled region, analytical solution of the momentum conservationequations gives the corresponding volumetric rates and conductivities: q s = c s h µ o L F s with c s = 12 − hw ∞ (cid:88) n =0 ( − n b n tanh (cid:18) b n w h (cid:19) , (4)and q p = c p h µ o L F p with c p = 13 − hw ∞ (cid:88) n =0 ( − n b n tanh (cid:18) b n w h (cid:19) . (5)Here, b n = ( n + 1 / π are the eigenvalues for each Fourier mode in the analytical solution of the momentumequation. For the nominal microgroove height and width in the experiments of Wexler et al. [43] we have c d = 0 . c s = 6 . × − , and c p = 4 . × − . Combining volume and momentum conservation lawsembodied in Eqs. (3)–(5) we arrive to a Lucas-Washburn-type (L-W) equation [43] dLdt = − c d µ o (cid:18) c s τ xy h − c p γh r min L (cid:19) . (6)This equation was derived in prior work by Wexler et al. [43] and predicts that for t → ∞ , for which dL/dt = 0, the system reaches a stationary or final length L ∞ = ( c p hγ ) / ( c s r min τ xy ). Introducing the finallength in Eq. (6) the equation for the displacement rate takes the simple form dL/dt = − U LW (1 − L ∞ /L ), where U LW = ( c s /c d )( τ xy h/µ ) determines the maximum displacement rate attained for L/L ∞ (cid:29)
1. Integrating thedisplacement rate dL/dt in Eq. (6) leads to an implicit expression for the column length: t = t S + L ∞ U LW (cid:20) log (cid:18) L ( t ) − L ∞ L ( t S ) − L ∞ (cid:19) + L ( t S ) − L ( t ) L ∞ (cid:21) , (7)where t S is the time after which stationary flow conditions are attained in the aqueous phase.A few comments are in order about the derivation of Eqs. (6)–(7). Predictions from Eqs. (6)–(7) are valid fora constant shear stress τ xy assuming Poiseuille flow in the aqueous phase, and thus t S (cid:39)
150 s in Eq. (7) is thefinite time required to reach steady state conditions in the outer phase (as discussed in Sec. II). The derivationassumes a contact line perimeter of length s = 2 h + w that is uniform and has a constant receding contact angle θ , which implies the assumption of a perfectly flat surface with constant and spatially homogeneous contactangle hysteresis. Nanoscale surface roughness and/or chemical heterogeneities induce spatial fluctuations of thecontact line position and local contact angle that are associated with “pinning” at localized surface defects.Thermally activated depinning becomes the dominant mechanism inducing contact line displacement as thesystem approaches the equilibrium length L → L ∞ where the effective driving force F d = − c s F s + c p F p → A. Surface heterogeneities and thermal motion
The L-W equation (Eq. (6)) describes a one-dimensional model of drainage dynamics characterized by a singlevariable L ( t ) when considering deterministic forces due to hydrodynamic and capillary effects on a macroscopicallysmooth surface. As shown in Fig. 2(a), 2D topographical imaging via AFM of a microscale section of the surfacereveals a random distribution of surface defects with a maximum (peak-to-peak) height of about 6 nm. Analysisof the surface topography reveals a nearly Gaussian probability distribution of defect heights h d (Fig. 2(b)) thatis commonly observed for random (non-patterned) surfaces. The surface height presents a small rms roughness h rms = 0.85 nm; the height distribution skewness is 0.3 and its kurtosis is 3.3, which are very close to the valuesexpected for a Gaussian distribution. The height autocorrelation is isotropic and presents a nearly Gaussiandecay (Fig. 2(c)) with the radial distance r and a radial correlation length r d = 26 . s d = √ r d (cid:39) . A d = πs d = 4 . × − µ m .As illustrated in Fig. 2(d), we will consider that the path x ( σ, t ) (0 ≤ σ ≤ s ) defined by the local streamwise po-sition of the contact line along its perimeter s is distorted by the surface defects detected in the AFM topographicimage (Fig. 2(a)). The average streamwise position of the contact line ¯ x ( t ) = (1 /s ) (cid:82) s x ( σ, t ) dη determines the(projected) surface area A = ( l − ¯ x ) s wetted by the liquid and thus the liquid column length L ( t ) = A/s . Hence,the wetting/dewetting of a single surface defect with (projected) surface area A d increases/reduces the liquidcolumn length by an amount λ = A d /s (see Fig. 2(d)). For simplicity we assume that the arclength s (cid:39) h + w of the contact line is approximately constant; assuming negligible variations of s implies neglecting contributionsto the system energy due to line tension [46]. We will further consider that surface defects with a finite height h d (cid:39) h rms > E in the energy E ( L ) required to vary theliquid column length L , as illustrated in Figs. 2(e)–(f). The energy fluctuation magnitude ∆ E is determined bycomplex morphological changes of the liquid-liquid and liquid-solid interfaces that are induced by surface defects.Moreover, adsorption of water or oil molecules at mesoscopic voids created by the substrate topography and in-terfacial phenomena induced by steric effects are likely to cause significant variations of the local surface energies.Given this complexity, the magnitude of the characteristic energy barrier ∆ E induced by surface defects will beconsidered as a model parameter that can be obtained by fitting experimental observations. Nevertheless, model-ing surface defects as cones with base area A d = πs d and height h d = h rms determined by AFM imaging we cananalytically estimate an energy barrier of magnitude ∆ E (cid:39) γs d h rms | − ( π/
2) cos θ | = 1 . × − J = 3 . k B T ;as illustrated in Fig. 2(f) the motion of the contact line over a modeled defect involves changes ∆ A wo = s d h rms inthe water-oil interfacial area and ∆ A od = ( π/ s d h rms in the surface area wetted by the oil phase. As expectedthe analytically estimated energy barrier vanishes for a perfectly flat surface with h rms = 0. 𝐴 𝑑 𝐴 𝑑 𝑬(𝑳) 𝛤 + 𝛤 − 𝑬 𝒐 𝑟 [nm] 𝐶 ( 𝑟 ) exp(−𝑟 /𝑟 𝑑 ) 𝑟 𝑑 ≃ 26.5nm AFM - (d)(a) (e)(b)(c) 𝑓 ( ℎ 𝑑 ) ℎ 𝑑 [nm] Gaussian
AFM 𝑠 𝑑 ≃ 2𝑟 𝑑 = 37.5nm ℎ 𝑟𝑚𝑠 = 0.85nm height(AFM) external flow 𝜆 𝜆 𝑥(𝜎, 𝑡) 𝐿 + 𝜆 𝐿 𝐿 − 𝜆 𝐴 𝑑 𝑥 𝜎 𝐿 + 𝜆
𝐿 𝐿 − 𝜆 (f)
𝐿 + 𝜆/2 𝐿 𝑥𝜎 𝑥(𝜎, 𝑡) Δ𝐸2 Δ𝐸2 Δ𝐸 ℎ 𝑟𝑚𝑠 𝑠 𝑑 Δ𝐴 𝑤𝑜 𝐿 − 𝜆/2 Δ𝐴 𝑜𝑑 FIG. 2. Nanoscale roughness and energy barriers. (a) Two-dimensional AFM image of a sample section of the microgroovesurface. (b) Local defect height h d distribution computed from AFM data, showing a nearly Gaussian distribution( h rms = 0 .
85 nm). (c) Autocorrelation function computed from AFM data (radial correlation length r d = 26 . s d (cid:39) . A d (cid:39) πs d .(e) Energy profiles E o ( L ) for h rms = 0 (dashed red line) and E ( L ) for h rms > E (cid:39) γs d h rms | − ( π/
2) cos θ | = 3 . k B T . In order to incorporate the effects of nanoscale surface defects and thermal fluctuations of the contact line wewill begin by considering L ( t ) as a generalized coordinate, or reaction coordinate, determined by the surface areawetted by the oil. Accordingly, we can recast Eq. (6) as dL/dt = − (1 /ξ )( dE o /dL ) where ξ = c d µ o s is an effectiveresistivity and E o ( L ) = s (cid:20) c s τ xy hL − c p γh r min log (cid:18) LL (cid:19)(cid:21) (8)is the energy required to change the liquid column length for the case of a smooth groove with h rms = 0 ( L is an arbitrary reference length, which results in the addition of an arbitrary constant in Eq. (8)). The energyprofile E o has a global minimum when the stationary length is reached and thus dE o /dL → L → L ∞ .For analytical simplicity, the effect of heterogeneities or localized surface defects will be modeled by adding asingle-mode perturbation to the smooth-surface energy E o so that the energy to vary the liquid column length is E ( L ) = E o ( L )+(∆ E/
2) sin(2 π ( L − L ∞ ) /λ + ϕ ); the arbitrary phase ϕ = − π/ L = L ∞ . Given that λ (cid:28) L , multiple local energy minima will exist at L o (cid:39) L ∞ ± nλ ( n is an integer) when the system is sufficiently close to equilibrium ( L → L ∞ ) where dE o /dL →
0. Therefore, for L → L ∞ the system exhibits multiple metastable configurations separated by different energy barriers ∆ E ± = E ( L o ± λ/ − E ( L o ) in the forward/backward (+ / − ) directions and thermal motion becomes the dominant effectinducing transitions between neighboring metastable states.To consider thermally activated processes, we incorporate in the L-W equation (Eq. (6)) for the column lengthdynamics a stochastic thermal force F th = √ k B T ξη ( t ), where η ( t ) is zero-mean and unit-variance Gaussiannoise; this thermal force F th is determined by means of the fluctuation-dissipation theorem. Including energyfluctuations caused by surface defects and stochastic forces induced by random thermal motion in Eq. (6) thedrainage dynamics is described by a Langevin-type equation dLdt = − ξ ddL (cid:20) E o + ∆ E (cid:18) πλ ( L − L ∞ ) − π (cid:19)(cid:21) + √ Dη ( t ) , (9)where D = k B T /ξ is the (long-time) diffusivity along the “reaction coordinate” defined by the liquid columnlength L . B. Near equilibrium dynamics
The smooth-surface energy in Eq. (8) has a global minimum at L = L ∞ and can be accurately approximatedby a second-order Taylor expansion E o ( L ) = ( d E o /dL ) | L = L ∞ × ( L − L ∞ ) for L − L ∞ < (3 / L ∞ . Hence for L/L ∞ < / E ( L ) = K L − L ∞ ) + ∆ E (cid:18) πλ ( L − L ∞ ) − π (cid:19) , (10)where K ≡ d E o dL (cid:12)(cid:12)(cid:12)(cid:12) L = L ∞ = c s τ xy r min sc p γ . (11)According to Eq. 9, as L → L ∞ and dE o /dL → L undergoes a random walk in a periodicpotential with multiple minima (i.e., metastable states) located at L o (cid:39) L ∞ ± nλ . Near equilibrium the columnlength L ( t ) will fluctuate around the local minima L o and will suddenly transition, or “hop”, to neighboringminima if crossing over the neighboring maxima at L ± = L o ± λ/ / − ) transition rates (cf. Fig. 2(e)) are given byΓ ± ( L ) = 12 πξ (cid:115) d E ( L o ) ∂L (cid:12)(cid:12)(cid:12)(cid:12) d E ( L ± ) ∂L (cid:12)(cid:12)(cid:12)(cid:12) exp (cid:20) − ( E ( L ± ) − E ( L o )) k B T (cid:21) (12)for | L − L o | < λ/
2. When “hopping” between metastable states at rates given by Eq. (12) the average drainagespeed can be estimated by a rate equation dL/dt = λ (Γ + − Γ − ) and thus we have [24] dLdt = − U H sinh (cid:18) L − L ∞ L H (cid:19) , (13)where the characteristic “hopping” velocity is U H = λ (cid:113) π/λ ) ∆ E − K πξ exp (cid:20) − (∆ E + Kλ / k B T (cid:21) , (14)and the “hopping” length is L H = 2 k B TKλ . (15)Integration of Eq. (13) leads to L ( t ) = L ∞ + L H arctanh (cid:20) exp (cid:18) − U H L H ( t − t o ) (cid:19)(cid:21) , (16)where t o is an initial time arising from the integration constant. Eq. (16) is valid for times t ≥ t c where t c is the crossover time after which the drainage dynamics is dominated by thermally activated processes. Aselaborated in the next section, one can analytically estimate a critical crossover length L c below which forcesresulting from surface heterogeneities and thermal motion are larger than forces due to hydrodynamic shear andcapillary pressure. Accordingly, the initial t o in Eq. (16) is determined to match the experimental condition L ( t c − t o ) = L c , where the crossover time t c in each experiment corresponds to the time elapsed to reach theanalytically estimated length L c . C. Regime crossover
Far from equilibrium conditions where the liquid column length is much larger than the equilibrium length L ( t ) (cid:29) L ∞ , the drainage dynamics is dominated by hydrodynamic shear and capillary forces, and can thus bedescribed with the L-W approach in Eqs. (6)–(7) [43]. As mechanical equilibrium is approached L → L ∞ and dE o /dL →
0, hydrodynamic and capillary forces balance out and the drainage of the microgrooves becomes athermally activated process described by Eqs. (13)–(16).Here, we aim to develop a criterion for predicting the crossover from shear-driven to thermally activateddrainage for different geometries and physical conditions. For this purpose we will analytically estimate a crit-ical column length L c below which the dynamics is dominated by random forces due to spatial fluctuationsof surface energy and thermal motion. For overdamped systems, the frictional force is equal to the sum (cid:80) F of all other (non-frictional) forces and thus ξ ( dL/dt ) = ( (cid:80) F ). While according to Eq. (9) the displacementrate is dL/dt = − (1 /ξ )( dE o /dL ) when hydrodynamic and capillary forces dominate, Eq. (13) determines that dL/dt = − U H sinh[( L − L ∞ ) /L H ] near equilibrium conditions where surface energy fluctuations and thermalmotion dominate. Hence, there must be critical column length L c for whichsinh (cid:18) L c − L ∞ L H (cid:19) = 1 ξU H dE o dL ( L c ) , (17)and forces resulting from random surface energy fluctuations and thermal motion are approximately equal to thesum of hydrodynamic and capillary forces. Once the critical length L c is obtained by solving Eq. (17) one canemploy Eq. (13) to determine a critical displacement rate magnitude U c = | U H sinh[( Lc − L ∞ ) /L H ] | below whichthe drainage process is thermally activated.It is worth remarking that the crossover between regimes is actually a gradual process and takes place overa range of lengths L ( t ) (cid:39) L c . For the sake of simplicity, however, we will assume the transition to thermallyactivated drainage occurs at a “crossover” point determined by the critical length L c implicitly defined byEq. (17). The integration constant in Eq. (16) will be determined to match the critical length L ( t c ) = L c that isexperimentally observed at a time t = t c for each studied condition, and thus t o = t c + ( L H /U H ) log { tanh[( L c − L ∞ ) / (2 L H )] } .In prior work [24] a simple explicit expression alternative to Eq. (17) was proposed to estimate the criticaldistance from equilibrium below which the final relaxation regime is dominated by thermally activated transitionsbetween metastable states. According to Eq. (10), metastable states induced by local energy minima where dE/dL = 0 can only exist for sufficiently small column lengths L ≤ L ∞ + ( π ∆ E ) / ( Kλ ). Hence, the approach toequilibrium is dominated by thermally activated transitions below a crossover length L c given by [24] L c − L ∞ L H = α π Ek B T (18)where α < α = 0.2–0.25. IV. RESULTS
The length of the wetted portion of a groove is determined by using automated image analysis on macroscalephotographs with a pixel size of 12.5 µ m. The pixel intensity is high in places that are wetted with oil (due tofluorescence) and low elsewhere. The upstream limit of the wetted length is determined by plotting the pixelintensity along the length of a groove, and finding the location where the slope changes most rapidly by applyinga third-order Savitzky-Golay filter with a window size of 50–70 pixels. These images are taken every 10 seconds,yielding a limit to the resolvable velocity of approximately 10 − m/s.Three different experimental conditions are studied where the outer flow rate and viscosity of the infused oil arevaried: (i) Q = 2 mL/min and µ o = 201 mPa-s (cf. Figs. 3(a)–(b)), (ii) Q = 2 mL/min and µ o = 42 . Q = 1 mL/min and µ o = 42 . dL/dt and time evolution of the column length L ( t ) measured experimentally are compared in Fig. 3 against analyticalpredictions from the L-W approach (Eqs. (6)–(7)) and the theory based on thermally activated transitions betweenmetastable states (Eqs. (13)–(16)).As discussed in Sec.II, a finite time t S = 150 s is employed in Eq. (7) to consider the time elapsed before steadyflow is attained in the aqueous phase; this is in agreement with experimental observations for the displacementrate magnitude reported in Fig. 3. For case (ii) where the highest volumetric rate ( Q = 2 mL/min) is employed 𝐿 [ 𝑚 ] ExperimentalEq. (16)Eq. (7) | 𝑑 𝐿 / 𝑑 𝑡 | [ 𝑚 / 𝑠 ] ExperimentalEq. (13)Eq. (6) 𝐿 𝑐 = 1.98𝐿 ∞ 𝐿 [ 𝑚 ] t [s] | 𝑑 𝐿 / 𝑑 𝑡 | [ 𝑚 / 𝑠 ] t [s] Experimental
Eq. (13)
Eq. (6) Experimental
Eq. (16)
Eq. (7) 𝐿 𝑐 = 1.82𝐿 ∞ (b)(a) (d)(c) (f) 𝐿 [ 𝑚 ] t [s] ExperimentalEq. (16)Eq. (7) (e) | 𝑑 𝐿 / 𝑑 𝑡 | [ 𝑚 / 𝑠 ] t [s] ExperimentalEq. (13)Eq. (6) crossovercrossover t [s] t [s] crossovercrossover 𝐿 𝑐 = 2.2𝐿 ∞ crossover crossover FIG. 3. Displacement rate magnitude | dL/dt | and column length L ( t ) versus time for three different experimental condi-tions. (a-b) Case (i): Q = 2 mL/min and µ o = 201 mPa-s. (c-d) Case (ii): Q = 2 mL/min and µ o = 42 . Q = 1 mL/min and µ o = 42 . t S = 150 s. Solid lines: analytical predictions for drainage dominated by thermallyactivated processes (Eqs. (13)–(16)) using λ = 0 .
15 nm and ∆ E = 3 . k B T ( T = 24 ◦ C ). Dashed-dotted (horizontal) lines:predictions from Eq. (17) for the crossover length L c . The initial time t o = t c + ( L H /U H ) log { tanh[( L c − L ∞ ) / (2 L H )] } inEq. (16) is determined to match the experimentally observed length at the crossover point L ( t c − t o ) = L c . and the liquid phase has the lowest viscosity ( µ o = 42 . τ xy = 4 .
04 Pa employed inEqs. (6)–(7) was 15% lower than predicted by assuming plane Poiseuille flow and a large viscosity ratio. For theother experimental conditions the shear stress employed in Eqs. (6)–(7) was the one predicted by assuming planePoiseuille flow; i.e., τ xy = 4 .
75 Pa for case (ii), and τ xy = 2 .
38 Pa for case (iii). After steady flow conditions areattained for t ≥ t S , there is good agreement between experimental observations and analytical predictions fromL-W equations (Eqs. (6)–(7)) during the initial stages of drainage where L ( t ) < L c and hydrodynamic shear andcapillary forces are expected to dominate.As the system approaches the final equilibrium length L ∞ there is a crossover to a slower drainage processpredicted by Eqs. (13)–(16), which are valid when the dynamics are dominated by thermally activated processes.In all studied cases, the period between metastable configurations λ = πs d / (2 h + w ) = 0 .
15 nm was determinedby the defect size s d (cid:39) . | 𝑑 𝐿 / 𝑑 𝑡 | / 𝑈 𝐻 (𝐿 − 𝐿 ∞ )/𝐿 𝐻 (a)(b) (c) Eq. (13)
Eq. (6)
Eq. (13)
Eq. (6)Eq. (13)
Eq. (6) crossover ( 𝛼 = 0.26 ) (d) ( 𝐿 − 𝐿 ∞ ) / 𝐿 𝐻 𝑄 = 2 mL/min 𝜇 = 42.7 Pa-s 𝑄 = 1 mL/min 𝜇 = 42.7 Pa-s 𝑄 = 2 mL/min 𝜇 = 201 Pa-s Eq. (16)
Eq. (7) crossover ( 𝛼 = 0.25 ) ( 𝑡 − 𝑡 𝑜 ) 𝑈 𝐻 /𝐿 𝐻 crossover ( 𝛼 = 0.25 )crossover ( 𝛼 = 0.21 ) | 𝑑 𝐿 / 𝑑 𝑡 | / 𝑈 𝐻 | 𝑑 𝐿 / 𝑑 𝑡 | / 𝑈 𝐻 FIG. 4. Approach to final equilibrium length L ∞ for three different experimental conditions: Case (i): Q = 2 mL/min, µ o = 201 mPa-s, L H /U H = 1087 s ( U H = 3 . × − m/s, L H = 3 . × − m). Case (ii): Q = 2 mL/min, µ o = 42 . L H /U H = 261 . U H = 1 . × − m/s, L H = 4 . × − m). Case (iii): Q = 1 mL/min, µ o = 42 . L H /U H = 1029 s ( U H = 1 . × − m/s, L H = 1 . × − m). (a) Normalized displacement ratemagnitude | dL/dt | /U H versus normalized distance from equilibrium length ( L ( t ) − L ∞ ) /L H . (b) Distance from equilibrium L ( t ) − L ∞ versus normalized time (length shown in logarithmic scale). A nearly exponential decay with a characteristictime T H = L H /U H is observed for all studied cases. Markers: experimental results for cases (i)–(iii). Solid lines: analyticalpredictions from Eqs. (13)–(16) using λ = 0 .
15 nm and ∆ E = 3 . k B T ( T = 24 ◦ C ). Dashed-dotted (horizontal) line:analytical estimation for the crossover length L c from Eq. (18). results reported in Figs. 3–4 an energy barrier magnitude ∆ E (cid:39) . k B T ( T = 24 ◦ C) is employed for all cases.Notably, the value of the energy barrier employed to fit experimental observations can be predicted via simplegeometric arguments (cf. Fig. 2) for the three studied conditions where the flow rate, viscosity, and surfacetension are varied. Moreover, the crossover criterion in Eq. (17) (see dashed-dotted horizontal lines in Fig. 3)can be used to estimate the critical lengths L c below which the drainage becomes a thermally activated processand L ( t ) is governed by Eq. (16). For the experimental conditions in case (i) (cf. Figs. 3(a)–(b)) the crossoverto thermally activated drainage occurs for t c (cid:39) L ( t c ) = 11 . L c = 1 . L ∞ ). In agreement with experimentalobservations for cases (ii) and (ii) (cf. Figs.( 3)(c)–(f)), Eq. (17) predicts an increase in the crossover lengthand an earlier transition to thermally activated drainage when the liquid viscosity is reduced. In particular, thecrossover criterion (Eq. (17)) indicates that for the lower flow rates employed in case (iii) (cf. Figs. 3(e)–(f))the crossover length is larger than the microgroove length and the entire drainage dynamics may be thermallyactivated.According to the theoretical model leading to Eqs. (13)–(16), all experimental observations near equilibriumconditions can be collapsed to a single curve when normalizing with the characteristic “hopping” velocity U H and length L H defined by Eq. (14) and Eq. (15), respectively. Indeed, Figs. 4(a)–(c) report that the displacementrate magnitude closely follows the single curve predicted by Eq. (13) for all studied cases (i)–(iii). Similarly, thedistance L ( t ) − L ∞ between the column length and the expected equilibrium length follows the single trajectorypredicted by Eq. (16) when normalized by the corresponding values of U H and L H for each case (Fig. 4(d)).The linear decay in the displacement rate magnitude for ( L − L ∞ ) /L H < L ( t ) − L ∞ ∝ exp( − t/T H ), near equilibrium conditions with a relaxation time T H = L H /U H varying from about200 to 1000 s (cf. Fig. 4(b)). In addition we observe that the simple crossover criterion in Eq. (18) can predictthe crossover length L c for scaling factors α (cid:39) V. CONCLUSIONS
The analysis and experimental observations in this work indicate that the interplay between nanoscale surfaceroughness and thermal motion needs to be carefully considered in order to describe the dynamics of drainage andimbibition in microscale capillaries. In the presence of significant energy barriers induced by nanoscale surfacedefects, the interface displacement is dominated by random thermally activated transitions between metastable1states. These random transitions give rise to a “kinetic” regime in the evolution of the surface area wetted by oneor other phase that cannot be described by conventional (continuum-based) wetting models (e.g., L-W equations)considering solely deterministic forces due to hydrodynamic and capillary effects. Therefore we have proposed astochastic Langevin equation that can be used to describe both the (far-from-equilibrium) dynamic and (near-equilibrium) kinetic regimes observed in the shear-driven drainage of microcapillaries infused with viscous liquid.The proposed model can be adopted to describe numerically diverse wetting processes, such as spreading ofmicrodroplets or colloidal particle adsorption, where thermal motion and nanoscale surface roughness give riseto the same fundamental phenomena considered in this work.To describe analytically the kinetic regime dominated by thermally-activated processes, we have employeda rate equation where transition rates are predicted by Kramers theory. Furthermore, we have considered anenergy profile exhibiting multiple metastable states with a characteristic period λ = 0 .
15 nm and separated bya characteristic energy barrier ∆ E (cid:39) . k B T . In the model proposed in this work, both the period and energybarrier are determined by nanoscale defects with characteristic size s d (cid:39) . h rms = 0 .
85 nmthat are observed in AFM topographic images. It is worth noticing that an energy barrier of magnitude 3 . k B T corresponds to the work of adhesion W a = γ (1 + cos θ ) A a on a molecular adsorption site of area A a = 0 .
32 nm .Thus, fitting experimental results by using an alternative wetting model such as MKT would have led us toinfer that the drainage dynamics near equilibrium is caused by surface defects of molecular dimensions s d (cid:39)√ A a = 0.6 nm. Notably, AFM imaging of the studied surfaces reported the presence nano- and mesoscaledefects with much larger dimensions ( s d >
10 nm) and areas ( A d > )). The model employed in this workdetermines that the very small separation between metastable states ( λ ∼ O (10 − m)) is given by the ratio ofthe surface defect area ( A d ∼ O (10 − m)) to the contact liner perimeter ( s ∼ O (10 − m)), i.e., it is not directlyprescribed by the physical distance between surface defects. The proper definition of model parameters made itpossible to predict both the crossover to the kinetic regime and the kinetic relaxation rate for all of the studiedexperimental conditions.The analysis in this work shows that it is feasible to characterize the nanoscale surface topography, usingAFM or alternative approaches, and then determine the system dimensions (e.g., capillary height and width)that will produce a desired drainage dynamics. While the final retention length L ∞ is prescribed by specificgeometric and physical parameters, the time to reach the final length can be significantly reduced/increased by(i) reducing/increasing the crossover length L c to the kinetic regime and (ii) decreasing/increasing the kineticrelaxation time T H = L H /U H , which varies exponentially with the energy barrier ∆ E prescribed by the surfacedefect area A d . The models employed in this work could aid the design of nanostructured surfaces to control thedynamics of drainage of capillaries as well as other wetting processes in microscale systems.We thank Dr. Chung-Chueh Chang at the Stony Brook University (SBU) ThInc for performing AFM imagingof the microgrooves samples. CEC acknowledges support from the SEED Grant Program by Brookhaven NationalLaboratory and SBU. JSW, YL, and HAS acknowledge support from ONR MURI Grants No. N00014-12-1- 0875and No. N00014-12-1-0962 (Program Manager Dr. Ki-Han Kim) [1] P. G. de Gennes, Rev. Mod. Phys. , 827 (1985).[2] D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley, Rev. Mod. Phys. , 739 (2009).[3] D. Qu´er´e, Annu. Rev. Mater. Res. , 71 (2008).[4] J. Joanny and P.-G. de Gennes, J. Chem. Phys. , 552 (1984).[5] M. O. Robbins and J.-F. Joanny, EPL , 729 (1987).[6] A. Prevost, E. Rolley, and C. Guthmann, Phys. Rev. B , 064517 (2002).[7] R. E. Johnson Jr and R. H. Dettre, J. Phys. Chem , 1744 (1964).[8] C. Huh and S. Mason, J. Colloid Interface Sci. , 11 (1977).[9] J. Oliver, C. Huh, and S. Mason, Colloids Surf. , 79 (1980).[10] C. Extrand and Y. Kumagai, J. Colloid Interface Sci. , 378 (1997).[11] S. Ramos, E. Charlaix, A. Benyagoub, and M. Toulemonde, Phys. Rev. E , 031604 (2003).[12] R. N. Wenzel, Ind. Eng. Chem. Res. , 988 (1936).[13] A. Cassie and S. Baxter, J. Chem. Soc. Faraday Trans. , 546 (1944).[14] G. McHale, Langmuir , 8200 (2007).[15] A. Marmur and E. Bittoun, Langmuir , 1277 (2009).[16] M. Ramiasa, J. Ralston, R. Fetzer, and R. Sedev, Adv. Colloid Interface Sci. , 275 (2014).[17] B. Cherry and C. Holmes, J. Colloid Interface Sci. , 174 (1969).[18] A. Marmur, Adv. Colloid Interface Sci. , 121 (1994).[19] E. Rolley and C. Guthmann, Phys. Rev. Lett. , 166105 (2007).[20] A. Prevost, E. Rolley, and C. Guthmann, Phys. Rev. Lett. , 348 (1999). [21] B. Davidovitch, E. Moro, and H. A. Stone, Phys. Rev. Lett. , 244505 (2005).[22] F. Restagno, L. Bocquet, T. Biben, and ´E. Charlaix, J. Phys. Condens. Matter , A419 (2000).[23] M. Ramiasa, J. Ralston, R. Fetzer, R. Sedev, D. M. Fopp-Spori, C. Morhard, C. Pacholski, and J. P. Spatz, J. Am.Chem. Soc. , 7159 (2013).[24] C. E. Colosqui, J. F. Morris, and J. Koplik, Phys. Rev. Lett. , 028302 (2013).[25] A. M. Rahmani, Y. Shao, M. Jupiterwala, and C. E. Colosqui, Phys. Fluids , 082004 (2015).[26] C. E. Colosqui, T. Teng, and A. M. Rahmani, Phys. Rev. Lett. , 154504 (2015).[27] T. Blake and J. Haynes, J. Colloid Interface Sci. , 421 (1969).[28] S. Semal, T. Blake, V. Geskin, M. J. De Ruijter, G. Castelein, and J. De Coninck, Langmuir , 8765 (1999).[29] M. J. De Ruijter, J. De Coninck, and G. Oshanin, Langmuir , 2209 (1999).[30] T. Blake and J. De Coninck, Adv. Colloid Interface Sci. , 21 (2002).[31] D. Duvivier, D. Seveno, R. Rioboo, T. Blake, and J. De Coninck, Langmuir , 13015 (2011).[32] D. M. Kaz, R. McGorty, M. Mani, M. P. Brenner, and V. N. Manoharan, Nat. Mater. , 138 (2012).[33] A. Wang, D. M. Kaz, R. McGorty, and V. N. Manoharan, AIP Conf. Proc. , 336 (2013).[34] P. Pieranski, Phys. Rev. Lett. , 569 (1980).[35] B. P. Binks and T. S. Horozov, Colloidal Particles at Liquid Interfaces (Cambridge University Press, 2006).[36] E. Rolley, C. Guthmann, and M. Pettersen, Phys. Rev. Lett. , 016101 (2009).[37] K. Davitt, M. S. Pettersen, and E. Rolley, Langmuir , 6884 (2013).[38] L. Du, H. Bodiguel, and A. Colin, Phys. Rev. E , 012402 (2014).[39] T. Blake and J. De Coninck, Eur. Phys. J. , 249 (2011).[40] S. Razavi, I. Kretzschmar, J. Koplik, and C. E. Colosqui, J. Chem. Phys. , 014904 (2014).[41] H. Kramers, Physica , 284 (1940).[42] P. Hanggi, J. Stat. Phys. , 105 (1986).[43] J. S. Wexler, I. Jacobi, and H. A. Stone, Phys. Rev. Lett. , 168301 (2015).[44] I. Jacobi, J. S. Wexler, and H. A. Stone, Phys. Fluids , 082101 (2015).[45] D. Bartolo, G. Degr´e, P. Nghe, and V. Studer, Lab Chip , 274 (2008).[46] A. Marmur, J. Colloid Interface Sci.186