Steiner triangular drop dynamics

TL;DR

This paper introduces the Steiner triangle as a minimal geometric model for liquid droplet dynamics, revealing complex behaviors such as bouncing, rocking, and quasiperiodic motions, and capturing symmetries of realistic models.

Contribution

It presents a novel minimal model based on Steiner's circumellipse for analyzing droplet dynamics, with detailed dynamical system analysis and identification of invariant manifolds.

Findings

Identified bouncing and rocking periodic motions.

Discovered nested quasiperiodic motions around equilibrium.

Model captures symmetries of more realistic droplet models.

Abstract

Steiner's circumellipse is the unique geometric regularization of any triangle to a circumscribed ellipse with the same centroid, a regularization that motivates our introduction of the Steiner triangle as a minimal model for liquid droplet dynamics. The Steiner drop is a deforming triangle with one side making sliding contact against a planar basal support. The center of mass of the triangle is governed by Newton's law. The resulting dynamical system lives in a four dimensional phase space and exhibits a rich one-parameter family of dynamics. Two invariant manifolds are identified with "bouncing" and "rocking" periodic motions; these intersect at the stable equilibrium and are surrounded by nested quasiperiodic motions. We study the inherently interesting dynamics and also find that this model, however minimal, can capture space-time symmetries of more realistic continuum drop models.

Equations

$$p=\sigma \kappa (\alpha ,\beta ;\ell ),$$

$${\mathbf{F}}_{\mathbf{surf}}$$

$${\mathbf{F}}_{\mathbf{pres}}$$

$$\rho V\frac{d^2\mathbf{x}}{d{t}^2}=\mathbf{F}={\mathbf{F}}_{\mathbf{surf}}+{\mathbf{F}}_{\mathbf{pres}}.$$

$$\kappa =k\left(\frac{1}{\sin \alpha }+\frac{1}{\sin \beta }+\frac{1}{\sin \gamma }\right).$$

$$k=\frac{4{\sin }^2{\alpha }_0}{c_0(4+\sec {\alpha }_0)}=\frac{2{\sin }^2{\alpha }_0\sqrt{\tan {\alpha }_0}}{\ell (4+\sec {\alpha }_0)}$$

$$q({\alpha }_0)\equiv \ell k=\frac{2{\sin }^2{\alpha }_0\sqrt{\tan {\alpha }_0}}{4+\sec {\alpha }_0}.$$

$$\cos \alpha =(x_C-x_A)/b,$$

$$\sin \alpha =y_C/b,$$

$$(x_A,x_B,x_C,y_C)=\left(-\frac{{\ell }^2}{3y},\frac{{\ell }^2}{3y},3x,3y\right),$$

$$a^2$$

$$b^2$$

$$c^2$$

$$\frac{F_x}{\sigma \ell }$$

$$\frac{F_y}{\sigma \ell }$$

$$[\ddot{x},\ddot{y}]$$

$$-3y\left(\frac{1}{a}+\frac{1}{b}\right)+\frac{q({\alpha }_0)}{3y}\left(\frac{2(a+b)}{3y}+ab\right)]$$

$$\equiv [f(x,y),h(x,y)]$$

$$a$$

$$b$$

$$\dot{x}$$

$$\dot{w}$$

$$\dot{y}$$

$$\dot{z}$$

$$\frac{\partial f}{\partial y}\ne \frac{\partial h}{\partial x}$$

$$f(-x,y)=-f(x,y),h(-x,y)=h(x,y).$$

$$G_1:\{t\mapsto -t,(x,w,y,z)\mapsto (-x,w,y,-z)\}$$

$$G_2:\{t\mapsto -t,(x,w,y,z)\mapsto (x,-w,y,-z)\}$$

$$S=G_2\circ G_1:(x,w,y,z)\mapsto (-x,-w,y,z)$$

$$(x,y)=\left(0,\frac{\sqrt{\tan \alpha }}{3}\right).$$

$$\ddot{x}$$

$$\ddot{y}$$

$$=-2\sin \alpha \left(1-\frac{q({\alpha }_0)}{q(\alpha )}\right).$$

$$q(\alpha )=q({\alpha }_0).$$

$${\mathbf{x}}_0=(x_0,y_0)=\left(0,\frac{\sqrt{\tan {\alpha }_0}}{3}\right).$$

$$\frac{486y^5}{\left(81y^4+1\right)\left(\sqrt{81y^4+1}+4\right)}=q({\alpha }_0)$$

