Steering droplets on substrates with plane-wave wettability patterns and deformations

Abstract

Motivated by strategies for targeted microfluidic transport of droplets, we investigate how sessile droplets can be steered toward a preferred direction using travelling waves in substrate wettability or deformations of the substrate. To perform our numerical study, we implement the boundary-element method to solve the governing Stokes equations for the fluid flow field inside the moving droplet. In both cases we find two distinct modes of droplet motion. For small wave speed the droplet surfs with a constant velocity on the wave, while beyond a critical wave speed a periodic wobbling motion occurs, the period of which diverges at the transition. These observation can be rationalized by the nonuniform oscillator model and the transition described by a SNIPER bifurcation. For the travelling waves in wettability the mean droplet velocity in the wobbling state decays with the inverse wave speed. In contrast, for travelling-wave deformations of the substrate it is proportional to the wave speed at large speed values since the droplet always has to move up and down. To rationalize this behavior, the nonuniform oscillator model has to be extended. Since the critical wave speed of the bifurcation depends on the droplet radius, this dependence can be used to sort droplets by size.

---

1 Introduction

Liquid droplets on a solid substrate occur naturally during, e.g., rain and condensation^1,2 or artificially in printing of ink³ or medical testing.⁴ Depending on their interaction with the substrate and gas phase, they either sit on a substrate, spread, contract, or even move laterally along the substrate.⁵ Lateral motion can be induced by a variety of mechanisms; through nonuniform electric fields (electrowetting),⁶ gravity,⁷ Marangoni advection,^8,9 gradients in wettability,^4,10,11 or deformations of the substrate itself.^12,13 Here, we focus on the latter two mechanisms.

In our previous work, we used the gradient in a wettability step to move a droplet forward using a feedback loop, which synchronizes step and droplet motion.¹⁴ However, due to surface imperfections at which droplets get pinned, such a motion can be difficult to realize in an experiment.⁵ Yet, pinning can be avoided by a large wettability gradient,¹⁵e.g., induced by structural changes of the substrate,¹⁶ or by using homogeneously clean and flat surfaces.¹⁷

Deformations of the substrate itself also move a droplet.¹⁸ In ref. 19 it has been demonstrated that capillary forces attract droplets towards regions with large inward curvature. An example of this are pipettes. It is possible to externally control such substrate deformations, e.g., using gels that reversibly swell and unswell in response to light with specific wavelengths.^13,20,21 Already more than three decades ago wetting of curved rigid substrates was investigated.²² In more recent work, elastic substrates were studied that deform in contact with a liquid droplet.^23–25 Internal tension builds up as a substrate deforms, which changes its surface tension what is known as Shuttleworth effect.^26,27 Here, we just consider curved rigid substrates, where this effect does not occur.

In this article we compare two mechanisms for inducing droplet motion by imposing the spatio-temporal pattern of a travelling wave onto the substrate; first in wettability and thereafter in the height profile of the substrate. This enables to steer the droplet into a prescribed direction; a typical task in lab-on-a-chip applications.²⁸ For possible experimental realizations of travelling waves in wettability and deformations, we refer the reader to Sections 4 and 5, respectively. Because the travelling wave is spatially periodic, the droplet is continuously driven out of equilibrium and cannot settle into a sessile state. Instead, the droplet can be considered as a driven nonuniform oscillator,²⁹ as we show, the driving force of which depends on its position relative to the travelling wave. In Fig. 1 we display two series of snapshots that illustrate the two possible modes of motion, which we term wobbling (left in Fig. 1) and surfing (right in Fig. 1), respectively. In ref. 14 and 30 we applied the boundary element method (BEM)³¹ to a wetting droplet on a plane substrate with a spatio-temporal wettability pattern. To describe deformation waves of the substrate, we needed to extend our approach to substrates with a spatio-temporal curvature pattern. However, we also needed to completely reformulate our formalism in order to guarantee conservation of the droplet volume.

Fig. 1 Time series (top to bottom) of snapshots for wettability waves travelling at different speeds: (left) wobbling motion at large speeds and (right) surfing motion of the droplet at small speeds. The greyscale shading indicates wettability and the red lines indicate the contact line of the droplet. The same two states are also observed for a travelling deformation wave.

First, in Section 2 we describe the theory of dynamic wetting for small droplets and our numerical approach based on the BEM. Next, in Section 3 we introduce the so-called nonuniform oscillator, which we use as an analytic model to interpret the observed droplet dynamics. We present and analyze the effect of travelling waves in wettability in Section 4 and travelling-wave deformations in Section 5. Finally, we comment on implications of our findings for the sorting of droplets by size in Section 6 and conclude with Section 7.

2 Boundary element method for dynamic wetting

We consider the motion of viscous droplets at small scales in the creeping flow regime, where viscous forces dominate, while inertia is negligible. Accordingly, fluid flow within droplets are described by the Stokes equations and the incompressibility condition,

η∇²v = ∇p and ∇·v = 0. (1)

Here, we use fluid velocity field v, pressure p, and shear viscosity η. In addition, appropriate boundary conditions have to be chosen, most importantly, describing forces acting at the fluid boundary. Notably, the flow fields determined by the Stokes equations adapt instantaneously to the boundary conditions because the equations do not contain any time derivatives.

In ref. 14 we formulated the boundary-element method for dynamic wetting. When extending this approach to droplets on deforming substrates, we encountered difficulties with the conservation of the droplet volume. This forced us to reformulate the boundary-element method in a more formalized way. As a result, for any point on the droplet surface we can write a force balance

where the friction force

K = −GṘ, (3)

is linear in the surface velocity field Ṙ and G is a linear operator that does not depend on v, since the Stokes equations are linear in the fluid velocity field v.^32,33 The force derives from the droplet free energy, as we explain in Section 2.2, and provides the Laplace pressure, Young's force at the contact line, and constraint forces. Finally, F governs the dynamics induced by deforming substrates.

After discretizing the droplet surface by a mesh, K, and F become vectors and G a friction matrix. In the following we systematically determine all these quantities, which will give us a systematic approach for calculating the motion of the droplet. We start with the mesh discretization in Section 2.1, then continue with the droplet free energy in Section 2.2 and discuss the relevant Lagrange parameters for the volume and rigidity constraint of droplet and deforming substrate, respectively, in Section 2.3. In Section 2.4 we analyze the boundary conditions at the droplet surface, which confirms the force balance of eqn (2), and we formulate the vectors K and F. In Section 2.5 we calculate the complete friction matrix by considering shear friction in the fluid using the boundary integral equation for the Stokes equation, as well as liquid–substrate and contact-line friction. In Section 2.6 we present the final dynamic equations, which we use in our simulations, and in Section 2.7 we collect the parameters used.

Before we proceed we add two comments. In our previous implementation of volume conservation we directly constrained the velocities of the droplet's mesh vertices following ref. 34. However, when the substrate deforms, this ad hoc constraint introduces spurious flows on the droplet's surface. Our current approach resolves this issue by properly introducing the volume constraint with the pressure as a Lagrange parameter in the droplet free energy as explained in Section 2.3. It constrains the restoring forces acting on the droplet's mesh vertices. This approach also allows us to directly incorporate an additional constraint that links the droplet to the height profile of the substrate. Thus, our reformulation of the boundary-element method puts an increased focus on forces rather than solely on velocities.

In parallel, we needed to modify our approach to contact-line friction. Whereas in ref. 14 we simply prescribed the velocity of the contact line using the Cox–Voinov law, now we first calculate the Young force acting on the contact line and then calculate its velocity with a Cox–Voinov type friction law (see Section 2.5.3).

2.1 Mesh discretization

