Why Capillary Flows in Slender Triangular Grooves Are So Stable Against Disturbances

TL;DR

This paper provides a theoretical analysis demonstrating that capillary flows in slender triangular grooves are highly stable against disturbances, explaining their widespread use in space and microfluidic applications.

Contribution

It introduces a novel theoretical framework showing the nonlinear and linear stability of capillary flows in open triangular channels, including Washburn dynamics.

Findings

Stationary interfaces are asymptotically nonlinearly stable.

Washburn dynamics are linearly stable to small perturbations.

Capillary flows are robust and stable in open grooved channels.

Abstract

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 that utilize surface tension to retain, channel and control flow in microgravity environments. Although it has been known for decades that capillary flows offer an ideal method of fuel management, there has 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 open triangular channels tend to be very stable against disturbances. By utilizing the gradient flow form of the governing fluid interface equation, we first prove that stationary interfaces in the presence of steady flow are asymptotically nonlinearly and exponentially stable in the Lyapunov sense. We then demonstrate that fluid interfaces exhibiting self-similar Washburn dynamics are transiently and asymptotically linearly stable to small perturbations. This second finding relies on a generalized non-modal stability analysis due to the non-normality of the governing disturbance operator. Taken together, these findings reveal the robust nature of transient and steady capillary flows in open grooved channels and likely explains the prevalent use of capillary flow management systems in many emerging technologies ranging from cubesats to point-of-care microfluidic diagnostic systems.

Tables

Table 1. Table 1: Characteristic scalings (lower case) and non-dimensional variables (uppercase) used to describe dimensionless system shown in Fig. 2 .

Quantity	Scaling	Rescaled
		variable
Slender parameter		$\epsilon =d/L\ll 1$
Coordinates	$x_c=d$	$X=x/x_c$
	$y_c=d$	$Y=y/y_c$
	$z_c=L$	$Z=z/z_c$
Velocity	$u_c={\epsilon }^2\gamma \text{Ca}/\mu$	$U=u/u_c$
	$v_c={\epsilon }^2\gamma \text{Ca}/\mu$	$V=v/v_c$
	$w_c=\epsilon \gamma \text{Ca}/\mu$	$W=w/w_c$
Pressure	$p_c=\gamma \text{Ca}/\left(\epsilon L\right)$	$P=p/p_c$
Time	$t_c=\mu L/\left(\epsilon \gamma \text{Ca}\right)$	$T=t/t_c$
		$\tau =\ln \left(T\right)$
Interface midline	$y_c=d$	$H=h(z,t)/y_c$
thickness
Interface shape	$y_c=d$	$\Sigma =\sigma (x,z,t)/y_c$
Interface radius	$y_c=d$	$R=r_c/y_c$
of curvature
Stationary state		$H_S\left(Z\right)$
midline thickness
Self-similar variable		$\eta =Z/\sqrt{T}$
Self-similar state		$S\left(\eta \right)$
midline thickness
Bond number		$\text{Bo}=\rho g{d}^2/\gamma$
Capillary number		$\text{Ca}=\mu w_c/\epsilon \gamma =\Phi$
Reynolds number		$\text{Re}=\rho w_{c}d/\mu$
		$=\left(\epsilon \rho \gamma d/{\mu }^2\right)\text{Ca}$

Equations

$$\nabla \cdot \vec{u}$$

$$\rho \left[\frac{\partial \vec{u}}{\partial t}+\left(\vec{u}\cdot \nabla \right)\vec{u}\right]$$

$$0$$

$${\epsilon }^3\text{Re}\frac{DU}{DT}$$

$${\epsilon }^3\text{Re}\frac{DV}{DT}$$

$$\epsilon \text{Re}\frac{DW}{DT}$$

$$\frac{D}{DT}$$

$$\Delta$$

$$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial Y}+\frac{\partial W}{\partial Z}$$

$$\frac{\partial P}{\partial X}=\frac{\partial P}{\partial Y}$$

$$\frac{\partial P}{\partial Z}$$

$${\left[\hat{n}\cdot (\tilde{\tau }-p\mathbb{I})\cdot \hat{n}+\gamma ({\nabla }_s\cdot \hat{n})\right]}_{y=\sigma (x,z,t)}$$

$${\left[{\hat{t}}_{i=1,2}\cdot \tilde{\tau }\cdot \hat{n}\right]}_{y=\sigma (x,z,t)}$$

$$\hat{N}$$

$${\hat{T}}_1$$

$${\hat{T}}_2$$

$$0$$

$$=-P-{\text{Ca}}^{-1}\text{K}(Z,T)+O({\epsilon }^2),$$

$$0$$

$$+2({\partial }_X\Sigma )(-{\partial }_{X}U+{\partial }_{Y}V)$$

$$-({\partial }_Z\Sigma )\left[{\partial }_{X}W+({\partial }_X\Sigma ){\partial }_{Y}W\right]+O({\epsilon }^2),$$

$$P(Z,T)=$$

$$\hat{R}(\alpha ,\theta )=\frac{\sin \alpha }{\cos \theta -\sin \alpha }.$$

$$\tilde{X}$$

$$\tilde{W}(\tilde{X},\tilde{Y})$$

$$\tilde{\Sigma }(\tilde{X})$$

$$Q(Z,T)$$

$$=-\frac{{\text{Ca}}^{-1}}{\hat{R}(\alpha ,\theta )}\Gamma (\alpha ,\theta )H^2\frac{\partial H}{\partial Z},$$

$$\hat{A}=\frac{\cos \theta \sin \alpha \cos (\alpha +\theta )-\left(\pi /2-\alpha -\theta \right){\sin }^2\alpha }{(\cos \theta -\sin \alpha {)}^2}$$

$$\hat{A}(\alpha ,\theta )\frac{\partial H^2}{\partial T}=\frac{{\text{Ca}}^{-1}}{\hat{R}(\alpha ,\theta )}\Gamma (\alpha ,\theta )\frac{\partial }{\partial Z}\left(H^2\frac{\partial H}{\partial Z}\right).$$

$$\Phi (\alpha ,\theta )=\frac{\Gamma (\alpha ,\theta )}{\hat{A}(\alpha ,\theta )\hat{R}(\alpha ,\theta )}.$$

$$\frac{\partial H^2}{\partial T}-\frac{\partial }{\partial Z}\left(H^2\frac{\partial H}{\partial Z}\right)=0,$$

$$Q_o$$

$$=-\hat{A}(\alpha ,\theta )H_S^2\frac{d{H}_S}{dZ}.$$

$$H_S={\left[\text{const}-\frac{3Q_{o}Z}{\hat{A}(\alpha ,\theta )}\right]}^{1/3}>0$$

$$H_S={\left[H_1^3-\frac{3Q_S}{\hat{A}(\alpha ,\theta )}(Z-Z_1)\right]}^{1/3}$$

$$H_S(Z)={\left[H_2^3\frac{Z-Z_1}{Z_2-Z_1}+H_1^3\frac{Z_2-Z}{Z_2-Z_1}\right]}^{1/3}.$$

$$V_S$$

$$=\frac{3^{5/3}{\hat{A}}^{1/3}}{5Q_S}\left[(C_S-Q_S{Z}_1{)}^{5/3}-(C_S-Q_S{Z}_2{)}^{5/3}\right].$$

Full text

I Introduction

It has been known for centuries that wetting liquids will rapidly and spontaneously creep along surfaces containing grooves, interior corners, crevices or roughened areas, a process known as wicking. Since the late 1960s, researchers have been incorporating this passive and reliable method of flow control in the design of novel propellant management devices able to store, channel and meter fuel resourcefully in microgravity environments Rollins et al. (1985); Jaekle (1991); Levine et al. (2015); Hartwig (2016, 2017). Such systems have significantly extended mission lifetimes of spacecraft and satellites, enabling future interplanetary explorations as well. Modern propellant management systems consist of combinations of sponges, traps, troughs, vanes, and wicks to channel propellant flow by capillary action, systems which have been investigated extensively Jaekle (1991, 1997); Weislogel (2001a); Weislogel and Collicott (2002); Collicott and Chen (2010); Darr et al. (2017). Shown in Fig. 1 are two common structures designed in such a way that the liquid film thickness is much smaller than the streamwise flow distance, a limiting ratio which as described in Section III leads to considerable simplification of the governing equations of motion. However, despite that by the 1920’s the mechanism describing internal capillary flow, whereby a liquid column spontaneously fills the interior of a slender capillary tube, was understood and the appropriate equations developed, the mechanism driving spontaneous capillary flow along an open grooved channel required a half-century more to be deduced. Ongoing efforts to miniaturize fluid management systems for many different applications continue to drive interest in the fundamentals of free surface capillary flow along structured substrates.