$${\alpha }_0^{\ast }\approx 1.391\approx 79.7^{\circ }.$$

Full text

Polygons are the simplest of planar shapes and triangles the simplest of all polygons. Here, we introduce a triangle with variable base length and height as a minimal model of a sliding and deforming liquid droplet, inspired by MJ Steiner’s (French geometer) classic theorem relating triangles to ellipses. The Steiner drop’s shape is related to Newtonian dynamics through a constitutive relationship that penalizes against sharp angles. The resulting dynamical system with a four dimensional phase space exhibits fixed points, periodic and quasiperiodic orbits, all related to two competing invariant manifolds. The symmetries of the dynamical system allow us to characterize the Steiner drop motions near equilibrium as “bouncing,” “rocking,” or quasiperiodic combinations of the two. By associating each triangle with its unique smallest circumscribed ellipse, the motions of the Steiner drops can be considered equivalently as motions of elliptical drops. To relate to more realistic drops, we also compute the center-of-mass motions of normal modes of perturbed spherical cap liquid drops and find that they also can be classified as bouncing or rocking. Thus the minimal model predicts the qualitative behavior of spherical-cap drop oscillations.

1 Introduction

We introduce the Steiner triangular model of a sliding and deforming liquid drop shaped by surface tension. In nature, of course, surface tension does not allow for corners. This deficiency is overlooked in favor of the simplicity of the triangle, whose deformations have just two degrees of freedom. The base can lengthen and shorten, consistent with a moving contact line, and the apex can move about in the plane. We seek a reduced order (ODE) model that captures the space-time dynamics of higher order (PDE) continuum descriptions. Our model is termed minimal since, on removing one degree of freedom, the model becomes trivial dynamically.

The triangular geometry can be regularized to an ellipse. Every triangle is uniquely associated with a circumscribed ellipse whose center coincides with the centroid of the triangle and whose area is minimal, the so-called ‘Steiner ellipse’ [28, 32]. Conversely, every ellipse has infinitely many maximal inscribed triangles, but only two of these have a side parallel to a given line, and these are reflections of each other. Thus each ellipse can be uniquely associated with a maximal inscribed triangle whose base is horizontal, which we refer to as its ‘Steiner triangle’ (not to be confused with the usage of the term ‘Steiner triangle’ as the Civian triangle of a Steiner point). This one-to-one mapping between Steiner triangles and ellipses yields an equivalence between triangular and elliptical drops. Motions can be equivalently viewed as those of triangles or of ellipses.

To complete the definition of the Steiner triangular drop, one needs to introduce a constitutive relationship between force and deformation that mimics surface tension. While we accept corners as natural to triangles, we do assess a penalty on the sharpness of the corners. Sharply acute angles are penalized by a pressure that blows up in the zero angle limit. This can be thought of as the analog of the Young-Laplace relationship, familiar from the continuum model for normal stress balance across an interface. For a minimal model, the exact nature of the force-deformation behavior is not important as long as there is a robust penalty for flattening the triangle.

Our interest is in the dynamical regime where liquid inertia is shaped by surface tension and dissipation can be ignored, the so-called capillary ballistic regime where the dynamics are conservative. Sources of dissipation that are neglected are bulk dissipation due to the inherent viscosity of the liquid and possibly contact line dissipation due to drag of the moving contact line [6, 7, 3]. Experiments with water drops on substrates fall into this regime [25, 24, 4], and PDE models have been successful at predicting observed frequencies of vibration [1, 16, 17, 10].

For conservative dynamics, the motion of the drop’s center of mass decouples from its deformation relative to the inertial frame. In this paper, we restrict our study to motions in the inertial frame, consistent with a focus on a drop’s natural vibrations. This eliminates the translational degree-of-freedom introduced by the sliding motion. Motions relative to the center of mass are solved for and motion of the contact points, $A$ and $B$ in Figure 1(a), relative to the substrate are preserved up to a (Galilean) translational velocity.

The relevant configurational variables are the coordinates of the triangle’s center of mass relative to its base. Newton’s law is applied to each of these, which leads to a first-order dynamical system of nonlinear ODEs of dimension four. Having a four dimensional phase-space means that, a priori, our minimal model can possibly exhibit long-term dynamics that includes fixed points, periodic, and quasi-periodic orbits. Surprisingly, this richness is actually realized.