As noted above, our aim is to reduce the droplet motion to the dynamics of its boundary consisting of the liquid–substrate and liquid–gas interfaces as well as the three-phase contact line, which requires special consideration. To treat the droplet dynamics numerically, we describe the entire droplet surface by a triangular mesh. Each triangle consists of three vertices and three edges, which it shares with its direct neighbour triangles. As long as the triangular mesh stays intact, the shape of the droplet is completely specified by the positions of the vertices r_i, which define the configuration vector R = (r₁, r₂, r₃, …).

In general, the part of the surface force field acting on the droplet at the liquid–solid and liquid–gas interactions derives from the functional derivative of a surface energy w.r.t. the configuration field R. However, since the droplet surface is discretized, the functional derivative becomes gradients w.r.t. the vertex positions. To continue, we construct an expression for the surface energy using surface integrals and discretize the integrals for our mesh in the following part.

2.2 Droplet free energy and discretization

The free energy of an incompressible liquid droplet at constant temperature can be written aswhere γ is the surface tension and the area A encompasses all the surfaces of the droplet. The bulk free energy is kept constant here. The surface integral can be split up into two parts: the liquid–gas (lg) interface and the solid–liquid interface (sl) with the respective surface tensions γ_lg and γ_sl. The solid–gas interface of the substrate not in contact with the droplet contributes to the free energy with surface tension γ_sg. Whenever part of this interface is covered by the droplet, the surface energy changes by γ_sl − γ_sg. Thus we obtainwhere A_lg and A_lg are the areas of the respective interfaces. Young's law for the equilibrium contact angle θ_eq is a force balance in the substrate plane that states that γ_sl − γ_sg = −γ_lg cos θ_eq.‡ Using it together with γ_lg = const., we obtainwhere θ_eq is determined by the local value of γ_sg − γ_sl. Because we want to study movable droplets on substrates with nonuniform wettability, the remaining integral cannot be further simplified here.

We now turn our attention toward the constraints on the shape of the droplet. They are two-fold: first, the volume of the droplet is conserved and, second, the normal forces on the liquid–solid interface sum up to zero since any fluid forces are compensated by the forces exerted by the rigid substrate that resists any elastic deformation. Using the method of Langrangian multipliers, each constraint gives rise to an additional term in the energy functional:

Here, the second term on the right-hand side constrains the droplet volume to the value V₀ and the pressure p₀ serves as the Langrange multiplier. The third term firmly connects the vertical position z of the liquid–solid interface to the height variations h(x) of the substrate written in Monge representation with x = (x, y). In the case of a flat substrate, h = 0, this term reduces towhich guarantees that the solid–liquid interface A_sl is defined by z = 0. Finally, note that the volume of the droplet can be expressed by the surface integralas one verifies immediately with Gauss's theorem.

Now, the functional derivative of is one part of the force balance of eqn (3). In our simulations we always calculate it numerically as detailed below. However, to illustrate its content, we perform the derivative and explicitely write it as

The last term on the right-hand side results from the rigidity constraint, which only affects the solid–liquid interface, and G(x) is the Jacobian determinant of the height profile h(x) of the curved substrate. The three contributions in eqn (10) are illustrated in Fig. 2, where the derivative of gives the acting forces due to surface tension and the Young force (see below). Their normal components acting on the substrate compensate the rigidity-constraint force and, as we demonstrate in Appendix D, the remaining terms of the functional derivative can be written asIn thermal equilibrium it is zero, while in non-equilibrium it drives the dynamics of the droplet. The first term contains the balance at the liquid–gas interface between static pressure p₀, which we introduced with the volume constraint, and the force due to surface tension. Here, κ is the local mean curvature of the liquid–gas interface, which is negative for the unit normal n pointing out of the droplet. The second term on the right-hand side of eqn (11) is the Young force acting locally on the contact line. It points in direction of unit vector t, which is perpendicular to the contact line and tangential to the substrate. In equilibrium the Young force vanishes (cos θ_eq = cos θ) and p₀ becomes the Laplace pressure p_L = −2κγ_lg = p₀.

Fig. 2 Schematic illustrating the three force contributions to in eqn (10) acting on the different interfaces of the droplet (black outline): static pressure p₀ (red), surface tension and Young force (blue), and rigidity-constraint force λGn at the solid–liquid interface (green). The rigid substrate is shaded in gray.

To discretize all the surface integrals of the free energy we use the triangular mesh of the droplet surface and sum up the contributions from each triangle so that

Here, we have to distinguish between functions f, the values of which are known for any position r and those functions, the values of which are only known at the vertices. In the first case, an example is cos θ_eq, we use Gaussian quadrature with 4 sample values of the prescribed profile to calculate the integral on each triangle. The sample values are taken at positions (1/3, 1/3, 1/3), (3/5, 1/5, 1/5), (1/5, 3/5, 1/5), and (1/5, 1/5, 3/5) in barycentric coordinates of the triangle and weighted by −27/48, 25/48, 25/48, and 25/48, respectively.³⁵ In the second case, an example is the calculation of the volume V, we linearly interpolate between the values at the vertices. If f has a constant value f_i on the area of the triangle, for example, in the case of γ_lg, we simply use A_if_i.

Note, after discretizing the free energy in terms of all positions of the vertices R, the functional derivative of becomes .

2.3 Method of lagrange multipliers

After having discretized all contributions to we obtain terms as a function of the configuration vector R, including all the constraint terms, which we write as λ_jg_j(R) = 0. For convenience, we label λ₀ = p₀ and assign the remaining multipliers to the rigidity constraint. Thus, eqn (7) becomes

We discretize the surface integral of the rigidity constraint in eqn (7),

where we subsume the area A_i surrounding each vertex into the λ_i, respectively, and we sum over all N vertices on the solid–liquid interface. The force due to this constraint is guaranteed to be normal to the substrate, since ∇_ig_i = n_i, as shown in Appendix B, and therefore equivalent to the desired rigidity constraint.§ This construction greatly simplifies the procedure below because ∇_Rg_i is easily accessible.

The purpose of the Lagrange multipliers λ_i is to constrain the force such that it does not have a component perpendicular to the manifold defined by the constraints g_i(R) = 0, otherwise the constraints would be violated.¶ Therefore, the condition

must hold. Using eqn (13) for in this condition withwe obtain a set of linear equations for the Lagrange multipliers λ_i (including λ₀),which we solve numerically. With this, is known up to a constant.||

2.4 Boundary condition and dynamic equation

So far we looked at the forces, which derive from the droplet free energy. They are part of the boundary condition of the Stokes equations, eqn (1), which balances all the forces acting at each point on the different interfaces of the droplet. There are three contributions to this force balance. First, we have the traction of the liquid on the interface, σn, where σ_ij = −pδ_ij + η(∂_iv_j + ∂_jv_i) is the stress tensor, n the unit normal vector, and δ_ij is the Kronecker symbol. For the following argumentation, we split the stress tensor into a part −p₀n, where p₀ is the spatially uniform static pressure due to the volume constraint already introduced in Section 2.2, and a spatially varying part [small sigma, Greek, tilde]n = σn + p₀n. Second, there is the traction of the substrate at the substrate–liquid interface. It consists of the constraint force f_rigid = λGn|_sl, already introduced in Section 2.2, with which the substrate resists deformations induced by the droplet, and f_deform, which acts on the liquid when the substrate deforms in time according to the prescribed height profile h. Finally and third, there is the traction due to the surface tension of both interfaces, which includes the Young force acting on the contact line and which we denoted in Section 2.2. Note that we neglect the vapour pressure of the gas phase surrounding the droplet, here and in the following.

In total, we have the force balance at each point of the droplet surface,

