Why Capillary Flows in Slender Triangular Grooves Are So Stable Against Disturbances
aa r X i v : . [ phy s i c s . f l u - dyn ] D ec Why Capillary Flows in Slender Triangular GroovesAre So Stable Against Disturbances
Nicholas C. White and Sandra M. Troian ∗ California Institute of Technology, T. J. Watson Sr. Laboratories of Applied Physics,1200 E. California Blvd., MC 128-95, Pasadena, CA 91125, USA (Dated: December 27, 2018)Ongoing development of fuel storage and delivery systems for space probes, interplanetary vehicles,satellites and orbital platforms continues to drive interest in propellant management systems thatutilize surface tension to retain, channel and control flow in microgravity environments. Althoughit has been known for decades that capillary flows offer an ideal method of fuel management, therehas been little research devoted to the general stability properties of such flows. In this work,we demonstrate theoretically why capillary flows which channel wetting liquids in slender opentriangular channels tend to be very stable against disturbances. By utilizing the gradient flow formof the governing fluid interface equation, we first prove that stationary interfaces in the presenceof steady flow are asymptotically nonlinearly and exponentially stable in the Lyapunov sense. Wethen demonstrate that fluid interfaces exhibiting self-similar Washburn dynamics are transiently andasymptotically linearly stable to small perturbations. This second finding relies on a generalizednon-modal stability analysis due to the non-normality of the governing disturbance operator. Takentogether, these findings reveal the robust nature of transient and steady capillary flows in opengrooved channels and likely explains the prevalent use of capillary flow management systems inmany emerging technologies ranging from cubesats to point-of-care microfluidic diagnostic systems.
I. INTRODUCTION
It has been known for centuries that wetting liquidswill rapidly and spontaneously creep along surfaces con-taining grooves, interior corners, crevices or roughenedareas, a process known as wicking. Since the late 1960s,researchers have been incorporating this passive and re-liable method of flow control in the design of novel pro-pellant management devices able to store, channel andmeter fuel resourcefully in microgravity environments [1–5]. Such systems have significantly extended missionlifetimes of spacecraft and satellites, enabling future in-terplanetary explorations as well. Modern propellantmanagement systems consist of combinations of sponges,traps, troughs, vanes, and wicks to channel propellantflow by capillary action, systems which have been inves-tigated extensively [2, 6–10]. Shown in Fig. 1 are twocommon structures designed in such a way that the liq-uid film thickness is much smaller than the streamwiseflow distance, a limiting ratio which as described in Sec-tion III leads to considerable simplification of the gov-erning equations of motion. However, despite that bythe 1920’s the mechanism describing internal capillaryflow, whereby a liquid column spontaneously fills the in-terior of a slender capillary tube, was understood and theappropriate equations developed, the mechanism drivingspontaneous capillary flow along an open grooved chan-nel required a half-century more to be deduced. Ongoingefforts to miniaturize fluid management systems for manydifferent applications continue to drive interest in the fun-damentals of free surface capillary flow along structuredsubstrates. ∗ (a)(c) (d)(b) © Edward Washburn appears to have provided the firsttheoretical analysis of the distance in time traversed bya Newtonian liquid flowing within a horizonal enclosedcylindrical tube under the action of capillary forces [11].The Washburn relation, as it is now known, is given bythe relation z = ( γ cos θ R t/ µ ) / , where z is the dis-tance traversed by the advancing (for wetting fluids) orreceding (for non-wetting fluids) meniscus, t is time, R is the cylinder radius, and µ , γ and θ , respectively, arethe liquid shear viscosity, surface tension and contact an-gle. Although this largely unidirectional flow is drivenby the capillary pressure drop across the curved menis-cus separating gas from liquid, the Washburn relationcan be regarded from a mathematical point of view asan example of diffusive driven interface growth since the“diffusion” coefficient preceding the time variable is de-scribed by units of length squared per unit time. Bytreating porous bodies as an assemblage of small cylin-drical capillaries, Washburn also went on to predict thatthe volume of liquid penetrating a porous medium shouldscale in time as ( γt/µ ) / , where geometric factors andwettability constants were incorporated into an overallproportionality constant. Washburn thereby revealed theessential physical mechanism governing the spontaneouscreep of liquids into enclosed spaces, a process that caneasily dominate and oppose the force of gravity for suffi-ciently small enclosures. The Washburn scaling has sincebeen used successfully to quantify capillary transport inmany different systems ranging from blood flow in mi-crovascular hemodynamics [12, 13] to oil extraction fromporous rocks [14] to water uptake in dried tree specimens[15], to name a few. Through modern day advances inmicrofabrication techniques, capillary action is now be-ing incorporated into the design of many microfluidic de-vices as well, some involving spontaneous flow throughchemically treated porous substrates, others relying ona combination of capillary action, positive displacementpumping and electrophoresis. The number of such ap-plications is multiplying rapidly with emphasis on dis-posable inexpensive platforms beneficial to global publichealth [16], drug discovery screening [17] and specializedfluid based logic circuitry [18, 19].While Washburn originally focused on the fluid dynam-ics underlying internal capillary flow, there soon camethe realization that capillary forces were also somehowresponsible for the spontaneous exterior flow of liquidalong open surfaces manifesting grooves, edges, cornersand roughened patches. This type of free surface flow istypically characterized by a solid boundary whose sur-face pattern is fixed and a gas/liquid interface whoseshape undergoes spatial and temporal variation as theflow evolves toward steady state conditions. Free sur-face wicking is now also being incorporated into manysmall-scale fluidic devices such as heat pipes for cool-ing of microelectronics [20–23], surface-enhanced Ramanspectroscopic devices [24] and small spacecraft fluid man-agement systems [25]. II. BACKGROUND
In one of the earliest studies of free surface wicking of aNewtonian liquid, Concus and Finn [26] in 1969 showedthat a liquid drop will penetrate the interior of an open αy xh ( z,t ) q (b) r c ( z ) d Ly xz (a) FIG. 2. (Color online) (a) Schematic of a wetting liquid film(0 ≤ θ < π/
2) flowing within a slender open triangular groovewith constant cross-section. The inlet midline film thicknessat the origin is denoted by d = h ( x = 0 , z = 0 , t ) and thechannel length by L where d/L ≪
1. (b) Cross-sectional viewof the flow geometry depicted in (a) where h ( z, t ) denotes thelocal midline film thickness, α is the groove interior half-angle, r c is the radius of curvature of the liquid interface and θ isthe contact angle of the liquid wetting the channel sidewalls,which is assumed constant. triangular channel with opening half-angle α so long as θ + α < π/
2. Shortly thereafter, Ayyaswamy et al. [27]derived and numerically solved for the stationary (i.e.steady state) streamwise velocity field for this geometry.Closed-form analytical solutions for the free surface tran-sient capillary flow in idealized containers with grooves orcorners were demonstrated two decades later by severalgroups working independently. Dong and Chatzis [28],Weislogel and Lichter [29, 30] and Romero and Yost [31]derived and solved for the interface equation h ( z, t ) de-scribing the thickness of a liquid film flowing within anopen slender triangular groove satisfying the condition d ≪ L , as sketched in Fig. 2. The mechanism responsi-ble for spontaneous capillary flow within an open grooveturned out to be more complex than its counterpart forsmall enclosures. While for interior flow the liquid col-umn is pulled forward by the net force of surface tensionacting along the triple line (i.e. the line defining the junc-tion between the gas, liquid and solid media), the freeboundary in an open capillary driven system experiencesa normal force across the entire interface whose shapetypically varies in space and time. By implementingtwo key assumptions, namely that the streamwise vari-ation in liquid height varies smoothly on a length scale L much larger than d , and that the flow is slow enoughsuch that capillary forces dominate inertial and viscousforces, these researchers were able to show that the cross-section of the liquid surface must describe a circular arcof radius R ( z ). Hence the pressure gradient driving theflow stems exclusively from the streamwise variation in R ( z ). These two assumptions therefore led to a non-linear diffusion equation describing the relaxation of thecross-sectional liquid area proportional to h ( z, t ). WhileRomero and Yost [31] and Dong and Chatzis [28] as-sumed unidirectional flow as well, Weislogel and Lichter[29, 30] were able to derive similar results without re-lying on this simplification by instead conducting a for-mal perturbation expansion based on the slender param-eter approximation ( ε ) = ( d/L ) ≪
1. These separateefforts led to the conclusion that there exist transientsolutions according to which the penetration distance z scales in time as ( D t ) / where D = γ d K ( θ, α ) /µ , d is acharacteristic liquid thickness (such as the groove depthor midline liquid thickness) and K ( θ, α ) is a geometricfactor. Chen, Weislogel, and Nardin [32] later general-ized this result to channels with constant cross-sectionranging from the nearly-rectangular to highly roundedgrooves. In all cases, it was confirmed that the penetra-tion distance retains the t / scaling. Numerous experi-mental studies [25, 29, 33–36] and references therein havesince confirmed their predictions.A few studies have gone beyond solution of the flowbase states to explore their stability spectrum. Weis-logel [37] examined the linear stability of an infinite,initially quiescent column of liquid supported within aslender triangular groove of constant cross-section us-ing conventional modal analysis. Other researchers havesince examined the free surface collapse and gas ingestionin inertial-capillary flows for rectangular [38], triangular[39] or asymmetrical [40] grooves, the stability of liquidmenisci attached to a solid edge [41, 42], and the dry-out instability for thermocapillary (i.e. non-isothermal)driven flow in triangular grooves [43]. Aside from thesespecialized studies, however, the majority of theoreticalstudies have been devoted exclusively to solutions govern-ing unperturbed base states i.e. solutions in the absenceof disturbances.In this work, we address for the first time the stabil-ity of stationary and transient interfaces associated withcapillary flows in slender, open triangular grooves of con-stant cross-section subject to the Concus-Finn inequality[26] ensuring liquid penetration ( θ + α < π/
2) and the re-striction that the liquid always wet the channel sidewalls.In Section III, we outline the derivation of the governinginterface equation [28–31]. In Section IV, we examinethe nonlinear stability of stationary states subject eitherto Dirichlet, Neumann or constant volume conditions.By invoking the gradient flow form of the interface dis-turbance equation and an associated Lyapunov energyfunctional governing the decay in time of arbitrary distur-bances, we demonstrate that stationary interface shapesare exponentially stable and therefore extremely robustto perturbations. In Section V, we examine the general-ized non-modal linear stability of self-similar states obey- ing Washburn scaling [44, 45] using both analytic argu-ments and direct numerical computation. Our resultsindicate that self-similar solutions characterized by thescaling z ∼ t / are both transiently and asymptoticallystable to infinitesimal perturbations. III. MODEL FOR FREE SURFACE CAPILLARYFLOW IN SLENDER OPEN TRIANGULARGROOVES
In this section, we outline the theoretical model de-scribing the capillary motion of a non-volatile, isother-mal Newtonian liquid film flowing within an open slen-der triangular groove in contact with an ambient pas-sive gas [28–31], as depicted in Fig. 2(a). The in-let midline film thickness at the origin is denoted by d = h ( x = 0 , z = 0 , t ) and the channel length by L wherethe slender limit is assumed, namely d/L ≪
1. Shown inFig. 2(b) is a cross-sectional view of the flow geometrywhere h ( z, t ) denotes the local midline film thickness, α is the groove interior half-angle, r c is the radius of cur-vature of the liquid interface and θ is the contact angleof the liquid wetting the channel sidewalls. The liquidfilm is assumed to maintain a constant contact angle in-dependent of location or flow speed. A. Slender limit form of the hydrodynamicequations
In order to conserve mass and momentum, an incom-pressible Newtonian liquid of constant density must sat-isfy the continuity and Navier-Stokes equation given bythe set of coupled equations ∇ · ~u = 0 , (1) ρ (cid:20) ∂~u∂t + ( ~u · ∇ ) ~u (cid:21) = −∇ p + µ ∇ ~u + ρ~g, (2)where the velocity field in Cartesian coordinates is rep-resented by ~u = ( u, v, w ), the fluid pressure by p ( x, y, z ),the gravitational acceleration by ~g and the constant fluiddensity by ρ . Assuming flow in a slender channel suchthat ε = d/L ≪ Bo , capillary Ca and Reynolds Re numbers in theslender limit. To order ε , the rescaled forms of Eqs. (1)Quantity Scaling RescaledvariableSlender parameter ε = d/L ≪ x c = d X = x/x c y c = d Y = y/y c z c = L Z = z/z c Velocity u c = ε γ Ca /µ U = u/u c v c = ε γ Ca /µ V = v/v c w c = εγ Ca /µ W = w/w c Pressure p c = γ Ca / ( εL ) P = p/p c Time t c = µL/ ( εγ Ca ) T = t/t c τ = ln( T )Interface midline y c = d H = h ( z, t ) /y c thicknessInterface shape y c = d Σ = σ ( x, z, t ) /y c Interface radius y c = d R = r c /y c of curvatureStationary state H S ( Z )midline thicknessSelf-similar variable η = Z/ √ T Self-similar state S ( η )midline thicknessBond number Bo = ρgd /γ Capillary number Ca = µw c /εγ = ΦReynolds number Re = ρw c d/µ = ( εργd/µ ) Ca TABLE I. Characteristic scalings (lower case) and non-dimensional variables (uppercase) used to describe dimension-less system shown in Fig. 2. and (2) are then given by0 = ∂U∂X + ∂V∂Y + ∂W∂Z , (3a) ε Re DUDT = BoCa G x − ∂P∂X + ε ∆ U, (3b) ε Re DVDT = BoCa G y − ∂P∂Y + ε ∆ V, (3c) ε Re DWDT = 1 ε BoCa G z − ∂P∂Z + ∆ W, (3d) where the substantial derivative D/DT and the Lapla-cian derivative ∆, respectively, are given by
DDT = ∂∂T + U ∂∂X + V ∂∂Y + W ∂∂Z , (4a)∆ = ∂∂X + ∂∂Y + ε ∂∂Z , (4b)and ~G = ~g/g . In the limits where ε ≪ ε Re ≪ Bo / Ca ≪ ε , the governing equations reduce to the form: ∂U∂X + ∂V∂Y + ∂W∂Z = 0 , (5a) ∂P∂X = ∂P∂Y = 0 , (5b) ∂P∂Z = ∂ W∂X + ∂ W∂Y . (5c)When subject to the slender limit, the flow is thereforeinertia-free and the fluid pressure is constant throughoutthe ( x, y ) plane. The pressure gradient driving the flow,which will stem solely from capillary forces, can thereforeonly vary along the streamwise axis and can only be coun-terbalanced the viscous force set in play by the no-slipboundary condition applied along the groove sidewalls,namely U = V = W = 0 at all liquid/solid interfaces. B. Boundary conditions at the liquid interface
The two (dimensional) boundary conditions specify-ing the jump in normal and shear stresses across thegas/liquid interface σ ( x, z, t ) are given by[ˆ n · ( e τ − p I ) · ˆ n + γ ( ∇ s · ˆ n )] y = σ ( x,z,t ) = 0 , (6a) (cid:2) ˆ t i =1 , · e τ · ˆ n (cid:3) y = σ ( x,z,t ) = 0 , (6b)where I is the 3 × e τ = µ [ ∇ ~u +( ∇ ~u ) T ] de-notes the shear stress tensor, ∇ s = ( ∇− ˆ n (ˆ n ·∇ )) denotesthe surface gradient operator and the triad (ˆ n, ˆ t , ˆ t ) de-notes the three unit vectors representing directions nor-mal and tangent to the moving interface with the conven-tion that ˆ n points away from the interface. In rescaledunits, these unit vectors are given byˆ N = 1[1 + ( ∂ X Σ) + ε ( ∂ Z Σ) ] / − ∂ X Σ1 − ε∂ Z Σ , (7a)ˆ T = 1[1 + ( ∂ X Σ) ] / ∂ X Σ0 , (7b)ˆ T = 1[1 + ε ( ∂ Z Σ) ] / ε∂ Z Σ1 , (7c)where Σ( X, Z, T ) = σ ( x, z, t ) /d denotes the non-dimensional interface function and subscripts denote dif-ferentiation with regard to the rescaled coordinates.Specifying a system for which the fluid pressure derivessolely from variations in the local interface curvature ofthe flowing liquid, the interfacial surface tension γ is ev-erywhere constant since the liquid is isothermal and con-tains no surfactant-like additives, and that the liquid re-mains in contact with a passive quiescent gas of negligibleviscosity and density with gauge pressure set to zero, thejump in normal stress is then strictly due to capillaryforces and the liquid interface is a surface of vanishingshear stress. To order O ( ε ) then, these dimensionlessboundary conditions reduce to the form0 = − P − Ca − (cid:18) ∂ X Σ[1 + ( ∂ X Σ) ] / (cid:19) + O ( ε ) (8a)= − P − Ca − K ( Z, T ) + O ( ε ) , (8b)0 = ∂ Y W − ( ∂ X Σ) ∂ X W + O ( ε ) , (8c)0 = (cid:2) − ( ∂ X Σ) (cid:3) ( ∂ Y U + ∂ X V )+ 2( ∂ X Σ)( − ∂ X U + ∂ Y V ) − ( ∂ Z Σ) [ ∂ X W + ( ∂ X Σ) ∂ Y W ] + O ( ε ) , (8d)where K(Z) represents the local curvature of the interfacefunction Σ (defined to be positive for a wetting liquid). Itis clear from Eq. (8b) that the curvature function K canin general only depend on ( Z, T ) since according to Eq.(5b), the pressure P is independent of ( X, Y ). This thenrequires that the cross-sectional shape of the liquid in-terface be described by a curve with constant curvature.This restriction limits the shape either to a flat interfaceor one described by a segment of a circle. Since the liq-uid must also satisfy a prescribed contact angle set bythe particulars of the liquid/solid interaction, a flat pro-file is disallowed and Σ(
Z, T ) must therefore trace out acircular arc of constant curvature. The non-dimensionalradius of curvature of the gas/liquid interface is thengiven by R ( Z, T ) = H ( Z, T ) sin α/ (cos θ − sin α ), or like-wise, the interface curvature is described by K ( Z, T ) =1 /R ( Z, T ) = (cos θ − sin α )(csc α ) H − ( Z, T ). (This re-lation differs slightly from that originally derived byRomero and Yost [31] where they adopted a sign con-vention for K opposite to ours and chose the referenceliquid thickness to be the height of the fluid intersectingthe groove wall, not the inlet midline film thickness.) Ac-cording to Eq. (8a), the capillary pressure is then givenby P ( Z, T ) = − Ca − b R ( α, θ ) H ( Z, T ) , (9)where b R ( α, θ ) = sin α cos θ − sin α . (10)The Concus-Finn condition α + θ < π/ b R ( α, θ ) >
0, or likewise P ( Z, T ) <
0, consis-tent with a liquid interface with positive curvature. The case α + θ > π/ C. Interface midline equation H ( Z, T ) for capillaryflow in slender open triangular grooves The fact that the interface shape can only be a segmentof a circle, and is therefore independent of the local coor-dinates (
X, Y ), leads to simplification of the expressionfor the streamwise volumetric flux Q ( Z, T ). The relevantvariables are then scaled by H according to e X = XH and e Y = YH , (11) f W ( e X, e Y ) = (cid:18) − Ca − b R ∂H∂Z (cid:19) − W, (12) e Σ( e X ) = 1 + b R − q b R − e X . (13)This rescaling allows for solution of f W independently ofthe local value of H . As a result, the geometric func-tion Γ( α, θ ) = RR f
W d e Xd e Y need only be computed nu-merically once. The dimensionless streamwise flux whichtraverses the local cross-sectional area A can then be re-expressed as Q ( Z, T ) =
Z Z A W dXdY (14)= − Ca − b R ( α, θ ) Γ( α, θ ) H ∂H∂Z , (15)where the integrated area A = b A ( α, θ ) H and b A = cos θ sin α cos( α + θ ) − ( π/ − α − θ ) sin α (cos θ − sin α ) (16)is a geometric factor [29, 31]. Since the local stream-wise gradient in liquid flux is directly related to thetime derivative of the local liquid cross-sectional area(see Appendix 2 of Ref. 46 for derivation) according to ∂A/∂T = − ∂Q/∂Z , the governing nonlinear diffusionequation for the midline height H ( Z, T ) is then given by b A ( α, θ ) ∂H ∂T = Ca − b R ( α, θ ) Γ( α, θ ) ∂∂Z (cid:18) H ∂H∂Z (cid:19) . (17)Without loss of generality and to recast this equationinto parameter-free form, the remaining scaling for thestreamwise velocity is chosen to be w c = ( εγ/µ )Φ, whereΦ( α, θ ) = Ca (see Table 1) is defined byΦ( α, θ ) = Γ( α, θ ) b A ( α, θ ) b R ( α, θ ) . (18)The resulting interface equation, whose stability proper-ties we examine in this work, is then given by ∂H ∂T − ∂∂Z (cid:18) H ∂H∂Z (cid:19) = 0 , (19)subject to the constraint that H ( Z, T ) is everywhere al-ways positive. In Fig. 3 are plotted the scaled functions b A ( α, θ ), b R ( α, θ ), Γ( α, θ ) and Φ( α, θ ) for four values of theliquid contact angle θ = 0 ◦ , 20 ◦ , 40 ◦ and 60 ◦ as a functionof increasing groove interior half-angle α . While b A ( α, θ ), b R ( α, θ ) and Γ( α, θ ) are of order O (10) or less, the valuesof Φ( α, θ ) are far smaller and tend toward O (10 − ) orless. Note too that since b R ( α, θ ) and Γ( α, θ ) are non-negative functions for systems obeying the Concus-Finncondition, the direction of the liquid flux specified by Eq.(14) is then strictly determined by the sign of the localinterface slope ∂H/∂Z . Interfaces with ∂H/∂Z < Q and vice versa. A vanishinglocal flux results whenever ∂H/∂Z = 0.The form of Eq. (19) falls within a class knownas the porous media equation generally given by ∂C ( Z, T ) /∂T = ∂ C m /∂Z , where C ( Z, T ) is a non-negative scalar function and m is a constant larger thanone [47–50]. As discussed in Ref. [50], this nonlineardiffusion equation describes the relaxation of the orderparameter C ( Z, T ) relevant to various phenomena whicharise in different branches of science and mathematicsand exhibiting properties such as scale invariance andself-similarity. To help track the energy flow associatedwith the evolution of C , Newman [47] outlined a gen-eral method for constructing Lyapunov functionals forsystems sustaining traveling wave solutions. He usedthat method to establish the actual rate of convergenceof an initial configuration to solutions exhibiting self-similarity. Ralston [48] complemented this work by show-ing that Newman’s choice of Lyapunov function allowedproof that initial conditions with finite mass converge toself-similar solutions asymptotically in time. When ap-plied to our system, these results indicate that self-similarstates which result from the spreading of an initial finitedrop within a slender open triangular groove are asymp-totically globally stable. However, despite longstandinginterest in the porous media equation and wicking phe-nomena in general, there has been no prior study of whichwe are aware of the general stability of stationary ortransient solutions corresponding to unconstrained vol-ume states. IV. STABILITY OF STATIONARY INTERFACESOLUTIONS H S ( Z ) In this section, we first derive the analytic form cor-responding to stationary solutions of the governing in-terface equation given by Eq. (19) with geometric fac-tors explicit. By invoking analogy to the porous me-dia equation and constructing an appropriate Lyapunov ˚ ˚ ˚ ˚ α ( degrees ) AR Γ θ = ˚ ˚ ˚ ˚ θ = ˚ ˚ ˚ ˚ θ = ˚ ˚ ˚ ˚ θ = Φ FIG. 3. (Color online) Geometric functions pertinent to cap-illary flow of a Newtonian liquid film with constant contactangle θ in a slender open triangular groove with half openingangle α satisfying the Concus-Finn condition α + θ < π/ b A ( θ, α ), (b) b R ( θ, α ), (c) Γ( θ, α )and (d) Φ( θ, α ) described in the text. functional, we then examine the nonlinear asymptoticstability of these states against arbitrary disturbances. A. Stationary states for time-independentDirichlet, Neumann and volume conditions
Stationary states of Eq. (19) correspond to those so-lutions for which ∂H /∂T = 0, which reduces the gov-erning equation to a second order, ordinary differentialequation requiring two boundary conditions. So long asthe interface slope does not vanish, such stationary in-terfaces H S ( Z ) are possible because the subsurface flowof liquid establishes a balance between the local capillaryand viscous stresses generates a constant flux Q o , whichaccording to Eq. (17) is given by Q o = − Ca − b R ( α, θ ) Γ( α, θ ) H S dH S dZ (20)= − b A ( α, θ ) H S dH S dZ . (21)The general form of these stationary solutions is thenrepresented by H S = " const − Q o Z b A ( α, θ ) / > Z / for positiveflux and power law increase for negative flux. Particularsolutions, of course, require specification of two bound-ary conditions for setting the values of const and Q o .For a triangular channel with fixed groove opening angleand liquid contact angle (i.e. constant value of b A ( α, θ ))extending between endpoints Z and Z , stationary so-lutions H S correspond to H S = " H − Q S b A ( α, θ ) ( Z − Z ) / (23)for Dirichlet and Neumann conditions H S ( Z ) = H and Q = Q S , respectively. Alternatively, Dirichlet condi-tions imposed at both endpoints, such that H S ( Z ) = H and H S ( Z ) = H , yield H S ( Z ) = (cid:20) H Z − Z Z − Z + H Z − ZZ − Z (cid:21) / . (24)Note from Eq. (9) and the scaling for the streamwise flowspeed w c that specification of the boundary film thick-ness is equivalent to specification of the boundary fluidpressure. Solutions H S to Eq. (22) can also be generatedsubject to constant flux Q S at one boundary ( Q = Q S or Q = Q S ) and conservation of volume V S , which setsthe constant value C S in the implicit relation : V S = b A Z Z Z H ( Z ) dZ (25)= 3 / b A / Q S h ( C S − Q S Z ) / − ( C S − Q S Z ) / i . (26)Likewise, constant volume and fixed liquid height at oneendpoint yield similar forms.Representative solutions H S ( Z ) are plotted in Fig. 4 Z = 0.01= 0.33= 0.67= 1.00 H S ( Z = 3 ) = 1.33 H S FIG. 4. (Color online) Representative stationary solutionssubject to Dirichlet conditions H S ( Z ) = H = 1 . H S ( Z ) = H = 0.01, 0.33, 0.67, 1.00 and 1.33 for the range0 ≤ ( Z , Z ) ≤ . for Dirichlet conditions which pin the inlet height to H = 1 . H to the five val-ues shown. From the expression for the flux given by Eq.(20), it is evident that the solution with H = 1 .
33, whichexhibits interface slope values dH/dZ which are every-where positive, describes a stationary solution with netstreamwise flux Q S < H ( Z ) = 1 . Z / . B. Nonlinear exponential stability of stationarysolutions
We next examine the time-asymptotic nonlinear sta-bility of stationary solutions by appealing to a Lyapunovanalysis. In particular, we construct a Lyapunov energyfunction (akin to a potential energy function in classi-cal mechanics) to determine whether and how rapidlyarbitrary disturbances equilibrate back to the stationarysolutions H S . A similar approach to asymptotic stabilityhas previously been used to characterize disturbance de-cay governed by nonlinear diffusion dynamics [51]. Herewe show that the capillary flow in a bounded domain withtime-independent boundary conditions has a uniquely de-termined stationary solution which is dynamically stable.This implies that any initial distribution which main-tains the condition H ( Z, T ) > globally stable - that is, any initialdistribution H ( Z, > H ( Z, T ) > ∂G ( Z, T ) ∂T − ∂∂Z (cid:26) M ( G ) ∂∂Z (cid:18) δ F δG (cid:19)(cid:27) = 0 , (27)where M ( G ) is the so-called mobility function and δ F /δG denotes the functional derivative of F ( H ) with respectto G . The quantity F ( H ) represents the Lyapunov freeenergy F which for stable stationary states must satisfythe relations δ F = 0 (i.e. extremal condition) and δ F > G ( H ) = H ( Z, T ) and letting M = 1 /
3, Eq.(19) is re-expressed as ∂G∂T = ∂∂Z (cid:20) ∂R∂Z (cid:21) (28)where the function R ( G ) is given by R ( G ) = ( G / − G / S ) = ( H − H S ) . (29)(The mapping from H to G is one-to-one since H isstrictly positive.) Comparison with Eq. (27) then yieldsthe correspondence δ F ( G ) /δG = R ( G ) subject to theconstraint that [ R∂R/∂Z ] Z ,Z = 0. This definition of R ( G ) ensures that stationary states are identified by theextrema of the free energy F ( H ) where δ F ( G ) /δG | G S = R ( G S ) = 0. Additionally, Dirichlet boundary conditionscorrespond to R | Z ,Z = 0 and Neumann boundary con-ditions to ( ∂R/∂Z ) | Z ,Z = 0. Multiplication of bothsides of Eq. (28) by R ( G ) followed by integration overthe finite domain [ Z , Z ] yields the relation Z Z Z R ∂G∂T dZ = Z Z Z R ∂∂Z (cid:20) ∂R∂Z (cid:21) dZ = (cid:20) R ∂R∂Z (cid:21) Z Z − Z Z Z (cid:20) ∂R∂Z (cid:21) dZ = − Z Z Z (cid:20) ∂R∂Z (cid:21) dZ ≤ . (30) The left term in Eq. (28) can also be expressed as Z Z Z R ∂G∂T dZ = Z Z Z ( G / − G / S ) ∂G∂T dZ = Z Z Z ∂∂T (cid:18) G / − GG / S (cid:19) dZ = ∂∂T Z Z Z (cid:18) G / − G / S G + 35 G / S (cid:19) dZ (31)where the last term in Eq. (31) has been added, withoutloss of generality, to ensure the integrand vanishes for G = G S , which essentially sets the value of the groundstate energy. Given that δ F /δG = R [ G ( Z )], it is evidentthat the integral defined in Eq. (31) is none other thanthe Lyapunov free energy, which can be re-expressed interms of the interface function H : F ( H ) = Z Z Z (cid:18) H − H S H + 35 H S (cid:19) dZ. (32)Equating Eq. (30) and Eq. (31) then yields the relation ∂ F ( H ) ∂T = − Z Z Z (cid:20) ∂∂Z ( H − H S ) (cid:21) dZ ≤ , (33)where the final equality holds only for initial states H ( Z, T ) exactly equal to the stationary state solutions H S ( Z ). This finding serves as proof that the flow de-scribed is asymptotically stable in the Lyapunov sensei.e. arbitrary initial distributions H ( Z, T ) decay in timetoward the stationary solution H S , characterized by avanishing first variation and positive second variation,namely δ F ( G ) | G S = Z Z Z δ F δG (cid:12)(cid:12)(cid:12)(cid:12) G S δG dZ (34)= Z Z Z R ( G ) | G S δG dZ = 0 and (35) δ F ( G ) (cid:12)(cid:12) G S = Z Z Z δ F δG (cid:12)(cid:12)(cid:12)(cid:12) G S ( δG ) dZ (36)= Z Z Z G / S ( δG ) dZ > . (37)To complete this part of the proof, we show that thestable stationary states H S are in fact also unique. Sup-pose then that in addition to the already specified steadystate solution H S there exists another steady solution H S satisfying the same boundary conditions. Since H S is time independent, it must satisfy the relation ∂ F ( H S ) /∂T = − (1 / R Z Z (cid:2) ∂ ( H S − H S ) /∂Z (cid:3) dZ = 0.Clearly this is only possible if H S = ( H S + C o ) / , wherethe constant C o must equal zero in order for H S tosatisfy the same boundary conditions as H S . This re-sult therefore establishes that the solution H S is indeedunique.The proof above requires that both the stationary H S ( Z ) and transient H ( Z, T ) solutions be everywherestrictly positive. For solutions satisfying Dirichlet bound-ary conditions, these requirements are easily met. Con-sider any positive initial state H ( Z, T = 0) which redis-tributes its height in time according to Eq. (41). Werethere a point within the domain where H ( Z, T ) →
0, thenthe local interface would have to satisfy ∂H/∂T < ∂ H/∂Z would have to become suffi-ciently negative to satisfy the balance of terms Eq. (41).The local interface would then have to develop a localprotrusion (negative curvature) and not a local dimpleleading to rupture (positive curvature), in contradictionto the assumption of positive curvature required at thatlocal minimum. Therefore, all positive initial states re-main positive in time. For Dirichlet conditions then, sta-tionary solutions H S ( Z ) are globally exponentially sta-ble. For the remainder stationary solutions subject to aNeumann boundary condition, it may become the casethat too high a flux condition leads to one of more pointswhere the local film thickness vanishes. In the vicinity ofsuch points, the local interface slope will become increas-ingly large and eventually violate the slender limit ap-proximation. Perhaps more importantly, disjoining pres-sure effects, known to modify the flow in ultrathin regionsof a liquid film, must then be incorporated into the modelinterface equation [54, 55]. In such cases then, the proofabove only establishes asymptotic stability, not globalasymptotic stability.Next, we demonstrate an even stronger statement re-garding the asymptotic stability of stationary states,namely that the stationary solutions H S represent expo-nentially stable equilibria of Eq. (19). A useful relationfor the function R ( Z ) defined in Eq. (29) can be obtainedby first applying the Poincar´e-Friedrichs inequality in onedimension [56] according to which Z Z Z (cid:18) dRdZ (cid:19) dZ ≥ Z − Z ) Z Z Z R ( Z ) dZ (38)subject to the constraint that there is some interior point Z ∗ in the bounded domain Z ≤ Z ∗ ≤ Z where R ( Z ∗ ) =0. This must always be the case, however, since H ( Z, T )must satisfy the same boundary conditions as H S . Giventhat H ( Z, T ) and H S ( Z ) are strictly positive throughoutthe domain, it is possible to construct an upper boundon ∂ F /∂T in Eq. (33) as follows: ∂ F ( H ) ∂T = − Z Z Z (cid:20) ∂∂Z ( H − H S ) (cid:21) dZ ≤ − Z − Z ) Z Z Z ( H − H S ) dZ< − C ( Z − Z ) Z Z Z ( H − H S ) × (cid:26) H + 4 H H S + 6 HH S + 3 H S H + HH S + H S ) (cid:27) dZ = − C ( Z − Z ) Z Z Z (cid:18) H − H S H + 35 H S (cid:19) dZ = − F ( H )3 C ( Z − Z ) , where C , which is strictly positive, is defined as C = min T,Z (cid:26) H + 4 H H S + 6 HH S + 3 H S H + HH S + H S ) (cid:27) . (39)(See Appendix A for further discussion regarding theboundedness property of C ). The final inequality soughtis then given by F ( H, T ) F ( H, T = 0) < exp (cid:26) − T C ( Z − Z ) (cid:27) , (40)which confirms that the Lyapunov free energy of anydisturbed state H ( Z, T ) satisfying the same boundaryconditions as the initial stationary state will decay backto the stationary state H S at least as fast as an expo-nential with a decay rate that scales with the square ofthe spatial domain size. This proof applies to all cate-gories of stationary solutions discussed in Section IV A.For disturbance functions subject either to two Dirichletconditions or one Dirichlet condition and one Neumann(constant flux) or constant volume condition, the proof istrivial since either R ( Z ) or R ( Z ) vanishes identically.For solutions H S that correspond to fixed flux Q S atone boundary and constant volume V S , the proof abovealso applies since there always exists an interior point Z ≤ Z ∗ ≤ Z where R ( Z ∗ ) = 0. This follows becausetwo solutions H and H S cannot have the same constantvolume unless the functions H ( Z, T ) and H S ( Z ) undergoat least one crossing point within the domain. V. STABILITY OF SELF-SIMILAR SOLUTIONS
We next seek unperturbed (i.e. base state) transientsolutions to Eq. (19) for non-conserved volume whichmanifest self-similarity. A global stability argument aspresented in Section IV B for stationary solutions provedunsuccessful since the self-similar form of the governinginterface equation cannot be converted into gradient flowform. Therefore we instead applied a generalized lin-0ear stability analysis [45] which helps determine whetherinfinitesimal non-modal disturbances undergo any tran-sient or asymptotic amplification.
A. Self-similar solutions with time independentDirichlet conditions
Previous studies have delineated the conditions lead-ing to existence and uniqueness of self-similar solutionsas well as the attraction of spatially confined initial dis-tributions toward self-similar base states [50]. Here wefocus on volume non-conserving positive states S ( η ) con-sistent with time-independent Dirichlet boundary condi-tions imposed at the domain endpoints, namely S (0) = 1and S ( η → η B ) = const where η B denotes a location fardownstream of the origin. The Dirichlet condition at theorigin can be set to unity without loss of generality sinceas evident from Eq. (46), a rescaling involving a mul-tiplicative factor of S (0) leaves the governing equationunchanged.In general, for self-similarity to hold, there can beno intrinsic length or time scale imposed on the flow,in contrast to the steady state solutions examined inthe previous section which depend on the groove length Z − Z . A simple scaling analysis of Eq. (19) re-veals that self-similar solutions may be possible when-ever T << L /H ∼ O ( L/ε ). To find such solutions, it isconvenient to expand and rewrite Eq. (19) in the form ∂H∂T − H ∂ H∂Z − (cid:18) ∂H∂Z (cid:19) = 0 . (41)The ansatz H sim ( η, T ) defined by H sim ( η ) = T β − S ( η ) where (42) η = ZT β for β > , (43)allows for a large class of self-similar solutions [30, 50]satisfying the general second order nonlinear differentialequation S S ηη + ( S η ) + βηS η + (1 − β ) S = 0 . (44)Inspection of the asymptotic behavior of S ( η ) as η → ∞ helps ascertain what range of exponents β are requiredfor bounded non-terminating (i.e. S >
0) states suchthat S ηη and S η asymptotically approach zero as η → ∞ .While the first two terms on the left side of Eq. (44) thenvanish identically, care must be taken with regard to thethird term which couples an increasingly large value of η with a diminishingly small term S η . Balancing the thirdand fourth terms yields the proper asymptotic scaling,namely dS/S ∼ [(2 β − /β ] dη/η , and hence S ( η → ∞ ) ∼ η (2 β − /β . Therefore, only the range 0 < β ≤ / β . For example, en-forcement of constant liquid volume V ≃ R Z Z H dZ = T β − R Z Z S dη is only consistent with β = 2 /
5. Accord-ing to Eq. (20), enforcement of a constant flux bound-ary condition Q = − b A ( α/θ ) H ( ∂H/∂Z ) = const = − b A ( α/θ ) T β − S S η is only consistent with β = 3 / β = 1 /
2, which accords with the Washburn relation andallows enforcement of time-independent Dirichlet bound-ary conditions. For this category of solutions, the non-dimensional flux defined in Eq. (14) is represented by Q ( η, T ) = − b A ( α, θ ) T / S ∂S∂η , (45)The self-similar solution S ( η ) then satisfies the equation[7]: SS ηη + ηS η + 2( S η ) = 0 . (46)To ascertain the interface shape of these solutions, wenumerically solved Eq. (46) by rewriting the second orderequation as a system of first order equations and usingthe ODE45 solver in Matlab [57]. Shown in Fig. 5 arerepresentative solutions for receding, uniform, advanc-ing and terminating states S ( η ) satisfying the far fieldDirichlet conditions shown. According to Eq. (45), so-lutions with S η > receding states , while solutionswith S η < advancing states. The solution for which Sη vanishes everywhere, which corresponds to a zero flux so-lution in self-similar coordinates, is designed a uniform state. It represents an exception in that it is the onlysolution which satisfies volume conservation. Solutionswhose advancing front are characterized by a vanishingvalue of S ( η ) are likewise designated terminating states.The numerical solutions indicate that solutions undergotermination only when the interface slope at the origin S η (0) & . B. Generalized non-modal linear stability ofinterface self-similar solutions
It is by now well understood that traditional modalstability analysis, wherein instability derives from ex-ponentially growing normal modes of the governing lin-earized autonomous disturbance operator, is an inappro-priate investigatory tool for examining transient stabilityof systems governed by non-normal operators [45]. Thisis because the modal spectrum, and in particular theeigenvalue with maximum real part, properly capturesthe system response to infinitesimal disturbances at all1 S ( η ) η
12 3 4 5 6 7 891011121314151617
Receding Advancing Uniform T e r m i n a t i n g FIG. 5. (Color online) Representative self-similar solutions S ( η ) for terminating (1-7), advancing (8-13), uniform (14)and receding (15-17) states. The computational domain usedin numerically solving for these solutions was [0 ≤ η ≤ ≤ η ≤ .
0] is shown in the figure since thedownstream behavior remains essentially unchanged beyondthat value.) All solutions satisfy the Dirichlet condition S ( η =0) = 1 at the origin. Solutions 1-7 exhibit interface slopes atthe origin given by S η (0) =-10, -2, -0.8, -0.5, -0.4, -0.36, and-0.3492. Solutions 8-13 satisfy far field Dirichlet conditionsgiven respectively by S (80) = 0.01, 0.17, 0.33, 0.5, 0.67 and0.83. (Resulting values for the interface slopes at the originare S η (0) = -0.3491, -0.3418, -0.3185, -0.2745, -0.2090, and-0.1185, respectively.) The line denoted by 14 represents auniform solution where S ( η ) = 1 .
0. Solutions 15-17 satisfy farfield Dirichlet conditions given respectively by S (80) = 1.17,1.33 and 1.5. (Resulting values for the interface slopes at theorigin are S η (0) = 0.1485, 0.3283 and 0.5451, respectively.) times only for linearized disturbance operators which arenormal; for disturbance operators which are non-normal,it describes the response only at infinite time. The case ofnon-normal operators therefore requires implementationof a generalized (non-modal) or transient growth stabil-ity analysis, as discussed in Refs. [44 and 45] and refer-ences therein. The important distinction is that while fornormal operators instability arises from a single, fastestgrowing, exponentially unstable eigenmode which dom-inates at all times, instability in systems governed bynon-normal operators derives from non-modal growth oftwo or more interacting nonorthogonal eigenmodes at fi-nite times while only asymptotically resolving to the dis-turbance function characterized by the eigenvalue withmaximum positive real part.Equations of motion which give rise to inhomogeneousbase states often yield linearized operators which are non-normal. Such non-normal operators commonly arise inmany thin film liquid systems [58, 59]. Given that thebase state solutions to Eq. (46) include advancing andreceding states which vary with the self-similar variable,we therefore appeal to a transient stability analysis toinvestigate their behavior. The stability of terminatingsolutions shown in Fig. 5, however, will not be examinedin this work due to the fact that the point of termination requires that additional terms be included in the gov-erning equation of motion to relieve the diverging stressknown to occur at a moving contact line [60–62].To begin, we linearize the general solution accordingto H ( η, τ ) = S ( η ) + δG ( η, τ ) , (47)where S ( η ) satisfies Eq. (46) with the given boundaryconditions, G ( η, τ ) represents the non-modal disturbancefunction and δ ≪ δ , namely: ∂G∂τ = L [ G ] , (48)where L = S ∂ ∂η + (cid:18) dSdη + η (cid:19) ∂∂η + 12 dSdη , L † = S ∂ ∂η − (cid:18) η dSdη (cid:19) ∂∂η − (cid:18)
12 + d Sdη (cid:19) . (49)Here, L † denotes the adjoint of the linear autonomousoperator L , τ = ln T and G ( η, τ ) = exp( L τ ) G ( η, L (Euclidean) norm, the adjoint L † representsthe unique linear operator where R η B v ( η ) L { w ( η ) } dη = R η B L † { v ( η ) } w ( η ) dη for all sufficiently smooth functions v and w which are L -integrable and satisfy homogeneousDirichlet boundary conditions at the boundary points(0 , η B ). For all non-uniform solutions S , LL † = L † L sothat L is non-normal. The transient amplification of dis-turbances G ( η, τ ) cannot therefore simply be ascertainedfrom the eigenspectrum of L . Disturbance amplificationover an interval in time τ will instead be quantified bythe function σ = k G ( η, τ ) kk G ( η, k = h G ( η, τ ) | G ( η, τ ) i / h G ( η, | G ( η, i / (50)= (cid:10) exp( L † τ ) exp( L τ ) G ( η, | G ( η, (cid:11) / h G ( η, | G ( η, i / (51) ≤ k exp( L τ ) k , (52)where the notation k · k represents the Euclidean norm h· | ·i / when applied to functions and the operator normwhen applied to operators, as in Eq. (52).To gain further insight into the growth or decay ofdisturbances at finite times, it is useful to consider thediscretized operator exp( L τ ) in diagonalized form suchthat exp( L τ ) = S exp(Λ τ ) S − (53)where exp(Λ τ ) denotes a diagonal matrix whose entriesare the eigenvalues of L arranged in decreasing order and2 S represents the matrix whose columns are the corre-sponding eigenfunctions. It can then be shown [45] thatexp( β max τ ) ≤ k exp( L τ ) k = kS exp(Λ τ ) S − k (54) ≤ kSkkSk − exp( β max τ ) , (55)where β max is the eigenvalue of L with maximum realpart, also known as the spectral abscissa of L . A matrixoperator is normal if and only if it is unitarily diago-nalizable such that k S kk S k − = 1, in which case themaximum amplification at time τ over all initial condi-tions G ( η,
0) is given by exp( β max τ ). For non-normaloperators, however, k S kk S k − ≥ τ can exceed theasymptotic value exp( β max τ ), sometimes significantly so.We note too from Eq. (52) that as τ → ∞ , the quan-tity d ln( k exp( L τ ) k ) /dτ → β max . While a finding that β max < asymptotically lin-early stable, it does not preclude the possibility of largepositive transient growth (i.e. linear instability) at finitetime.We now determine the upper bound on the instanta-neous rate of disturbance growth at time τ given by ω = 1 σ ∂σ∂τ = ∂ ln σ∂τ (56)= 1 k G ( η, τ ) k ∂∂τ h G ( η, τ ) | G ( η, τ ) i / (57)= hL G ( η, τ ) | G ( η, τ ) i + h G ( η, τ ) |L G ( η, τ ) i k G ( η, τ ) k (58)= 1 k G ( η, τ ) k (cid:28) L + L † G ( η, τ ) | G ( η, τ ) (cid:29) (59) ≤ max (cid:26) λ (cid:18) L + L † (cid:19)(cid:27) (60) ≡ ω max . (61)The quantity ω max , also known as the numerical abscissa,represents the maximum eigenvalue of the operator sum,( L + L † ) /
2. Since this combined operator is self-adjoint,its eigenvalues are strictly real. Therefore, while β max —the eigenvalue of L with maximum real part—governs theasymptotic stability as τ → ∞ , it is the value max { λ } which governs the transient stability of the linearized sys-tem at finite times. Since physical systems undergo flowinstabilities at finite time, it is this latter value which is ofmost interest. According to this definition then, if it canbe shown that ω max < τ , then instanta-neous disturbance growth will always be suppressed andthe system can be deemed linearly stable to infinitesi-mal perturbations. Below we establish the generalizedlinear stability of self-similar solutions S ( η ) > d confined to the interior corner of an open triangu-lar channel. He showed that axial sinusoidal perturba-tions of any wavelength ℓ all undergo rapid decay witha dimensional time constant proportional to µℓ /γd . Inour study, the governing dimensionless linear operatorcorresponding to a stationary uniform state is given by L = ( d/ ∂ /∂Z , an operator which is self-adjointand therefore normal, whose eigenvalues are strictly real,positive and given by ( d/ nπ/L ) for n = 1, 2, etc.According to Eq. 48 (when expressed in the originallaboratory frame coordinates ( Z, T )), the least stabledisturbance δG ( Z, T ) will evolve in time according toexp[ − ( d/ π/L ) T ] and therefore decay away, restoringthe system back to the initial uniform state. This ob-servation confirms Weislogel’s result that uniform quies-cent liquid filaments confined to a slender open triangularchannel are linearly stable to arbitrary infinitesimal dis-turbances of any wavenumber. This result holds for anyof the four sets of boundary conditions discussed earlier,namely two Dirichlet conditions, one Dirichlet and oneNeumann condition, one Dirichlet condition and constantvolume or one Neumann condition and constant volume. C. Generalized stability of volume non-conservingself-similar solutions
In what follows, we examine the generalized linear sta-bility of self-similar non-terminating states on the finitedomain 0 ≤ η ≤ η B . To proceed, we first note that thelinear autonomous operator ( L + L † ) / L + L † ∂∂η (cid:18) S ( η )2 ∂∂η (cid:19) − (cid:18) ∂ S∂η (cid:19) (62)= S ∂ ∂η + 12 ∂S∂η ∂∂η − (cid:18) ∂ S∂η (cid:19) , (63)satisfies a Sturm-Liouville eigenvalue equation where (cid:18) L + L † (cid:19) G ( η, τ ) = − λ G ( η, τ ) . (64)provided S ( η ) > S , dS/dη and d S/dS are con-tinuous within the interval 0 ≤ η ≤ ∞ . As a result,the eigenvalues λ are strictly real and the correspondingeigenfunctions form a complete orthogonal set of basisfunctions. According to the definitions in Eq. (61) then,if it can be shown that all the eigenvalues of L + L † / ω max < S ( η ) = 1 is easy to confirm. From thedefinitions given by Eqs. (63) and (64), the relevant op-erator sum reduces to (1 /
2) ( ∂ /∂η ) − (1 / / nπ/L ) +1 / n = 1, 2, etc. According to Eq. (64), the leaststable disturbance δG ( η, τ ) will evolve in time accord-ing to exp − (1 / nπ/L ) + 1 / T ] τ and therefore decayaway, restoring the system back to the initial uniformself-similar state.It is often the case that the Rayleigh quotient [63] canbe used to establish bounds on growth or decay by help-ing determine the overall sign of eigenvalues correspond-ing to Sturm-Liouville operators. In our case, however,this approach yields no useful information. We thereforeinstead appeal to the following lemma: Lemma 1.
Given a second order Sturm-Liouville eigen-value equation of the form ( P G η ) η + Q G = − λG where P > , G = 0 , λ is real, and G (0) = G ( η B ) = 0 , then λ min > ( − Q ) min , where subscripts denote differentiation with respect tothe domain coordinate η . In order for G to satisfy ho-mogeneous conditions at the boundary points, it must bethe case that there exists an extremum within the domain[0 , η B ] where G η ( η ∗ ) = 0. Without loss of generality, thefirst such extremum from the origin may be assumed tobe a local maximum in G for G >
0. (Note that G canalways be made positive at this point by multiplying thelocal value by -1 and still remain an eigenfunction). Be-cause P G η must change sign on either side of this max-imum, then at some point in its vicinity, ( P G η ) η < λ > − Q ( η ∗ ), which therefore im-plies λ min > − Q ( η ∗ ) ≥ ( − Q ) min , where ( − Q ) min denotesthe smallest value of − Q within the domain.
1. Transient and asymptotic stability of advancingsolutions S ( η ) We first examine the stability of advancing self-similarsolutions shown in Fig. 5. These solutions are all charac-terized by a downstream boundary value smaller than theinlet value, namely S ( η B ) = const < S (0) = 1, and aninlet slope which is always negative, namely S η (0) < η B represents the coordinatewhere the downstream boundary condition is applied.Establishing that such states are linearly stable re-quires a finding that the smallest eigenvalue of ( L + L † ) / λ min , is positive, which then requires λ min > ( − Q ) min = (1 + S min ηη ) / > , (65)which in turn requires that S min ηη > −
1. This minimumvalue can occur either at some interior point η int or at the boundary points η = 0 , η B , which shall be examinedseparately. When the minimum value S min ηη occurs at thedownstream boundary point η B , the inequality is easilysatisfied for all solutions S , whether advancing or reced-ing, since as evident from Fig. 5, all the self-similar solu-tions asymptote to a uniform thickness, which thereforeyields values S ηη ( η → η B ) → > − S min ηη occurs at an interior point η int , it will thenbe the case that S min ηηη ( η int ) = 0. Differentiation of Eq.(46) then yields the corresponding third order equationevaluated at the point η int (cid:20) (cid:18) S η S (cid:19) (cid:2) − S + ( η + 2 S η ) ( η + 5 S η ) (cid:3)(cid:21) ( η int ) = 0 (66)whose solutions are given by S η ( η int ) = 0 or (67)= − η int ± p η int + 40 S ( η int )20 , (68)which when substituted into Eq. (46) provides the valueof the curvature at that interior point, namely S ηη ( η int ) = 0 or (69)= −
15 + 3 η int ± η int p η int + 40 S ( η int )50 S ( η int ) . (70)The minimal value of this relation is attained for thenegative root, which as shown is always less than -1: S min ηη ( η int ) = −
15 + 3 η int − η int p η int + 40 S ( η int )50 S ( η int )= −
15 + 350 η int p S ( η int ) s η int S ( η int ) − s η int S ( η int ) + 409 ! ≥ −
15 + 350 × (cid:18) − (cid:19) (71) > − . (72)The last inequality derives from the general relation η ( p η − p η + k ) ≥ − k/ η and k real and posi-tive.When the minimum value of S ηη occurs at the origin,evaluation of Eq. (46) subject to the Dirichlet condi-tion S (0) = 1 yields the criterion for stability, namely S ηη (0) = − S η (0) > −
1. Since advancing solutionsalways exhibit negative slopes at the origin, this cri-terion can be rewritten as S η (0) < − / √
2. Our nu-merical results show that in order for advancing solu-tions to remain strictly positive throughout the domain[0 , Lη B ] and not yield a termination point, it must be thecase that S η (0) . − . − / √
2. The following geometric argument supports thisconclusion as well.Since S ηη (0) = − S η (0) < S ηη ( η B ) = 0, theremust occur at least one inflection point in the domain4 S ( η ) η + η p S η (0) η p S ηη ( η p )=0 FIG. 6. (Color online) Representative self-similar solution S ( η ) (blue solid line) for advancing state. The (red) dot des-ignates the point η P where there occurs the first inflectionpoint away from the origin where S ηη ( η p ) = 0. The (pur-ple) dashed line represents a linear extension of the Taylorexpansion of S ( η ) about the origin to the point η P . [0 , η B ]. Consider the first such inflection point at position η p where S ηη ( η p ) = 0. As shown in Fig. 6, since S ( η p )lies below the extended line 1 + η p S η (0), then S ( η p ) < η p S η (0) = 1 − S η ( η p ) S η (0) < − S η (0). Theserelations are found by noting that S η ( η p ) = − η p / S η ( η p ) < S η (0) where both slopesare negative. This then establishes that S η (0) > − / √
2. Transient and asymptotic stability of receding solutions S ( η ) We next examine the stability of receding solutionswhich are characterized by a downstream boundary valuelarger than the inlet condition, namely S ( η B ) = const >S (0) = 1, and inlet slopes which are always posi-tive, namely S η (0) >
0. In order to make use of thelemma above, we first apply a change of variable toEq. (64) such that G ( η ) is replaced by the product S ( η ) G ( η ) >
0. Straightforward calculation yields yet an-other Sturm-Liouville eigenvalue equation with weightingfactor S ( η ): (cid:18) M + M † (cid:19) G ( η, τ ) = − λS ( η ) G ( η, τ ) , (73)where (cid:18) M + M † (cid:19) = ∂∂η (cid:18) S ∂∂η (cid:19) − S ηS η + S ) . (74)The lemma then applies here, with an adjustment for theweighting factor: λ min > ( − Q/S ) min . This then yields the result λ min > min (cid:26) S ( ηS η + S ) (cid:27) or likewise (75) − λ min < max (cid:26) − S ( ηS η + S ) (cid:27) . (76)Since for receding solutions all terms in the expression( ηS η + S ) /S are strictly positive, then λ min is strictly pos-itive and therefore ω max is strictly negative. This demon-stration establishes that receding solutions are tran-siently and asymptotically linearly stable to any infinites-imal disturbances G satisfying homogeneous Dirichletboundary conditions. D. Numerical results and comparison to analyticbounds
Our numerical results for the operator exponential gov-erning transient growth are next compared to analyticbounds derived above. Derivatives were constructed us-ing second-order finite differences. (In discrete form,these operators are simply square matrices.) Homoge-neous Dirichlet boundary conditions were enforced bydirectly reducing the operator matrices. The numeri-cal abscissa ω max was obtained by identifying the small-est eigenvalue of the corresponding reduced operator ma-trix ( L + L † ) / ω max = − λ min . Thespectral abscissa β max was obtained by identifying thesmallest eigenvalue of the reduced matrix L . The quan-tity ln( k exp( L τ ) k ), representing the maximum instanta-neous disturbance amplification, was obtained from theoperator norm k · k given by the matrix maximum sin-gular value. Shown in Fig. 7 are the numerical resultsfor these three quantities plotted alongside the analyticbounds for − λ min given by max {− ( ηS η + S ) / (4 S ) } for re-ceding solutions and max {− (1+ S ηη ) / } = − (1+ S min ηη ) / k exp( L τ ) k for all times τ fall intermediate between the analytic upper bound,namely max {− ( ηS η + S ) / (4 S ) } for receding solutions or − (1 + S min ηη ) / − β max . The numerical resultsalso confirm that as τ → ∞ , the slope of the functionln k exp( L τ ) k exactly equals the value − β max , as mustbe the case since the asymptotic decay rate of distur-bances is dictated by the eigenvalue of L with maximumreal part.The numerical results in Fig. 7 (and similar studiesnot shown of other receding and advancing states) revealthat there is no transient nor asymptotic growth of dis-turbances. Furthermore, the suppression of disturbancesas quantified by ln k exp( L τ ) k asymptotes to the eigen-value − β max already at early times τ < .
0. The suddendrop off in the value ln k exp( L τ ) k observed near τ = 8in Fig. 7(b) simply reflects onset of increased dampingof disturbances incurred in the thin precursor film. All5 τ
10 15 20-25-20-15-10-50 50 S ( ) = 1.33 τ
10 15 20-25-20-15-10-50 50 w max τ : Eq (62) ln || exp ( L τ ) || : Eq (54) − β max τ : Eq (56) max { − S ( η S η + S )} τ /4 : Eq (77) S ( ) = 0.01 w max τ : Eq (62) ln || exp ( L τ ) || : Eq (54) − β max τ : Eq (56)Receding soln Advancing soln − ( S ηη ) τ /4 : Eq (66) min FIG. 7. (Color online) Analytic bounds and numerical resultsconfirming transient and asymptotic stability of representa-tive receding and advancing self-similar states shown in Fig.5. (a) Results for receding state subject to Dirichlet condi-tion S ( η = 80) = 1 .
33. (b) Results for advancing state withDirichlet condition S ( η = 80) = 0 . results presented were computed using spatial domains0 ≤ η ≤ η B where 5 ≤ η B ≤
80 with fixed mesh sizeranging from 0.05 to 0.005. As shown in Fig. 8, a do-main length of 80 was sufficiently long to ensure numer-ical convergence irrespective of mesh size.
VI. CONCLUSION
In this work, we have examined the global and linearstability of solutions describing inertia-free flow of a thinwetting Newtonian film within an open slender triangu-lar groove of constant cross-section with constant liquidcontact angle. The study is limited to flow states forwhich the interface function representing the local filmthickness is always strictly positive thereby ruling outrupture states. The slender approximation (also knownas the lubrication or long wavelength approximation) en-sures that to order ε the local value of the fluid pressuredepends only on the local film thickness whose interface - 0.50- 0.45- 0.40- 0.35- 0.30- 0.25- 0.20 Domain length η B Analytic bounds w max S ( η B ) AdvAdvAdvRec
FIG. 8. (Color online) Results of convergence studies as afunction of the domain length η B . The labels Rec and
Adv refer respectively to receding and advancing self-similar solu-tions shown in Fig. 5. The quantity ω max is defined in Eq.(62). The analytic bound for the receding solution is givenby max − ( ηS η + S ) / (4 S ) defined in Eq. (77). The analyticbounds for the advancing solutions are given by (1 + S min ηη ) / shape in the plane normal to the streamwise flow repre-sents a segment of a circle. This approximation simplifiesthe governing equation of motion, which is described bya second order nonlinear diffusion equation.The relevant base states for stability analysis includestationary interface states H S ( Z ) and self-similar states S ( η ) which adopt the (dimensionless) Washburn scaling η = Z/T / . The stationary states allow for variousboundary conditions imposed at the ends of the chan-nel including Dirichlet-Dirichlet, Dirichlet-Neumann (i.e.flux condition), Neumann-constant volume or Dirichlet-constant volume. The self-similar states are volumenon-conserving and result from application of Dirichletboundary conditions. By exploiting an analogy withother gradient flow equations and examining the asymp-totic behavior of the corresponding Lyapunov function, itwas shown that (strictly positive) stationary interface so-lutions represent exponentially stable equilibrium points.Disturbances of any type therefore decay away at leastexponentially fast to restore the system back to the ini-tial stationary state. The second important result isthat advancing, uniform and receding self-similar statessatisfying Washburn dynamics are both transiently andasymptotically linearly stable to infinitesimal perturba-tions. This finding required implementation of a general-ized non-modal linear stability analysis since the inhomo-geneous nature of the self-similar states naturally givesrise to non-normal disturbance operators.These two results highlight the reasons why the inter-face states describing either transient flow (self-similarstates) or asymptotic steady state flow (stationary states)within slender open triangular grooves tend to be so sta-6ble against perturbations so long as the liquid always wetsthe channel sidewalls. This result, however, requires thatthe slender channel approximation be everywhere satis-fied, which in turn sets an upper bound on the interfaceslope, namely ( ∂h/∂z ) ∼ ( d/L ) ∼ O ( ε ). Alterna-tively, this constraint can be viewed as a low-pass filter-ing requirement which prevents low-frequency modes ordisturbances from ever generating high-frequency modes.In particular, neither the base states nor disturbances ap-plied to these base states can trigger streamwise capillarywaves. Consequently, the results obtained in this studycannot address effects such as interface leveling in thinfilms [64]. As noted by Yang and Homsy [43], flow causedby streamwise curvature can also not be neglected in thelimit ( θ + α ) → π/ θ is a constantequal to the equilibrium contact angle, an assumptionknown to break down for sufficiently large flow speeds. Afully numerical model incorporating a velocity dependentcontact angle [65, 66] can certainly be developed, butsuch an extension precludes analytic solutions for basestates and therefore development of analytic bounds per-tinent to stability. We are currently extending this modelin a new promising direction by including effects due tosubstrate axial curvature. We hope that the stabilityfindings reported in this work as extended to curved sub-strates will ultimately provide a full analytic and numer-ical treatment particularly useful to the design of novelpropellant management systems for various space appli-cations. ACKNOWLEDGMENTS
The authors gratefully acknowledge financial supportfrom a 2016 NASA/Jet Propulsion Laboratory Presi-dent’s and Director’s Fund and a 2017 NASA Space Tech-nology Research Fellowship (NCW).
Appendix A: Boundedness of constant C in Eq. (39) We first note that the constant C given by Eq. (39),namely C = min T,Z (cid:26) H + 4 H H S + 6 HH S + 3 H S H + HH S + H S ) (cid:27) . (A1)satisfies the condition C ≥ H ( Z, T ) > H S ( Z ) >
0. We can eliminate the troublesome possi-bility C = 0, however, as follows. Assuming H S ( Z ) and H ( Z, T ) remain bounded, then the function in bracketsis observed to decrease monotonically with H and H S since the derivatives ∂∂H (cid:18) H + 4 H H S + 6 HH S + 3 H S H + HH S + H S ) (cid:19) = (A2) − H (cid:0) H + 3 H H S + 6 HH S + 5 H S (cid:1) H + HH S + H S ) < ∂∂H S (cid:18) H + 4 H H S + 6 HH S + 3 H S H + HH S + H S ) (cid:19) = (A4) − H S (cid:0) H + 3 HH S + H S (cid:1) H + HH S + H S ) < . (A5)The smallest possible value of the bracketed term in C is attained when that ratio is evaluated at the maximaof H S ( Z ) and H ( Z, T ) (which don’t necessarily occur atthe same point Z ). It then follows that C ≥ (cid:18) H + 4 H H S + 6 HH S + 3 H S H + HH S + H S ) (cid:19) , (A6)where H and H S are evaluated at their respective max-ima. Therefore, so long as H ( Z, T ) < ∞ , C > H ( Z, T ) might potentially diverge at finitetimes depending on the boundary conditions imposed.We eliminate the possibility H ( Z, T ) → ∞ by appealingto the governing equation of motion, which essentiallydescribes the intrinsic diffusive behavior of the H .Consider first the behavior of Eq.(41) given by ∂H∂T = (cid:18) ∂H∂Z (cid:19) + H ∂ H∂Z (A7)in the vicinity of local extrema. Any local maximumof H will satisfy ∂H/∂Z = 0, ∂ H/∂Z ≤ ∂H/∂T ≤
0, which leads to a diminishment in heightdue to the diffusive nature of the underlying equation ofmotion. Similarly, a local minimum of H cannot furtherdecrease in value. Therefore, any new extrema beyondthose present in the initial condition can only be reachedby possible extrema at the boundaries. For those solu-tions satisfying Dirichlet conditions at both endpoints ofthe domain, H can never therefore attain a new maxi-mum above the maximum value of the initial conditionnor attain a new minimum below the minimum value ofthe initial condition. This behavior then guarantees thatso long as the initial condition satisfies the constraint H ( Z, T = 0) >
0, then C will always remain boundedfrom below.A Neumann, or equivalently, a flux boundary conditionposes a potential problem since too large a flux for a giveninitial condition could potentially lead to a local drainage7spot where the film thickness vanishes. To assess thiscase, we recast Eq. (19) not in terms of the local variable H ( Z, T ) but the local flux Q ( Z, T ) = − H ∂H/∂Z = − (1 // ∂H /∂Z by first multiplying Eq. (19) by − H/ Z , which yields ∂Q∂T = ∂∂Z (cid:18) H ∂Q∂Z (cid:19) = H ∂ Q∂Z + ∂H∂Z ∂Q∂Z . (A8)Adopting a similar argument as above, for any case re-quiring constant liquid volume, the prescribed fluxes atthe endpoint must be equal, thereby imposing Dirichletconditions on the flux Q . It then follows that Q ( Z, T ) cannever exceed the extrema present in the initial condition Q ( Z, T = 0), which therefore precludes H from becom-ing arbitrarily large; hence C is again bounded below atall times.For cases subject to one Dirichlet and one flux bound-ary condition, the only way in which H might approachinfinity is if the flux boundary condition is made verylarge. However, this would give rise either to an arbi-trarily large internal flux extremum, which is disallowedby the previous argument, or an arbitrarily large flux ofliquid at the boundary subject to the Dirichlet condition. For the latter case, the divergence in H at one endpointwould be driven by the flux condition at other endpoint,such that the system would approach infinite volume andthe Lyapunov functional would therefore not decrease asrequired. Hence, such a situation cannot occur, and H is therefore bounded above at all times and C boundedbelow at all times.Finally, for the case of a Dirichlet boundary conditioncoupled with a requirement of constant volume, were H to approach infinity at one boundary, the flux there wouldalso diverge and be matched by identical divergence atthe other end point in order to satisfy constant volume.But according to Eq. (A8), an interior flux minimumwould have to undergo increase thus increasing the min-imum slope of H , which would eventually violate theconstraint of constant volume. Hence, in this case aswell, H cannot diverge at any point in time and C there-fore remains bounded below at all times.In conclusion, the demonstration above therefore es-tablishes that all stationary states, irrespective of theboundary conditions imposed above, are exponentiallystable in the Lyapunov sense. [1] J. Rollins, R. Grove, and D. Jaekle, in AIAA 21st JointPropulsion Conference , edited by A. I. of Aeronautics andAstronautics (1985) pp. AIAA–85–1199.[2] D. E. Jaekle, in
AIAA/SAE/ASME/ASEE 27th JointPropulsion Conference, June 24-26, 1991, Sacramento,CA (American Institute of Aeronautics and Astronau-tics, 1991) pp. AIAA–91–2172.[3] D. Levine, B. Wise, R. Schulman, H. Gutierrez, D. Kirk,N. Turlesque, W. Tam, M. Bhatia, and D. Jaekle, J.Propul. Power , 429 (2015).[4] J. W. Hartwig, in , AIAA Propulsion and Energy Forum(2016) pp. AIAA–16–4772.[5] J. W. Hartwig, J. Spacecr. Rockets , 1 (2017).[6] D. E. Jaekle, in (American Institute of Aeronautics and Astro-nautics, 1997) pp. AIAA–97–2811.[7] M. M. Weislogel, AIAA Journal , 2320 (2001).[8] M. M. Weislogel and S. H. Collicott, in (2002).[9] S. H. Collicott and Y. Chen, Microgravity Sci. Technol. , 487 (2010).[10] S. R. Darr, C. F. Camarotti, J. W. Hartwig, and J. N.Chung, Phys. Fluids , 017101 (2017).[11] E. W. Washburn, Phy. Rev. XVII , 273 (1921).[12] R. T. Jones, Annu. Rev. Fluid Mech. , 223 (1969).[13] R. Skalak, N. Ozkaya, and T. C. Skalak, Annu. Rev.Fluid Mech. , 167 (1989).[14] R. A. Wooding and H. J. Morel-Seytoux, Annu. Rev.Fluid Mech. , 233 (1976).[15] J. Johansson and G. Kifetew, Eur. J. Wood Prod. , 77(2010). [16] P. Yager, T. Edwards, E. Fu, K. Helton, K. Nelson, M. R.Tam, and B. H. Weigl, Nature , 412 (2006).[17] P. S. Dittrich and A. Manz, Nat. Rev. Drug. Discov. ,210 (2006).[18] T. Thorsen, S. J. Maerkl, and S. R. Quake, Science ,580 (2002).[19] M. Prakash and N. Gershenfeld, Science , 832 (2007).[20] T. P. Cotter, Principles and prospects for micro heatpipes , Library Without Walls Project (Los Alamos Na-tional Lab., NM (USA), 1984).[21] A. K. Mallik, G. P. Peterson, and M. H. Weichold, J.Electron. Packaging , 436 (1992).[22] G. P. Peterson, A. B. Duncan, and M. H. Weichold, J.Heat Trans. T. ASME , 751 (1993).[23] J. Qu, H. Wu, P. Cheng, Q. Wang, and Q. Sun, Int J.Heat Mass Trans. , 294 (2017).[24] B. D. Piorek, S. J. Lee, J. G. Santiago, M. Moskovits,S. Banerjee, and C. D. Meinhart, PNAS , 18898(2007).[25] Y. Chen, L. S. Melvin, S. Rodriguez, D. Bell, and M. M.Weislogel, Microelectron. Eng. , 1317 (2009).[26] P. Concus and R. Finn, PNAS , 292 (1969).[27] P. S. Ayyaswamy, I. Catton, and D. K. Edwards, ASMEJ. Appl. Mech. , 248 (1974).[28] M. Dong and I. Chatzis, J. Colloid Interface Sci. ,278 (1995).[29] M. M. Weislogel, Capillary Flow in an Interior Corner ,Ph.D. thesis, Northwestern University (1996).[30] M. M. Weislogel and S. Lichter, J. Fluid Mech. , 349(1998).[31] L. A. Romero and F. G. Yost, J. Fluid Mech. , 109(1996).[32] Y. Chen, M. M. Weislogel, and C. L. Nardin, J. Fluid Mech. , 235 (2006).[33] R. R. Rye, J. A. Mann, and F. G. Yost, Langmuir ,555 (1996).[34] R. R. Rye, F. G. Yost, and E. J. O’Toole, Langmuir ,3937 (1998).[35] D. Deng, Y. Tang, J. Zeng, S. Yang, and H. Shao, Int J.Heat Mass Trans. , 311 (2014).[36] J. Berthier, K. A. Brakke, E. P. Furlani, I. H. Karam-pelas, V. Poher, D. Gosselin, M. Cubizolles, andP. Pouteau, Sensor Actuat. B Chem. , 258 (2015).[37] M. M. Weislogel, Phys. Fluids , 3101 (2001).[38] D. Haake, J. Klatte, A. Grah, and M. E. Dreyer, Micro-gravity Sci. Tec. , 129 (2010).[39] Y. Wei, X. Chen, and Y. Huang, Science China Physics,Mechanics & Astronomy , 1551 (2013).[40] Y. Tang, X. Chen, and Y. Huang, Chinese J. Aeronaut. , 720 (2015).[41] D. Langbein, J. Fluid Mech. , 251 (1990).[42] L. Yang and G. M. Homsy, Phys. Fluids (2007).[43] L. Yang and G. M. Homsy, Phys. Fluids (2006).[44] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A.Driscoll, Science , 30 (1993).[45] B. F. Farrell and P. J. Ioannou, J. Atmos. Sci. , 2025(1996).[46] R. Lenormand and C. Zarcone, in , Society of Petroleum Engineer-ing (American Institute of Mining, Metallurgical andPetroleum Engineers, Inc., 1984).[47] W. I. Newman, J. Math. Phys. , 3120 (1984).[48] J. Ralston, J. Math. Phys. , 3124 (1984).[49] F. Otto, Commun. Part. Diff. Eq. , 101 (2001).[50] J. L. V´azquez, The Porous Medium Equation: Mathe-matical Theory , Oxford Science Publications (ClarendonPress, 2007).[51] W. Kern and B. U. Felderhof, Zeitschrift f¨ur Physik BCondensed Matter and Quanta , 129 (1977).[52] J. W. Cahn and J. E. Taylor, Acta Metall. Mater. , 1045 (1994).[53] L. Giacomelli and F. Otto, Calc. Var. Partial Dif. , 377(2001).[54] G. F. Teletzke, H. T. Davis, and L. Scriven, Chem. Eng.Comm. , 41 (1987).[55] L. W. Schwartz, R. V. Roy, R. R. Eley, and S. Petrash,J. Colloid Interface Sci. , 363 (2001).[56] F. Smith, T. Fearn, and S. Bullett, Advanced TechniquesIn Applied Mathematics , Ltcc Advanced Mathematics Se-ries (World Scientific Publishing Company, 2016).[57] MATLAB and Statistics Toolbox Release 2015a, TheMathWorks, Inc., Natick, MA, USA (2015).[58] J. M. Davis and S. M. Troian, Phys. Fluids , 1344(2003).[59] J. M. Davis, B. J. Fischer, and S. M. Troian, in In-terfacial Fluid Dynamics and Transport Properties , Lec-ture Notes in Physics, edited by R. R. Narayanan andD. Schwabe (Springer-Verlag, Berlin, 2003) pp. 79–105.[60] J. M. Davis and S. M. Troian, Phys. Rev. E , 016308(2003).[61] J. M. Davis and S. M. Troian, Phys. Rev. E , 046309(2004).[62] J. M. Davis, D. E. Kataoka, and S. M. Troian, Phys.Fluids , 092101 (2006).[63] R. Haberman, Applied Partial Differential Equationswith Fourier Series and Boundary Value Problems , 4thed. (Pearson Prentice Hall, 2004).[64] S. E. Orchard, Appl. Sci. Res. , 451 (1963).[65] M. Bracke, F. De Voeght, and P. Joos, The kinetics ofwetting: the dynamic contact angle , edited by P. Both-orel and E. Dufourc, Trends in Colloid and Interface Sci-ence III, Progress in Colloid & Polymer Science, Vol. 79(Steinkopff, 1989).[66] P. Joos, P. van Remoortere, and M. Bracke, J. ColloidInterface Sci.136