Liquid drops that deform and translate along a solid support are found widely in nature and in application [13, 12]. The contact line represents the locus where liquid, solid and gas phases meet and is characteristic of partially wetting liquids such as water on a silicon wafer. Drop translation requires contact lines to move, and rapidly moving contact lines remain a modeling challenge in the continuum description [26, 29]. In our $2D$ model, contact lines become contact points, the points $A$ and $B$, Figure 1(a). These points move independently according to the capillary forces acting on them and the Newtonian inertial dynamics of the center of mass.

From an application viewpoint, there is a growing need to simulate large populations of droplets that are translating, deforming and colliding, splitting up (atomizing) or perhaps joining up (coalescing) [30, 18, 34]. Having a minimal model that captures contact line motions with some fidelity provides a computationally economical path to studying the dynamics of large populations.

Vibrated drops of water exhibit capillary ballistic motions associated with a series of mode shapes and resonant frequencies [5]. The shapes and frequencies, or spectrum, of the free liquid sphere modeled as a continuum, were first reported by Rayleigh [20]. By symmetry, vibrations the hemispherical drop maintaining a $\pi/2$ contact angle against a planar substrate are part of Rayleigh’s family of solutions. Studies of the dynamical response of the hemispherical sessile drop subject to other contacting conditions are more recent [19, 24, 33]. Analyses have been extended to spherical cap drops both with pinned and moving contact lines[3]. The space-time symmetries represented are non-trivial. In this paper using the center-of-mass metric, we show that all these responses map into one of three classes of motions: steady, bouncing or rocking. Our Steiner drops faithfully predict these classes.

This paper is organized as follows. In Sections 2-3 we introduce the model, derive the equations of motion, and convert to center-of-mass coordinates. The resulting dynamical system and its symmetries are described in Section 4. In Sections 5-6 we identify the equilibria of the system and analyze their stability and bifurcation. In Section 7 we use the symmetries of the system to identify two-dimensional invariant manifolds filled with “bouncing” and “rocking” periodic orbits. Section 8 describes the torus trajectories in the region surrounding the stable equilibrium. Finally, in Section 9, we show that for small perturbations, the center-of-mass motion of sessile spherical-cap drops can be characterized as “bouncing,” “rocking,” or a combination of those. Thus the triangular drop model qualitatively predicts motions of fluid drops.

2 The model

The triangle $ABC$ sits with its base on the $x$ axis which models the planar support. The triangle has vertices $A,B,C$, with corresponding interior angles $\alpha,\beta,\gamma$ and opposite edges $a,b,c$, respectively. Side lengths are denoted by $a,b,c$: Figure 1 (a). We shall speak of a constant volume $V$ and a length scale $\ell\equiv V^{1/3}$ even though, for our Steiner triangles, it is an area $\ell^{2}$ which is preserved to model liquid incompressibility. Alternatively, one can think of the Steiner drops as prismatic with a uniform depth of dimension $\ell$. In either case, one ends up with a description of a deforming $2D$ object with physical units that correspond to a $3D$ drop, for convenience in formulation, and without loss of generality.

The instantaneous configuration of the triangle is determined by two angles and the length scale $\ell$. In order to apply Newton’s law, the instantaneous forces acting on the object must be identified.

This is done using the free-body diagram method, Figure 1 (b), where we imagine a planar section just above the support which cuts the interface and is ‘shrink-wrapped’ around the rest of the triangle. The net force on the triangle thus has two contributions, one from the surface tension $\sigma$, acting tangent to the liquid/gas interface at the contact line, $\mathbf{F_{surf}}$, and one from the pressure $p$ acting plane-normal to the support, $\mathbf{F_{pres}}$. The surface tension contribution is of constant magnitude, with varying components determined by the contact angle, $\alpha$ or $\beta$ via the Young-Dupre equation at the points $A$ and $B$. The pressure contribution is also proportional to surface tension but depends on the global configuration in a way that needs to be modeled,

[TABLE]

as described below. These contributions are analogous to those exerted on a small supported drop in a full three dimensional continuum description [3]. Thus we have

[TABLE]

The net force gives the acceleration relative to the inertial frame of the center of mass of the triangle, $\mathbf{x}$, according to Newton’s law:

[TABLE]

Inertia, proportional to $\rho V$, and capillarity, proportional to $\sigma$, compete to determine the dynamics. Note that the support is fixed in the lab frame and that, relative to the lab frame, the net force is non-zero, in both directions. In contrast, the inertial frame is attached to the droplet center of mass, a frame in which the net force is zero and with respect to which there is no relative motion.