Edward Washburn appears to have provided the first theoretical analysis of the distance in time traversed by a Newtonian liquid flowing within a horizonal enclosed cylindrical tube under the action of capillary forces Washburn (1921). The Washburn relation, as it is now known, is given by the relation $z=(\gamma\,\cos\theta\,R\,t/2\mu)^{1/2}$, where $z$ is the distance traversed by the advancing (for wetting fluids) or receding (for non-wetting fluids) meniscus, $t$ is time, $R$ is the cylinder radius, and $\mu$, $\gamma$ and $\theta$, respectively, are the liquid shear viscosity, surface tension and contact angle. Although this largely unidirectional flow is driven by the capillary pressure drop across the curved meniscus separating gas from liquid, the Washburn relation can be regarded from a mathematical point of view as an example of diffusive driven interface growth since the “diffusion” coefficient preceding the time variable is described by units of length squared per unit time. By treating porous bodies as an assemblage of small cylindrical capillaries, Washburn also went on to predict that the volume of liquid penetrating a porous medium should scale in time as $(\gamma t/\mu)^{1/2}$, where geometric factors and wettability constants were incorporated into an overall proportionality constant. Washburn thereby revealed the essential physical mechanism governing the spontaneous creep of liquids into enclosed spaces, a process that can easily dominate and oppose the force of gravity for sufficiently small enclosures. The Washburn scaling has since been used successfully to quantify capillary transport in many different systems ranging from blood flow in microvascular hemodynamics Jones (1969); Skalak et al. (1989) to oil extraction from porous rocks Wooding and Morel-Seytoux (1976) to water uptake in dried tree specimens Johansson and Kifetew (2010), to name a few. Through modern day advances in microfabrication techniques, capillary action is now being incorporated into the design of many microfluidic devices as well, some involving spontaneous flow through chemically treated porous substrates, others relying on a combination of capillary action, positive displacement pumping and electrophoresis. The number of such applications is multiplying rapidly with emphasis on disposable inexpensive platforms beneficial to global public health Yager et al. (2006), drug discovery screening Dittrich and Manz (2006) and specialized fluid based logic circuitry Thorsen et al. (2002); Prakash and Gershenfeld (2007).

While Washburn originally focused on the fluid dynamics underlying internal capillary flow, there soon came the realization that capillary forces were also somehow responsible for the spontaneous exterior flow of liquid along open surfaces manifesting grooves, edges, corners and roughened patches. This type of free surface flow is typically characterized by a solid boundary whose surface pattern is fixed and a gas/liquid interface whose shape undergoes spatial and temporal variation as the flow evolves toward steady state conditions. Free surface wicking is now also being incorporated into many small-scale fluidic devices such as heat pipes for cooling of microelectronics Cotter (1984); Mallik et al. (1992); Peterson et al. (1993); Qu et al. (2017), surface-enhanced Raman spectroscopic devices Piorek et al. (2007) and small spacecraft fluid management systems Chen et al. (2009).

II Background

In one of the earliest studies of free surface wicking of a Newtonian liquid, Concus and Finn Concus and Finn (1969) in 1969 showed that a liquid drop will penetrate the interior of an open triangular channel with opening half-angle $\alpha$ so long as $\theta+\alpha<\pi/2$. Shortly thereafter, Ayyaswamy et al. Ayyaswamy et al. (1974) 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 transient capillary flow in idealized containers with grooves or corners were demonstrated two decades later by several groups working independently. Dong and Chatzis Dong and Chatzis (1995), Weislogel and Lichter Weislogel (1996); Weislogel and Lichter (1998) and Romero and Yost Romero and Yost (1996) derived and solved for the interface equation $h(z,t)$ describing the thickness of a liquid film flowing within an open slender triangular groove satisfying the condition $d\ll L$, as sketched in Fig. 2. The mechanism responsible for spontaneous capillary flow within an open groove turned out to be more complex than its counterpart for small enclosures. While for interior flow the liquid column is pulled forward by the net force of surface tension acting along the triple line (i.e. the line defining the junction between the gas, liquid and solid media), the free boundary in an open capillary driven system experiences a normal force across the entire interface whose shape typically varies in space and time. By implementing two key assumptions, namely that the streamwise variation in liquid height varies smoothly on a length scale $L$ much larger than $d$, and that the flow is slow enough such that capillary forces dominate inertial and viscous forces, these researchers were able to show that the cross-section of the liquid surface must describe a circular arc of radius $R(z)$. Hence the pressure gradient driving the flow stems exclusively from the streamwise variation in $R(z)$. These two assumptions therefore led to a nonlinear diffusion equation describing the relaxation of the cross-sectional liquid area proportional to $h^{2}(z,t)$. While Romero and Yost Romero and Yost (1996) and Dong and Chatzis Dong and Chatzis (1995) assumed unidirectional flow as well, Weislogel and Lichter Weislogel (1996); Weislogel and Lichter (1998) were able to derive similar results without relying on this simplification by instead conducting a formal perturbation expansion based on the slender parameter approximation $(\varepsilon)^{2}=(d/L)^{2}\ll 1$. These separate efforts led to the conclusion that there exist transient solutions according to which the penetration distance $z$ scales in time as $(D\,t)^{1/2}$ where $D=\gamma\,d\,K(\theta,\alpha)/\mu$, $d$ is a characteristic liquid thickness (such as the groove depth or midline liquid thickness) and $K(\theta,\alpha)$ is a geometric factor. Chen, Weislogel, and Nardin Chen et al. (2006) later generalized this result to channels with constant cross-section ranging from the nearly-rectangular to highly rounded grooves. In all cases, it was confirmed that the penetration distance retains the $t^{1/2}$ scaling. Numerous experimental studies Rye et al. (1996, 1998); Weislogel (1996); Chen et al. (2009); Deng et al. (2014); Berthier et al. (2015) and references therein have since confirmed their predictions.

A few studies have gone beyond solution of the flow base states to explore their stability spectrum. Weislogel Weislogel (2001b) examined the linear stability of an infinite, initially quiescent column of liquid supported within a slender triangular groove of constant cross-section using conventional modal analysis. Other researchers have since examined the free surface collapse and gas ingestion in inertial-capillary flows for rectangular Haake et al. (2010), triangular Wei et al. (2013) or asymmetrical Tang et al. (2015) grooves, the stability of liquid menisci attached to a solid edge Langbein (1990); Yang and Homsy (2007), and the dryout instability for thermocapillary (i.e. non-isothermal) driven flow in triangular grooves Yang and Homsy (2006). Aside from these specialized studies, however, the majority of theoretical studies have been devoted exclusively to solutions governing unperturbed base states i.e. solutions in the absence of disturbances.

In this work, we address for the first time the stability of stationary and transient interfaces associated with capillary flows in slender, open triangular grooves of constant cross-section subject to the Concus-Finn inequality Concus and Finn (1969) ensuring liquid penetration ($\theta+\alpha<\pi/2$) and the restriction that the liquid always wet the channel sidewalls. In Section III, we outline the derivation of the governing interface equation Romero and Yost (1996); Dong and Chatzis (1995); Weislogel (1996); Weislogel and Lichter (1998). In Section IV, we examine the nonlinear stability of stationary states subject either to Dirichlet, Neumann or constant volume conditions. By invoking the gradient flow form of the interface disturbance equation and an associated Lyapunov energy functional governing the decay in time of arbitrary disturbances, we demonstrate that stationary interface shapes are exponentially stable and therefore extremely robust to perturbations. In Section V, we examine the generalized non-modal linear stability of self-similar states obeying Washburn scaling Trefethen et al. (1993); Farrell and Ioannou (1996) using both analytic arguments and direct numerical computation. Our results indicate that self-similar solutions characterized by the scaling $z\sim t^{1/2}$ are both transiently and asymptotically stable to infinitesimal perturbations.

III Model for free surface capillary flow in slender open triangular grooves

In this section, we outline the theoretical model describing the capillary motion of a non-volatile, isothermal Newtonian liquid film flowing within an open slender triangular groove in contact with an ambient passive gas Dong and Chatzis (1995); Romero and Yost (1996); Weislogel (1996); Weislogel and Lichter (1998), as depicted in Fig. 2(a). The inlet midline film thickness at the origin is denoted by $d=h(x=0,z=0,t)$ and the channel length by $L$ where the slender limit is assumed, namely $d/L\ll 1$. Shown in Fig. 2(b) is a cross-sectional view of the flow geometry where $h(z,t)$ denotes the local midline film thickness, $\alpha$ is the groove interior half-angle, $r_{c}$ is the radius of curvature of the liquid interface and $\theta$ is the contact angle of the liquid wetting the channel sidewalls. The liquid film is assumed to maintain a constant contact angle independent of location or flow speed.

