Stick-Slip Contact Line Motion on Kelvin-Voigt Model Substrates

The capillary traction of a liquid contact line causes highly localized deformations in soft solids, tremendously slowing down wetting and dewetting dynamics by viscoelastic braking. Enforcing nonetheless large velocities leads to the so-called stick-slip instability, during which the contact line periodically depins from its own wetting ridge. The mechanism of this periodic motion and, especially, the role of the dynamics in the fluid have remained elusive, partly because a theoretical description of the unsteady soft wetting problem is not available so far. Here we present the first numerical simulations of the full unsteady soft wetting problem, with a full coupling between the liquid and the solid dynamics. We observe three regimes of soft wetting dynamics: steady viscoelastic braking at slow speeds, stick-slip motion at intermediate speeds, followed by a region of viscoelastic braking where stick-slip is suppressed by liquid damping, which ultimately gives way to classical wetting dynamics, dominated by liquid dissipation.

Introduction.-The capillary interaction of liquids with soft solids is a ubiquitous situation in natural or technological settings [1][2][3][4], ranging from droplets that interact with epithelia, for instance in the human eye [5], or epithelial cells governed by capillary physics [6], to ink-jet printing on flexible materials [7]. The capillary tractions, exerted by the liquid onto their soft support, cause strong deformations if the substrate is soft or the considered length scale is sufficiently small [8,9]. The typical scale below which capillarity deforms solids is given by the elastocapillary length, = γ/G 0 , the ratio of surface tension γ and static shear modulus G 0 . At three-phase contact lines, the length scale of the exerted traction lies in the molecular domain, deforming the solid into a sharp-tipped wetting ridge [10]. As a liquid spreads over a soft surface, the traction moves relative to the material points of the substrate. The necessary rearrangement of the solid deformation leads to strong viscoelastic dissipation which counteracts the motion, a phenomenon called viscoelastic braking [11][12][13][14][15][16][17][18][19]. At small speeds, the motion remains steady [11,12], whereas at large speeds, unsteady motion, frequently termed stick-slip, has been observed [13,[20][21][22][23][24]. In this mode, the contact line velocity and the apparent contact angle undergo strong, periodic oscillations. On paraffin gels, Kajiya et al. observed stick-slip motion only in an intermediate velocity range, returning to continuous motion if the speed was increased even further [21].
The origin of this stick-slip motion remains debated in literature. It is clear that the pinning and depinning is not associated with permanent surface features, but rather with the dynamics of the wetting ridge itself: the solid deformation cannot follow the fast contact line motion of the depinned (slip) phase of a stick-slip cycle [23,24]. Unclear, however, are the conditions upon which a contact line may escape from its ridge, thus eliminating the viscoelastic braking force. The depinning of a contact line from a sharp-tipped feature on a surface is governed by the Gibbs inequality [23,25]. Van Gorcum et al. [23] postulated a dynamical solid surface tension, which would alter the local force balance and thus allow the contact line to slide down the slope of the ridge. Still, the physicochemical origin of such dynamic solid surface tension remains elusive. Roche et al. [15] postulated the existence of a point force due to bulk viscoelasticity, but the shearthinning nature of typical soft polymeric materials would prevent such a singularity at the strain rates encountered in soft wetting [24,26]. Unclear as well is the role of the fluid phase during the cyclic motion, mainly because a comprehensive multi-physics model for the unsteady soft wetting problem is not available to date.
Here we present the first fully unsteady numerical simulations of dynamical soft wetting, fully accounting for liquid and solid mechanics, and for the capillarity of the interfaces, by which we reveal the life cycle of stick-slip motion. We derive phase diagrams of steady and unsteady contact line motion by tuning parameters over large ranges, recov-  ering stick-slip behavior at intermediate speeds. At small and large speeds, we observe steady motion, in quantitative agreement with an analytical model.
Setup.- Figure 1 (a) shows the geometric setup of the numerical simulations. A hollow cylinder (undeformed inner radius R), made of a soft viscoelastic material (gray, thickness h s R), with a fixed (rigid) outer surface, is partially filled with a liquid (blue) and an ambient fluid phase (transparent). To keep physics conceivable, we use a minimal model and apply the Stokes limit and identical viscosities for fluid and ambient, and assume constant and equal surface tensions γ = γ s on all interfaces (liquid-ambient, solid-liquid, and solid-ambient). The inner surface of the soft viscoelastic cylinder wall is deformed into a wetting ridge due to the capillary action of the liquid meniscus (cf. panel (b)). We use an incompressible Kelvin-Voigt constitutive model, characterized by a frequency dependent complex modulus G * = G 0 +i η s ω, with static shear modulus G 0 and effective substrate viscosity η s , obtaining an elastocapillary length = γ/G 0 h s , a characteristic time scale τ = η s /G 0 , and a characteristic velocity v = /τ = γ/η s .
Our numerical model is formulated with cylindrical symmetry, implementing the two-phase fluid by a phasefield approach. Thus the liquid-ambient interface has a finite thickness h s , and the capillary traction of the meniscus onto the solid is distributed over this characteristic width. The solid is modelled by a finite element approach, with a sharp interface toward the fluid. The grid size is about 5% of the elastocapillary length at the liquid-ambient interface, and typically about 20% outside of the interface region. All phases are fully coupled to each other by kinematic and stress boundary conditions. The material parameters are listed in table 1, details on the numerical model are given in [27] and the supplementary material [28]. The liquid meniscus is forced to move by imposing the fluxes Φ on either end of the cylinder, but can freely change its shape (curvature) in response to the the fluid flow. Thus the instantaneous contact line velocity v is not imposed, but rather it's long-term mean v = Φ/(π R 2 ). All simulations are started at t = 0 with a flat substrate, a flat meniscus, and a constant imposed flux at the boundaries, and run until a steady state or limit cycle has been reached.
We compare our simulation results for the solid deformation to the analytical plane-strain model from [13], imposing a constant contact line velocity v = v and replacing the δ-shaped traction by which can be derived for the phase field model in equilibrium (see supplementary material for details [28]). Since h s R, the substrate deformation is well approximated by plane-strain conditions. Importantly, since , our analytical and numerical results do not significantly depend on the actual value of . Figure 1 (b) shows the quasistationary substrate deformation for several imposed velocities, comparing simulation results (markers) to the analytical model (lines), in excellent agreement. Note, that this comparison is only possible as steady ridge shapes are observed for the chosen velocities. In an intermediate velocity range we find unsteady cyclic shape dynamics, as detailed below, which cannot be grasped by the analytical model.
Modes of contact line motion.-The dynamics of the contact line motion are characterized by the timedependent rotation φ = θ−θ eq of the liquid interface at the triple line (cf. Fig. 1 (a)). Figure 2 shows φ as a function of contact line position for the three characteristic regimes that we find in our simulations. At small speed (v v , blue), after some initial transient the contact line moves steadily, with a constant dynamic contact angle. Here, the relation between v and φ is permanently dominated by viscoelastic braking [11][12][13]. Once the forcing velocity exceed a critical value, the motion becomes unsteady, finding a limit cycle after an initial transient (red): the liquid interface rotation φ shows large oscillations, of peak-to-peak amplitude ∆φ on the order of the mean rotation φ, with a non-trivial waveform, as the contact line advances. This behavior is not captured by the simple analytical model. For larger speeds (yellow), we observe again a constant φ after an initial transient. We note here that the motion in this regime is very sensitive to discretization artifacts and requires rather fine grid resolutions to give consistent results. Movies illustrating contact line motion and substrate dynamics for the three modes in Figure 2 can be found in the supplementary data. Figure 3 shows a phase portrait of the contact line motion i.e., in terms of the physically relevant variables φ and v: φ = θ − θ eq is a measure of the imbalance in Young's equation, and thus a measure of the total dissipative force (liquid and solid). Multiplied with the instantaneous velocity, one obtains the total dissipated power per unit length of contact line, since our equations of motion are overdamped. For slow forcing speeds (blue), we observe a continuous, steady contact line motion, up to the scale of grid artifacts (mind the logarithmic scales). For intermediate forcing speeds (red), we observe a limitcycle: As the liquid rotation exceeds a well-defined maximum ( 1 ○ in Fig. 3), the contact line accelerates. In this phase, it surfs down it's own wetting ridge, releasing energy stored in the meniscus curvature, rate-limited partly by liquid dissipation ( 2 ○ in Fig. 3). It thus decelerates, and a new wetting ridge starts to grow, opposing the contact line motion further ( 3 ○ in Fig. 3) until the next cycle starts. For larger forcing speeds, the region covered by the limit cycle decreases until it virtually vanishes (yellow), up to grid artifacts. This is caused by the growing importance of liquid dissipation, which effectively limits, and finally prevents, the large-speed excursions during the slip phases.
Regimes of contact line motion.-We characterize the contact line motion by φ max ( 1 ○ in Fig. 3), φ min ( 3 ○ in Fig. 3), and φ, the minimum, maximum, and mean values p-3 of φ in the stationary/limit cycle regime. Figure 4 shows these values as a function of the imposed (long-term mean) v. For small speeds, v = v, and the simulated φ (symbols) coincides with the result of the analytical model (black line), indicating the stability of steady contact line motion. The onset of stick-slip motion (red region) aligns with the maximum of φ observed in the analytical model where v is imposed instead of v. This was stipulated in [13] since the rotation φ is a measure for the dissipative (viscoelastic braking) force: A dissipative force that decreases with speed causes acceleration, and thus an unstable motion. The maximum braking force is observed at v = v = γ/(G 0 τ ), the elastocapillary velocity: The finite width of the traction distribution regularizes the dissipation singularity at the scale h s . Thus the contact line motion excites a dominant frequency ∼ v/ in the solid, corresponding to a dynamical elastocapillary length v ∼ γ s /G * (v/ ) = γ s /(η s v). Resonance is expected at ∼ v i.e., v ∼ γ s /η s , independent of the choice of , which is confirmed in our analytical model [28]. φ max remains approximately constant upon entering the stickslip regime, indicating a well-defined upper limit of the viscoelastic braking force also in unsteady situations (cf. location 1 ○ in Fig. 3). However, this force periodically drops to much smaller values, as indicated by the much smaller values of φ min . In these surfing phases ( 2 ○ in Fig. 3), liquid dissipation and the finite capillary energy stored in the curved meniscus are the rate-limiting factors.
As the imposed v is increased further, the amplitude of the oscillation ∆φ shrinks, reaching virtually zero (indicated by the fading red region). In this regime, the reduced viscoelastic braking force (φ max ) limits the build-up of capillary energy in the meniscus, while liquid dissipation prevents its fast release. Thus the oscillatory motion is effectively damped out by liquid dissipation, while the overall motion is still governed by viscoelastic braking: φ ∼ v −1 closely follows the result from the analytical model. The increased mean liquid rotation for the largest velocity ∼ 0.28v , relative to the prediction of the analytical model, is caused by liquid dissipation. This can be rationalized by a comparison with the Cox-Voinov law for moving contact lines on rigid surfaces, which, given the capillary number ∼ 10 −2 , predicts rotations on this order of magnitude [29][30][31]. In this hydrodynamic regime, one returns to the classical wetting physics on rigid surfaces. In Figure 5, we summarize the dynamical wetting behavior in terms of these three modes, as a function of the imposed mean speed and the solid parameters. Steady small-speed, stick-slip, and steady high-speed modes are indicated by blue, red, and yellow discs, respectively. On panel (a), we vary on the vertical axis the solid and liquid surface tensions and thus the elastocapillary number α s = /h s , while keeping the Neumann angles of static wetting constant. The onset of stick-slip is located near v = v , given by the maximum of φ vs. v. This maximum is independent of α s , up to a small correction due to the finite thickness of the fluid interface, as can be shown by the analytical model (see Fig. 1 of the supplementary material [28]). In physical units, however, the onset of stick-slip is inversely proportional to γ s since v ∼ γ s . The transition to the fast continuous mode is, in scaled units, nearly independent of the surface tension. Consequentially, the stick-slip mode disappears at very low γ.
Similarly, the solid viscosity η s (panel (b)) has no measurable impact on the critical v/v for the transition to stick slip, but the physical critical velocity is proportional to η s since v ∼ η −1 s . Thus, for small solid viscosities, the damping effect of liquid dissipation and thus the transition back to steady motion becomes noticeable already at smaller v/v , such that the stick-slip region ultimately disappears at very low η s . In any case, at very large speeds, liquid dissipation will take over, leading to wetting dynamics equivalent to those on rigid surfaces.
Discussion.-In this Letter, we provide a comprehensive numerical analysis of dynamical soft wetting, including the physics of all relevant elements, the liquid, the solid, and the interfaces. For each element, we used the minimal required level of complexity, to keep physics intact and conceivable: The Stokes limit for the fluid, a Kelvin-Voigt constitutive relation for the soft solid, regularized at a constant scale , and constant and equal solid surface tensions on all three interfaces. This simple model already requires a complex, strongly coupled multi-physics modelling approach, and exhibits rich behaviors.
Our numerical experiments cover a wide range of system parameters, and we reveal three regimes in which the dominant physical mechanisms differ: (i) a slow regime, in which the contact line motion is entirely dominated by the the dissipation in the solid. This regime is observed as long as the viscoelastic braking force increases with speed. (ii) an intermediate regime, in which the dominant ratelimiting mechanism periodically switches from solid to liquid dissipation. This regime starts where the viscoelastic braking force exhibits a maximum with respect to the contact line velocity. This maximum is caused by a resonance effect, due to the regularization of a singular dissipation at some finite (constant) length scale. Other mechanisms, like dynamic solid surface tensions (surface constitutive relations [23,[32][33][34][35][36]) or a constitutive relation that exhibits resonance (e.g., standard-linear solid [13]) would lead to the same phenomenology. (iii) a large-v-regime with continuous motion, yet governed by viscoelastic braking, in which liquid dissipation prevents strong oscillations of the meniscus. Since the viscoelastic braking force, in contrast to liquid dissipation, does not increase with velocity, we ultimately find again liquid dissipation dominating the contact line motion, and one recovers the wetting physics of rigid surfaces.
With this first comprehensive overview of soft wetting physics scenarios, we provide a strong basis for interpreting different phenomenology observed in experiments, ranging from paraffins [21] over microelectronic sealants [8,23] to biology [5,6,37], and motivate experiments in the so-far little explored large-v regimes. * * * SK and SA acknowledge funding by the German Research Foundation (DFG project no. KA4747/2-1 to SK and AL1705/5-1 to SA). Simulations were performed at the Center for Information Services and High Performance Computing (ZIH) at TU Dresden.