Our guide to pressure closure is the Laplace pressure across a curved surface, which penalizes against small radii of curvature. Analogously, we choose a functional form for the curvature $\kappa$, and hence pressure (1), that penalizes against small internal angles:

[TABLE]

Here $k$ is constant, determined as now described. For the fully $3D$ droplet, equilibrium requires that the integral of the horizontal force component around the closed curve representing the contact line vanishes. The corresponding condition here is that the sum of horizontal force components at vertices $A$ and $B$ be zero. This implies that, at equilibrium, $\alpha=\beta\equiv\alpha_{0}$ and $c\equiv c_{0}$, the length of side $c$. The equilibrium contact angle $\alpha_{0}$ is readily measured in real drops, and makes a good point of comparison for the model.

All equilibria, $(\alpha,\beta,c)=(\alpha_{0},\alpha_{0},c_{0})$ are evidently isosceles triangles. Using the identity $\alpha+\beta+\gamma=\pi$ to eliminate $\gamma$ in (5) and equating the vertical force component to zero in (4) as required by equilibrium, one finds,

[TABLE]

where the second expression uses $c_{0}=2\ell\sqrt{\cot\alpha_{0}}$, by the area constraint. Thus we obtain a dimensionless parameter

[TABLE]

which is shown in Figure 2. Notice that each value of $q$ is attained at two values of $\alpha_{0}$, except its maximum value, which we denote by $q(\alpha_{0}^{*})$.

3 Center of mass description

To make Newton’s law into a dynamical system, we express the net force as a function of $(x,y)$, the Cartesian coordinates of the triangle’s center of mass. Let the vertices have Cartesian coordinates $(x_{A},0)$, $(x_{B},0)$, and $(x_{C},y_{C})$ respectively. Then we have

[TABLE]

The origin is put at the center of the base, stationary in the lab frame. Let $d$ be the half-width of the base, so $x_{A}=-d$ and $x_{B}=d$. Trigonometry and volume constraint lead to

[TABLE]

and the lengths of the three sides in terms of $x$ and $y$,

[TABLE]

The components of the net force, $\mathbf{F}=(F_{x},F_{y})$, can now be written:

[TABLE]

A dimensionless form of Newton’s law (4) follows from (13) and (14) by scaling lengths by $\ell$, and time by the inertial time, $t=\hat{t}\sqrt{\rho V/\sigma}$.

4 Dynamical System

Rewriting (4) using (13) and (14) yields a dimensionless second-order ODE for $\hat{x},\hat{y}$ in $\hat{t}$. We drop all hats and let dots represent time derivatives to record the resulting dynamical system,

[TABLE]

where the non-dimensionalized side lengths are

[TABLE]

and $q(\alpha_{0})$ is defined in Equation (7).

It is useful to convert the two second-order equations to four first-order equations:

[TABLE]

The system is not Hamiltonian using the natural generalized coordinates $\mathbf{q}=(x,y)$, $\mathbf{p}=(w,z)$, since

[TABLE]

However, notice that $f(x,y)$ is odd and $h(x,y)$ is even with respect to $x$:

[TABLE]

This means that the system is reversible, and in particular it is invariant [21] under the time-reversing phase space involutions

[TABLE]

and

[TABLE]

and hence also under the symmetry [14]

[TABLE]

Although Hamiltonian systems need not be reversible or equivariant, and vice versa, [14] these classes of systems share many dynamical properties. In particular, there are reversible analogues of the Lyapunov center theorem [8] and KAM theory. [22, 23]

These results allow us to characterize the orbits of our system near the neutrally stable equilibrium, discussed below. Our system is time-reversible and hence Lyapunov stability is the relevant definition; below we shall simply refer to Lyapunov stable states, as ‘stable’.

5 Equilibria

In analyzing the equilibria of the drop, it is easiest to write the components of the net force in terms of the contact angles. By the symmetry of the system, any equilibrium configuration must be an isosceles triangle, i.e. $\alpha=\beta$. By trigonometry, the scaled (dimensionless) center of mass of such a triangle is

[TABLE]

Substituting these values into Equation (15) gives

[TABLE]

Thus any equilibrium configuration is an isosceles triangle with contact angle $\alpha$, where

[TABLE]

Since each value of $q(\alpha_{0})$, except the maximum, is attained at two values of $\alpha_{0}$, this means that in general there are two equilibria.

The equilibrium with contact angle $\alpha_{0}$ is given by

[TABLE]