III.1 Slender limit form of the hydrodynamic equations

In order to conserve mass and momentum, an incompressible Newtonian liquid of constant density must satisfy the continuity and Navier-Stokes equation given by the set of coupled equations

[TABLE]

where the velocity field in Cartesian coordinates is represented by $\vec{u}=(u,v,w)$, the fluid pressure by $p(x,y,z)$, the gravitational acceleration by $\vec{g}$ and the constant fluid density by $\rho$. Assuming flow in a slender channel such that $\varepsilon=d/L\ll 1$ and a dominant balance between the pressure gradient and the viscous force per unit volume leads to the characteristic scalings and non-dimensional variables listed in Table 1 along with the corresponding Bond Bo, capillary Ca and Reynolds Re numbers in the slender limit.

To order $\varepsilon^{2}$, the rescaled forms of Eqs. (1) and (2) are then given by

[TABLE]

where the substantial derivative $D/DT$ and the Laplacian derivative $\Delta$, respectively, are given by

[TABLE]

and $\vec{G}=\vec{g}/g$. In the limits where $\varepsilon^{2}\ll 1$, $\varepsilon\,\textsf{Re}\ll 1$ and $\textsf{Bo}/\textsf{Ca}\ll\varepsilon$, the governing equations reduce to the form:

[TABLE]

When subject to the slender limit, the flow is therefore inertia-free and the fluid pressure is constant throughout the $(x,y)$ plane. The pressure gradient driving the flow, which will stem solely from capillary forces, can therefore only vary along the streamwise axis and can only be counterbalanced the viscous force set in play by the no-slip boundary condition applied along the groove sidewalls, namely $U=V=W=0$ at all liquid/solid interfaces.

III.2 Boundary conditions at the liquid interface

The two (dimensional) boundary conditions specifying the jump in normal and shear stresses across the gas/liquid interface $\sigma(x,z,t)$ are given by

[TABLE]

where $\mathbb{I}$ is the 3$\times$3 identity matrix, $\widetilde{\tau}=\mu[\nabla\vec{u}+(\nabla\vec{u})^{T}]$ denotes the shear stress tensor, $\nabla_{s}=(\nabla-\hat{n}(\hat{n}\cdot\nabla))$ denotes the surface gradient operator and the triad ($\hat{n},\hat{t}_{1},\hat{t}_{2}$) denotes the three unit vectors representing directions normal and tangent to the moving interface with the convention that $\hat{n}$ points away from the interface. In rescaled units, these unit vectors are given by

[TABLE]

where $\Sigma(X,Z,T)=\sigma(x,z,t)/d$ denotes the non-dimensional interface function and subscripts denote differentiation with regard to the rescaled coordinates. Specifying a system for which the fluid pressure derives solely from variations in the local interface curvature of the flowing liquid, the interfacial surface tension $\gamma$ is everywhere constant since the liquid is isothermal and contains no surfactant-like additives, and that the liquid remains in contact with a passive quiescent gas of negligible viscosity and density with gauge pressure set to zero, the jump in normal stress is then strictly due to capillary forces and the liquid interface is a surface of vanishing shear stress. To order $O(\varepsilon^{2})$ then, these dimensionless boundary conditions reduce to the form

[TABLE]

where K(Z) represents the local curvature of the interface function $\Sigma$ (defined to be positive for a wetting liquid). It is clear from Eq. (8b) that the curvature function K can in general only depend on $(Z,T)$ since according to Eq. (5b), the pressure $P$ is independent of $(X,Y)$. This then requires that the cross-sectional shape of the liquid interface be described by a curve with constant curvature. This restriction limits the shape either to a flat interface or one described by a segment of a circle. Since the liquid must also satisfy a prescribed contact angle set by the particulars of the liquid/solid interaction, a flat profile is disallowed and $\Sigma(Z,T)$ must therefore trace out a circular arc of constant curvature. The non-dimensional radius of curvature of the gas/liquid interface is then given by $R(Z,T)=H(Z,T)\,\sin\alpha/(\cos\theta-\sin\alpha)$, or likewise, the interface curvature is described by $K(Z,T)=1/R(Z,T)=(\cos\theta-\sin\alpha)(\csc\alpha)H^{-1}(Z,T)$. (This relation differs slightly from that originally derived by Romero and Yost Romero and Yost (1996) where they adopted a sign convention for K opposite to ours and chose the reference liquid thickness to be the height of the fluid intersecting the groove wall, not the inlet midline film thickness.) According to Eq. (8a), the capillary pressure is then given by

[TABLE]

where

[TABLE]

The Concus-Finn condition $\alpha+\theta<\pi/2$ for liquid imbibition yields $\widehat{R}(\alpha,\theta)>0$, or likewise $P(Z,T)<0$, consistent with a liquid interface with positive curvature. The case $\alpha+\theta>\pi/2$ is not relevant to our study since it ultimately leads to dewetting configurations resulting in a cascade instability resembling a linear array of primary, secondary, and tertiary droplets Yang and Homsy (2007)).

III.3 Interface midline equation $H(Z,T)$ for capillary flow in slender open triangular grooves

The fact that the interface shape can only be a segment of a circle, and is therefore independent of the local coordinates $(X,Y)$, leads to simplification of the expression for the streamwise volumetric flux $Q(Z,T)$. The relevant variables are then scaled by $H$ according to

[TABLE]

This rescaling allows for solution of ${\widetilde{W}}$ independently of the local value of $H$. As a result, the geometric function $\Gamma(\alpha,\theta)=\iint{\widetilde{W}}d{\widetilde{X}}d{\widetilde{Y}}$ need only be computed numerically once. The dimensionless streamwise flux which traverses the local cross-sectional area $A$ can then be re-expressed as

[TABLE]

where the integrated area $A=\widehat{A}(\alpha,\theta)H^{2}$ and

[TABLE]

is a geometric factor Weislogel (1996); Romero and Yost (1996). Since the local streamwise gradient in liquid flux is directly related to the time derivative of the local liquid cross-sectional area (see Appendix 2 of Ref. Lenormand and Zarcone, 1984 for derivation) according to $\partial A/\partial T=-\partial Q/\partial Z$, the governing nonlinear diffusion equation for the midline height $H(Z,T)$ is then given by

[TABLE]

Without loss of generality and to recast this equation into parameter-free form, the remaining scaling for the streamwise velocity is chosen to be $w_{c}=(\varepsilon\gamma/\mu)\Phi$, where $\Phi(\alpha,\theta)=\textsf{Ca}$ (see Table 1) is defined by

[TABLE]

The resulting interface equation, whose stability properties we examine in this work, is then given by

[TABLE]

subject to the constraint that $H(Z,T)$ is everywhere always positive. In Fig. 3 are plotted the scaled functions $\widehat{A}(\alpha,\theta)$, $\widehat{R}(\alpha,\theta)$, $\Gamma(\alpha,\theta)$ and $\Phi(\alpha,\theta)$ for four values of the liquid contact angle $\theta=0^{\circ}$, $20^{\circ}$, $40^{\circ}$ and $60^{\circ}$ as a function of increasing groove interior half-angle $\alpha$. While $\widehat{A}(\alpha,\theta)$, $\widehat{R}(\alpha,\theta)$ and $\Gamma(\alpha,\theta)$ are of order $O(10)$ or less, the values of $\Phi(\alpha,\theta)$ are far smaller and tend toward $O(10^{-2})$ or less. Note too that since $\widehat{R}(\alpha,\theta)$ and $\Gamma(\alpha,\theta)$ are non-negative functions for systems obeying the Concus-Finn condition, the direction of the liquid flux specified by Eq. (14) is then strictly determined by the sign of the local interface slope $\partial H/\partial Z$. Interfaces with $\partial H/\partial Z<0$ engender a positive local flux $Q$ and vice versa. A vanishing local flux results whenever $\partial H/\partial Z=0$.

The form of Eq. (19) falls within a class known as the porous media equation generally given by $\partial C(Z,T)/\partial T=\partial^{2}C^{m}/\partial Z^{2}$, where $C(Z,T)$ is a non-negative scalar function and $m$ is a constant larger than one Newman (1984); Ralston (1984); Otto (2001); Vázquez (2007). As discussed in Ref. [Vázquez, 2007], this nonlinear diffusion equation describes the relaxation of the order parameter $C(Z,T)$ relevant to various phenomena which arise in different branches of science and mathematics and exhibiting properties such as scale invariance and self-similarity. To help track the energy flow associated with the evolution of $C$, Newman Newman (1984) outlined a general method for constructing Lyapunov functionals for systems sustaining traveling wave solutions. He used that method to establish the actual rate of convergence of an initial configuration to solutions exhibiting self-similarity. Ralston Ralston (1984) complemented this work by showing that Newman’s choice of Lyapunov function allowed proof that initial conditions with finite mass converge to self-similar solutions asymptotically in time. When applied to our system, these results indicate that self-similar states which result from the spreading of an initial finite drop within a slender open triangular groove are asymptotically globally stable. However, despite longstanding interest in the porous media equation and wicking phenomena in general, there has been no prior study of which we are aware of the general stability of stationary or transient solutions corresponding to unconstrained volume states.