where f_rigid + f_deform is zero at the gas–liquid interface. Note that f_deform is due to the boundary conditions at the liquid–substrate interface, which include continuity of the normal velocity component and a slip boundary condition of the tangential component, as we will specify in Section 2.5.2. At the liquid–gas interface eqn (18) gives zero tangential stresses or traction forces, while the normal forces obey the usual Neumann boundary condition, n·[small sigma, Greek, tilde]n − p₀ = 2κγ_lg, where we have used eqn (11). Finally, using eqn (10) the first two terms on the r.h.s. of eqn (18) together with p₀n can be rewritten as and we obtainThe equilibrium case with zero dynamic stresses, [small sigma, Greek, tilde]n = 0, and f_deform = 0 gives and was already discussed in Section 2.2.

We now rephrase boundary condition (19) for the discretized droplet surface. At each vertex i the traction [small sigma, Greek, tilde]_in_i from the fluid is multiplied by area A_i and collected in the friction vector

K = (−[small sigma, Greek, tilde]n₁A₁, −[small sigma, Greek, tilde]n₂A₂, −[small sigma, Greek, tilde]n₃A₃, …), (20)

Similarly, we collect the substrate traction f_deform in a force vector

F = (f⁽¹⁾_deformA₁, f⁽²⁾_deformA₂, f⁽³⁾_deformA₃, …). (21)

Our method to calculate F is described in detail in Section 2.5. As usual, the gradient of describes a force which takes the place of the functional derivative . Thus, condition (19) is rewritten as which proves the force balance of eqn (2), with which we started. In the next section we will show that K = −GṘ and thereby calculate the friction matrix. This completes the reformulation of our boundary element approach.

2.5 Friction matrix

Now, we specify the friction matrix G such that all sources of friction related to the droplet are accounted for. Three types of friction are relevant in our setup. First, there is shear friction or viscous shear stresses in the droplet fluid, when neighboring fluid layers move relative to each other. Second, there is slip friction between the substrate and the directly adjacent fluid layer. And finally, there is friction of the moving contact line. The latter should naturally arise from the first and second mechanisms. However, so far this has proved difficult to show analytically⁵ and also to implement numerically due to the flow singularity that arises at the contact line.³⁶ Therefore, we use a hydrodynamically motivated friction law for the contact line. It is closely related to the Cox–Voinov law,³⁷ which has successfully been matched to experiments, e.g., in ref. 38.

2.5.1 Shear friction. We start with shear friction in the droplet fluid. The solution of the Stokes equations within a given domain is equivalent to a boundary-integral equation that relates flow velocity v(r) to surface forces σ(r)n and force/source dipoles at the domain surface,³¹where α is the inward solid angle, O the Oseen tensor, and T the associated stress tensor. For the surface velocities relevant here, α = 2π at the fluid interface and α = 2θ at the contact line, while inside the droplet α = 4π. We note that σn in eqn (22) can be replaced by [small sigma, Greek, tilde]n since a constant pressure does not contribute to the integral as shown in Appendix A. To discretize the boundary-integral equation, we proceed as in our previous work.¹⁴ We assign to each surface vertex parts of the surrounding triangles such that the resulting polygonal cells with area A_j do not overlap and cover the whole region of integration. The result is a set of linear equations for the vertex velocities v_j,Recalling the definition of K in eqn (20) and introducing block matrices, which we specify in Appendix C,we are able to rewrite eqn (23) as

CṘ = −XK − YṘ. (26)

or as

K_sf = −X⁻¹(C + Y)Ṙ, (27)

where X⁻¹(C + Y) is the shear contribution to the friction matrix.

2.5.2 Liquid–substrate friction. We now turn to friction between liquid and substrate due to a non-zero slip velocity. The conventional slip boundary condition, also called Navier condition, is l_sP_tσn = ηP_tv, where l_s is the slip length and P_t = 1 − n⊗n the projection operator on the tangential plane. Note that the surface normal n points out of the fluid according to ref. 39 and that P_tσn = P_t[small sigma, Greek, tilde]n. However, since we are interested in moving substrates, the fluid motion relative to the local substrate velocity v_subs and the deforming traction f_deform become relevant and we have l_sP_t([small sigma, Greek, tilde]n − f_deform) = ηP_t(v − v_subs). On the discretized solid/fluid interface of the substrate using eqn (20) and (21), the boundary condition becomes

P_‖(K_ls + F) = −QP_‖(Ṙ − V_subs). (28)

Here, is the projection operator from all the n droplet vertices onto the tangential space of the substrate with m vertices. For consistency, we defined V_subs = (v⁽¹⁾_subs, v⁽²⁾_subs,…) as the vector of the substrate velocities at the vertices, which also includes the vertices on the liquid–gas interface with v⁽ⁱ⁾_subs = 0. Finally,is the slip contribution to the friction matrix, where δ_ij is the Kronecker symbol and i, j refer to vertices on the substrate. Here we consider that, according to eqn (3), K must vanish if Ṙ = 0, and similarly, the parallel components of vector F in eqn (2) must vanish if there is no substrate deformation, i.e., if V_subs = 0. We therefore conclude that P_‖F = QP_‖V_subs which satisfies both conditions.

So far, we have not addressed the continuity condition for the normal velocity component, v·n = v_subs·n, at the liquid–substrate interface, which we already mentioned in Section 2.4. Similar to P_‖, it is convenient to introduce the projection operator onto the normal space of the substrate vertices . The boundary condition can then be restated for our discretized mesh as P_⊥Ṙ = P_⊥V_subs. It directly determines the normal components of the vertex velocities at the substrate, which are contained in Ṙ. However, the condition also means that the substrate pushes or pulls on the liquid as it deforms its height profile h and the liquid resists with friction. We have already introduced the traction forces due to the substrate deformations. Their normal components are contained in the vector F_⊥ = P^T_⊥P_⊥F, which has to be calculated numerically at the same time as those components of Ṙ which are still unknown.

For completeness we give the normal velocity component of the substrate velocity for a deforming substrate with height profile z = h(x, t). It amounts to v_subs(x)·n = ∂_th(x,t)e_z·n.

2.5.3 Contact-line friction. Finally, we address the friction of the contact line. In 1964 Moffatt solved the Stokes flow at the contact line of a moving wedge filled with a viscous liquid and found an analytic expression for the traction σn along the liquid–gas interface,⁴⁰Here, H is the distance of the interface from the substrate, θ the contact angle, and v_cl the tip velocity of the wedge along the substrate. Voinov relied on eqn (30) to derive the well-known Cov–Voinov law,³⁷ which we have previously used to prescribe contact line velocity in ref. 14. However, since in our reformulated approach the Young forces appear explicitely through the droplet free energy, we derive a direct relationship between Young forces and contact line velocity based on eqn (30). As in the Cox–Voinov law we do not take into acount contact angle hysteresis and pinning, which are not our focus here.

Again, if the substrate moves, we have v_cl = v − v_subs and its deforming traction f_deform, which modify eqn (30). The velocity lies in the tangential plane of the substrate and is normal to the contact line. To proceed, we integrate the traction over height H along the liquid/gas interface from a microscopic length l_s, i.e., the slip length, up to a mesoscopic length ζ

In the discretized representation of the droplet surface, we choose for the integral the surface force σ_in_i assigned to the vertex i on the contact line, which is contained in the friction vector K, i.e.,For the contact line, we therefore have the additional contribution

P_cl(K_cl + F) = −MP_cl(Ṙ − V_subs), (33)

where is the projection operator on the direction normal to the contact line and tangential to the substrate, where k is the number of vertices on the contact line. Obviously, these vertices contribute withto the friction matrix. Here we use the Kronecker symbol δ_ij, which means that M is diagonal. As in Section 2.5.2, we identify P_clF = MP_clV_subs as contribution to F in eqn (2). Note that the instantaneous local contact angle θ_i is determined by the shape of the droplet mesh, i.e., the value of R and so M_ij is not a function of the contact line velocity v_cl. However, θ_i also appears in the driving Young force [see second line of eqn (11)] and thereby determines v_cl.