We will denote the other equilibrium by $\mathbf{x}_{1}=(0,y_{1})$. Using Equations (7) and (28) in Equation (31), we see that $y_{1}$ is the second real, positive solution of

[TABLE]

There is no closed-form expression for $y_{1}$, but we can find it numerically for a given value of $\alpha_{0}$.

The parameter $q$ attains its maximum at the critical value

[TABLE]

For values of $\alpha_{0}$ below $\alpha_{0}^{*}$, $\mathbf{x}_{0}$ is the lower equilibrium. At $\alpha_{0}=\alpha_{0}^{*}$ the two equilibria coincide, and for $\alpha_{0}$ above $\alpha_{0}^{*}$, $\mathbf{x}_{1}$ is the lower equilibrium. See Figure 3.

6 Stability of equilibria

The Jacobian of the system (19)-(22) is given by

[TABLE]

At $x=0$, we find that $f_{y}=h_{x}=0$, so the eigenvalues at either equilibrium are

[TABLE]

As predicted by the reversibility of the system, [14] for all values of $\alpha_{0}$ the eigenvalues come in pairs $\{\lambda,-\lambda\}$.

Using Equation (33), the partial derivatives can be written in terms of the $y$-coordinate of the equilibrium as

[TABLE]

where $y=y_{0}$ or $y_{1}$.

At $\mathbf{x}_{0}$, we have in terms of $\alpha_{0}$

[TABLE]

For $\alpha_{0}<\alpha_{0}^{*}$, all four eigenvalues of the Jacobian at $\mathbf{x}_{0}$ are pure imaginary, so the equilibrium of the linearized system is (Lyapunov) stable. By the symmetry of the system, $\mathbf{x}_{0}$ is a nonlinear center.

At $\alpha_{0}=\alpha_{0}^{*}$, $\lambda_{3}$ and $\lambda_{4}$ pass through 0 and become real, so the equilibrium becomes a saddle. See Figure 4.

This change of stability occurs in a transcritical (exchange of stability) bifurcation with $\mathbf{x}_{1}$. See Figure 5, in which the $y$ coordinate of the equilibria is given by Equation (33).

Notice that invariance of the dynamics under $S$ implies that the two-dimensional manifold $x=w=0$ is invariant. On this manifold, we have the reduced system

[TABLE]

Taking $x=0$, we can plot a phase portrait of the reduced system Equation (41)-(42) in the $(z,y)$ plane. Below $\alpha_{0}^{*}$, the prescribed equilibrium $\mathbf{x}_{0}$ is stable and surrounded by an island of periodic orbits, while $\mathbf{x}_{1}$ is a saddle; at $\alpha_{0}^{*}$, there is a single equilibrium, a degenerate saddle; and above $\alpha_{0}^{*}$, $\mathbf{x}_{0}$ is a saddle, while $\mathbf{x}_{1}$ is stable. Figure 6 shows the phase portrait for $\alpha_{0}=1.2$. Note that, by the phase flow, initial conditions outside the island escape to infinity. Escape may anticipate the droplet jumping that occurs for real droplets [9, 15]. For $y=y_{0}$ there are critical velocities and for $z=0$, critical displacements, beyond which escape occurs. ‘Escape displacements’ correspond to sufficiently flattened as well as stretched Steiner drops. That these can escape was surprising to us. See Supplementary Material for an animation.

Physically, this bifurcation shows that above a critical contact angle, the stable drop configuration does not adopt the prescribed contact angle at all, but a lower one.

7 Symmetry and invariant manifolds

At the stable equilibrium, the linearized system has two pairs of pure imaginary eigenvalues. The reversible Lyapunov center theorem [8, 14] asserts that for each of these pairs, there is a one-parameter family of nested symmetric periodic orbits. Therefore, we look for two 2-dimensional invariant manifolds containing these families.

We have already observed the first such manifold: as described in the previous section, there is an invariant manifold $x=w=0$. This is easily seen, since $f(x,y)$ is odd in $x$, so taking $x=0$ reduces the dynamical system to Equations (41)-(42). Note that this reduced system is Hamiltonian: if $U(y)$ satisfies $U^{\prime}(y)=-h(0,y)$ then we have the Hamiltonian

[TABLE]

which satisfies

[TABLE]

and is constant in time:

[TABLE]

Therefore a region around the center $\mathbf{x}_{0}$ in this invariant manifold is filled with periodic orbits, corresponding to a vertically oscillating “bouncing” mode of the drop. We will refer to the manifold $x=w=0$ as the bouncing manifold for this reason.