IV Stability of stationary interface solutions $H_{S}(Z)$

In this section, we first derive the analytic form corresponding to stationary solutions of the governing interface equation given by Eq. (19) with geometric factors explicit. By invoking analogy to the porous media equation and constructing an appropriate Lyapunov functional, we then examine the nonlinear asymptotic stability of these states against arbitrary disturbances.

IV.1 Stationary states for time-independent Dirichlet, Neumann and volume conditions

Stationary states of Eq. (19) correspond to those solutions for which $\partial H^{2}/\partial T=0$, which reduces the governing equation to a second order, ordinary differential equation requiring two boundary conditions. So long as the interface slope does not vanish, such stationary interfaces $H_{S}(Z)$ are possible because the subsurface flow of liquid establishes a balance between the local capillary and viscous stresses generates a constant flux $Q_{o}$, which according to Eq. (17) is given by

[TABLE]

The general form of these stationary solutions is then represented by

[TABLE]

described by a power law decrease of $Z^{1/3}$ for positive flux and power law increase for negative flux. Particular solutions, of course, require specification of two boundary conditions for setting the values of const and $Q_{o}$. For a triangular channel with fixed groove opening angle and liquid contact angle (i.e. constant value of $\widehat{A}(\alpha,\theta)$) extending between endpoints $Z_{1}$ and $Z_{2}$, stationary solutions $H_{S}$ correspond to

[TABLE]

for Dirichlet and Neumann conditions $H_{S}(Z_{1})=H_{1}$ and $Q_{2}=Q_{S}$, respectively. Alternatively, Dirichlet conditions imposed at both endpoints, such that $H_{S}(Z_{1})=H_{1}$ and $H_{S}(Z_{2})=H_{2}$, yield

[TABLE]

Note from Eq. (9) and the scaling for the streamwise flow speed $w_{c}$ that specification of the boundary film thickness is equivalent to specification of the boundary fluid pressure. Solutions $H_{S}$ to Eq. (22) can also be generated subject to constant flux $Q_{S}$ at one boundary ($Q_{1}=Q_{S}$ or $Q_{2}=Q_{S}$) and conservation of volume $V_{S}$, which sets the constant value $C_{S}$ in the implicit relation :

[TABLE]

Likewise, constant volume and fixed liquid height at one endpoint yield similar forms.

Representative solutions $H_{S}(Z)$ are plotted in Fig. 4 for Dirichlet conditions which pin the inlet height to $H_{1}=1.0$ and pin the outlet height $H_{2}$ to the five values shown. From the expression for the flux given by Eq. (20), it is evident that the solution with $H_{2}=1.33$, which exhibits interface slope values $dH/dZ$ which are everywhere positive, describes a stationary solution with net streamwise flux $Q_{S}<0$ i.e. net flow directed from right to left. The uniform solution $H_{2}(Z)=1.0$ clearly then represents a case with no flux and no subsurface flow i.e. a quiesence liquid filament. The remaining curves with negative interface slopes throughout the domain correspond to stationary solutions with positive flux i.e. net flow directed from left to right. Other boundary conditions yield similar shapes with characteristic scaling $Z^{1/3}$.

IV.2 Nonlinear exponential stability of stationary solutions

We next examine the time-asymptotic nonlinear stability of stationary solutions by appealing to a Lyapunov analysis. In particular, we construct a Lyapunov energy function (akin to a potential energy function in classical mechanics) to determine whether and how rapidly arbitrary disturbances equilibrate back to the stationary solutions $H_{S}$. A similar approach to asymptotic stability has previously been used to characterize disturbance decay governed by nonlinear diffusion dynamics Kern and Felderhof (1977). Here we show that the capillary flow in a bounded domain with time-independent boundary conditions has a uniquely determined stationary solution which is dynamically stable. This implies that any initial distribution which maintains the condition $H(Z,T)>0$ for all time will ultimately evolve toward the stationary state. Furthermore, it is shown that initial states satisfying Dirichlet boundary conditions are globally stable - that is, any initial distribution $H(Z,0)>0$ will converge to the stationary state. For all sets of boundary conditions posed, we show that the convergence to the stationary state is exponential in time so long as it satisfies the positivity condition $H(Z,T)>0$.

To begin, it proves convenient to recast Eq. (19) into gradient flow form Cahn and Taylor (1994); Giacomelli and Otto (2001)

[TABLE]

where $M(G)$ is the so-called mobility function and $\delta\mathfrak{F}/\delta G$ denotes the functional derivative of $\mathfrak{F}(H)$ with respect to $G$. The quantity $\mathfrak{F}(H)$ represents the Lyapunov free energy $\mathfrak{F}$ which for stable stationary states must satisfy the relations $\delta\mathfrak{F}=0$ (i.e. extremal condition) and $\delta^{2}\mathfrak{F}>0$ (positivity condition on Hessian for stable flow). By introducing $G(H)=H^{2}(Z,T)$ and letting $M=1/3$, Eq. (19) is re-expressed as

[TABLE]

where the function $R(G)$ is given by

[TABLE]

(The mapping from $H$ to $G$ is one-to-one since $H$ is strictly positive.) Comparison with Eq. (27) then yields the correspondence $\delta\mathfrak{F}(G)/\delta G=R(G)$ subject to the constraint that $[R\partial R/\partial Z]_{Z_{1},Z_{2}}=0$. This definition of $R(G)$ ensures that stationary states are identified by the extrema of the free energy $\mathfrak{F}(H)$ where $\left.\delta\mathfrak{F}(G)/\delta G\right|_{G_{S}}=R(G_{S})=0$. Additionally, Dirichlet boundary conditions correspond to $R|_{Z_{1},Z_{2}}=0$ and Neumann boundary conditions to $(\partial R/\partial Z)|_{Z_{1},Z_{2}}=0$. Multiplication of both sides of Eq. (28) by $R(G)$ followed by integration over the finite domain $[Z_{1},Z_{2}]$ yields the relation

[TABLE]

The left term in Eq. (28) can also be expressed as

[TABLE]

where the last term in Eq. (31) has been added, without loss of generality, to ensure the integrand vanishes for $G=G_{S}$, which essentially sets the value of the ground state energy. Given that $\delta\mathfrak{F}/\delta G=R[G(Z)]$, it is evident that the integral defined in Eq. (31) is none other than the Lyapunov free energy, which can be re-expressed in terms of the interface function $H$:

[TABLE]

Equating Eq. (30) and Eq. (31) then yields the relation

[TABLE]

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 described is asymptotically stable in the Lyapunov sense i.e. arbitrary initial distributions $H(Z,T)$ decay in time toward the stationary solution $H_{S}$, characterized by a vanishing first variation and positive second variation, namely

[TABLE]

To complete this part of the proof, we show that the stable stationary states $H_{S}$ are in fact also unique. Suppose then that in addition to the already specified steady state solution $H_{S}$ there exists another steady solution $H_{S2}$ satisfying the same boundary conditions. Since $H_{S2}$ is time independent, it must satisfy the relation $\partial\mathfrak{F}(H_{S2})/\partial T=-(1/3)\int_{Z_{1}}^{Z_{2}}\left[\partial(H_{S2}^{3}-H_{S}^{3})/\partial Z\right]^{2}dZ=0$. Clearly this is only possible if $H_{S2}=(H^{3}_{S}+C_{o})^{1/3}$, where the constant $C_{o}$ must equal zero in order for $H_{S2}$ to satisfy the same boundary conditions as $H_{S2}$. This result therefore establishes that the solution $H_{S}$ is indeed unique.