We summarize here our result for the substrate traction force F. While the components tangential to the substrate are determined by the tangential velocity of the substrate, the normal component needs to be determined. Thus we have

F = (P^T_‖QP_‖ + P^T_clMP_cl)V_subs + F_⊥ (35)

where T indicates matrix transposition.

2.6 Final dynamic equations

Finally, we are able to formulate the full set of dynamic equations. Collecting all the contributions to the friction force K linear in Ṙ from eqn (27), (28) and (33), the total friction matrix G becomes

G = X⁻¹(C + Y) +P^T_‖QP_‖ + P^T_clMP_cl. (36)

Thus the force balance of eqn (2) is completely determined and readswhere he have used the specific expression for F from eqn (35). To solve this equation for the velocities Ṙ for a static substrate (V_subs = 0 and F = 0), one just needs to invert the friction matrix G. To test our reformulated implementation of the boundary-element method, in Appendix E we simulated the relaxation of a droplet towards its equilibrium shape. The relaxation of the dynamic contact angle is in very good agreement with experimental measurements.

For substrates deforming in time, the right-hand side of eqn (37) contains the unknown non-zero components of F_⊥, while on the left-hand side in Ṙ the normal velocity components at the solid–liquid interface are given by the normal substrate velocity components, i.e., P_⊥Ṙ = P_⊥V_subs. By introducing the projector that projects Ṙ onto the subspace of its unknown entries, we can reorder eqn (37) so that on the right-hand side only known quantities appear

where we used P_⊥Ṙ = P_⊥V_subs as well as P_uV_subs = P_uF_⊥ = 0. Eqn (36) further simplifies the right-hand side so thatNow, one first solves for P_uṘ by inverting −P_uGP^T_u and then uses P_uṘ in the first row to determine P_⊥F_⊥ if needed. This is a procedure, which we implemented in our computer code.

2.7 Nondimensionalization and parameters

We choose characteristic units, which derive from the motion of the droplet surfaces. First, as a unit of length l, we choose the initial radius R₀ of the liquid–solid interface of the droplet. The dominant time scale for the deformation of the droplet is determined by the motion of the contact line. A good estimate for the speed of the contact line was derived by Voinov in ref. 37 using eqn (30) and (18). He ultimately arrives the Cox–Voinov lawFrom the prefactor of the Cox–Voinov law and R₀, we obtain an inherent time scale

τ = 9γ_lg⁻¹ηR₀ ln(ζ/l_s). (41)

for the motion of the contact line. Furthermore, a characteristic force follows from combining the liquid–gas surface tension γ_lg and R₀ as f_c = γ_lgR₀. Now, writing eqn (1) in nondimensional form using these characteristic units, gives rise to a dimensionless viscosity [small eta, Greek, tilde] = [9 ln(ζ/l_s)]⁻¹. Thus, our results become independent from specific values of γ_lg and R₀. They only depend on the ratios [small eta, Greek, tilde], ζ/R₀ and l_s/R₀, for which we choose appropriate values.

As in ref. 14 the parameters of our system are matched to droplets from a 90%-glycerol/10%-water mixture for which de Ruijter et al.³⁸ measured in experiments ln(ζ/l_s) = 44, η = 209 mPa s, l_s = 1 nm, γ_lg = 65.3 mN m⁻¹, and mass density ρ = 1.24 g ml⁻¹. To determine a value for ζ needed in the friction matrix of eqn (34), we fitted the experimental data for the droplet relaxation in ref. 38 in our simulations and obtained ζ = 0.17R₀. All simulations are performed using a surface mesh of 1199 vertices which, in their initial configuration, have an average distance of 0.094R₀ from the nearest neighbor. In the following, we use R₀ = 100 μm.

From the characteristic parameters, we determine the Reynolds number as

which matches the assumption of negligible inertia inherent in eqn (1).

3 Driven droplets as nonuniform oscillators

Our droplets are driven by travelling waves either in wettability or in deformations of the substrate. For small wave speeds they are strongly connected to and surf on the wave, while for large speeds they can no longer follow and perform a wobbling motion. Before we describe the phenomenology of both substrate types, we introduce and study a simplified model for the motion of a droplet under a spatially periodic external stimulus. Neglecting inertia, we start with a simplified force balance for the center of mass y of the droplet, −Γẏ + F(y, t) = 0, where Γ is a friction constant and F is the effective force exerted by the substrate, on which the droplet is moving. The travelling wave pattern is modeled by the effective forcewith wavelength λ and wave speed v_wave. It models that there are specific positions in the travelling wave, y − v_wavet = nπ with n = 0, 1, …, where no force is excerted on the droplet. For our specific wave patterns, these are the minima and maxima. In the co-moving reference frame of the substrate pattern, y_com = y − v_wavet, the droplet's velocity ẏ_com = ẏ − v_wave generally oscillates around −v_wave, i.e.,with amplitude v_c = F₀/Γ. Depending on its position on the pattern, the droplet's speed is increased or reduced. In particular, for y_com = nλ/2 the force on the droplet is zero, which means ẏ = 0 or ẏ_com = −v_wave. Although ẏ_com oscillates around −v_wave, that is not its time average. Instead −v_wave is the spatial average of ẏ_com.

We can map eqn (44) onto the so-called ideal nonuniform oscillator²⁹ for phase variable φ with constants a and c by identifying

This results into the expression

[small phi (variant), Greek, dot above] = a − c sin φ (46)

also called the Adler equation.⁴¹ Unlike the uniform (harmonic) oscillator, the phase of which has a constant rate of change, φ determines its own rate of change. In particular, the oscillation stops for |c| ≥ |a| when the two terms on the r.h.s. balance each other.⁴²

Strogatz derives the period T of the nonuniform oscillator in ref. 29 and following his derivation we find with our parameters

To predict the long-time behavior of the droplet, for example, its position, the time averaged droplet speed v̄ is useful, which we derive now. Generally, in steady state the droplet position oscillates relative to the substrate pattern with a periodwhere the difference between v_wave and v̄ gives the droplet's speed relative to the substrate pattern with wavelength λ. Differently speaking, T is the period of the forcing experienced by the droplet through the travelling wave. Setting eqn (47) equal to eqn (48), we obtainwhich relates v̄ to v_wave with only a single free parameter, v_c. Consequently, in this model v_c contains all material properties of the driving mechanism and droplet-substrate interaction. The limiting case for large v_wave is v̄ ≈ v_c²/(2v_wave) and thereforewhich is useful to find v_c in the simulations but also in experiments.

However, v_c is not just a fitting parameter because it also marks a change in the dynamics of the droplet. For v_wave ≤ v_c the oscillation stops, as predicted by the non-uniform oscillator, which means ẏ_com = 0 or, equivalently, ẏ = v_wave = v̄. In our following analysis, we will refer to this stationary dynamics as ‘surfing’ and to the oscillatory dynamics as ‘wobbling’.

The main limitation of the nonuniform oscillator model is that the constant v_c is a priori unknown and has to be determined, e.g., from experimental or (in our case) simulation data. From its definition v_c = F₀/Γ we see that it is specific for the periodic forcing, eqn (43), which is applied to the droplet. Consider the case of deforming substrates with a travelling-wave height profile. On a deforming substrate, the amplitude of the driving force F₀ comes from the substrate curvature. However, the curvature is determined by the amplitude of the travelling wave and its wavelength λ. Therefore, v_c is not a independently tunable parameter but also a function of λ. Similarly, in the case of travelling waves in wettability, λ determines the local gradient in wettability and thereby F₀. Because v_c generally depends on λ this means that the relations expressed in eqn (47), (49) and (50) are only useful for predictions as long as λ (and therefore v_c) is kept constant.

4 Travelling waves in wettability

We first study how travelling waves in wettability can move droplets on a completely flat substrate, similar to our earlier study of moving steps in wettability in ref. 14. However, there are clear differences, which we explain in the following.