The second type of invariant 2-dimensional manifold is less obvious. We find that there are invariant manifolds containing the equilibria, where $y=g(x,w)$. For most values of $\alpha_{0}$, $g(x,w)$ is relatively flat, so the trajectories in this manifold exhibit a near-horizontally oscillating “rocking” motion. The reversible KAM theory[22] predicts that a region around the center in this manifold is filled with periodic orbits. We will call these orbits the rocking mode of the drop, and this type of invariant manifold the rocking manifold. See Supplementary Material for animations of bouncing and rocking Steiner drops.

To find this invariant manifold, we let $y=g(x,w)$, so the dynamical system reduces to

[TABLE]

By the chain rule, we have

[TABLE]

and

[TABLE]

where $f$ and its derivatives are evaluated at $(x,g(x,w))$.

By hypothesis,

[TABLE]

We express $g$ as a power series,

[TABLE]

and expand $f(x,g(x,w))$ and $h(x,g(x,w))$ as series in $x$ and $w$. Then, we set Equations (7) and (52) equal and collect like terms, and thus solve for the coefficients $a_{i}$.

The linear and quadratic coefficients are

[TABLE]

See Figure 7. Notice that these have singularities at $\alpha_{0}=\alpha_{0}^{*}$, which is to be expected, and at $\alpha_{0}^{\dagger}\approx 0.870\approx 49.8^{\circ}$, which is surprising.

Computing the coefficients of higher-degree terms, we find that the cubic terms vanish, as do the $xw^{3}$ and $x^{3}w$ terms. The coefficients of the $x^{4}$, $x^{2}w^{2}$, and $w^{4}$ terms are cumbersome but can be computed numerically. Like the quadratic coefficients, the quartic coefficients have singularities at $\alpha_{0}^{*}$ and $\alpha_{0}^{\dagger}$; there is also a singular value $\alpha_{0}\approx 0.517\approx 29.6^{\circ}$. See Figure 8.

It makes sense that $g$ should be even with respect to both variables. The time-reversing involutions $G_{1}$ and $G_{2}$ applied to the reduced system give

[TABLE]

and

[TABLE]

respectively. If $g$ is even with respect to both variables, the system constrained to the invariant manifold inherits the full system’s invariance under $G_{1}$ and $G_{2}$.

Away from the singular values of $\alpha_{0}$, the symmetry of the reduced system Equations (48)-(49) lying in the manifold $y=g(x,w)$ shows that a region around the origin in this manifold is filled with periodic orbits.

We use this same method to find the invariant manifold containing the second equilibrium $\mathbf{x}_{1}$, for values of $\alpha_{0}$ above the transcritical bifurcation at $\alpha_{0}=\alpha_{0}^{*}$. Note that since there is no closed-form expression for $\mathbf{x}_{1}$, we must find the invariant manifold numerically for a specified value of $\alpha_{0}$. For example, the result when $\alpha_{0}=1.45$ is

[TABLE]

8 Numerical solutions: Torus trajectories

The reversible KAM theory [22] says that in the $2n$-dimensional phase space of a smooth dynamical system that is reversible with respect to an involution which fixes an $n$-dimensional submanifold, for each $0\leq m\leq n$ there is an $m$-parameter family of invariant $m$-tori. The union of these tori has positive $2m$-dimensional measure.

In particular, in our 4-dimensional phase space with two involutions which each fix a 2-dimensional subspace, there is a 0-dimensional manifold of 0-tori (the stable equilibrium); 2-dimensional manifolds of 1-tori (the periodic orbits discussed in the previous section) and a 4-dimensional region of 2-tori. The quasiperiodic frequencies of trajectories on these 2-tori are strongly incommensurable, [22] meaning that each such trajectory is dense in its invariant torus.

Physically, these torus trajectories combine both vertical and horizontal motion. They can be visualized in three dimensions as lying on tori where one direction of revolution is in the $y$-$z$ plane, and the other is in a traveling $x$-$w$ plane. See Figure 12.

To show the trajectories in this way, we parameterize them via

[TABLE]

Several torus trajectories are shown in the Appendix. See Supplementary Material for an animation of this.

Figures 16, where $\alpha_{0}=\pi/4$, and 18, where $\alpha_{0}=2\pi/5$, show typical sets of trajectories. In all of these, the initial condition is perturbed from $(x_{0},y_{0})$ by a distance of 0.05. The figures show a range of angles of the perturbation, from vertical (angle $\phi=0$), lying in the bouncing manifold, to lying in the $y=g(x,w)$ rocking manifold).