The proof above requires that both the stationary $H_{S}(Z)$ and transient $H(Z,T)$ solutions be everywhere strictly positive. For solutions satisfying Dirichlet boundary conditions, these requirements are easily met. Consider any positive initial state $H(Z,T=0)$ which redistributes its height in time according to Eq. (41). Were there a point within the domain where $H(Z,T)\to 0$, then the local interface would have to satisfy $\partial H/\partial T<0$ and the local curvature $\partial^{2}H/\partial Z^{2}$ would have to become sufficiently negative to satisfy the balance of terms Eq. (41). The local interface would then have to develop a local protrusion (negative curvature) and not a local dimple leading to rupture (positive curvature), in contradiction to the assumption of positive curvature required at that local minimum. Therefore, all positive initial states remain positive in time. For Dirichlet conditions then, stationary solutions $H_{S}(Z)$ are globally exponentially stable. For the remainder stationary solutions subject to a Neumann boundary condition, it may become the case that too high a flux condition leads to one of more points where the local film thickness vanishes. In the vicinity of such points, the local interface slope will become increasingly large and eventually violate the slender limit approximation. Perhaps more importantly, disjoining pressure effects, known to modify the flow in ultrathin regions of a liquid film, must then be incorporated into the model interface equation Teletzke et al. (1987); Schwartz et al. (2001). In such cases then, the proof above only establishes asymptotic stability, not global asymptotic stability.

Next, we demonstrate an even stronger statement regarding the asymptotic stability of stationary states, namely that the stationary solutions $H_{S}$ represent exponentially stable equilibria of Eq. (19). A useful relation for the function $R(Z)$ defined in Eq. (29) can be obtained by first applying the Poincaré-Friedrichs inequality in one dimension Smith et al. (2016) according to which

[TABLE]

subject to the constraint that there is some interior point $Z^{*}$ in the bounded domain $Z_{1}\leq Z^{*}\leq Z_{2}$ where $R(Z^{*})=0$. This must always be the case, however, since $H(Z,T)$ must satisfy the same boundary conditions as $H_{S}$. Given that $H(Z,T)$ and $H_{S}(Z)$ are strictly positive throughout the domain, it is possible to construct an upper bound on $\partial\mathfrak{F}/\partial T$ in Eq. (33) as follows:

[TABLE]

where $\mathcal{C}$, which is strictly positive, is defined as

[TABLE]

(See Appendix A for further discussion regarding the boundedness property of $\mathcal{C}$). The final inequality sought is then given by

[TABLE]

which confirms that the Lyapunov free energy of any disturbed state $H(Z,T)$ satisfying the same boundary conditions as the initial stationary state will decay back to the stationary state $H_{S}$ at least as fast as an exponential with a decay rate that scales with the square of the spatial domain size. This proof applies to all categories of stationary solutions discussed in Section IV.1. For disturbance functions subject either to two Dirichlet conditions or one Dirichlet condition and one Neumann (constant flux) or constant volume condition, the proof is trivial since either $R(Z_{1})$ or $R(Z_{2})$ vanishes identically. For solutions $H_{S}$ that correspond to fixed flux $Q_{S}$ at one boundary and constant volume $V_{S}$, the proof above also applies since there always exists an interior point $Z_{1}\leq Z^{*}\leq Z_{2}$ where $R\left(Z^{*}\right)=0$. This follows because two solutions $H$ and $H_{S}$ cannot have the same constant volume unless the functions $H(Z,T)$ and $H_{S}(Z)$ undergo at least one crossing point within the domain.

V Stability of self-similar solutions

We next seek unperturbed (i.e. base state) transient solutions to Eq. (19) for non-conserved volume which manifest self-similarity. A global stability argument as presented in Section IV.2 for stationary solutions proved unsuccessful since the self-similar form of the governing interface equation cannot be converted into gradient flow form. Therefore we instead applied a generalized linear stability analysis Farrell and Ioannou (1996) which helps determine whether infinitesimal non-modal disturbances undergo any transient or asymptotic amplification.

V.1 Self-similar solutions with time independent Dirichlet conditions

Previous studies have delineated the conditions leading to existence and uniqueness of self-similar solutions as well as the attraction of spatially confined initial distributions toward self-similar base states Vázquez (2007). Here we focus on volume non-conserving positive states $S(\eta)$ consistent with time-independent Dirichlet boundary conditions imposed at the domain endpoints, namely $S(0)=1$ and $S(\eta\rightarrow\eta_{B})=\textrm{const}$ where $\eta_{B}$ denotes a location far downstream of the origin. The Dirichlet condition at the origin can be set to unity without loss of generality since as evident from Eq. (46), a rescaling involving a multiplicative factor of $S(0)$ leaves the governing equation unchanged.

In general, for self-similarity to hold, there can be no intrinsic length or time scale imposed on the flow, in contrast to the steady state solutions examined in the previous section which depend on the groove length $Z_{2}-Z_{1}$ . A simple scaling analysis of Eq. (19) reveals that self-similar solutions may be possible whenever $T<<L^{2}/H\sim O(L/\varepsilon)$. To find such solutions, it is convenient to expand and rewrite Eq. (19) in the form

[TABLE]

The ansatz $H_{\textrm{sim}}(\eta,T)$ defined by

[TABLE]

allows for a large class of self-similar solutions Vázquez (2007); Weislogel and Lichter (1998) satisfying the general second order nonlinear differential equation

[TABLE]

Inspection of the asymptotic behavior of $S(\eta)$ as $\eta\to\infty$ helps ascertain what range of exponents $\beta$ are required for bounded non-terminating (i.e. $S>0$) states such that $S_{\eta\eta}$ and $S_{\eta}$ asymptotically approach zero as $\eta\to\infty$ . While the first two terms on the left side of Eq. (44) then vanish identically, care must be taken with regard to the third term which couples an increasingly large value of $\eta$ with a diminishingly small term $S_{\eta}$. Balancing the third and fourth terms yields the proper asymptotic scaling, namely $dS/S\sim[(2\beta-1)/\beta]d\eta/\eta$, and hence $S(\eta\to\infty)\sim\eta^{(2\beta-1)/\beta}$. Therefore, only the range $0<\beta\leq 1/2$ yields bounded non-terminating self similar states.

Boundary conditions also impose constraints on the allowable values of the exponent $\beta$. For example, enforcement of constant liquid volume $V\simeq\int_{Z_{1}}^{Z_{2}}H^{2}dZ=T^{5\beta-2}\int_{Z_{1}}^{Z_{2}}S^{2}d\eta$ is only consistent with $\beta=2/5$. According to Eq. (20), enforcement of a constant flux boundary condition $Q=-\widehat{A}(\alpha/\theta)H^{2}(\partial H/\partial Z)=\textsf{const}=-\widehat{A}(\alpha/\theta)T^{3\beta-2}S^{2}S_{\eta}$ is only consistent with $\beta=3/5$. Clearly then, a constant flux boundary condition (Neumann condition) is therefore inconsistent with bounded solutions.

In what follows, we restrict attention to the value $\beta=1/2$, which accords with the Washburn relation and allows enforcement of time-independent Dirichlet boundary conditions. For this category of solutions, the non-dimensional flux defined in Eq. (14) is represented by

[TABLE]

The self-similar solution $S(\eta)$ then satisfies the equation Weislogel (2001a):

[TABLE]

To ascertain the interface shape of these solutions, we numerically solved Eq. (46) by rewriting the second order equation as a system of first order equations and using the ODE45 solver in Matlab Mat (2015). Shown in Fig. 5 are representative solutions for receding, uniform, advancing and terminating states $S(\eta)$ satisfying the far field Dirichlet conditions shown. According to Eq. (45), solutions with $S_{\eta}>0$ correspond to states with net liquid flux to the left, designated receding states, while solutions with $S_{\eta}<0$ correspond to a net liquid flux to the right, designated advancing states. The solution for which $S\eta$ vanishes everywhere, which corresponds to a zero flux solution in self-similar coordinates, is designed a uniform state. It represents an exception in that it is the only solution which satisfies volume conservation. Solutions whose advancing front are characterized by a vanishing value of $S(\eta)$ are likewise designated terminating states. The numerical solutions indicate that solutions undergo termination only when the interface slope at the origin $S_{\eta}(0)\gtrsim 0.349$.

V.2 Generalized non-modal linear stability of interface self-similar solutions

It is by now well understood that traditional modal stability analysis, wherein instability derives from exponentially growing normal modes of the governing linearized autonomous disturbance operator, is an inappropriate investigatory tool for examining transient stability of systems governed by non-normal operators Farrell and Ioannou (1996). This is because the modal spectrum, and in particular the eigenvalue with maximum real part, properly captures the system response to infinitesimal disturbances at all times only for linearized disturbance operators which are normal; for disturbance operators which are non-normal, it describes the response only at infinite time. The case of non-normal operators therefore requires implementation of a generalized (non-modal) or transient growth stability analysis, as discussed in Refs. [Trefethen et al., 1993; Farrell and Ioannou, 1996] and references therein. The important distinction is that while for normal operators instability arises from a single, fastest growing, exponentially unstable eigenmode which dominates at all times, instability in systems governed by non-normal operators derives from non-modal growth of two or more interacting nonorthogonal eigenmodes at finite times while only asymptotically resolving to the disturbance function characterized by the eigenvalue with maximum positive real part.