4.1 Substrate dynamic

Here, we impose a travelling wave in wettability on a flat substrate such that the prescribed equilibrium contact angle varies according towith wavelength λ = 2R₀, v_wave the speed of the wave in y-direction, the Heaviside or Θ-function, the amplitude Δθ = 30°, and θ₀ = 90°.**

In experiment such a travelling wave can be generated, for example, by first exposing an azobenzene-coated substrate (cmp. ref. 11) to structured blue light with a plane-wave intensity profile, then UV light of uniform intensity, then again structured blue light with a slightly shifted intensity profile, and so on. By alternating between these UV and blue light patterns, one produces a wettability profile which slowly changes in time.

The droplet moves toward regions of high wettability, which means the valleys of the travelling wave in θ_eq. Because these valleys are steadily moving forward in y-direction, the droplet is biased to also move in that direction. However, as noted in Section 3, the droplet can perform different kinds of motion under the influence of this driving mechanism, namely surfing and wobbling motion (see Fig. 1), which we will analyse in detail now.

4.2 Surfing and wobbling

We first review the basic phenomenology. A surfing droplet moves at the same speed as the travelling wave. Consequently, its shape remains constant and a stationary flow field exists inside the droplet (see Movie M1 in the ESI†). This stationary state is identical to the surfing dynamics we observed in ref. 14.

A wobbling droplet moves at a slower speed than the travelling wave pattern. It oscillates in speed and shape as it repeatedly passes over the peaks and valleys of the wave (see Movie M2 in the ESI†). The wobbling dynamics was not present in our previous study of moving steps in wettability¹⁴ because it requires a periodic wettability pattern.

According to the nonuniform oscillator model, we outlined in Section 3, wobbling occurs above a critical speed of the travelling wave pattern. As the speed of the travelling wave pattern decreases, the period at which the droplet shape oscillates grows and it diverges exactly at the transition. So, the resulting surfing state can be considered as a wobbling which is frozen in time. This bifurcation is analyzed in detail in Section 4.4.

4.3 Quantitative analysis

One characteristic quantity of the surfing and wobbling dynamics is the droplet speed and we ask how it relates to the speed of the wave pattern. The droplet speed needs to be defined carefully, especially for wobbling, because it varies periodically over time. It is appropriate to average the droplet speed over one period of oscillation T. The period of oscillation refers to the shortest duration between two identical shapes of the droplet, albeit not identical positions in the lab reference frame. The average droplet speed is then formally given bywhere the starting time t is set after the decay of any transients and the number of periods n is chosen according to the available simulation data. Note that, for the case of a surfing droplet, we have y(t) = v_wavet + y₀ with position y₀ at t = 0. Here, the time interval T is arbitrary and the definition above reduces to v̄ = v_wave.

In Fig. 3(a) we plot the speed ratio v̄/v_waveversus v_wave as determined from our BEM simulations. The speed ratio is 1 in the surfing state and it decays as v_wave⁻² in the wobbling state (dashed line), as predicted from the nonuniform oscillator model and eqn (49) in the limit v_wave ≫ v_c. This means v̄ is inversely proportional to v_wave for wobbling and equal to v_wave for surfing. The time-averaged speed of the droplet decays to zero which also matches our observation from an earlier article that droplets are less susceptible to rapid changes in wettability.³⁰ For a full quantitative comparison with our model, we need to determine the velocity v_c. Following eqn (50), we plot versus v_wave in Fig. 3(b) and read off v_c = 2.65R₀τ⁻¹ at large wave speeds (dash-dotted red line). The dashed line in Fig. 3(a) is the prediction from eqn (49) using the determined value for v_c.

Fig. 3 Diagrams of (a) the mean droplet speed v̄, (b) the helper quantity and (c) the wobbling oscillation period T plotted versus v_wave for a travelling wave in wettability with wavelength λ = 2R₀. Blue dots indicate our results from BEM simulation, the dashed black lines indicate predictions from the nonuniform oscillator model and the dash-dotted red line indicates the critical wave speed v_c, which is the sole fitting parameter of the model.

The transition from surfing to wobbling appears continuous, meaning there is no jump in droplet speed, as predicted by the nonuniform oscillator model. However, the transient dynamics decays slowly close to the transition so that we lack data close to v_wave = v_c. If the transition is indeed continuous, the droplet moves fastest exactly at the critical wave speed v_c. The period T of wobbling oscillations, displayed in Fig. 3(c), grows inversely proportional to v_wave − v_c as v_wave approaches v_c from above. These observations for the wobbling state again quantitatively match the prediction of the nonuniform oscillator model.

4.4 Comparison with travelling steps in wettability

In ref. 14 we observed a droplet surfing along a substrate on a wettability step that moves with velocity v_step and a sigmoidal profile σ(y). Similar to eqn (44) we can describe the droplet velocity in the co-moving frame of the step withwhere for the sigmoidal function we chose σ(y) = (1 + e^−y)⁻¹, and δ is the step width. Most importantly, for v_step < v_c two surfing states exist. They sit outside the center of the step and are stable and unstable fixed points, which eventually meet in a saddle-node bifurcation at v_step = v_c. Furthermore, there is a sessile steady state where the droplet sits outside of the step. It occurs for any value of v_step in the limit |y_com/δ| ≫ 0 so that the contribution from σ′ vanishes. The saddle node bifurcation with the stable and unstable fixed points as well as the stable fixed point of the sessile droplet are schematically illustrated in Fig. 4(a) using an appropriate shape variable.

Fig. 4 Schematic bifurcation diagrams for an appropriate droplet-shape variable plotted versus pattern velocity. (a) The droplet dynamics in the step wettability pattern (cf. ref. 14) corresponds to a saddle-node bifurcation at the critical velocity v_c. (b) The travelling wave pattern generates an additional limit cycle that collides with the saddle and node, giving rise to a saddle-node-infinite-period (SNIPER) bifurcation.

Comparing the moving wettability step treated in ref. 14 to the travelling wettability wave considered here and modeled by eqn (44), we note that a sessile state is impossible due to the periodicity of the wave pattern. Instead, the wobbling dynamic takes the place of the sessile state for v_wave > v_c. In Fig. 4(b) the resulting bifurcation diagram is sketched. The saddle-node bifurcation of the moving wettability step is still present but for travelling waves it collides with the limit cycle associated with wobbling. This collision comprises a global bifurcation which is called Saddle-Node-Infinite-PERiod (SNIPER) bifurcation.²⁹

A SNIPER bifurcation was previously observed for droplets exposed to an oscillating external force on a static heterogeneous substrate⁴¹ and for a droplet subject to a constant external force on a rotating cylinder^43,44 with the force acting perpendicular to the cylinder's axis of rotation. In both cases, the bifurcation marks the transition from a pinned to a sliding droplet.

5 Travelling-wave deformations

We now turn to the second driving mechanism introduced in Section 2. The substrate deforms such that its height profile follows a travelling wave. The wettability is kept uniform and we solely rely on the travelling height profile to move the droplet. In experiment such a travelling wave can be generated, for example, by light-responsive liquid-crystal elastomers exposed to structured light patterns.²⁰ The influence of gravity in such a deformation wave is negligible because the capillary length is much larger than the droplet radius. This raises the question of the driving mechanism that moves the droplet forward.