Figure 17 shows a set of trajectories where $\alpha_{0}=\alpha_{0}^{\dagger}\approx 0.870\approx 49.8^{\circ}$ is the critical value at which the invariant manifold has a singularity. Again, the initial condition is perturbed from $(x_{0},y_{0})$ by a distance of 0.05. We see that no periodic orbit corresponding to $y=g(x,w)$ can be found.

Figure 19 shows a set of trajectories where $\alpha_{0}=1.45>\alpha_{0}^{*}$. In this case the stable equilibrium is $\mathbf{x}_{1}$, for which we have no closed-form expression. Thus the invariant manifold must be computed numerically. Again, the initial condition is perturbed from $(x^{*},y^{*})$ by a distance of 0.05.

9 Comparison to continuum spherical-cap drop model predictions

Continuum models consider deformations with infinite degrees of freedom. Studies of these models typically solve, in some limit, the Navier-Stokes equations with boundary conditions that accommodate moving contact lines [29]. Bostwick and Steen [2] perform a linear stability analysis of a 2-parameter family of spherical-cap base states to find the resonant mode shapes and frequencies in the Euler equation limit. These continuum predictions have guided Chang et al.[4] to the laboratory discovery of the first 35 modes. The space-time symmetries represented by these various solutions are non-trivial [5, 27]. In this section we post process the mode shapes presented by Bostwick and Steen [2] to find their center-of-mass dynamics in order to compare with the bouncing and rocking modes reported above.

A static spherical-cap drop with equilbrium contact angle $\alpha$ is scaled so that the contact line has radius 1; thus the radius of the unperturbed drop is $\csc(\alpha)$. Normal modes for perturbations of the surface are parameterized by azimuthal wavenumber $l$ and frequency $\Omega$, which depends on the polar wavenumber $k$:

[TABLE]

where $s$ is the polar angle. See Figure 13. This ansatz holds for all $(k,l)$ except $(1,1)$, in which case it is inappropriate since there is instability [2]. To compute the normal modes, Bostwick and Steen reduce the dynamic pressure balance to an operator equation on a function space satisfying Laplace’s equation and a no-penetration boundary condition. The normal modes are the eigenfunctions of this operator equation, given as sums of the harmonic basis functions[2].

To calculate the center of mass of the perturbed drop, we integrate over two regions, taking the origin to be the center of the unperturbed drop: the cone whose base is the wetted disk and apex is the origin, and the sector with $0\leq s\leq\alpha$.

Note that if $\alpha<\pi/2$ the cone volume will be negative. Also, no domain perturbation is necessary for $\epsilon\ll 1$, the errors introduced will be higher order, that is, $O(\epsilon^{2})$.

We compute

[TABLE]

where $M$ is the volume of the drop

[TABLE]

and the integrals of $x$, $y$ and $z$ over the drop are similarly taken in spherical coordinates.

Evaluating the inner integral, we find that there are terms whose denominators vanish at $l=0$ and $l=1$. We treat these cases separately. As the integrand is a bounded function of $l$, we use the dominated convergence theorem to justify taking limits inside the integral using L’Hôpital’s rule.

9.1 L=0 modes

We find that for $l=0$

[TABLE]

but the $O(\epsilon)$ term integrates to 0 when $\xi(s)$ is one of the eigenfunctions computed by Bostwick and Steen[2]. The unscaled coordinates of the center of mass are

[TABLE]

which tells us that the center-of-mass motion of the $l=0$ modes is vertical, corresponding to the bouncing mode of the triangular drop. These integrals can be computed numerically, taking $\xi(s)$ to be an eigenfunction. For example, see Figure 14, and compare to Figure 11 (a).

9.2 L=1 modes

For $l=1$ (and $k\neq 1$) the volume and unscaled center-of-mass coordinates are

[TABLE]

which shows that the center-of-mass motion of the $l=1$ modes is mostly horizontal, and confined to a vertical plane determined by the phase of the $\phi$ term. We can draw a qualitative correspondence to the rocking mode of the triangular drop. A numerically integrated trajectory, taking $\xi(s)$ to be an eigenfunction, is shown in Figure 15. Compare to Figure 11 (b).

9.3 L$>$1 modes

For all other whole-number values of $l$ there are no $O(\epsilon)$ terms in the center-of-mass motion:

[TABLE]

Thus for $l>1$ the center of mass of the sessile drop is stationary to $O(\epsilon)$.