Equations of motion which give rise to inhomogeneous base states often yield linearized operators which are non-normal. Such non-normal operators commonly arise in many thin film liquid systems Davis and Troian (2003a); Davis et al. (2003). Given that the base state solutions to Eq. (46) include advancing and receding states which vary with the self-similar variable, we therefore appeal to a transient stability analysis to investigate their behavior. The stability of terminating solutions shown in Fig. 5, however, will not be examined in this work due to the fact that the point of termination requires that additional terms be included in the governing equation of motion to relieve the diverging stress known to occur at a moving contact line Davis and Troian (2003b, 2004); Davis et al. (2006).

To begin, we linearize the general solution according to

[TABLE]

where $S(\eta)$ satisfies Eq. (46) with the given boundary conditions, $G(\eta,\tau)$ represents the non-modal disturbance function and $\delta\ll 1$ denotes the small expansion parameter for linearization. Substituting this form into Eq. (41) yields the governing linear disturbance equation to order $\delta$, namely:

[TABLE]

where

[TABLE]

Here, $\mathcal{L}^{\dagger}$ denotes the adjoint of the linear autonomous operator $\mathcal{L}$, $\tau=\ln T$ and $G(\eta,\tau)=\exp(\mathcal{L}\tau)G(\eta,0)$. Under the $L_{2}$ (Euclidean) norm, the adjoint $\mathcal{L}^{\dagger}$ represents the unique linear operator where $\int_{0}^{\eta_{B}}v(\eta)\mathcal{L}\left\{w(\eta)\right\}d\eta=\int_{0}^{\eta_{B}}\mathcal{L}^{\dagger}\left\{v(\eta)\right\}w(\eta)d\eta$ for all sufficiently smooth functions $v$ and $w$ which are $L_{2}$-integrable and satisfy homogeneous Dirichlet boundary conditions at the boundary points $(0,{\eta_{B}})$. For all non-uniform solutions $S$, $\mathcal{L}\mathcal{L}^{\dagger}\neq\mathcal{L}^{\dagger}\mathcal{L}$ so that $\mathcal{L}$ is non-normal. The transient amplification of disturbances $G(\eta,\tau)$ cannot therefore simply be ascertained from the eigenspectrum of $\mathcal{L}$. Disturbance amplification over an interval in time $\tau$ will instead be quantified by the function

[TABLE]

where the notation $\|\cdot\|^{2}$ represents the Euclidean norm $\left\langle\cdot\,|\,\cdot\right\rangle^{1/2}$ when applied to functions and the operator norm when applied to operators, as in Eq. (52).

To gain further insight into the growth or decay of disturbances at finite times, it is useful to consider the discretized operator $\exp(\mathcal{L}\tau)$ in diagonalized form such that

[TABLE]

where $\exp(\Lambda\tau)$ denotes a diagonal matrix whose entries are the eigenvalues of $\mathcal{L}$ arranged in decreasing order and $\mathcal{S}$ represents the matrix whose columns are the corresponding eigenfunctions. It can then be shown Farrell and Ioannou (1996) that

[TABLE]

where $\beta_{\textrm{max}}$ is the eigenvalue of $\mathcal{L}$ with maximum real part, also known as the spectral abscissa of $\mathcal{L}$. A matrix operator is normal if and only if it is unitarily diagonalizable such that $\|S\|\|S\|^{-1}=1$, in which case the maximum amplification at time $\tau$ over all initial conditions $G(\eta,0)$ is given by $\exp(\beta_{\textrm{max}}\tau)$. For non-normal operators, however, $\|S\|\|S\|^{-1}\geq 1$ and so the maximum disturbance amplification at finite time $\tau$ can exceed the asymptotic value $\exp(\beta_{\textrm{max}}\tau)$, sometimes significantly so. We note too from Eq. (52) that as $\tau\rightarrow\infty$, the quantity $d\ln(\|\exp(\mathcal{L}\tau)\|)/d\tau\rightarrow\beta_{\textrm{max}}$. While a finding that $\beta_{\textrm{max}}<0$ establishes that a system is asymptotically linearly stable, it does not preclude the possibility of large positive transient growth (i.e. linear instability) at finite time.

We now determine the upper bound on the instantaneous rate of disturbance growth at time $\tau$ given by

[TABLE]

The quantity $\omega_{\textrm{max}}$, also known as the numerical abscissa, represents the maximum eigenvalue of the operator sum, $(\mathcal{L}+\mathcal{L}^{\dagger})/2$. Since this combined operator is self-adjoint, its eigenvalues are strictly real. Therefore, while $\beta_{\textrm{max}}$—the eigenvalue of $\mathcal{L}$ with maximum real part—governs the asymptotic stability as $\tau\rightarrow\infty$, it is the value $\textrm{max}\,\{\lambda\}$ which governs the transient stability of the linearized system at finite times. Since physical systems undergo flow instabilities at finite time, it is this latter value which is of most interest. According to this definition then, if it can be shown that $\omega_{\textrm{max}}<0$ for all times $\tau$, then instantaneous disturbance growth will always be suppressed and the system can be deemed linearly stable to infinitesimal perturbations. Below we establish the generalized linear stability of self-similar solutions $S(\eta)>0$ subject to time-independent Dirichlet conditions for volume non-conserving base states describing inertia-less flow in a slender open triangular channel. We first compute analytic bounds on instantaneous disturbance growth for advancing and receding films followed by results obtained from direct numerical computation.

Having introduced the main concepts underlying transient growth analysis, we conclude this section with a few words about the stability of uniform stationary solutions to Eq. (41). Weislogel Weislogel (2001b) previously conducted a modal linear stability analysis for an infinitely long, slender, quiescent liquid filament of uniform thickness $d$ confined to the interior corner of an open triangular channel. He showed that axial sinusoidal perturbations of any wavelength $\ell$ all undergo rapid decay with a dimensional time constant proportional to $\mu\ell^{2}/\gamma d$. In our study, the governing dimensionless linear operator corresponding to a stationary uniform state is given by $\mathcal{L}=(d/2)\,\partial^{2}/\partial Z^{2}$, an operator which is self-adjoint and therefore normal, whose eigenvalues are strictly real, positive and given by $(d/2)(n\pi/L)^{2}$ for $n=1$, 2, etc. According to Eq. 48 (when expressed in the original laboratory frame coordinates $(Z,T))$, the least stable disturbance $\delta G(Z,T)$ will evolve in time according to $\exp[-(d/2)(\pi/L)^{2}T]$ and therefore decay away, restoring the system back to the initial uniform state. This observation confirms Weislogel’s result that uniform quiescent liquid filaments confined to a slender open triangular channel are linearly stable to arbitrary infinitesimal disturbances of any wavenumber. This result holds for any of the four sets of boundary conditions discussed earlier, namely two Dirichlet conditions, one Dirichlet and one Neumann condition, one Dirichlet condition and constant volume or one Neumann condition and constant volume.

V.3 Generalized stability of volume non-conserving self-similar solutions

In what follows, we examine the generalized linear stability of self-similar non-terminating states on the finite domain $0\leq\eta\leq\eta_{B}$. To proceed, we first note that the linear autonomous operator $(\mathcal{L}+\mathcal{L}^{\dagger})/2$ defined in Eq. (49) (which is also self-adjoint and therefore normal) and given by

[TABLE]

satisfies a Sturm-Liouville eigenvalue equation where

[TABLE]

provided $S(\eta)>0$ and $S$, $dS/d\eta$ and $d^{2}S/dS^{2}$ are continuous within the interval $0\leq\eta\leq\infty$. As a result, the eigenvalues $\lambda$ are strictly real and the corresponding eigenfunctions form a complete orthogonal set of basis functions. According to the definitions in Eq. (61) then, if it can be shown that all the eigenvalues of $\mathcal{L}+\mathcal{L}^{\dagger}/2$ are strictly positive, then $\omega_{max}<0$ and disturbance growth is suppressed.

Before proceeding with the behavior of advancing or receding states, we first note that the linear stability of the uniform state $S(\eta)=1$ is easy to confirm. From the definitions given by Eqs. (63) and (64), the relevant operator sum reduces to $(1/2)\,(\partial^{2}/\partial\eta^{2})-(1/4)$, which is self-adjoint and therefore normal and whose eigenvalues are strictly real, positive and given by $(1/2)(n\pi/L)^{2}+1/4$ for $n=1$, 2, etc. According to Eq. (64), the least stable disturbance $\delta G(\eta,\tau)$ will evolve in time according to $\exp{-(1/2)[(n\pi/L)^{2}+1/4T]\,\tau}$ and therefore decay away, restoring the system back to the initial uniform self-similar state.