For an equilibrium contact angle of θ_eq = 90°, which according to Young's law occurs at surface tensions γ_sl = γ_sg, only the liquid–gas interface is relevant and the droplet deforms to minimize the area of this interface. Hence, on a curved surface the droplet moves away from peaks and toward valleys because there the liquid–gas interface occupies less area compared to a droplet with the same volume sitting on a flat part or a peak of the substrate. For θ_eq < 90° the droplet is even more attracted to valleys in the substrate, as it is energetically advantageous to cover more of the solid with liquid, while for θ_eq > 90° the situation is less clear. In ref. 19 Lv et al. demonstrate how a droplet moves on the inside and outside of a conical substrate without any additional external stimulus. On the outside of the cone the liquid–substrate curvature is negative and on the inside of the cone it is positive. Notably, for all cases studied with θ_eq ranging from 30° to 120°, the droplets move in the direction of increasing curvature. The same behavior was recently confirmed for droplets on cylinders of varying size.⁴⁵ Indeed, when we place a droplet onto a substrate peak with θ_eq = 120° in our BEM simulations and then perturb the droplet, the position is unstable. Lv et al. term this mechanism curvi-propulsion and in this section we investigate if it generates the same wobbling and surfing dynamics we found for the wettability waves in Section 4 and if there are any differences between the two cases.

5.1 Surfing and wobbling

Instead of a wettability profile, we now impose a travelling height profile of the substrate with amplitude h₀,where ξ determines the initial phase shift or position of the wave relative to the droplet at y = 0. Especially, in the surfing dynamic we use a value of ξ = 0.3R₀ to place the droplet close to its surfing position on the wave profile, while ξ is irrelevant for the wobbling dynamics. We introduce the factorandThe purpose of the factor q(t) is to allow a continuous evolution of the surface deformation from a flat to the sinusoidal shape to which the droplet can adjust gradually starting at a specific position relative to the wave. Then, the deformation wave starts to travel at t = t₀, where we choose t₀ in the interval 0.1τ ≤ t₀ ≤ 0.2τ. To be specific, we set λ = 4R₀ for the wavelength, h₀ = 0.1R₀ for the amplitude, and an equilibrium contact angle of θ_eq = 90°.†† In Section 6, we will explore travelling deformation waves for a range of wavelengths.

With time the droplet again settles into either surfing or wobbling motion depending on the wave speed v_wave. Because it is attracted to the valleys in the substrate, it tries to follow them as the travelling deformation wave progresses. If the travelling wave moves sufficiently slowly, the droplet can surf on a position behind the valley (see Movie M3 in the ESI†). If the travelling wave moves faster than a critical speed, the droplet is always overtaken by the peak behind it, relaxes into the closest valley, and falls back again giving rise to the wobbling motion (see Movie M4 in the ESI†).

5.2 Quantitative analysis

In analogy to the travelling wettability pattern, we analyze the time-averaged speed of the droplet as defined in eqn (52). In Fig. 5(a) we again plot v̄/v_waveversus v_wave. In contrast to the wettability wave, the speed of the wobbling droplet now grows linearly with v_wave at large v_wave. Due to this clear difference, we realize that our model needs to be adjusted.

Fig. 5 Diagrams of (a) the mean droplet speed v̄, (b) the helper quantity v*, and (c) the wobbling oscillation period T as functions of v_wave for a travelling deformation wave with wavelength λ = 4R₀. Blue dots indicate our results from BEM simulation, the dashed black line indicates the revised nonuniform oscillator model, eqn (59), the dash-dotted red line indicates the critical wave speed v_c, and the dotted red line indicates the parameter ε. The latter two are the fitting parameters of the model.

To start, we re-examine the nonuniform oscillator model of eqn (44), where a spatially periodic force drives the droplet with a constant amplitude v_c. For the deforming substrate we introduced a traction force F_⊥ in eqn (37), which imposes the deformation of the substrate onto the droplet shape. However, this force depends on v_wave: it is zero for v_wave ≤ v_c, where the droplet's shape is constant and beyond v_c it grows with v_wave as the droplet is deformed more rapidly. As before, the force is also spatially periodic since it originates from the wave pattern of the substrate's height profile. Therefore, to include the dependence of the amplitude of the driving force on v_wave, we extend eqn (44) and write

The dimensionless coefficient ε ≪ 1 quantifies the relative strength of the wave-speed depending part of the amplitude. In analogy to eqn (47), we find for the time period,which together with eqn (48) gives the mean droplet speed,Here, we have used in eqn (58). The small parameter ε means that v̄ increases proportional to v_wave for large v_wave. Note that ε does not affect the scaling of the wobbling period T with the inverse wave speed well above v_c.

Now, our procedure to determine ε and v_c, the central parameters of the model, is as follows. First, from eqn (59) we determine the limiting value of v̄/v_wave = ε for large v_wave, which we can directly read off from Fig. 5(a) as ε = 3.565 × 10⁻³. Then, we subtract εv_wave from both sides of eqn (59),

and expand the r.h.s. for large v_wave to obtainSolving for v_c yieldswhich for ε = 0 indeed simplifies to eqn (50). Accordingly, in Fig. 5(b) we plot v* as defined in eqn (62). Because ε appears as coefficient of v_wave², v* is sensitive to small errors in ε at large v_wave. However, around v_wave/R₀τ⁻¹ = 2 it takes on a constant value, from which we determine v_c = 0.185R₀τ⁻¹. The dashed black line labelled model in Fig. 5(a) shows that with these parameters our revised model quantitatively matches the BEM data (blue dots). Similarly, it also reproduces the wobbling period displayed in Fig. 5(c). Because around v_wave = v_c, the original and the revised model behave the same, the latter also predicts a SNIPER bifurcation at v_wave = v_c for droplets driven by travelling substrate deformations.

Our revised model fits the BEM data well and thus provides an interpretation of the asymptotic behavior of v̄ for large v_wave. To capture the asymptote, we introduced an amplitude of the driving force, which depends on v_wave. It accounts for the fact that a travelling substrate deformation always needs to displace the droplet, which requires a larger driving force with increasing v_wave. This is the clear difference to travelling wettability patterns treated in Section 4, where the pattern can just pass beneath the droplet without the need to lift it up.

6 A note on particle sorting

In Sections 4 and 5 we described that droplets move in response to spatial gradients in wettability or curvature of the substrate, respectively. Both, gradient and curvature, change in magnitude when we adjust the wavelength λ of the travelling wave pattern. Therefore, we can expect that the speed of the droplet v̄ also depends on λ. However, in the surfing state the droplet speed is identical to the wave speed v_wave, regardless of λ. Therefore, we only need to consider a wobbling droplet.

In Fig. 6 we display for both driving mechanisms the mean droplet velocity v̄ as a function of R₀/λ for v_wave = 50R₀τ⁻¹, which places the droplet far into the regime of wobbling dynamics. For travelling deformation waves (green crosses in Fig. 6) no clear trend is discernible because mean droplet speed v̄ has a minimum near R₀ = 0.17λ. In contrast, for travelling waves in wettability (blue dots in Fig. 6) we observe that v̄ grows like (R₀/λ)². Further decreasing λ, one would expect v̄ to decrease again. However, for larger ratios R₀/λ the creeping flow approximation does not hold anymore for deforming substrates and one enters the regime of unsteady Stokes flow due to the increased oscillation frequency of the height profile.

Fig. 6 Diagram of mean droplet speed v̄ as a function of R₀/λ, i.e., the ratio between droplet radius and wavelength, for travelling waves with speed v_wave = 50R₀τ⁻¹. Blue dots indicate results from BEM simulation with waves in wettability (right axis), and green crosses indicate results from BEM simulation with deformation waves (left axis) with h₀ = 0.1R₀. Connecting dotted and dashed lines are added to guide the eye. The solid black line indicates a quadratic scaling.

Now consider a poly-disperse collection of droplets with R₀ in the range of 0.1λ to 0.5λ on a single substrate with a travelling wave in wettability with fixed λ. According to Fig. 6 there is a unique value of v̄ for each droplet size. This means, after we let the droplets run for some time (and with sufficient distance between them to prevent mergers), the droplets will separate according to their size. Because small droplets lag behind larger droplets, this naturally leads to a sorting mechanism for such a poly-disperse collection of droplets.