We conclude that the triangular drop model in fact qualitatively predicts all of the $O(\epsilon)$ pure-mode motions of the fluid drop’s center of mass. In general, the motion of a fluid drop is a linear combination of modes, each with its own frequency. A combination of one $l=0$ and one $l=1$ mode may be identified with the torus trajectories of the triangular drop.

10 Conclusion

We have introduced the sliding Steiner triangular drop in the spirit of a minimal model. We are unaware of any prior studies of triangular drop models. The dynamical system is surprisingly rich, organized by two 1-parameter families of fixed points, with the stable equilibrium at the intersection of two invariant manifolds. As the parameter varies, the fixed points exchange stabilities at a critical value where the “rocking” manifold folds up into the “bouncing” manifold. Evidence presented pre- and post-collision suggests a complicated transition, whose exact nature remains an open question. Off the manifolds, but in a nested structure, exist 2-tori where quasi-periodic dynamics occur.

To our surprise, the model captures space-time symmetry-breaking motions observed in real drops and may anticipate nonlinear behavior like torus dynamics and droplet jumping, currently under active study for real drops [31].

The dynamical system resides in a four dimensional phase space. Although not Hamiltonian, the system is invariant under two time-reversing phase-space involutions and by these spatial symmetries the system falls under a theorem of Sevryuk that predicts 2 families of nested $m$-tori up to $2m=4$, where $2m$ is the dimension of the full phase space and $m$ is the dimension of the subspace fixed by each involution. That is, for our $4D$ phase space, nested 0-tori (fixed points), 1-tori (periodic orbits), and 2-tori are predicted. These are identified. An open question is whether chaotic orbits also exist: to our knowledge, the reversible KAM theory does not rule them out[23] but we have not observed them.

The dynamical system depends on a single parameter, the rest state contact angle $\alpha_{0}$. For every stable fixed point in this 1-parameter family of isosceles triangles, there is an unstable (saddle) fixed point. In a bifurcation diagram, these two families of fixed points cross transcritically at $\alpha_{0}=\alpha_{0}^{*}$. We identify two 2-dimensional invariant manifolds $\{x=w=0\}$ and $\{y=g(x,w),z=dg(x,w)/dt\}$ where the stable fixed points are surrounded by periodic orbits. These manifolds are flat and nearly flat and constitute the bouncing and rocking modes of the drop, respectively. A region around these periodic orbits in phase space is filled with nested invariant tori, which correspond to quasiperiodic orbits. We observe several of these numerically.

The rocking manifold is apparently not smooth for all values of $\alpha_{0}$ as higher-order approximations reveal singular points where the power-series approximation fails. The nature of this manifold near these points remains an open question.

Here we remind the reader of the spirit of our minimal model. Most engineering studies seek models that can make detailed predictions of observables. We are not seeking such a model, rather a model with the fewest degrees of freedom that can predict qualitative features of the motions of droplets with moving contact lines. Our Steiner triangular drop predicts bouncing and rocking motions. Our model is minimal in the sense that, if any one degree of freedom is removed, our Steiner drop no longer functions as a drop. For example, pinning either basal vertex confines the apex to a hyperbola while pinning the apex fixes the base length owing to the area constraint. In both cases, the dynamics become trivial. Like removing one leg of a three-legged stool, removing one degree of freedom from our model leaves it non-functional.

As a minimal model, adding features goes against its minimalist nature. If one chooses nonetheless to do so, extensions in any number of directions are possible. Contact line damping or bulk viscous dissipation could be added. Extensions to three dimensions, say, to polyhedral ‘drops’, would add deformational degrees of freedom. Tetrahedral drops might be related to the Steiner ellipsoid, if that aspect of the model is of interest. The pressure constitutive relationship we introduce is purposely ad hoc and one could postulate other such relationships. To study homogeneous (same $\alpha_{0}$) populations of interacting Steiner drops, one might add coalescence, jumping or ejection behavior depending on the phenomenon of interest. Alternatively, for heterogeneous (different $\alpha_{0}$) populations, disparities in size and/or chemical make-up might be included.

Supplementary material

See the supplementary material for videos of Steiner triangular drops: bouncing, rocking, torus trajectory, and unbounded trajectory.

Acknowledgments

ENW acknowledges support from Cornell University as the inaugural CAM multidisciplinary postdoctoral fellow. PHS acknowledges support by NSF Grant CBET-1637960.

Appendix: Torus trajectory figures