It is often the case that the Rayleigh quotient Haberman (2004) can be used to establish bounds on growth or decay by helping determine the overall sign of eigenvalues corresponding to Sturm-Liouville operators. In our case, however, this approach yields no useful information. We therefore instead appeal to the following lemma:

Lemma 1**.**

Given a second order Sturm-Liouville eigenvalue equation of the form $(P\,G_{\eta})_{\eta}+Q\,G=-\lambda G$ where $P>0$, $G\neq 0$, $\lambda$ is real, and $G(0)=G(\eta_{B})=0$, then $\lambda_{\textrm{min}}>(-Q)_{\textrm{min}}$,

where subscripts denote differentiation with respect to the domain coordinate $\eta$. In order for $G$ to satisfy homogeneous conditions at the boundary points, it must be the case that there exists an extremum within the domain $[0,\eta_{B}]$ where $G_{\eta}(\eta^{*})=0$. Without loss of generality, the first such extremum from the origin may be assumed to be a local maximum in $G$ for $G>0$. (Note that $G$ can always be made positive at this point by multiplying the local value by -1 and still remain an eigenfunction). Because $PG_{\eta}$ must change sign on either side of this maximum, then at some point in its vicinity, $(PG_{\eta})_{\eta}<0$. This yields the relation $\lambda>-Q(\eta^{*})$, which therefore implies $\lambda_{\min}>-Q(\eta^{*})\geq(-Q)_{\min}$, where $(-Q)_{\min}$ denotes the smallest value of $-Q$ within the domain.

V.3.1 Transient and asymptotic stability of advancing solutions $S(\eta)$

We first examine the stability of advancing self-similar solutions shown in Fig. 5. These solutions are all characterized by a downstream boundary value smaller than the inlet value, namely $S(\eta_{B})=\textrm{const}<S(0)=1$, and an inlet slope which is always negative, namely $S_{\eta}(0)<0$. Here and in what follows, $\eta_{B}$ represents the coordinate where the downstream boundary condition is applied.

Establishing that such states are linearly stable requires a finding that the smallest eigenvalue of $(\mathcal{L}+\mathcal{L}^{\dagger})/2$, namely $\lambda_{\textrm{min}}$, is positive, which then requires

[TABLE]

which in turn requires that $S^{\textrm{min}}_{\eta\eta}>-1$. This minimum value can occur either at some interior point $\eta_{\textsf{int}}$ or at the boundary points $\eta=0,\eta_{B}$, which shall be examined separately. When the minimum value $S^{\textrm{min}}_{\eta\eta}$ occurs at the downstream boundary point $\eta_{B}$, the inequality is easily satisfied for all solutions $S$, whether advancing or receding, since as evident from Fig. 5, all the self-similar solutions asymptote to a uniform thickness, which therefore yields values $S_{\eta\eta}(\eta\rightarrow\eta_{B})\to 0>-1$.

When $S^{\textrm{min}}_{\eta\eta}$ occurs at an interior point $\eta_{\textsf{int}}$, it will then be the case that $S^{\textrm{min}}_{\eta\eta\eta}(\eta_{\textsf{int}})=0$. Differentiation of Eq. (46) then yields the corresponding third order equation evaluated at the point $\eta_{\textsf{int}}$

[TABLE]

whose solutions are given by

[TABLE]

which when substituted into Eq. (46) provides the value of the curvature at that interior point, namely

[TABLE]

The minimal value of this relation is attained for the negative root, which as shown is always less than -1:

[TABLE]

The last inequality derives from the general relation $\eta(\sqrt{\eta^{2}}-\sqrt{\eta^{2}+k})\geq-k/2$ for $\eta$ and $k$ real and positive.

When the minimum value of $S_{\eta\eta}$ occurs at the origin, evaluation of Eq. (46) subject to the Dirichlet condition $S(0)=1$ yields the criterion for stability, namely $S_{\eta\eta}(0)=-2S_{\eta}^{2}(0)>-1$. Since advancing solutions always exhibit negative slopes at the origin, this criterion can be rewritten as $S_{\eta}(0)<-1/\sqrt{2}$. Our numerical results show that in order for advancing solutions to remain strictly positive throughout the domain $[0,L\eta_{B}]$ and not yield a termination point, it must be the case that $S_{\eta}(0)\lesssim-0.349$, which is clearly greater than $-1/\sqrt{2}$. The following geometric argument supports this conclusion as well.

Since $S_{\eta\eta}(0)=-2S_{\eta}^{2}(0)<0$ while $S_{\eta\eta}(\eta_{B})=0$, there must occur at least one inflection point in the domain $[0,\eta_{B}]$. Consider the first such inflection point at position $\eta_{p}$ where $S_{\eta\eta}(\eta_{p})=0$. As shown in Fig. 6, since $S(\eta_{p})$ lies below the extended line $1+\eta_{p}S_{\eta}(0)$, then $S(\eta_{p})<1+\eta_{p}S_{\eta}(0)=1-2S_{\eta}(\eta_{p})S_{\eta}(0)<1-2S_{\eta}^{2}(0)$. These relations are found by noting that $S_{\eta}(\eta_{p})=-\eta_{p}/2$ from Eq. (46) and that $S_{\eta}(\eta_{p})<S_{\eta}(0)$ where both slopes are negative. This then establishes that $S_{\eta}(0)>-1/\sqrt{2}$, which is proof that advancing solutions are linearly stable at all times.

V.3.2 Transient and asymptotic stability of receding solutions $S(\eta)$

We next examine the stability of receding solutions which are characterized by a downstream boundary value larger than the inlet condition, namely $S(\eta_{B})=\textrm{const}>S(0)=1$, and inlet slopes which are always positive, namely $S_{\eta}(0)>0$. In order to make use of the lemma above, we first apply a change of variable to Eq. (64) such that $G(\eta)$ is replaced by the product $S(\eta)G(\eta)>0$. Straightforward calculation yields yet another Sturm-Liouville eigenvalue equation with weighting factor $S^{2}(\eta)$:

[TABLE]

where

[TABLE]

The lemma then applies here, with an adjustment for the weighting factor: $\lambda_{\textrm{min}}>(-Q/S^{2})_{\textrm{min}}$. This then yields the result

[TABLE]

Since for receding solutions all terms in the expression $(\eta S_{\eta}+S)/S$ are strictly positive, then $\lambda_{\textrm{min}}$ is strictly positive and therefore $\omega_{\textsf{max}}$ is strictly negative. This demonstration establishes that receding solutions are transiently and asymptotically linearly stable to any infinitesimal disturbances $G$ satisfying homogeneous Dirichlet boundary conditions.

V.4 Numerical results and comparison to analytic bounds

Our numerical results for the operator exponential governing transient growth are next compared to analytic bounds derived above. Derivatives were constructed using second-order finite differences. (In discrete form, these operators are simply square matrices.) Homogeneous Dirichlet boundary conditions were enforced by directly reducing the operator matrices. The numerical abscissa $\omega_{\textrm{max}}$ was obtained by identifying the smallest eigenvalue of the corresponding reduced operator matrix $(\mathcal{L}+\mathcal{L}^{\dagger})/2$ according to which $\omega_{\textrm{max}}=-\lambda_{\textrm{min}}$. The spectral abscissa $\beta_{\textrm{max}}$ was obtained by identifying the smallest eigenvalue of the reduced matrix $\mathcal{L}$. The quantity $\ln(\|\exp(\mathcal{L}\tau)\|)$, representing the maximum instantaneous disturbance amplification, was obtained from the operator norm $\|\cdot\|$ given by the matrix maximum singular value. Shown in Fig. 7 are the numerical results for these three quantities plotted alongside the analytic bounds for $-\lambda_{\textrm{min}}$ given by $\max\{-(\eta S_{\eta}+S)/(4S)\}$ for receding solutions and $\max\{-(1+S_{\eta\eta})/4\}=-(1+S^{\textrm{min}}_{\eta\eta})/4$ for advancing solutions.

As predicted, the values of $\ln\|\exp(\mathcal{L}\tau)\|$ for all times $\tau$ fall intermediate between the analytic upper bound, namely $\max\{-(\eta S_{\eta}+S)/(4S)\}$ for receding solutions or $-(1+S^{\textrm{min}}_{\eta\eta})/4$ for advancing solutions, and the analytic lower bound given by $-\beta_{\textrm{max}}$. The numerical results also confirm that as $\tau\rightarrow\infty$, the slope of the function $\ln\|\exp(\mathcal{L}\tau)\|$ exactly equals the value $-\beta_{\textrm{max}}$, as must be the case since the asymptotic decay rate of disturbances is dictated by the eigenvalue of $\mathcal{L}$ with maximum real part.