7 Conclusions

We have studied liquid-droplet motion on substrates using two different mechanisms: travelling waves in wettability and travelling-wave deformations of the substrate. To that end, we first implemented the boundary element method (BEM) for the fluid velocity field inside the droplet starting from the free energy of the droplet, the balance of forces at the droplet interface, and the relevant friction contributions.

To interpret our numerical findings from applying the BEM, we first introduced an analytic model for a particle driven by a travelling-wave force, which maps onto the so-called nonuniform oscillator. For increasing wave velocity, the analytic model exhibits a bifurcation between two states, which we call surfing and wobbling and which directly connect to the phenomenology observed in our BEM simulations.

For travelling waves in wettability the results from BEM simulations for wobbling period and mean droplet speed are consistent with the analytic model and indeed confirm the bifurcation in the full dynamical system. This further implies that the motion of the droplet's center of mass can be characterized solely by the critical wave speed for the bifurcation, which is determined by the material properties of the droplet and substrate. For a specific set of parameters we determined the critical wave speed by varying the wave speed in the BEM simulations.

We also investigated travelling-wave deformations of the substrate, which give rise to what is called curvi-propulsion in ref. 19. Again, we could identify a transition from surfing to wobbling with increasing wave speed but, interestingly, the simulation results reveal that the droplet speed does not approach zero for large wave speeds. By including a speed-dependent amplitude of the driving force, the modified analytic model could account for this new behavior, as we demonstrated for a specific set of parameters in the BEM simulations.

Finally, we also performed simulations where we varied the wavelength λ of the travelling wave for both driving mechanisms in the wobbling regime. For waves in wettability the droplet speed decays as λ⁻², which implies that the critical wave speed is proportional to droplet radius over wavelength. This suggests a potential sorting mechanism for a collection of droplets by droplet size. For travelling-wave deformations of the substrate, we could not observe a clear trend in droplet speed w.r.t. wavelength.

Our findings compare and contrast two different mechanisms for steering droplets toward a preferred direction, one based on switchable wettability and one based on externally-induced deformations of the substrate. They demonstrate that the relevant quantities such as droplet speed and wobbling period can be reproduced and interpreted with an analytic model using only one or two fitting parameters, respectively. Moving droplets is relevant for lab-on-a-chip applications. Our work suggest two methods for realizing this. Another attractive possibility to explore are substrates with switchable softness⁴⁶ and the occurence of durotaxis since our reformulated approach is already well-adapted to deforming substrates. In particular, the elasticity of a soft substrate can be represented by an additional contribution to the free energy.

Conflicts of interest

There are no conflicts to declare.

Appendices

A Invariance of the boundary-integral equation

In eqn (22) one can replace σn by [small sigma, Greek, tilde]n = σn + p₀n, where the constant pressure p₀ does not contribute to the surface integral. To show this, we write for an arbitrary force vector f₀,The closed surface integral of the Stokeslet flow field, u_Stokeslet(r) = O(r − r₀)f₀, vanishes, since u_Stokeslet(r) is incompressible by construction and after applying Gauss's theorem.

B Normal vector in Monge representation

We parameterize the substrate according to Monge representation by r̃ = (x,y,h(x,y)). The normal n then points perpendicular to both tangential vectors of the surface which are given by the derivatives of r̃ w.r.t. x and y, respectively, so thatNote that n, as written here, is generally not of unit length. However, as noted in Section 2.3, the normalization constant is subsumed into the Lagrange multiplier λ.

C Block matrices

The matrices in eqn (23) are defined asandwhere the surface integral is performed over the polygonal cell C_j around vertex j and r_i is the position of vertex i. Furthermore, in eqn (23), if vertex i is on the contact line, c_i = θ_i/2π with contact angle θ_i and otherwise c_i = 1/2.

D Functional derivative and rigidity constraint

The functional derivative of w.r.t. shape deformation iswhich contributes to the r.h.s. of eqn (10). The first and second terms are due to surface tension within the liquid–gas and liquid–solid interfaces of the droplet, respectively, and κ is the local curvature of the respective interface.⁴⁷ The third and fourth terms on the right-hand side of eqn (67) comprise the Young force acting locally on the contact line, where δ_cl is the Dirac δ-function centered at the location of the contact line. The Young force has two components: one along unit vector t tangential to the plane of the substrate and one along unit vector n_sl normal to the substrate.

In eqn (10) the last term is a rigidity force with which the substrate acts to keep its shape against the interior pressure p₀, against the curvature force at the solid–liquid interface in the second term on the r.h.s. of eqn (67), and the normal component of the Young force in the fourth term. Because the rigidity constraint with the Lagrange multiplier λ needs to balance these forces, the rigidity-constraint force f_rigid = G(x)λ(x)n|_sl introduced in Section 2.4 reads

f_rigid = (2κγ_lg cos θ_eq + γ_lg sin(θ)δ_cl|_cl − p₀)n|_sl. (68)

Together with eqn (10) and (67) it results in eqn (11).

E Validation

To test our method we compare its result to an experiment reported in ref. 38 and to a fit with our previous simulation approach in ref. 14. In the experiment, a glycerol droplet (with the material properties written in Section 2.7) is placed on a substrate with equilibrium contact angle 64.3° starting from a contact angle above this value. During the relaxation of the droplet, its dynamic contact angle is recorded over time. To transfer the parameters of the experiment to BEM, we simulate a droplet which is initially in equilibrium with a contact angle close to the first recorded angle in the experiment on a substrate with equilibrium contact angle 64.3°. As in the experiment, we measure the dynamic contact angle θ in our simulation by fitting a spherical cap shape to the discretized liquid–gas interface at each time step and then calculate θ from the cap's tangent at the contact line.

To compare our simulation to the experiment, we need to convert our timesteps to SI units by explicitely calculating the value of τ for the experimental data as specified in eqn (41). In the experiment, the volume of the droplet is reported to be at least V = 5 mm³ from which we calculate . With this value for R₀ and the material properties written in Section 2.7, we calculate τ = 1.7 s. We display both the experimental data and the result of our simulation in Fig. 7 where we can see good agreement, especially in the initial relaxation. The long-time relaxation of the experimental data appears to be a bit slower than our simulation data which may be due to defects on the substrate which we do not consider in our method. The test demonstrates that our method is able to reproduce the relaxation of a real droplet as it is observed in experiment. Notably, because we have carefully nondimensionalized our method, we do not need to alter any simulation parameters to simulate smaller or larger droplets; any change in droplet volume merely affects the conversion of our timescale from τ to seconds.

Fig. 7 Diagram of the dynamic contact angle θ over time t. Red dots indicate the experimental data visually extracted from Fig. 2(a) in ref. 38 and the solid blue line shows our simulation data which has been converted to ESI,† units with τ = 1.7 s, the time scale implied by the parameters of the experiment. The dash-dotted green line shows data from our previous simulation approach in ref. 14 and the dashed black line indicates the equilibrium contact angle θ_eq = 64.3°.

Acknowledgements

We thank Uwe Thiele and Akash Choudhary for fruitful discussions and acknowledge financial support from German Research Foundation (DFG) as part of priority program 2171 (project number 505839720).

Notes and references