The numerical results in Fig. 7 (and similar studies not shown of other receding and advancing states) reveal that there is no transient nor asymptotic growth of disturbances. Furthermore, the suppression of disturbances as quantified by $\ln\|\exp(\mathcal{L}\tau)\|$ asymptotes to the eigenvalue $-\beta_{\textrm{max}}$ already at early times $\tau<9.0$. The sudden drop off in the value $\ln\|\exp(\mathcal{L}\tau)\|$ observed near $\tau=8$ in Fig. 7(b) simply reflects onset of increased damping of disturbances incurred in the thin precursor film. All results presented were computed using spatial domains $0\leq\eta\leq\eta_{B}$ where $5\leq\eta_{B}\leq 80$ with fixed mesh size ranging from 0.05 to 0.005. As shown in Fig. 8, a domain length of 80 was sufficiently long to ensure numerical convergence irrespective of mesh size.

VI Conclusion

In this work, we have examined the global and linear stability of solutions describing inertia-free flow of a thin wetting Newtonian film within an open slender triangular groove of constant cross-section with constant liquid contact angle. The study is limited to flow states for which the interface function representing the local film thickness is always strictly positive thereby ruling out rupture states. The slender approximation (also known as the lubrication or long wavelength approximation) ensures that to order $\varepsilon^{2}$ the local value of the fluid pressure depends only on the local film thickness whose interface shape in the plane normal to the streamwise flow represents a segment of a circle. This approximation simplifies the governing equation of motion, which is described by a second order nonlinear diffusion equation.

The relevant base states for stability analysis include stationary interface states $H_{S}(Z)$ and self-similar states $S(\eta)$ which adopt the (dimensionless) Washburn scaling $\eta=Z/T^{1/2}$. The stationary states allow for various boundary conditions imposed at the ends of the channel including Dirichlet-Dirichlet, Dirichlet-Neumann (i.e. flux condition), Neumann-constant volume or Dirichlet-constant volume. The self-similar states are volume non-conserving and result from application of Dirichlet boundary conditions. By exploiting an analogy with other gradient flow equations and examining the asymptotic behavior of the corresponding Lyapunov function, it was shown that (strictly positive) stationary interface solutions represent exponentially stable equilibrium points. Disturbances of any type therefore decay away at least exponentially fast to restore the system back to the initial stationary state. The second important result is that advancing, uniform and receding self-similar states satisfying Washburn dynamics are both transiently and asymptotically linearly stable to infinitesimal perturbations. This finding required implementation of a generalized non-modal linear stability analysis since the inhomogeneous nature of the self-similar states naturally gives rise to non-normal disturbance operators.

These two results highlight the reasons why the interface states describing either transient flow (self-similar states) or asymptotic steady state flow (stationary states) within slender open triangular grooves tend to be so stable against perturbations so long as the liquid always wets the channel sidewalls. This result, however, requires that the slender channel approximation be everywhere satisfied, which in turn sets an upper bound on the interface slope, namely $(\partial h/\partial z)^{2}\sim(d/L)^{2}\sim O(\varepsilon^{2})$. Alternatively, this constraint can be viewed as a low-pass filtering requirement which prevents low-frequency modes or disturbances from ever generating high-frequency modes. In particular, neither the base states nor disturbances applied to these base states can trigger streamwise capillary waves. Consequently, the results obtained in this study cannot address effects such as interface leveling in thin films Orchard (1963). As noted by Yang and Homsy Yang and Homsy (2006), flow caused by streamwise curvature can also not be neglected in the limit $(\theta+\alpha)\to\pi/2$.

The model investigated in this work also relies on the assumption that the liquid contact angle $\theta$ is a constant equal to the equilibrium contact angle, an assumption known to break down for sufficiently large flow speeds. A fully numerical model incorporating a velocity dependent contact angle Bracke et al. (1989); Joos et al. (1990) can certainly be developed, but such an extension precludes analytic solutions for base states and therefore development of analytic bounds pertinent to stability. We are currently extending this model in a new promising direction by including effects due to substrate axial curvature. We hope that the stability findings reported in this work as extended to curved substrates will ultimately provide a full analytic and numerical treatment particularly useful to the design of novel propellant management systems for various space applications.

Acknowledgements.

The authors gratefully acknowledge financial support from a 2016 NASA/Jet Propulsion Laboratory President’s and Director’s Fund and a 2017 NASA Space Technology Research Fellowship (NCW).

Appendix A Boundedness of constant $\mathcal{C}$ in Eq. (39)

We first note that the constant $\mathcal{C}$ given by Eq. (39), namely

[TABLE]

satisfies the condition $\mathcal{C}\geq 0$ for $H(Z,T)>0$ and $H_{S}(Z)>0$. We can eliminate the troublesome possibility $\mathcal{C}=0$, however, as follows. Assuming $H_{S}(Z)$ and $H(Z,T)$ remain bounded, then the function in brackets is observed to decrease monotonically with $H$ and $H_{S}$ since the derivatives

[TABLE]

and

[TABLE]

The smallest possible value of the bracketed term in $\mathcal{C}$ is attained when that ratio is evaluated at the maxima of $H_{S}(Z)$ and $H(Z,T)$ (which don’t necessarily occur at the same point $Z$). It then follows that

[TABLE]

where $H$ and $H_{S}$ are evaluated at their respective maxima. Therefore, so long as $H(Z,T)<\infty$, $\mathcal{C}>0$.

A related issue arises in the context of the stability arguments presented in Section IV.2. It could be the case that although the integrated value of the free energy defined in Eq. (32) remains bounded at all times, the transient function $H(Z,T)$ might potentially diverge at finite times depending on the boundary conditions imposed. We eliminate the possibility $H(Z,T)\to\infty$ by appealing to the governing equation of motion, which essentially describes the intrinsic diffusive behavior of the $H^{2}$.

Consider first the behavior of Eq.(41) given by

[TABLE]

in the vicinity of local extrema. Any local maximum of $H$ will satisfy $\partial H/\partial Z=0$, $\partial^{2}H/\partial Z^{2}\leq 0$ and thus $\partial H/\partial T\leq 0$, which leads to a diminishment in height due to the diffusive nature of the underlying equation of motion. Similarly, a local minimum of $H$ cannot further decrease in value. Therefore, any new extrema beyond those present in the initial condition can only be reached by possible extrema at the boundaries. For those solutions satisfying Dirichlet conditions at both endpoints of the domain, $H$ can never therefore attain a new maximum above the maximum value of the initial condition nor attain a new minimum below the minimum value of the initial condition. This behavior then guarantees that so long as the initial condition satisfies the constraint $H(Z,T=0)>0$, then $\mathcal{C}$ will always remain bounded from below.

A Neumann, or equivalently, a flux boundary condition poses a potential problem since too large a flux for a given initial condition could potentially lead to a local drainage spot where the film thickness vanishes. To assess this case, we recast Eq. (19) not in terms of the local variable $H(Z,T)$ but the local flux $Q(Z,T)=-H^{2}\partial H/\partial Z=-(1//3)\partial H^{3}/\partial Z^{3}$ by first multiplying Eq. (19) by $-H/2$ followed by differentiation by $Z$, which yields

[TABLE]

Adopting a similar argument as above, for any case requiring constant liquid volume, the prescribed fluxes at the endpoint must be equal, thereby imposing Dirichlet conditions on the flux $Q$. It then follows that $Q(Z,T)$ can never exceed the extrema present in the initial condition $Q(Z,T=0)$, which therefore precludes $H$ from becoming arbitrarily large; hence $\mathcal{C}$ is again bounded below at all times.

For cases subject to one Dirichlet and one flux boundary condition, the only way in which $H$ might approach infinity is if the flux boundary condition is made very large. However, this would give rise either to an arbitrarily large internal flux extremum, which is disallowed by the previous argument, or an arbitrarily large flux of liquid at the boundary subject to the Dirichlet condition. For the latter case, the divergence in $H$ at one endpoint would be driven by the flux condition at other endpoint, such that the system would approach infinite volume and the Lyapunov functional would therefore not decrease as required. Hence, such a situation cannot occur, and $H$ is therefore bounded above at all times and $\mathcal{C}$ bounded below at all times.

Finally, for the case of a Dirichlet boundary condition coupled with a requirement of constant volume, were $H$ to approach infinity at one boundary, the flux there would also diverge and be matched by identical divergence at the other end point in order to satisfy constant volume. But according to Eq. (84), an interior flux minimum would have to undergo increase thus increasing the minimum slope of $H^{3}$, which would eventually violate the constraint of constant volume. Hence, in this case as well, $H$ cannot diverge at any point in time and $\mathcal{C}$ therefore remains bounded below at all times.

In conclusion, the demonstration above therefore establishes that all stationary states, irrespective of the boundary conditions imposed above, are exponentially stable in the Lyapunov sense.