1. M. Alwazzan, K. Egab, B. Peng, J. Khan and C. Li, Int. J. Heat Mass Transf., 2017, 112, 991.
2. M. Edalatpour, L. Liu, A. Jacobi, K. Eid and A. Sommers, Appl. Energy, 2018, 222, 967.
3. S. N. Varanakkottu, M. Anyfantakis, M. Morel, S. Rudiuk and D. Baigl, Nano Lett., 2016, 16, 644.
4. L. Qi, Y. Niu, C. Ruck and Y. Zhao, Lab Chip, 2019, 19, 223–232.
5. D. Bonn, J. Eggers, J. Indekeu, J. Meunier and E. Rolley, Rev. Mod. Phys., 2009, 81, 739.
6. P. Teng, D. Tian, H. Fu and S. Wang, Mater. Chem. Front., 2020, 4, 140–154.
7. U. Thiele, K. Neuffer, M. Bestehorn, Y. Pomeau and M. G. Velarde, Colloids Surf., A, 2002, 206, 87.
8. D. Baigl, Lab Chip, 2012, 12, 3637.
9. A. Venancio-Marques and D. Baigl, Langmuir, 2014, 30, 4207–4212.
10. M. K. Chaudhury and G. M. Whitesides, Science, 1992, 256, 1539.
11. K. Ichimura, S. Oh and M. Nakagawa, Science, 2000, 288, 1624.
12. J. Vialetto, M. Hayakawa, N. Kavokine, M. Takinoue, S. N. Varanakkottu, S. Rudiuk, M. Anyfantakis, M. Morel and D. Baigl, Angew. Chem., Int. Ed., 2017, 56, 16565–16570.
13. J. ter Schiphorst, J. Saez, D. Diamond, F. Benito-Lopez and A. P. H. J. Schenning, Lab Chip, 2018, 18, 699–709.
14. J. Grawitter and H. Stark, Soft Matter, 2021, 17, 2454.
15. H. B. Eral, D. J. C. M. t Mannetje and J. M. Oh, Colloid Polym. Sci., 2013, 291, 247.
16. H. S. Lim, J. T. Han, D. Kwak, M. Jin and K. Cho, J. Am. Chem. Soc., 2006, 128, 14458.
17. H.-J. Butt, J. Liu, K. Koynov, B. Straub, C. Hinduja, I. Roismann, R. Berger, X. Li, D. Vollmer, W. Steffen and M. Kappl, Curr. Opin. Colloid Interface Sci., 2022, 59, 101574.
18. S. Karpitschka, A. Pandey, L. A. Lubbers, J. H. Weijs, L. Botto, S. Das, B. Andreotti and J. H. Snoeijer, Proc. Natl. Acad. Sci. U. S. A., 2016, 113, 7403–7407.
19. C. Lv, C. Chen, Y.-C. Chuang, F.-G. Tseng, Y. Yin, F. Grey and Q. Zheng, Phys. Rev. Lett., 2014, 113, 026101.
20. S. Palagi, A. G. Mark, S. Y. Reigh, K. Melde, T. Qiu, H. Zeng, C. Parmeggiani, D. Martella, A. Sanchez-Castillo, N. Kapernaum, F. Giesselmann, D. S. Wiersma, E. Lauga and P. Fischer, Nat. Mater., 2016, 15, 647–653.
21. I. Rehor, C. Maslen, P. G. Moerman, B. G. van Ravensteijn, R. van Alst, J. Groenewold, H. B. Eral and W. K. Kegel, Soft Robot., 2021, 8, 10–18.
22. M. P. Gelfand and R. Lipowsky, Phys. Rev. B: Condens. Matter Mater. Phys., 1987, 36, 8725–8735.
23. C. Duprat and H. A. Stone, Fluid-Structure Interactions in Low-Reynolds-Number Flows, The Royal Society of Chemistry, 2016, pp. 193–246.
24. J. Bico, É. Reyssat and B. Roman, Annu. Rev. Fluid Mech., 2018, 50, 629–659.
25. S. Aland and D. Mokbel, Int. J. Numer. Methods Eng., 2021, 122, 903–918.
26. R. Shuttleworth, Proc. Phys. Soc. Sect. A, 1950, 63, 444–457.
27. B. Andreotti and J. H. Snoeijer, Annu. Rev. Fluid Mech., 2020, 52, 285–308.
28. T. M. Squires and S. R. Quake, Rev. Mod. Phys., 2005, 77, 977.
29. S. Strogatz, Nonlinear dynamics and Chaos, Westview Press, 1994.
30. J. Grawitter and H. Stark, Soft Matter, 2021, 17, 9469.
31. C. Pozrikidis, Boundary integral and singularity methods for linearized viscous flow, Cambridge University Press, 1992.
32. M. Doi, J. Phys.: Condens. Matter, 2011, 23, 284118.
33. A. Zafferi, D. Peschka and M. Thomas, Z. Angew. Math. Mech., 2023, 103, e202100254.
34. E. Alinovi and A. Bottaro, J. Comput. Phys., 2018, 356, 261–281.
35. in The Boundary Element Method for Engineers and Scientists, ed. J. T. Katsikadelis, Academic Press, Oxford, 2nd edn, 2016, pp. 403–421.
36. C. Huh and L. E. Scriven, J. Colloid Interface Sci., 1971, 35, 85.
37. O. V. Voinov, Fluid Dyn., 1976, 11, 714.
38. M. J. de Ruijter, J. De Coninck, T. D. Blake, A. Clarke and A. Rankin, Langmuir, 1997, 13, 7293.
39. S. J. Bolaños and B. Vernescu, Phys. Fluids, 2017, 29, 057103.
40. H. K. Moffatt, J. Fluid Mech., 1964, 18, 1–18.
41. U. Thiele and E. Knobloch, Phys. Rev. Lett., 2006, 97, 204501.
42. R. Adler, Proc. IRE, 1946, 34, 351.
43. U. Thiele, J. Fluid Mech., 2011, 671, 121.
44. T.-S. Lin, S. Rogers, D. Tseluiko and U. Thiele, Phys. Fluids, 2016, 28, 082102.
45. J. Liu, Z. Feng, W. Ouyang, L. Shui and Z. Liu, ACS Omega, 2022, 7, 20975–20982.
46. N. Nekoonam, G. Vera, A. Goralczyk, F. Mayoussi, P. Zhu, D. Böcherer, A. Shakeel and D. Helmer, ACS Appl. Mater. Interfaces, 2023, 15, 27234–27242.
47. K. Deckelnick, G. Dziuk and C. M. Elliott, Acta Numerica, 2005, 14, 139–232.

---

Footnotes

† Electronic supplementary information (ESI) available: Consists of four animations, M1–M4, corresponding to specific BEM simuluations: M1 and M2 correspond to simulations with v_wave = 1R₀τ⁻¹ (surfing) and v_wave = 10R₀τ⁻¹ (wobbling) in Fig. 3, respectively. M3 corresponds to a simulation as described in Section 5 with ξ = 1R₀, t₀ = 0.2τ and v_wave = 0.18R₀τ⁻¹ (surfing). M4 corresponds to the simulation with v_wave = 70R₀τ⁻¹ (wobbling) in Fig. 5. See DOI: https://doi.org/10.1039/d4sm00213j

‡ In fact, Young's law follows from minimizing under a volume constraint.

§ In the case of a flat substrate placed at a height h = 0, we have n_i = e_z and so the constraint reduces to λ_ig_i = λ_iz_i, which is aquivalent to eqn (8).

¶ This is commonly referred to as the principle of virtual work.

|| Note that we can reuse the gradients ∇_Rg_i from above in calculating in eqn (16). Note also that in our previous study ref. 14 when we studied a droplet on a rigid flat substrate there was no need to calculate the constraint forces due to the rigid substrate because in that case we calculated the force on each vertex individually and imposed the volume constraint only on the resulting velocity field. That approach proved inapplicable for our current study.

** Note that in ref. 14 we studied a moving wettability step and explored the influence of different contact angles θ₀ and step widths Δθ on the surfing speed of a droplet. One important finding was that θ₀ = 90° does not result in a qualitatively different result than other values of θ₀.

†† Note that θ_eq = 90° is not a special case here because curvi-propulsion acts in the same direction regardless of wettability, i.e., regardless whether a substrate is hydrophobic or hydrophilic.¹⁹

---