# Breakup Dynamics of Capillary Bridges on Hydrophobic Stripes

Maximilian Hartmann  Nano- and Microfluidics Group, TU Darmstadt, Alarich-
Weiss-Straße 10, 64287 Darmstadt, Germany Mathis Fricke Corresponding auhor.
Mathematical Modeling and Analysis Group, TU Darmstadt, Alarich-Weiss-Straße
10, 64287 Darmstadt, Germany  Lukas Weimar  Nano- and Microfluidics Group, TU
Darmstadt, Alarich-Weiss-Straße 10, 64287 Darmstadt, Germany Dirk Gründing
Mathematical Modeling and Analysis Group, TU Darmstadt, Alarich-Weiss-Straße
10, 64287 Darmstadt, Germany  Tomislav Marić  Mathematical Modeling and
Analysis Group, TU Darmstadt, Alarich-Weiss-Straße 10, 64287 Darmstadt,
Germany  Dieter Bothe  Mathematical Modeling and Analysis Group, TU Darmstadt,
Alarich-Weiss-Straße 10, 64287 Darmstadt, Germany  Steffen Hardt  Nano- and
Microfluidics Group, TU Darmstadt, Alarich-Weiss-Straße 10, 64287 Darmstadt,
Germany

###### Abstract

The breakup dynamics of a capillary bridge on a hydrophobic stripe between two
hydrophilic stripes is studied experimentally and numerically using direct
numerical simulations. The capillary bridge is formed from an evaporating
water droplet wetting three neighboring stripes of a chemically patterned
surface. By considering the breakup process in a phase space representation,
the breakup dynamics can be evaluated without the uncertainty in determining
the precise breakup time. The simulations are based on the Volume-of-Fluid
(VOF) method implemented in Free Surface 3D (FS3D). In order to construct
physically realistic initial data for the VOF simulation, Surface Evolver is
employed to calculate an initial configuration consistent with experiments.
Numerical instabilities at the contact line are reduced by a novel
discretization of the Navier-slip boundary condition on staggered grids. The
breakup of the capillary bridge cannot be characterized by a unique scaling
relationship. Instead, at different stages of the breakup process different
scaling exponents apply, and the structure of the bridge undergoes a
qualitative change. In the final stage of breakup, the capillary bridge forms
a liquid thread that breaks up consistently with the Rayleigh-Plateau
instability.

This preprint was accepted for publication in the _International Journal of
Multiphase Flow_. When citing this work, please refer to the journal article:
DOI: 10.1016/j.ijmultiphaseflow.2021.103582.

Keywords: Structured surface, Capillary bridge, Volume-of-Fluid, Breakup
dynamics,
Rayleigh-Plateau instability ††footnotetext: E-Mail addresses:
<EMAIL_ADDRESS>(M. Hartmann<EMAIL_ADDRESS>(M.
Fricke<EMAIL_ADDRESS>(L. Weimar<EMAIL_ADDRESS>darmstadt.de (D. Gründing<EMAIL_ADDRESS>(T. Marić),
<EMAIL_ADDRESS>(D. Bothe<EMAIL_ADDRESS>(S. Hardt)

##  1 Introduction

Wetting of patterned surfaces is omni-present in nature. The Lotus effect [3]
or the fog harvesting of the Stenocara desert beetle in the Namib Desert [45]
are only two examples. Also when it comes to engineering applications like
ink-jet printing [58] or water harvesting [61], the understanding of wetting
behavior on (chemically) patterned surfaces is of crucial importance.

In the present study, the focus lies on chemically patterned striped surfaces
that are wetted by water droplets with radii of the order of magnitude of the
stripe width. For such kind of systems, static wetting behavior has been
studied by using energy minimization techniques [12, 42, 32], lattice
Boltzmann simulations [34, 35] and experiments [4]. Also inertia-driven
spreading [35] and the splitting [38, 51, 62] of impinging droplets on such
kind of surfaces have been investigated in some detail. When it comes to
evaporation, simulations using the phase field [59] method have been performed
more recently. Hartmann and Hardt [30] showed that an evaporating droplet
wetting two hydrophilic stripes, with a hydrophobic one in between them, is
stable as long as the pressure inside the liquid forming the capillary bridge
above the hydrophobic stripe can be balanced in the liquid above the two
hydrophilic stripes. This is the case until a certain width of the capillary
bridge is reached. When the width of the bridge decreases further due to
evaporation, it breaks up. While in the latter publication the focus is on the
statics of an evaporating droplet, the subject of the present article is the
dynamics of the breakup process itself.

Qualitatively, this process shows similarities to the collapse of a soap film
between two circular rings, as it was investigated by Chen and Steen [11]
using numerical calculations. The authors found two breakup regimes, a
$d\sim\tau^{2/3}$-, and a $d\sim\tau^{2/5}$-regime (with $d$ the minimum
bridge width and $\tau$ the time before breakup), which are both dominated by
the balance between capillary and inertial forces. The change in regime is due
to a geometric transition during breakup. The $d\sim\tau^{2/3}$-scaling is a
classical result for the evolution of the minimum width $d$ of a “free liquid
bridge” (i.e. without contact to a substrate) obtained from dimensional
analysis [36]

$\displaystyle d(\tau)=C\left(\frac{\sigma\tau^{2}}{\rho}\right)^{1/3},$ (1)

where $\sigma$ is the surface tension, $\rho$ is the fluid density and
$\tau=t_{0}-t$ is the time $t$ before the bridge breaks up at $t_{0}$. The
prefactor $C$ has been believed to be a universal constant with a value close
to $1.4$ for quite some time (see, e.g., Eggers and Fontelos [19]). Assuming
the prefactor to be universal with $C=0.9\pm 0.01$, Hauner et al. [31]
concluded that water exhibits a dynamic surface tension, characterizing
freshly created surface sections, of $\sigma\approx 90\,\text{mN}/\text{m}$,
which is significantly larger than the equilibrium value of
$72\,\text{mN}/\text{m}$. However, more recent studies show that $C$ is indeed
not universal. Instead, the experimentally observable value of the prefactor
may depend on the fluid and system parameters [14].

According to Li & Sprittles [40], the appropriate dimensionless number
characterizing the breakup of a free capillary bridge of a liquid is the
Ohnesorge number

$\displaystyle\operatorname{Oh}=\frac{\eta}{\sqrt{\rho\sigma R}}\,,$ (2)

which is the square root of the ratio of the viscous length scale
$\eta^{2}/(\rho\sigma)$ (with $\eta$ being the dynamic viscosity) and a
characteristic length scale $R$. We choose the width of the hydrophobic stripe
$250\operatorname{{\mu}m}\leq w_{\text{phob}}\leq 750\operatorname{{\mu}m}$ as
the characteristic length scale i.e., $R=w_{\mathrm{phob}}$. With this
definition, the Ohnesorge number for water in air for the relevant hydrophobic
stripe widths lies within the range

$\displaystyle 4.3\cdot 10^{-3}\leq\operatorname{Oh}\leq 7.5\cdot 10^{-3}.$
(3)

Clearly, the limiting cases $\operatorname{Oh}\rightarrow 0$ and
$\operatorname{Oh}\rightarrow\infty$ characterize the inviscid and viscous
breakup regime, respectively. The inviscid breakup regime is already mentioned
above for the soap film. Within the viscous regime, viscous forces play the
dominant role and $d$ scales linearly with time $\tau$, $d\sim K\tau$, with
$K$ being a constant. In the same publication [40], a phase diagram is
presented, in which, besides the above-mentioned two regimes, also a third
regime, the viscous inertial regime for intermediate values of
$\operatorname{Oh}$, can be identified in the parameter space that is spanned
by $d$ and $\operatorname{Oh}$. Here, as well as in the viscous regime, $d\sim
K\tau$ is valid but with a different constant $K$. Li and Sprittles [40]
performed numerical simulations for a free capillary bridge with a similar
value of $\operatorname{Oh}$ compared to equation (3) and observed dynamic
transitions into the inviscid regime and from the inviscid into the viscous
regime as the minimum diameter of the bridge decreases approaching breakup.

Note that Chen and Steen [11] as well as Li and Sprittles [40] deal with a
free capillary bridge without substrate contact. This is not the case in the
present investigation, where the capillary bridge is in contact with the
substrate. For liquids in contact with a substrate, Bostwick and Steen [5]
provided a review, where the focus lies on the stability of constrained
capillary surfaces in general. The same authors theoretically investigated the
instability of static rivulets [6] and considered varicose (symmetric) and
sinuous (anti-symmetric) modes for pinned and free contact lines. For
symmetric modes, they confirmed the results by Davis [13]: For all contact
angles, breakup can happen in the case of a free, i.e. unpinned, contact line,
while for pinned contact lines, contact angles must be greater than
$90^{\circ}$. Other studies exist for unpinned [16, 9] and pinned [52] contact
lines, as well. In the case of anti-symmetric modes, a static rivulet is
always stable in the case of a pinned contact line and unstable only for
contact angles larger than $90^{\circ}$.

The above-mentioned articles deal with the stability of rivulets. To the best
knowledge of the authors, no investigations of the breakup and dynamics of
capillary bridges in contact with surfaces have been published so far. The aim
of this study is to examine the corresponding breakup dynamics and to compare
it to a “free” liquid bridge. A collection of the reported data is published
in an open research data repository [29].

##  2 Experimental Methods

Substrates are prepared by creating a stripe pattern of a positive photo
resist (AZ 9260, MicroChemicals GmbH, Germany) on a borofloat33 glass wafer
(Siegert Wafer GmbH, Germany) using standard photolithography steps. This
allows to produce covered and uncovered glass areas. After oxygen plasma
treatment, these wafers are silanized in a low-pressure chemical vapor
deposition process, for which 1H,1H,2H,2H-Perfluorodecyltrichlorosilane
(PFDTS, CAS: 78560-44-8, abcr GmbH, Germany) is used. In the following step,
the photo resist is rinsed off using acetone and isopropanol. Subsequently,
the substrate is dried in a nitrogen stream. This leads to a pattern of
hydrophilic and hydrophobic stripes. The substrates are then stored until
experiments are performed. Different hydrophilic stripe widths
$w_{\mathrm{phil}}$ and different ratios $\alpha$ of hydrophobic and
hydrophilic stripe widths are used, where

$\alpha=w_{\mathrm{phob}}/w_{\mathrm{phil}}.$

Experiments are performed by placing de-ionized water droplets (Milli-Q
device; specific resistance 18.2 M$\Omega\cdot$cm at 25 ∘C) with a pipette
onto the stripe pattern. The volume is chosen in a way that two hydrophilic
stripes with one hydrophobic stripe in between them are wetted. The water then
evaporates in the lab environment (temperature $\sim$25 ∘C) until a critical
width of the bridge on the hydrophobic stripe is reached. The critical width
marks the transition to a configuration beyond which no stationary wetting
state of the droplet exists (see Hartmann and Hardt [30] for a discussion of
the stability of the liquid bridge). We did not control the lab humidity since
the final breakup process occurs at a small time scale which we expect to be
of the order of the capillary time scale. It can be shown that the capillary
time scale is about 5 orders of magnitude smaller than the evaporation time
scale that therefore can be neglected in the present case (see Appendix C for
more details).

The final time span ($\sim$ 0.13 s) before breakup, i.e. the time span in
which the capillary bridge on the hydrophobic stripe decays, is recorded with
a high-speed camera (Photron FASTCAM SA-1.1) in top view mode. In order to
obtain the values necessary for fixing the initial condition for the numerical
simulations, e.g. the wetted length and the contact angle on the hydrophilic
stripe in the moment of breakup, a second high-speed camera (Photron FASTCAM
SA-X2) is used, which records the droplet from one side. The cameras are
synchronized via dedicated software (Photron FASTCAM Viewer) and triggered
externally. A frame rate of 75,000 frames per second (fps) is used. Each
camera is connected to a macro objective (Navitar-12X). Illumination is
performed via backlight in side view and co-axially in top view using two
cold-light sources (VOLPI intra LED 5). A silicon wafer that acts as a mirror
is placed below the glass substrate in order to achieve good illumination
conditions and record the breakup process at a sufficient frame rate and
magnification. We will also report results on the very last stages of bridge
breakup. In these cases, the capillary bridge is recorded in bottom view mode
with the Photron camera attached to a microscopy body. A more detailed
description of this experimental setup, used materials and substrate
preparation steps can be found in a previously published article [30]. Contact
angles on the silanized glass wafers are measured using the evaporation
method, in which the receding contact angle can be measured in the constant-
contact angle mode during evaporation [44].

##  3 Numerical Methods

###  3.1 Stationary States from Energy Minimization

To be able to simulate the process of droplet breakup, a precise geometric
description of the liquid distribution at the beginning of the simulation is
indispensable. Hence, the interface geometry has to be extracted from the
experimental observations and transformed into a phase fraction field that can
be processed by Free Surface 3D (FS3D). For this purpose, a triangulated
surface is produced using Surface Evolver [8], a tool for calculating minimal
surfaces by minimizing the free energy functional that, in the present case,
can be written as

$\mathcal{F}=\sigma_{\mathrm{lg}}\cdot|A_{\mathrm{lg}}|+\sum_{i\in\\{\mathrm{phil,phob}\\}}(\sigma_{\mathrm{sl,}i}-\sigma_{\mathrm{sg,}i})\cdot|A_{\mathrm{sl,}i}|$
(4)

subject to a prescribed volume $V$, contact angle $\Theta$ and wetted length
on the hydrophilic stripe $l_{\mathrm{phil}}$. In the present case, the
interfacial tension between liquid and gas $\sigma_{\mathrm{lg}}$ is constant,
while the specific interfacial areas $|A_{\mathrm{j,i}}|$ can change to
minimize the free energy $\mathcal{F}$. The subscripts sl, sg, and lg denote
solid-liquid, solid-gas and gas-liquid, respectively. Since the droplet wets
alternating hydrophilic and hydrophobic (abbreviated with phil and phob)
stripes, both must be taken into account for calculating $\mathcal{F}$. Note
that $\sigma_{\mathrm{sl,}i}-\sigma_{\mathrm{sg,}i}$ can be substituted by
Young’s equation, which is useful in the present case, because the contact
angles on $i\in\\{\mathrm{phil,phob}\\}$ can be determined experimentally.

FIG. 1: Initial geometry for the numerical simulations obtained from Surface
Evolver (blue). Due to symmetry, only 1/4th of the droplet is calculated. Red
areas denote hydrophobic, green areas hydrophilic stripes. The wetted length
$l_{\mathrm{phil}}$, the contact angle on the hydrophilic stripe
$\Theta_{\mathrm{phil}}$, the hydrophilic stripe width $w_{\mathrm{phil}}$ and
$\alpha$, the ratio between the hydrophobic and the hydrophilic stripe width,
are held constant during the calculation. Used values can be found in Table 1.

As input data, the wetted length of the hydrophilic stripe $l_{\mathrm{phil}}$
as well as the contact angle $\Theta_{\mathrm{phil}}$ on it, as observed in
the experiments, are used. These values are summarized in Table 1.

$\alpha\leavevmode\nobreak\ (-)$ | $\Theta_{\mathrm{phil}}\leavevmode\nobreak\ (^{\circ})$ | $l_{\mathrm{phil}}\leavevmode\nobreak\ (\operatorname{{\mu}m})$
---|---|---
0.5 | 21 | 1966
1 | 28 | 2638
1.5 | 31 | 2936
Table 1: Contact angle and wetted length on the hydrophilic stripe as observed
in the experiments and used in the Surface Evolver calculations.

Furthermore, the contact angle on the hydrophobic stripe
$\Theta_{\mathrm{phob}}$ is assumed to be the receding contact angle on a
silanized glass wafer. The corresponding experimental value is 102∘. All of
these three parameters are kept constant during the whole minimization
procedure. Note that at the boundary between the hydrophilic and the
hydrophobic stripe, the contact angle may lie between $\Theta_{\mathrm{phil}}$
and $\Theta_{\mathrm{phob}}$. Due to symmetry considerations and to reduce
calculation time, only a quarter of a droplet is calculated, as it is
indicated with blue color in Fig. 1. For that purpose, two mirror planes are
considered. Both are oriented perpendicular to each other and to the surface
itself. While one mirror plane is located parallel to the stripes in the
middle of the hydrophobic one, the other plane is oriented perpendicular to
the stripes, cutting the hydrophobic bridge in the middle along its axis. At
these mirror planes, the contact angle is set to 90∘. This means that the
liquid–gas interface meets the symmetry plane at an angle of 90∘, which is a
symmetry boundary condition.

The goal of the Surface Evolver calculations is to find the equilibrium
solution of the capillary surface that obeys the above described boundary
conditions (i.e. wets both the hydrophilic and the hydrophobic stripes) and
belongs to the smallest volume of liquid, which corresponds to the minimum
stable liquid bridge configuration. For that purpose an algorithm is
implemented as described in the following:

  1. 1.

Start calculation with a certain volume $V$. If this is the initial step of
the droplet evolution, the volume is $V=V_{0}$.

  2. 2.

Evolve surface and calculate

$c_{v,n}=\frac{s_{n}}{\overline{\mathcal{F}}}=\frac{\frac{1}{\sqrt{n(n-1)}}\sqrt{\sum_{i}^{i-n}(\overline{\mathcal{F}}-\mathcal{F}_{i})^{2}}}{\frac{1}{n}\sum_{i}^{i-n}\mathcal{F}_{i}},$
(5)

the coefficient of variation of the energy of the $n$ last iterations.
$c_{v,n}$ is the ratio between the standard deviation $s_{n}$ and the mean
value of the free energy $\overline{\mathcal{F}}$ at iteration $i$. Note that
in the present work $n=7$.

  3. 3.

Calculate bridge width $d$.

  4. 4.

Check for convergence. Convergence is reached when $c_{v,7}<10^{-7}$. If the
calculation has converged, save $V$ as $V_{\mathrm{stable}}$, and decrease
volume by $\Delta V=V_{0}/2$ if convergence has been reached for the first
time, or by $\Delta V=\Delta V/2$ if convergence has been reached at least
once, so that $V=V-\Delta V$. Then check if the criterion of step 6 is
fulfilled. If not, continue with step 5.

  5. 5.

Check for breakup. The capillary bridge is considered to be broken up if
$d<0.1\cdot w_{\mathrm{phil}}$. This value is well below the critical bridge
width according to Hartmann and Hardt [30]. If this criterion is fulfilled,
increase volume by $\Delta V=\Delta V/2$ so that $V=V+\Delta V$, and check if
the criterion of step 6 is fulfilled. If not, repeat steps 1 - 5.

  6. 6.

Check if $\Delta V$ is smaller than 0.5 % of the last stable volume
$V_{\mathrm{stable}}$. If yes, stop the iteration.

During the evolution of the surface, the mesh is successively refined until
the maximum length of each cell edge is smaller than
$1/100\cdot(w_{\mathrm{phil}}+w_{\mathrm{phob}})$. Convergence can only be
achieved if this is true. The overall procedure of volume variation, as
described above, follows the principle of nested intervals and is similar to a
binary search algorithm. From the stopping criterion $\Delta V=0.005\cdot
V_{\mathrm{stable}}$ it is expected that the last found critical volume is
close enough to the real value. Note that it is impossible to calculate
arbitrarily close to the real critical volume. The result of the algorithm is
the triangulated surface mesh with the smallest stable volume. This geometry
is exported as a triangulated surface mesh in STL (Standard Triangulation
Language) format and subsequently transformed into a volume-fraction field for
FS3D (see Section 3.3). Note that Surface Evolver is a tool that is usually
used to calculate stable minimal surfaces. In the present case, it is used to
calculate an equilibrium shape close to the onset of instability. Moreover,
from experiments it is known that before the capillary bridge becomes unstable
(during evaporation), the contact angle on the hydrophobic stripe changes
within the range that is given by the advancing and receding contact angle. In
Surface Evolver, the contact angles are set to fixed values, i.e. contact
angle hysteresis is not accounted for, although it is observed in the
experiments. Therefore, the initial configuration differs from the one
observed in the experiments (see Fig. 19 in Appendix B).

###  3.2 Continuum Mechanical Model

FIG. 2: Illustration of employed mathematical notations.

The continuum mechanical model is based on the incompressible two-phase Navier
Stokes equations in the sharp interface formulation. The conservation
equations for momentum and mass in the bulk phases read

$\displaystyle\rho\frac{D\vec{v}}{Dt}-\nabla\cdot S+\nabla p=\rho\vec{g},\
\nabla\cdot\vec{v}=0\ \quad$ $\displaystyle\text{in}\
\Omega\setminus\Sigma(t),$ (6)

where $p$ is the pressure and $S=\eta(\nabla\vec{v}+{\nabla\vec{v}}^{\sf T})$
is the viscous stress tensor. Note that $\vec{g}$ is the gravitational
acceleration and the Lagrangian time-derivative is defined as

$\frac{D}{Dt}=\partial_{t}+\vec{v}\cdot\nabla\,.$

To simplify the model, it is assumed that no mass is transferred across the
gas-liquid interface. This can be expected to be a good approximation on the
timescale of the breakup process, which is much smaller than the evaporation
timescale (see estimation in Appendix C). Together with the assumption of no
slip at the gas-liquid interface, we obtain continuity of the velocity field
$v$, i.e.

$\displaystyle\left\llbracket\vec{v}\right\rrbracket=0\quad\text{on}\
\Sigma(t),$ (7)

where $\left\llbracket\psi\right\rrbracket(t,x)=\lim_{h\rightarrow
0^{+}}\psi(t,x+hn_{\Sigma})-\psi(t,x-hn_{\Sigma})$ denotes the jump of a given
(discontinuous) quantity across the interface. The interfacial transmission
condition for momentum in case of constant surface tension reads

$\displaystyle\left\llbracket
p\mathds{1}-S\right\rrbracket\,n_{\Sigma}=\sigma\kappa n_{\Sigma}\quad$
$\displaystyle\text{on}\ \Sigma(t).$ (8)

The right-hand side is the surface tension force with $\sigma>0$ the surface
tension and $\kappa=-\operatorname{div_{\Sigma}}n_{\Sigma}$ the mean
curvature. The motion of the interface is coupled to the bulk flow through the
kinematic boundary condition

$\displaystyle V_{\Sigma}=\vec{v}\cdot n_{\Sigma}\quad$
$\displaystyle\text{on}\ \Sigma(t),$

where $V_{\Sigma}$ is the interface normal velocity. Formally, the latter
condition can be reformulated as the advection equation

$\displaystyle\partial_{t}\chi+\vec{v}\cdot\nabla\chi=0,$ (9)

where

$\displaystyle\chi(t,x)=\begin{cases}1&\text{if}\quad x\in\Omega^{-}(t),\\\
0&\text{if}\quad x\notin\Omega^{-}(t)\end{cases}$

is the indicator function for the liquid phase. The transport equation (9)
forms the basis for the Volume-of-Fluid interface capturing method (see
Section 3.3). In order to regularize††Note that the Navier slip condition with
$L>0$ modifies the singularity at the moving contact line such that it becomes
integrable and the viscous dissipation rate becomes finite. On the other hand,
it has been shown mathematically that the solution cannot be completely
regular if the slip length is finite at the contact line [23]. the moving
contact line singularity [33], the Navier slip boundary condition is applied
for the velocity at the solid wall, i.e.

$\displaystyle\vec{v}\cdot
n_{\partial\Omega}=0,\quad\beta\vec{v}_{\parallel}+(Sn_{\partial\Omega})_{\parallel}=0\quad$
$\displaystyle\text{on}\ \partial\Omega\setminus\Gamma(t),$ (10)

where $L=\eta/\beta>0$ is the so-called slip length. For simplicity, we assume
that the slip length is the same in the liquid and in the gas phase, i.e.
$L=L^{\pm}$. The wettability of the solid is modeled through the contact angle
boundary condition

$\displaystyle n_{\Sigma}\cdot n_{\partial\Omega}=-\cos\theta\quad\text{on}\
\Gamma(t),$ (11)

where the contact angle may be a function of position, modeling regions of
different wettability. Here, for simplicity, we apply a fixed contact angle in
both the hydrophilic part of the boundary $\partial\Omega_{\text{phil}}$ and
the hydrophobic part of the boundary $\partial\Omega_{\text{phob}}$, i.e.

$\displaystyle n_{\Sigma}\cdot
n_{\partial\Omega}=\begin{cases}-\cos\operatorname{\theta_{phil}}&\text{on}\quad\Gamma(t)\cap\partial\Omega_{\text{phil}}\\\
-\cos\operatorname{\theta_{phob}}&\text{on}\quad\Gamma(t)\cap\partial\Omega_{\text{phob}}\end{cases},$

where $0<\operatorname{\theta_{phil}}<\operatorname{\theta_{phob}}<\pi$ are
constants. Note that this is a simplified model, neglecting, e.g., contact
angle hysteresis.

###  3.3 Volume-of-Fluid Discretization

In the present study, the two-phase flow solver _Free Surface 3D_ (FS3D),
originally developed by Rieber and Frohn [49, 50], is employed to solve the
incompressible two-phase Navier Stokes equations. Here we only discuss the
most important aspects of the discretization and refer to the literature [22,
49, 20, 41, 24] for details.

The two-phase Navier Stokes equations in the Continuum Surface Force (CSF)
formulation [7] are discretized using the finite volume approach on a fixed
Cartesian grid. Within the CSF formulation, the effect of surface tension is
modeled as a singular source term in the Navier Stokes equations and the
equations (6) and (8) are replaced by

$\displaystyle\rho\frac{D\vec{v}}{Dt}-\nabla\cdot S+\nabla p$
$\displaystyle=\rho\vec{g}+\sigma\kappa
n_{\Sigma}\delta_{\Sigma},\quad\nabla\cdot\vec{v}=0\quad\text{in}\quad\Omega,$
(12)

where $\delta_{\Sigma}$ denotes the surface delta distribution on $\Sigma(t)$.
The phase indicator function $\chi$ is replaced by the discrete volume
fraction

$\displaystyle f=\frac{1}{|V_{0}|}\int_{V_{0}}\chi\,dV$ (13)

in each computational cell, which is used to track the location of the
interface. Integration of (9) over a control volume $V_{0}$ yields the
transport equation for the volume fraction, i.e.

$\displaystyle\frac{df}{dt}=-\frac{1}{|V_{0}|}\int_{\partial
V_{0}}\chi\,\vec{v}\cdot\vec{n}\,\text{dA}.$ (14)

Within the geometrical Volume-of-Fluid (VOF) method, the flux on the right-
hand-side of (14) is approximated by means of geometrical methods. Based on a
reconstructed interface geometry, the numerical fluxes for the volume fraction
are computed using an operator-splitting method [54], which decomposes the
transport problem in a series of one-dimensional transport problems along the
coordinate axes. The interface is locally reconstructed as a plane in each
cell (also known as piecewise linear interface calculation - PLIC) [48], where
the Youngs method [60] is used to estimate the interface normal vector based
on the volume fraction field. At the contact line, we employ a three-
dimensional variant of the _Boundary Youngs reconstruction method_ [24].

Within the CSF formulation, the two-phase flow is treated as a single fluid
where the density $\rho$ is volume averaged according to

$\rho=f\rho_{l}+(1-f)\rho_{g}.$

Here the subscripts $l$ and $g$ refer to the liquid and gas phases,
respectively. The averaging of the viscosity in interface cells involves both
arithmetic and harmonic averaging [56, 21].

The time integration is based on an explicit Euler method, where the pressure-
velocity coupling is realized by Chorin’s projection method which leads to an
elliptic equation for the pressure to perform a projection onto the space of
solenoidal velocity fields. Note that the grids for velocity and pressure are
staggered to enhance the stability of the method [28].

The surface tension force is discretized with the _balanced_ CSF method
introduced by Popinet [46]. A height-function representation of the interface
is constructed in order to approximate the mean curvature. It has been
demonstrated that this method is able to significantly reduce spurious
currents at the interface away from the boundary. Following the approach by
Afkhami and Bussmann [1, 2], the height function is also used to indirectly
enforce the contact angle boundary condition. The idea is to linearly
extrapolate the height function at the contact line into a ghost cell layer,
where the slope of the extrapolated interface is determined by the prescribed
contact angle. As a result, the approximated value of the mean curvature is
altered, leading to a “numerical force” that drives the interface towards the
desired contact angle. A drawback of this method is that it may create
spurious currents at the contact line.

In order to ensure the stability of the numerical method, the time step is
chosen according to the stability criterion

$\Delta t=\min\\{(\Delta t)_{\sigma},(\Delta t)_{\eta},(\Delta
t)_{\boldsymbol{v}}\\}$

with the timescales (see Tryggvason et al. [57] for a similar criterion) given
by

$(\Delta t)_{\sigma}=\sqrt{\frac{(\rho_{l}+\rho_{g})(\Delta
x)^{3}}{4\pi\sigma}},\quad(\Delta t)_{\eta}=\frac{\rho_{l}(\Delta
x)^{2}}{6\eta_{l}},\quad(\Delta t)_{\boldsymbol{v}}=\frac{\Delta
x}{\lVert\vec{v}\rVert}_{\infty}.$

##### No-slip and Navier Slip boundary conditions:

To allow for a motion of the contact line, one can either make use of the
_numerical slip_ inherent to the method or prescribe the (staggered) Navier
slip condition. The numerical slip is a property of the advection algorithm
since the latter uses face-centered values to transport the volume fraction
field. In fact, the velocity boundary condition is indirectly enforced within
the finite volume method using the concept of “ghost cells”.

(a) No-slip (“Numerical slip”).

(b) Standard Navier slip.

(c) _Staggered_ Navier slip.

FIG. 3: Ghost-cell based numerical realization of the no-slip, Navier slip and
_staggered_ Navier slip boundary conditions in FS3D (for (a) and (b) see also
Gründing et al. [26]).

  1. (i)

In the classical approach, the no-slip boundary condition is enforced by
extrapolating the tangential velocity field into the layer of ghost cells next
to the physical boundary, such that the tangential velocity interpolates to
zero exactly at the boundary (see Fig. 3(a)), i.e.

$v_{\text{ghost}}=-v.$

The velocity in the ghost cells modifies the discrete viscous forces at the
solid boundary such that the motion of the contact line is inhibited. This
approach enforces the no-slip condition only in the limit of the mesh size
going to zero. Therefore, the face-centered velocity at the boundary cell
layer may be non-zero leading to a motion of the contact line. However, the
extent of numerical slip decreases with increasing mesh resolution, typically
leading to a significant mesh dependence of the solution (see, e.g, [26]).
This numerical effect has been first described in the context of VOF methods
by Renardy et al. [47].

  2. (ii)

The Navier slip boundary condition is implemented following the same approach
by setting a different velocity in the ghost-cell layer, such that the
tangential velocity interpolates to zero at a fixed distance $L>0$ from the
physical boundary (see Fig. 3(b)), i.e.

$\displaystyle v_{\text{ghost}}=v\ \frac{2L-\Delta x}{2L+\Delta x}.$ (15)

It has been demonstrated in the literature [53, 26] that this may lead to
mesh-convergent results if the resolution of the computational mesh is well
below the slip length $L$. However, the physically expected slip length is on
the scale of nanometers [43]. Since the macroscopic length scale of the
problem is one the mm scale, we cannot resolve the nanometer scale with the
present numerical method without massive computational costs. Moreover, note
that for $L>0$ the magnitude of the “counter velocity” $v_{\text{ghost}}$ in
the ghost cell is always less than or equal to $|v_{\text{ghost}}|$ in the
case of no-slip.

##### “Staggered” slip:

The goal of the “staggered slip” boundary condition, which is introduced
below, is to reduce the amount of artificial numerical slip compared to the
standard implementation of the Navier slip condition. The idea is to apply a
given slip length with respect to a _virtual_ “staggered” boundary located in
between the physical boundary and the location of the face-centered velocity
nodes (see Fig. 3(c)). Mathematically, the distance $\tilde{L}$ to the virtual
boundary can be expressed as

$\displaystyle\tilde{L}=L+\frac{\omega\Delta x}{2},$ (16)

where $L$ is the distance to the physical boundary. Consequently, the applied
ghost velocity is

$\displaystyle v_{\text{ghost}}=v\ \frac{2\tilde{L}-\Delta
x(1+\omega)}{2\tilde{L}+\Delta x(1-\omega)},\quad\omega\in[0,1].$ (17)

Note that here we keep $\tilde{L}$ fixed at the (physically) prescribed value
when refining the mesh. The virtual boundary allows to apply larger “counter
velocities” in the ghost cell layer than the standard Navier slip or no-slip
conditions. Through (discrete) viscous forces this leads to a stronger damping
or numerical dissipation. The quantity $\omega\in[0,1]$ is an empirical
parameter that allows to control the amount of damping (see Appendix A). For a
staggered slip length $\tilde{L}$ smaller than $\omega\Delta x/2$, the
staggered slip can be interpreted as standard Navier slip with a _negative_
slip length, while $\tilde{L}=\omega\Delta x/2$ corresponds to $L=0$
(numerical slip). Note that a small staggered slip length may lead to very
small numerical timesteps since the velocity $v_{\text{ghost}}$ in the ghost
cells may become very large. Therefore, one cannot choose arbitrary small
$\tilde{L}$ in practice.

In the numerical simulations of the breakup process, we observed spurious
velocities at the contact line on the hydrophilic stripe leading to capillary
waves propagating along the droplet. The staggered slip condition allows to
significantly dampen these unphysical velocities at the contact line.
Moreover, it is found that the breakup dynamics of the capillary bridge is
unaffected by the choice of $\omega$; see Appendix A for more details. Hence,
the staggered slip condition only serves as a stabilization method in the
present study. Unless stated otherwise, we choose $\omega=2/3$ and
$\tilde{L}=500\,\text{nm}$ in the numerical simulations. See Appendix A and
[22] for more information on the staggered slip condition and its effect on
dynamic wetting simulations.

$f$
---
$0.0$
---
$1.0$
---
$0.2$
---
$0.4$
---
$0.6$
---
$0.8$
---

$\tilde{\Sigma}$
---

FIG. 4: Initialization of the volume fraction field from STL-data. (a)
Overview of the STL surface from Surface Evolver from Figure 1 (blue),
together with the background mesh required for FS3D (greyscale). The enlarged
image shows the STL surface with the resulting volume fraction in greyscale.
(b) 2D illustration of the intersection between the background mesh (squares)
and a volume bounded by the discrete surface $\tilde{\Sigma}$ (light grey)
that may lead to non-convex intersection volumes (dark grey).

##### Computing volume fractions from the surface mesh given by Surface
Evolver

As described in Section 3.1, the gas-liquid interface as shown in Figure 1 has
been computed as a minimal surface using Surface Evolver. This tool provides
the interface in the form of a triangulated surface (an STL-file) such as the
blue triangulated surface shown in Figure 4. The dynamic breakup simulation
using the VOF method requires the initial condition for the volume fraction
defined in (13). The computation of the initial volume fraction field $f$ uses
the triangulated surface generated by Surface Evolver as indicated by the
greyscale in the Figure 4. Detailed information on the algorithm used for
initializing volume fractions from surface meshes is available in [55], here
we only briefly outline the algorithm. The calculation of volume fractions
from a surface mesh relies on the intersection of the surface mesh with each
cell of the background mesh, as illustrated for the 2D case in Figure 4. Here,
the background mesh points are indicated by circles, while squares indicate
the mesh points of the STL surface. Arrows indicate the orientation of the STL
surface. Each intersected cell is decomposed into tetrahedra to simplify the
geometrical intersection, shown as triangles in Figure 4 in the cell with the
phase-specific volume indicated with a darker gray color. A sphere is defined
for every tetrahedron in an intersected cell that encompasses the tetrahedron,
schematically shown with the dashed line. The tetrahedron centroid is its
center, and the maximal Euclidean distance from the tetrahedron centroid to
its corner points is the sphere radius. The tetrahedron sphere identifies a
subset of triangles from the surface that intersects the tetrahedron. The
phase-specific volume contained within the tetrahedron (black triangle in
Figure 4), for which $\chi=1$, is found by geometrically intersecting the
tetrahedron with the triangulated surface subset. Therefore, the total phase-
specific volume of a cell is the sum of the phase-specific volumes from all
the tetrahedrons of the cell (total shaded volume in the cell in Figure 4).
The volume fraction is then the ratio of the total phase-specific volume in
the cell and the total cell volume. This algorithm returns a highly accurate
geometrical value of ${0<\chi<1}$ for the intersected cells and its accuracy
depends on the resolution of the surface mesh generated by the Surface
Evolver.

FIG. 5: Setup for FS3D using symmetry conditions.

##### Numerical setup for FS3D:

The setup for FS3D is shown in Figure 5. In order to save computational
resources, we make use of the symmetry of the problem and simulate only a
quarter of the droplet, while applying symmetry boundary conditions at the
respective symmetry boundaries ($x=0$ and $y=0$). At the outer boundaries of
the domain (i.e., $x=x_{\text{max}}$, $y=y_{\text{max}}$ and
$z=z_{\text{max}}$ ), we apply no-slip for the velocity, which turns out to be
irrelevant for the breakup dynamics (compared to, e.g., fixed pressure outflow
boundary conditions). The staggered Navier slip condition is applied at the
solid boundary ($z=0$). The computational domain has a size of
$1000\operatorname{{\mu}m}\times 2000\operatorname{{\mu}m}\times
500\operatorname{{\mu}m}$ in the $x$-, $y$\- and $z$-direction and is
subdivided into an equidistant mesh of $2N\times 4N\times N$ cells, where $N$
is varied between $48$ and $128$. Therefore, the maximum resolution of the
mesh is $500\operatorname{{\mu}m}/128\approx 3.9\operatorname{{\mu}m}$ in each
direction. Unfortunately, adaptive mesh-refinement is currently not available
for our version of FS3D. The physical parameters used in the numerical
simulation are listed in Table 2.

$\rho_{l}$ [$\text{kg}/\text{m}^{3}$] |  $\rho_{g}$ [$\text{kg}/\text{m}^{3}$] |  $\eta_{l}$ [$10^{-6}\frac{\text{kg}}{\text{m}\cdot\text{s}}$] |  $\eta_{g}$ [$10^{-6}\frac{\text{kg}}{\text{m}\cdot\text{s}}$] |  $\sigma$ [$\frac{\text{mN}}{\text{m}}$]
---|---|---|---|---
$997.05$ | $1.17$ | $890.45$ | $18.5$ | $71.96$
Table 2: Physical parameters for the numerical simulations (literature values
for $T=298\,\text{K},\,p=100\,\text{kPa}$).

##  4 Results and Discussion

###  4.1 Qualitative Comparison between Experiments and Numerics

Since the focus of the present study is on the dynamics of the capillary
bridge on the hydrophobic stripe, this part of the droplet is shown in more
detail in Fig. 6. It shows the bridge at three different instants in time
$\tau$ before the breakup event, i.e. $\tau=t_{0}-t$, where $t_{0}$ is the
breakup time, for both, experiment and simulation. The geometrical parameters
are $w_{\mathrm{phil}}=500\operatorname{{\mu}m}$ and $\alpha=1$. In the top
row, the black and dark gray regions represent the liquid (with a reflection
of the light in the middle of the capillary bridge), while the liquid is white
in the bottom row. The images in both rows also show some liquid that is
wetting the hydrophilic stripes in the top and bottom part of each frame.

FIG. 6: Qualitative comparison between experiment and simulation for
$w_{\mathrm{phil}}=500\operatorname{{\mu}m}$ and $\alpha=1$ at the same
instants in time before breakup. The scale bar is valid for both experiment
and simulation.

Qualitatively, in both experiments and simulations, the capillary bridge
develops from a catenoid type to a narrower shape. Then a liquid thread is
formed that is getting constricted at two points. At $\tau$ = 0.4267 ms, which
is close to the initial condition for the simulations, the minimal width $d$
of the capillary bridge in both simulation and experiment is approximately the
same. This is also true for smaller $\tau$ until a liquid thread forms in the
final instants before breakup at $\tau$ = 0.02667 ms, when the minimum width
of the bridge, as well as the bulge that forms in the middle, is larger in the
simulations. Note that $d$ is initially located in the middle of the capillary
bridge. After the elongated thread has formed, the position of the minimum
width moves away from the center, as indicated at $\tau$ = 0.02667 ms in Fig.
6. Below we will show that this is the starting point of a second dynamic
regime. Besides the bulge and the minimal width $d$, also the length of the
thread differs between experiment and simulation. The width of the capillary
bridge at the instance when it becomes unstable is of the order of 200
$\operatorname{{\mu}m}$. On this length scale, unavoidable contaminations can
lead to significantly smaller contact angles than the measured receding
contact angle (see e.g. Park et al. [44] and the supplementary material of
Hartmann and Hardt [30]). These phenomena can neither be captured in the
simulations nor in the calculation of the initial droplet configuration with
Surface Evolver. This might explain the deviations in the initial
configuration between experiments and simulations, for example, the slightly
bigger curvature of the capillary bridge in the plane parallel to the surface
which is obvious at $\tau$ = 0.4267 ms. Also the absolute value of the
parameter $R$, which is the maximum width of the capillary bridge at the
boundary between the hydrophobic and the hydrophilic stripe, differs slightly
between experiment and simulation, though it stays constant in both.
Because not all of the influencing factors relevant in the experiments could
be incorporated in the simulations, quantitative agreement between experiments
and simulations is not expected in all aspects of the static and dynamic
wetting phenomena observed. For quantitative comparison, the focus, therefore,
lies on one specific aspect, namely the breakup dynamics of the liquid bridge,
i.e. the time evolution of the narrowest section of the bridge.

###  4.2 Phase Space Representation of the Breakup Dynamics

Following the literature, the breakup dynamics is usually described via the
minimum width $d$ as a function of the time $\tau$ before the breakup event,
i.e. $\tau=t_{0}-t$, where $t_{0}$ is the breakup time. For the inviscid
breakup of a free capillary bridge, it can be shown by means of an asymptotic
analysis that this function follows the power law (1). However, the precise
time of the breakup event is hard to determine both in experiments (due to
finite spatial and temporal resolution) and in the simulation. Note that the
choice of $t_{0}$ can have a large effect on the effective exponent that is
extracted from the data. It has been reported that the same set of data
appears to be represented by power laws $d(\tau)\propto\tau^{\nu}$ with an
exponent $\nu$ ranging from $0.6$ to $0.8$ depending on the choice of $t_{0}$
[10]. Both Li and Sprittles [40] and Deblais et al. [14] eliminated the
dependency on the breakup time by rewriting the power law (1) according to

$\displaystyle
d=C\left(\frac{\sigma}{\rho}\right)^{1/3}\,(t_{0}-t)^{2/3}\quad\Leftrightarrow\quad
d^{3/2}=C^{3/2}\left(\frac{\sigma}{\rho}\right)^{1/2}(t_{0}-t).$ (18)

Hence, relation (1) holds if and only if the quantity $d^{3/2}$ is linear in
time and, in this case, the value of $C$ can be found from the slope (i.e. the
time derivative) of $d^{3/2}$.

(a) $t_{0}=0.632\,\text{ms}$, $\nu=2/3$.

(b) $t_{0}=0.622\,\text{ms}$, $\nu=0.4755$.

FIG. 7: Breakup dynamics for different choices of the breakup time $t_{0}$ for
the numerical data.

In the Volume-of-Fluid simulation, the actual breakup is usually mesh-
dependent since it is ultimately performed by the interface reconstruction
algorithm. Moreover, since the breakup process involves very small length
scales, it cannot be fully resolved by the numerics. Therefore, the numerical
results can only be considered meaningful down to a certain length scale
determined by the computational mesh. In the present study, this length scale
is approximately $20\,\mu m$. Consequently, the breakup time cannot be
extracted from the numerics in a meaningful way without extrapolating the
data.

The experimental value of the breakup time is determined from the pictures
taken by the high-speed camera. The first image where the capillary bridge is
pinched off defines the time $t^{\star}$ which is always larger than the real
breakup time $t_{0}$. Clearly, the breakup must occur between $t^{\star}$ and
the time associated with the previous image. We estimate the breakup time to
be $t_{0}=t^{\star}-\Delta t/2$, where $\Delta t=1.33\cdot 10^{-2}$ ms follows
from the frame rate of the high-speed camera which is 75,000 fps. To account
for the uncertainty in breakup time in the case of experiments, horizontal
error bars with a total length of $\Delta t$ are drawn. In the vertical
direction, the error bar represents the standard deviation obtained from 5
experiments.

A concrete example for the sensitivity of the simulation results with respect
to the choice of the breakup time is given in Fig. 7. Besides the experimental
results, Fig. 7(a) shows the simulation data for the choice
$t_{0}=0.632\,\text{ms}$. The latter value is found by fitting the data in the
final regime $d\lesssim 50\operatorname{{\mu}m}$ with the exponent $\nu=2/3$.
Fig. 7(b) shows the same numerical data set for a slightly smaller breakup
time $t_{0}=0.622\,\text{ms}$. In this case, the numerical data agrees well
with a power law with a smaller exponent $\nu\approx 0.48$, which has been
obtained from the systematic phase space analysis described below. Note that
the difference between the two choices $\Delta t_{0}=10\,\mu\text{s}$ is
smaller than the inverse frame rate of the high-speed camera which is
approximately $\Delta t=13\,\mu\text{s}$. This example clearly shows the need
for a systematic method that does not rely on the breakup time.

##### Phase space representation.

We apply a different approach to describe the breakup dynamics which is
_independent_ of the choice of the breakup time and allows to identify
different dynamic regimes in a systematic way. The basic idea is to use the
variables $(d,-\dot{d})$ rather than $(\tau,d)$ to describe the evolution of
the minimum width. This approach was first applied by Li and Sprittles [40]
and is systematically elaborated in the following.

Given the minimum width of the capillary bridge as a function of physical time
$t$, we consider the _breakup speed_ , i.e. the time derivative

$V=-\dot{d}$

as a function of the _minimum width_ itself, i.e. we formally define

$\displaystyle
V(\tilde{d}):=-\dot{d}(\tilde{t}),\quad\text{where}\quad\tilde{t}:=d^{-1}(\tilde{d}).$
(19)

The latter quantity is well-defined since the minimum width $d$ is a
monotonically decreasing function (over the time scale of the breakup
process). Obviously, the function $V(d)$ is invariant with respect to shifts
in the time coordinate. Moreover, it can be shown mathematically that the
remaining information of the function $d(\tau)$ is still contained in $V(d)$.
Indeed, the function $d(\tau)$ can be reconstructed from $V(d)$ up to a shift
in time; see [22] for details.

For the power law (1) describing the inviscid regime we have

$d(t)=C\left(\frac{\sigma(t_{0}-t)^{2}}{\rho}\right)^{1/3}\quad\Rightarrow\quad\dot{d}(t)=-\frac{2}{3}C\left(\frac{\sigma}{\rho(t_{0}-t)}\right)^{1/3}=-\frac{2}{3}C^{3/2}\left(\frac{\sigma}{\rho\,d(t)}\right)^{1/2}.$

Hence the power law (1) translates to

$\displaystyle V(d)=\frac{2}{3}C^{3/2}\left(\frac{\sigma}{\rho
d}\right)^{1/2}.$ (20)

More generally, one can easily show the relation

$\displaystyle d(t)=c(t_{0}-t)^{\nu}\quad\Rightarrow\quad V(d)=\nu
c^{1/\nu}\,d^{1-1/\nu}=\tilde{c}\,d^{\tilde{\nu}},$ (21)

which is valid for an arbitrary power law ($\nu,\,c>0$). So the exponent
$\tilde{v}$ and the prefactor $\tilde{c}$ obtained from the phase space
diagram can be transformed into the standard representation via

$\displaystyle\nu=\frac{1}{1-\tilde{\nu}},\quad
c=\left(\frac{\tilde{c}}{\nu}\right)^{\nu}=[(1-\tilde{\nu})\,\tilde{c}]^{1/(1-\tilde{\nu})}.$
(22)

For the material parameters of water in air (temperature $T=298\,\text{K}$,
pressure $p=100\,\text{kPa}$, see Table 2), the relations (1) and (20) for a
free liquid bridge take the form

$\displaystyle d(\tau)=C\cdot
416\,\hat{\tau}^{2/3}\,\operatorname{{\mu}m},\quad V(d)=C^{3/2}\cdot
5657\,\hat{d}^{-1/2}\,\frac{\operatorname{{\mu}m}}{\text{ms}},$

where $\tau=\hat{\tau}\,\text{ms}$ and $d=\hat{d}\operatorname{{\mu}m}$. In
the present paper, these relations will be referred to as “inviscid theory”.

Note that the above method requires to differentiate potentially noisy data
with respect to time. This issue has been addressed via filtering out high-
frequency oscillations in the experimental values of $V(d)$ by locally fitting
a straight line to the data (using six neighboring points). Despite this
difficulty, the method allows studying the breakup dynamics in detail without
the uncertainty in choosing $t_{0}$. In the following, we will report both
$V(d)$ and $d(t_{0}-t)$ for completeness.

###  4.3 Quantitative Comparison between Experiments and Simulations

####  4.3.1 The case $\alpha=1$:

FIG. 8: Minimum bridge width as a function of time before breakup, experiment
vs. simulation ($\alpha=1$).

We first consider the case $\alpha=1$. Figure 8 compares the experimental data
with numerical simulations for $\operatorname{\theta_{phob}}=102^{\circ}$ and
$\operatorname{\theta_{phil}}=28^{\circ}$. Before the onset of breakup process
($d\gtrsim 180\operatorname{{\mu}m}$), the bridge width in the experiments is
decreasing due to evaporation. This continues until the critical width is
reached, indicated by the starting point of the simulation curve in Fig. 8
corresponding to the configuration computed by Surface Evolver.
In the region between about $100\operatorname{{\mu}m}$ and
$180\operatorname{{\mu}m}$, the numerical results are in good quantitative
agreement with the experimental data and the inviscid theory for $C\approx
0.75$. However, a careful inspection shows that the data in this region do not
exactly lie on a straight line in the double logarithmic diagram. Indeed there
is a non-vanishing curvature visible. Below approximately
$100\operatorname{{\mu}m}$, there is a transition to a second dynamic regime
that follows a power law with an exponent $\nu$ close to 0.5. This value is
significantly smaller than the one for a free capillary bridge in the inviscid
regime. Note that the exponent is obtained from fitting the numerical data in
the phase space diagram, see below.

FIG. 9: Breakup speed ($\alpha=1$) as a function of the bridge width compared
to 3 different experiments. The pictures inside of the graph show the
(closest) images for simulations (sim) and experiments (exp) to the minimum
bridge width indicated in the top row.

##### Phase space representation:

Figure 9 shows the same data set in the phase space diagram. This
representation of the data allows a much more detailed study of the breakup
process independent of the breakup time. Note that the experimental data refer
to three _individual_ repetitions of the same experiment. By considering
individual experiments rather than averaged quantities, we are able to
visualize detailed features of the dynamics that might otherwise be averaged
out. Moreover, note that the noise in the experimental data for the bridge
width translates to a noise in the calculated breakup speed. This explains the
oscillations in the experimental data which are not present in the numerical
data.

The phase space diagram in Fig. 9 shows that the dynamics in the region
$100\operatorname{{\mu}m}$ and $180\operatorname{{\mu}m}$ agrees with the
inviscid theory for $C\approx 0.75$ only in an average sense. In fact, there
is no clear scaling relation in the latter region. At approximately
$70\operatorname{{\mu}m}$, the speed in the experiment reaches a local maximum
of approximately $570\operatorname{{\mu}m}/\text{ms}$ on average (indicated by
the dashed horizontal line in Fig. 9) before it starts to decrease further
towards a local minimum. This behavior is also clearly visible in the
numerical data. We note that the position of the local minimum differs for
different repetitions of the experiment (see Fig. 9). Note that the simulation
predicts a smaller local maximum of approximately
$500\operatorname{{\mu}m}/\text{ms}$.

FIG. 10: Breakup speed $V$ and position of the local minima $x_{\text{min}}$
(measured relative to the center of the hydrophobic stripe) in the numerical
simulation ($\alpha=1$).

Moreover, the analysis of the numerical data shows that the local minimum in
the breakup speed corresponds to the time instant in the process when the
liquid filament has formed and the position of the minimum width
$x_{\text{min}}$ starts to move symmetrically††The symmetry is enforced
explicitly in the simulation but also observed in the experiment; see Fig. 6.
from the middle of the liquid bridge ($x=0$) towards the points of the final
pinch-off; see Fig. 10. So, in fact, this can be understood as the starting
point of a second dynamic regime. In the experiments, this final breakup
behaves slightly differently. Though a filament has formed at the time when
the breakup speed reaches a local minimum, the width of this filament first
further decreases, while the filament has approximately the same width along
its axis without any constriction point. Compared to the simulations, the
location where the bridge has its minimum starts to move in a symmetrical
manner to the constriction points at a later instance in time, i.e. for a
smaller value of $\tau$.
As can be seen from both the experimental and the numerical data in Fig. 9,
the second dynamic process is not consistent with the inviscid theory for a
free capillary bridge (1). In fact, by fitting the numerical data, the breakup
velocity in the second breakup regime follows the power law

$V(d)=23879\,\frac{\operatorname{{\mu}m}}{\text{ms}}\cdot\left(\frac{d}{\operatorname{{\mu}m}}\right)^{-1.1029}$
(23)

which corresponds to

$d(\tau)=171.95\,\operatorname{{\mu}m}\cdot\left(\frac{\tau}{\text{ms}}\right)^{0.4755};$
(24)

see Fig. 8. Note that the last data point ($d\approx 21\operatorname{{\mu}m}$)
from the numerical simulation in Fig. 9 has been excluded from the fit since
the breakup speed appears to reach another local maximum in this region.
Hence, the numerical data suggest that another regime transition may happen at
smaller scales that cannot be resolved with the current approach. We attribute
the observed discrepancy to the inviscid scaling law to the complex geometric
configurations found in the second dynamic regime. A simple scaling law such
as (1) is only expected to work if the configuration of a system is
essentially rescaled during its time evolution. We do not expect that the
scaling law to capture more complex scenarios such as a shift of the position
of the minimum width of the filament. In fact, the self-similarity of the
interface shape, which is an essential assumption for the derivation of (1)
(see [36]), is lost due to the presence of the solid substrate. Moreover, as
it has been pointed out by Deblais et al., “the approach to asymptotic power
laws can be slow and may pass through one or several transient regimes before
the final universal regime is reached” [14]. Hence the exponent $\nu$ may
change on smaller length scales below the resolution of the present
experiments and simulations.

(a) Breakup speed for different densities.

(b) Dimensionless form of the breakup speed.

FIG. 11: Variation of the liquid density and gas density for $\alpha=1$,
$\operatorname{\theta_{phob}}=102^{\circ}$ and
$\operatorname{\theta_{phil}}=28^{\circ}$.

Influence of inertial forces: We vary both the liquid density $\rho_{l}$ and
the gas density $\rho_{g}$ in the simulation to study how inertial forces
modify the breakup dynamics. Fig. 11(a) shows the breakup speed in the phase
space diagram for four different values of the liquid and gas densities, i.e.

$\rho_{l}=K\cdot 997.05\ \text{kg}/\text{m}^{3},\quad\rho_{g}=K\cdot 1.17\
\text{kg}/\text{m}^{3}\quad\text{where}\quad K\in\\{0.75,1,1.25,1.5\\}.$

As expected, the breakup speed decreases with increasing density. Remarkably,
the data can be collapsed onto a single master curve by an appropriate
nondimensionalization. Fixing a length scale $L$ in (1) yields

$\displaystyle\frac{d(\tau)}{L}=C\left(\frac{\sigma\tau^{2}}{\rho
L^{3}}\right)^{1/3}.$

Hence, with the choice

$T=\sqrt{\frac{\rho L^{3}}{\sigma}}$

for the time scale, equation (1) takes the form

$\displaystyle\hat{d}(\hat{\tau})=C\,\hat{\tau}^{2/3},$ (25)

where $\hat{d}=d/L$ and $\hat{\tau}=\tau/T$. The corresponding velocity scale
is

$\displaystyle V_{\text{scale}}=\frac{L}{T}=\sqrt{\frac{\sigma}{\rho
L}}\propto\frac{1}{\sqrt{\rho}}.$ (26)

Indeed, plotting the non-dimensional breakup speed
$\hat{V}=V/V_{\text{scale}}$ for

$L=w_{\text{phob}}=500\mu\text{m}$

collapses the data onto a single curve; see Fig. 11(b). This indicates that
the observed process is indeed dominated by a balance of inertial and
capillary forces despite the fact that the dynamics does _not_ simply follow
the power law (1).

Mesh study: To quantify the influence of the numerical discretization, Fig.
12(a) shows numerical data at the breakup speed for different mesh resolutions
($8N^{3}$ computational cells where $N\in\\{64,96,128\\}$). Besides some
oscillations in the breakup speed on the coarse mesh, the results appear to be
reasonably mesh-independent for $d$ larger than approximately
$40\operatorname{{\mu}m}$, whereas the simulation on the coarsest mesh
($N=48,\,\Delta x\approx 10.4\operatorname{{\mu}m}$) appears to be under-
resolved. The finest mesh ($N=128,\,\Delta x\approx 3.9\operatorname{{\mu}m}$)
delivers reasonable results for $d\gtrsim 10\operatorname{{\mu}m}$. This
corresponds to only $2.5$ computational cells within the bridge width. Since
mesh convergence cannot be assured on this scale, the numerical data below
approximately $20\operatorname{{\mu}m}$ should not be used to draw a
quantitative conclusion.

(a) Mesh dependence study ($\tilde{L}=500\operatorname{nm}$).

(b) Slip-length dependence study ($N=128$).

FIG. 12: Mesh and slip-length dependence study for $\alpha=1$,
$\operatorname{\theta_{phob}}=102^{\circ}$ and
$\operatorname{\theta_{phil}}=28^{\circ}$.

Influence of slip and viscous effects: As shown in Fig. 12(b), where the
staggered slip length $\tilde{L}$ for $\omega=2/3$ is varied from
$100\,\text{nm}$ to $500\,\text{nm}$, the choice of the (staggered) slip
length shows only a minor influence on the breakup dynamics. This behavior is
to be expected in the inviscid regime (low $\operatorname{Oh}$, here
$\operatorname{Oh}=4.7\cdot 10^{-3}$), where viscous dissipation due to slip
is irrelevant for the breakup process. A slight monotonic decrease of the
breakup speed with the slip length is visible when the breakup speed
approaches the local minimum at approximately $45\operatorname{{\mu}m}$. In
all subsequent numerical simulations, we apply a staggered slip length of
$500\,\text{nm}$ and $\omega=2/3$. In Appendix A it is shown that the choice
of $\omega$ has no relevant influence on the breakup dynamics. The influence
of the dynamic viscosity on the dynamics in the simulation is studied in Fig.
13(a). The green curve is computed with an $50\%$ increased dynamic viscosity
in both the liquid and the gas phase. Despite this large increase in the
viscosity, there is only a small decrease in the breakup speed of
approximately $4\%$. In particular, the overall shape is unchanged including
the position of the local minimum and the dynamic exponents. From this, we
conclude that viscous effects including boundary slip (see Fig. 12(b)) play no
significant role in the considered range of parameters. Instead, the process
is controlled mainly by a balance of inertial and capillary forces (see also
Fig. 11).

(a) Viscosity variation.

(b) Contact angle variation.

FIG. 13: Influence of the dynamic viscosity and the hydrophobic contact angle.

Influence of the wetting conditions: The influence of the wetting conditions
through the contact angle
$\operatorname{\theta_{phob}}\in\\{90^{\circ},102^{\circ},110^{\circ}\\}$ on
the hydrophobic stripe is studied in Fig. 13(b). Note that we employed Surface
Evolver to compute individual initial liquid surfaces for each case. The
results show that the dynamics stays qualitatively similar, while the whole
graph is shifted to the right with decreasing $\operatorname{\theta_{phob}}$.
The colored lines show fits to the second dynamic regime. Note that the
difference in the exponents is very small. The numerical values for the
exponent $\nu$ and the constant $c$ obtained from the phase space analysis are
given by

$\displaystyle\operatorname{\theta_{phob}}=90^{\circ}:\nu=0.47596,\quad c$
$\displaystyle=179.3147\operatorname{{\mu}m}$
$\displaystyle\operatorname{\theta_{phob}}=102^{\circ}:\nu=0.47553,\quad c$
$\displaystyle=171.9538\operatorname{{\mu}m}$
$\displaystyle\operatorname{\theta_{phob}}=110^{\circ}:\nu=0.47623,\quad c$
$\displaystyle=162.5795\operatorname{{\mu}m}.$

These values have been obtained from a least-squares fit of the data starting
from the local minimum in the breakup speed. Note that the last data point in
Fig. 13(b) has been excluded from the fit for
$\operatorname{\theta_{phob}}=102^{\circ}$ and
$\operatorname{\theta_{phob}}=110^{\circ}$. Similarly, the last two data
points have been excluded from the fit for
$\operatorname{\theta_{phob}}=90^{\circ}$. In summary, it can be concluded
that the wetting condition, i.e. the contact angle on the hydrophobic stripe,
has no major influence on the breakup dynamics in the second dynamic regime.

####  4.3.2 The case $\alpha=0.5$

FIG. 14: Breakup dynamics for $\alpha=0.5$. Left: Minimum bridge width as a
function of time. Right: Breakup speed as a function of the bridge width.

The breakup dynamics for $\alpha=0.5$ is shown in Fig. 14. The overall
dynamics is qualitatively similar to the case $\alpha=1$. After a short
initial phase, a good quantitative agreement between the experimental data and
the numerical data is found. In the region $90\operatorname{{\mu}m}\lesssim
d\lesssim 160\operatorname{{\mu}m}$, the dynamics is consistent to the
prediction of the inviscid theory with $C\approx 0.77$ in an averaged sense.
At a bridge width of about $60\operatorname{{\mu}m}$, the experimental breakup
speed reaches a maximum at approximately $630\operatorname{{\mu}m}/\text{ms}$,
which is captured accurately by the numerics. For smaller bridge widths the
velocity decreases until a local minimum is reached. In the numerical data,
the local minimum lies at approximately 32 $\operatorname{{\mu}m}$, while the
minimum in the experiment is shifted to smaller widths (except for Exp3).
After the minimum is reached, a second dynamic regime is initiated which is
similar to the one obtained for $\alpha=1$. A least-squares fit of the
numerical data in the phase space representation yields the exponent
$\nu\approx 0.48$ for the second dynamic regime.

####  4.3.3 The case $\alpha=1.5$

FIG. 15: Breakup dynamics for $\alpha=1.5$. Left: Minimum bridge width as a
function of time. Right: Breakup speed as a function of the bridge width.

The results for $\alpha=1.5$ are displayed in Fig. 15. The dynamics of the
breakup process in the simulation is qualitatively similar to the cases
$\alpha=1$ and $\alpha=0.5$. Shortly after the onset of the instability at
approximately $280\operatorname{{\mu}m}$, the dynamics in the numerics is (on
average) consistent with the inviscid theory for $C\approx 0.75$. At
approximately $90\operatorname{{\mu}m}$ the breakup speed reaches a maximum
before it decreases toward a local minimum at approximately
$56\operatorname{{\mu}m}$. Below this length scale, a second dynamic regime
occurs. The observed exponent $\nu\approx 0.47$ is very close to the values
observed for $\alpha=1$ and $\alpha=0.5$. Interestingly, the numerical data
for the breakup speed show another local maximum at $d\approx
30\operatorname{{\mu}m}$.

In contrast to the cases $\alpha=1$ and $\alpha=0.5$, there is a significant
discrepancy between experiment and simulation at the onset of the breakup
process. While the breakup process in the simulation starts at approximately
$280\operatorname{{\mu}m}$, the onset of instability in the experiments occurs
at a width of the capillary bridge of approximately
$200\operatorname{{\mu}m}$. It appears that the whole process is shifted
towards smaller values of the bridge width compared to the numerical data. In
particular, in the simulation the position of the local minimum is
significantly shifted towards larger values of the minimum width. In fact, the
initial shape of the droplet delivered by Surface Evolver shows a significant
deviation from the experimentally observed shape. Notably, there is a
different curvature of the interface in the plane of the substrate for all
investigated values of $\alpha$; see Fig.19 in Appendix B. In the experiments,
the volume decreases due to evaporation before the critical bridge width is
reached. Consequently, the minimum width decreases and the pressure within the
capillary bridge increases. Shortly before the point of instability is
reached, the capillary system responds to the increased pressure within the
bridge by increasing the curvature, and therefore the pressure, above the
hydrophilic stripes. As a consequence, the contact angle
$\Theta_{\mathrm{phil}}$ increases, since liquid is pumped from the region of
highest pressure (location of minimum bridge width) towards the liquid
reservoirs above the hydrophilic stripes (see [30]). This process cannot be
captured by the Surface Evolver simulations and, as a consequence, the initial
curvature in the numerical simulations deviates from the one in the
experiments when the capillary structure reaches the point of instability.
This deviation might explain the shift in the dynamics.

###  4.4 Rayleigh-Plateau Instability

In Fig. 16a, the final moments of bridge breakup are shown for
$w_{\mathrm{phil}}=500\operatorname{{\mu}m}$ and $\alpha=1$. Emerging from a
catenoid type structure, as already shown in Fig. 6, a liquid thread (a
cylinder cut by a plane) with diameter $D$ is formed between the tips of two
cones with opening angle $\beta$ ($\tau=0.04$ ms). The thread then gets
pinched at the tips of the cones ($\tau=0.013$ ms) and finally breaks up,
leaving a primary droplet in the middle of the hydrophobic stripe
($\tau=-0.013$ ms). Besides the primary droplet, some smaller secondary
droplets can be seen in the experiments. These droplets have their origin in a
self-similar breakup process. This phenomenon is also found in the breakup of
a soap film that is spanned between two circular rings [11], or in the breakup
of a liquid jet [37, 17].

FIG. 16: Analysis of the final breakup behavior. In (a), the capillary bridge
has evolved into a liquid cylinder between two cones. The final breakup
happens at two pinch-off points leaving behind a primary droplet with smaller
secondary droplets in the middle of the hydrophobic stripe. In (b), the cone’s
opening angles $\beta$ are evaluated for different $\alpha$’s and hydrophilic
stripe widths. In (c), the diameter $D$ of the cylindrical segment is plotted
versus the measured wavelength $\lambda$ for
$w_{\mathrm{phil}}=500\operatorname{{\mu}m}$ and different $\alpha$. Each
point represents one experiment. These are condensed in the ellipses,
representing the standard deviation around the mean value for each $\alpha$.
The dashed line corresponds to equation (27).

The described breakup process exhibits two correlated features: First of all,
cones, which in the two-dimensional projection appear as triangles, are formed
in each experiment (at all observed values of $\alpha$) at a certain point in
time. This is close to the point where a liquid thread is formed. As can be
seen in Fig. 16b, these cones have unique opening angles $\beta$ for each
$\alpha$. Their structure is independent of the hydrophilic stripe width. The
value of $\beta$ is highly reproducible between different experiments.
The liquid thread breaks up in a self-similar process. Fig. 16c shows all data
of experiments performed with a hydrophilic stripe width of
$500\leavevmode\nobreak\ \operatorname{{\mu}m}$ for $\alpha=0.5-1.5$. The
diameter $D$ of the cylindrical thread is plotted versus the parameter
$\lambda$, which is the distance between the two pinch-off points and
simultaneously the distance between the two cone tips in Fig. 16a. Note that
these data points, as well as the ones depicted in figure 16b, originate from
different experiments that were performed using a microscope. This setup has
the advantage of a higher optical resolution. As a drawback, it does not allow
side view imaging of the droplet. However, the side view data is needed for
the comparison of the experiments with the simulations. Different colors are
used for different $\alpha$’s, and the mean values and standard deviations for
each $\alpha$ are represented by ellipses around black crosses. The values are
compared with the classical Rayleigh-Plateau instability, which predicts that
a cylindrical liquid jet of diameter $D$ becomes unstable since a
decomposition into droplets is energetically more favorable. $D$ is connected
to the wavelength $\lambda$ that belongs to the fastest growing mode according
to the formula (see e.g. [15, 18, 19])

$\frac{\pi D}{\lambda}\approx 0.7.$ (27)

From Fig. 16c it can be seen that, within the experimental error, the final
breakup process on the substrate follows the classical Rayleigh-Plateau
instability. The measured receding contact angle is 102∘, which is slightly
higher than 90∘. For the case of 90∘ a liquid jet with unpinned contact lines
becomes unstable following the classical Rayleigh-Plateau result (see, e.g.,
Bostwick and Steen [5, 6]). This also follows from symmetry considerations.
Remarkably, the wavelength of the fastest growing mode corresponds to the
distance between the tips of the cones and the distance between the secondary
satellite droplets. The reasonable agreement with the classical Rayleigh-
Plateau theory, which is based on an inviscid fluid, also indicates that the
influence of viscous stresses that shows up in the final stages of breakup, is
not significant for the decay of the thread. The scattering of the data points
of each experiment might be due to the uncertainty in the receding contact
angle. This is already discussed in Section 4.1. Furthermore, a part of the
data scatter visible in Fig. 16c is due to the spatial resolution of the
camera. Since the final breakup process occurs at a relatively small time
scale ($\approx$0.04 ms in case of $w_{\mathrm{phil}}$=500 $\mu$m and $\alpha$
= 1, see Fig. 16a), a high frame rate is needed for image acquisition.
Additionally, the region of interest is smaller than 1 mm. In order to achieve
a good image quality with the given light source, image acquisition is limited
to the given pixel resolution.

##  5 Conclusion

The breakup dynamics of a capillary bridge on a hydrophobic stripe, which
forms during evaporation of a droplet, was studied experimentally and
numerically. The droplet wets two hydrophilic stripes, separated by a
hydrophobic stripe. Different ratios of the hydrophobic and hydrophilic stripe
width $\alpha$ were considered. By performing experiments with two
synchronized high-speed cameras for $\alpha=$ 0.5, 1 and 1.5, the breakup
dynamics of the capillary bridge could be observed simultaneously to the
contact angle change on the hydrophilic stripe. Consequently, the contact
angle data could be used to calculate physically realistic initial conditions
for the continuum mechanical simulations using Surface Evolver.

The data import into the numerical simulation tool Free Surface 3D (FS3D) was
achieved using an algorithm within OpenFOAM that converts the surface mesh
from Surface Evolver into a volume fraction field. The geometric Volume-of-
Fluid method was employed to solve the three-dimensional two-phase Navier
Stokes equations. In order to dampen spurious currents at the contact line,
the Boundary Youngs interface reconstruction algorithm [24] was adapted to
three dimensions and a modification of the Navier Slip boundary condition was
introduced (“staggered slip”).

The phase space picture of the breakup dynamics, that allows for a detailed
and systematic study of the dynamics of the process, is employed. In
particular, the ambiguity related to the choice of the breakup time $t_{0}$ is
removed. Two distinct regimes were observed both in the experiments and the
numerical simulations which are connected by a transition region. The initial
regime follows approximately (in an average sense) the well-known inviscid
relation for a _free_ capillary bridge, i.e. $d(\tau)\propto
C(\sigma\tau^{2}/\rho)^{1/3}$, where $C$ is found to be $0.75-0.77$ in both
experiments and simulations. The regime transition is characterized by a
maximum in the breakup speed followed by a local minimum which corresponds to
the onset of the second dynamic regime. Within the second regime, a different
exponent $\nu$ can be determined ranging from $0.47$ to $0.48$ in the
investigated range of scales and parameters. The latter values differ
significantly from those obtained from dimensional analysis for a free
capillary bridge in the inviscid regime. This means that either the presence
of the moving contact line affects the exponent or that the final universal
regime is not yet reached. Nevertheless, it is shown that the complex dynamics
is governed by a balance of inertial and capillary forces. Viscous forces and
boundary slip show only a minor influence on the time evolution of the minimum
width. Experiments and simulations show a good quantitative agreement in terms
of the breakup speed down to a bridge width of approximately
$20\operatorname{{\mu}m}$.

In both experiments and simulations, a liquid thread (a cylinder cut by a
plane) is formed between two cones. These have a well-defined opening angle
that decreases with increasing $\alpha$. From the experiments, it could be
inferred that the final liquid thread breaks up in a Rayleigh-Plateau type
instability. Remarkably, the wavelength of the breakup belonging to the
fastest growing mode corresponds to the distance between the tips of the
cones.

##### Acknowledgments:

We kindly acknowledge the financial support by the German Research Foundation
(DFG) within the Collaborative Research Centre 1194 “Interaction of Transport
and Wetting Processes” – Project-ID 265191195, subprojects A02b and B01, B02
and Z-INF. Calculations for this research were conducted on the Lichtenberg
high performance computer of the TU Darmstadt.

Maximilian Hartmann and Mathis Fricke contributed equally to this work.

##  Appendix A On the Staggered Slip Boundary Condition

The staggered slip condition introduced in Section 3.3 reduces the amount of
artificial numerical slip and, thereby, is able to reduce spurious velocities
at the domain boundary and at the contact line. The method is conceptually
simple, can be implemented in a single line of code and is surprisingly
effective (see below). The purpose of this section is to give some more
details on the approach and its effect on dynamic wetting simulations. A more
comprehensive discussion can be found in Chapter 11 of [22].

The motivation for the adaptation of the Navier slip condition is the
observation that the contact line speed is usually over-predicted on coarse
meshes. There are (at least) two reasons for this numerical phenomenon: The
first reason is that the interface in a cell located at the boundary is
advected by a face-centered velocity located at a distance $\Delta x/2$ above
the physical boundary. The latter velocity is typically much larger than the
expected velocity right at the boundary leading to an over-prediction of the
contact line speed. The second reason is that the loss of energy due to
dissipation in the region close to the contact line is not properly resolved
on coarse meshes. This “missing dissipation” is mimicked by the increased
discrete viscous dissipation in a boundary cell caused by the increased
opposite velocity in the ghost cell resulting from the staggered slip
condition. Clearly, the proposed method is somewhat ad hoc and the results
could probably be further improved by a more sophisticated numerical modeling
of the missing viscous dissipation. Nevertheless, the method is able to
improve the convergence in a dynamic wetting simulation on coarse meshes when
the slip length is not resolved by the computational mesh (see below).

(a) $\omega=0$ (b) $\omega=1$

(c) $\omega=1/2$

(d) Convergence in the maximum norm.

FIG. 17: Capillary rise dynamics using the standard ($\omega=0$) and the
staggered ($\omega=1/2,1$) slip boundary condition compared to reference data
by Gründing et al. [26, 27] for $L=R/50$ and $\Omega=1$.

##### Capillary rise benchmark:

We examine the effect of the staggered slip boundary condition for the
numerical simulation of wetting flows based on the capillary rise problem as
described in detail by Gründing et al. [26, 27]. The latter study provides a
comprehensive numerical dataset on the dynamics of the two-dimensional
capillary rise in different regimes ranging from a highly oscillatory behavior
to a purely monotonic rise towards the stationary state. Following the work by
Fries and Dreyer [25] and Gründing et al. [26], the physical parameters of the
systems are varied such that different values of the non-dimensional parameter
(with $R$ the radius of the capillary)

$\Omega=\sqrt{\frac{9\sigma\cos\theta\eta^{2}}{\rho^{3}g^{2}R^{5}}}$

are realized. The numerical results for different values of $\omega$ and the
reference data [27] in terms of the rise height††The rise height is defined as
the height of the meniscus relative to the liquid bath. over time for the
setup and parameters specified by Gründing et al. [26] for a slip length
$L=R/50$ and $\Omega=1$ are shown in Fig. 17. The reference data are computed
on a fine mesh resolving the slip length and are confirmed by four different
numerical methods††The study [26] compares the present geometrical VOF method
FS3D with the algebraic VOF method interFoam, the Arbitrary-Lagrangian-
Eulerian method interTrackFoam and the level-set based extended discontinuous
Galerkin method BoSSS; see [26] for more details. including the present
implementation of FS3D with the standard Navier slip condition (15). For the
standard Navier slip condition ($\omega=0$) it is found that the reference
solution is approached from _above_ , i.e. the contact line velocity is over-
estimated on coarse grids. Conversely, the reference solution is approached
from _below_ for $\omega=1$, i.e. the contact line velocity is under-estimated
on coarse meshes. The error in terms of the maximum deviation of the rise
height over time is similar for $\omega=0$ and $\omega=1$; see Fig. 17(d). The
error is reduced by approximately one order of magnitude for the choice
$\omega=1/2$.

##### Single-phase channel flow:

Another instructive example for the validation of the boundary condition is a
single-phase flow in a two-dimensional channel driven by a prescribed pressure
gradient $G=-\partial_{x}p$. An analytical solution for the Navier Stokes
equations with the Navier slip boundary condition is available in this case
serving as a reference for validation. In particular, one can show that the
mass transport rate across the channel is proportional to $1+6L/H$, where $H$
is the height of the channel. The channel flow problem is studied with the
staggered slip implementation in [22, p. 156]. It is found that the order of
convergence for the mass transport rate drops from two to one when the
staggered slip condition is applied (unless $\omega$ is chosen to be very
small). Hence, the standard Navier slip condition is more accurate in the
single-phase case; see [22] for more details.

##### Effect on the break dynamics:

It is important to note that in the present study of the breakup dynamics of a
liquid bridge, the staggered slip condition only serves the purpose of damping
spurious currents at the contact line. The choice of $\omega$ shows only a
weak influence of the breakup dynamics in the simulation; see Fig. 18. The
dynamic exponents remain unchanged and a difference is visible only when the
local minimum in the velocity is approached. This behavior is consistent with
the observation that viscous effects show only a minor influence on the
breakup process on the considered length scales.

FIG. 18: Influence of the choice of $\omega$ in (17) on the breakup dynamics
($\alpha=1$).

In summary, the staggered slip condition is easy to implement and appears to
effectively model the missing viscous dissipation in the contact line region
which is not resolved by the grid. However, the optimal value for the free
parameter $\omega$ cannot be estimated a priori. Further research is necessary
to understand the underlying mechanisms at the moving contact line more
quantitatively.

##  Appendix B Initial Condition

Figure 19 shows the initial condition for the numerical simulations generated
with Surface Evolver compared to the corresponding state of three different
experiments for different $\alpha$. The experimental images shown are chosen
such that the minimum width $d$ of the capillary bridge matches the initial
condition. Obviously, the curvatures in the plane parallel to the substrate
deviate between the experiments and the deviation increases for increasing
$\alpha$. This is especially obvious in the quantitative comparison between
experiment and simulation for $\alpha=1.5$. In the experiments, the capillary
surface adopts its shape shortly before the final point of instability.
However, Surface Evolver is not able to calculate this transient process,
since it can only be used to calculate equilibrium shapes. This final time
span in the experiments, even though it is very short, might lead to these
small deviations. Nevertheless, employing Surface Evolver appears to be the
only method to compute the initial condition with reasonable effort.

FIG. 19: Initial condition for the numerical simulations compared to the
corresponding image from the experiments. In the first column, a scale bar is
depicted, which is valid for each respective row. In the column for the
simulations, the liquid is depicted in blue, in the experimental part in black
color. In both cases, the substrate is visible in gray color.

##  Appendix C Evaporation and Capillary Time Scales

We expect the breakup process to occur on the capillary time scale $T$, see
equation (1). Without going into detail, $T$ can be estimated as
$T=\sqrt{\frac{\rho}{\sigma}w_{\mathrm{phob}}^{3}}$, with $\rho$ being the
density of the liquid, $\sigma$ being the surface tension and
$w_{\mathrm{phob}}$ being the hydrophobic stripe width. The evaporation time
scale of a sessile droplet can be estimated as $T_{\mathrm{evap}}=\frac{\rho
LH}{D_{\mathrm{w,a}}\Delta c_{\mathrm{w}}}$, with $L$ and $H$ being two
different length scales, $D_{\mathrm{w,a}}=2.4\cdot
10^{-5}\leavevmode\nobreak\ \mathrm{m}^{2}/\mathrm{s}$ being the diffusion
coefficient of water vapor in air and $\Delta c_{\mathrm{w}}\approx 1\cdot
10^{-2}\leavevmode\nobreak\ \mathrm{kg}/\mathrm{m}^{3}$ being the
concentration difference of a saturated atmosphere to the surrounding air
[39]. For typical widths of the capillary bridge, which are in the order of
100 $\mu$m, $\frac{T_{\mathrm{evap}}}{T}\approx 10^{5}$ and therefore
evaporation can be neglected during breakup.

## References

  * [1] S. Afkhami and M. Bussmann.  Height functions for applying contact angles to 2D VOF simulations.  International Journal for Numerical Methods in Fluids, 57(4):453–472, 2008.  doi:10.1002/fld.1651.
  * [2] S. Afkhami and M. Bussmann.  Height functions for applying contact angles to 3D VOF simulations.  International Journal for Numerical Methods in Fluids, 61(8):827–847, 2009.  doi:10.1002/fld.1974.
  * [3] W. Barthlott and C. Neinhuis.  Purity of the sacred lotus, or escape from contamination in biological surfaces.  Planta, 202(1):1–8, 1997.  doi:10.1007/s004250050096.
  * [4] O. Bliznyuk, E. Vereshchagina, E. S. Kooij, and B. Poelsema.  Scaling of anisotropic droplet shapes on chemically stripe-patterned surfaces.  Physical Review E, 79(4):041601, 2009.  doi:10.1103/PhysRevE.79.041601.
  * [5] J. B. Bostwick and P. H. Steen.  Stability of Constrained Capillary Surfaces.  Annual Review of Fluid Mechanics, 47(1):539–568, 2015.  doi:10.1146/annurev-fluid-010814-013626.
  * [6] J. B. Bostwick and P. H. Steen.  Static rivulet instabilities: Varicose and sinuous modes.  Journal of Fluid Mechanics, 837:819–838, 2018.  doi:10.1017/jfm.2017.876.
  * [7] J. U. Brackbill, D. B. Kothe, and C. Zemach.  A continuum method for modeling surface tension.  Journal of Computational Physics, 100(2):335–354, 1992.  doi:10.1016/0021-9991(92)90240-Y.
  * [8] K. A. Brakke.  The Surface Evolver.  Experimental Mathematics, 1(2):141–165, 1992.  doi:10.1080/10586458.1992.10504253.
  * [9] M. Brinkmann, J. Kierfeld, and R. Lipowsky.  Stability of liquid channels or filaments in the presence of line tension.  Journal of Physics: Condensed Matter, 17(15):2349–2364, 2005.  doi:10.1088/0953-8984/17/15/008.
  * [10] J. C. Burton, J. E. Rutledge, and P. Taborek.  Fluid pinch-off in superfluid and normal ${}^{4}\mathrm{He}$.  Phys. Rev. E, 75:036311, 2007.  doi:10.1103/PhysRevE.75.036311.
  * [11] Y.-J. Chen and P. H. Steen.  Dynamics of inviscid capillary breakup: collapse and pinchoff of a film bridge.  Journal of Fluid Mechanics, 341(1997):245–267, 1997.  doi:10.1017/S002211209700548X.
  * [12] R. David and A. W. Neumann.  Anisotropic drop shapes on chemically striped surfaces.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 393:32–36, 2012.  doi:10.1016/j.colsurfa.2011.10.020.
  * [13] S. H. Davis.  Moving contact lines and rivulet instabilities. Part 1. The static rivulet.  Journal of Fluid Mechanics, 98(02):225, 1980.  doi:10.1017/S0022112080000110.
  * [14] A. Deblais, M. A. Herrada, I. Hauner, K. P. Velikov, T. van Roon, H. Kellay, J. Eggers, and D. Bonn.  Viscous Effects on Inertial Drop Formation.  Physical Review Letters, 121(25):254501, 2018.  doi:10.1103/PhysRevLett.121.254501.
  * [15] P. G. Drazin.  Introduction to Hydrodynamic Stability.  Cambridge University Press, sep 2002.  doi:10.1017/CBO9780511809064.
  * [16] A. Dziedzic, M. Nakrani, B. Ezra, M. Syed, S. Popinet, and S. Afkhami.  Breakup of finite-size liquid filaments: Transition from no-breakup to breakup including substrate effects.  The European Physical Journal E, 42(2):18, 2019.  doi:10.1140/epje/i2019-11785-y.
  * [17] J. Eggers.  Nonlinear dynamics and breakup of free-surface flows.  Rev. Mod. Phys., 69:865–930, Jul 1997.  doi:10.1103/RevModPhys.69.865.
  * [18] J. Eggers and E. Villermaux.  Physics of liquid jets.  Reports on Progress in Physics, 71(3), 2008.  doi:10.1088/0034-4885/71/3/036601.
  * [19] J. G. Eggers and M. A. Fontelos.  Singularities : Formation, Structure, and Propagation.  Cambridge, 2015.  doi:10.1017/CBO9781316161692.
  * [20] A. Fath and D. Bothe.  Direct numerical simulations of thermocapillary migration of a droplet attached to a solid wall.  International Journal of Multiphase Flow, 77:209–221, 2015.  doi:10.1016/j.ijmultiphaseflow.2015.08.018.
  * [21] C. Focke and D. Bothe.  Direct numerical simulation of binary off-center collisions of shear thinning droplets at high Weber numbers.  Physics of Fluids, 24(7):073105, 2012.  doi:10.1063/1.4737582.
  * [22] M. Fricke.  Mathematical modeling and Volume-of-Fluid based simulation of dynamic wetting.  PhD thesis, Darmstadt, 2021.  doi:10.12921/tuprints-00014274.
  * [23] M. Fricke, M. Köhne, and D. Bothe.  A kinematic evolution equation for the dynamic contact angle and some consequences.  Physica D: Nonlinear Phenomena, 394:26–43, 2019.  doi:10.1016/j.physd.2019.01.008.
  * [24] M. Fricke, T. Marić, and D. Bothe.  Contact Line Advection using the geometrical Volume-of-Fluid Method.  Journal of Computational Physics, 407:109221, 2020.  doi:10.1016/j.jcp.2019.109221.
  * [25] N. Fries and M. Dreyer.  Dimensionless scaling methods for capillary rise.  Journal of Colloid and Interface Science, 338(2):514–518, Oct 2009\.  doi:10.1016/j.jcis.2009.06.036.
  * [26] D. Gründing, M. Smuda, T. Antritter, M. Fricke, D. Rettenmaier, F. Kummer, P. Stephan, H. Marschall, and D. Bothe.  A comparative study of transient capillary rise using direct numerical simulations.  Applied Mathematical Modelling, 86:142–165, 2020.  doi:10.1016/j.apm.2020.04.020.
  * [27] D. Gründing, M. Smuda, T. Antritter, M. Fricke, D. Rettenmaier, F. Kummer, P. Stephan, H. Marschall, and D. Bothe.  A comparative study of transient capillary rise using direct numerical simulations: Benchmark Data, 2020.  doi:10.25534/tudatalib-173.
  * [28] F. H. Harlow and J. E. Welch.  Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface.  Physics of Fluids, 8(12):2182, 1965.  doi:10.1063/1.1761178.
  * [29] M. Hartmann, M. Fricke, L. Weimar, D. Gründing, T. Marić, D. Bothe, and S. Hardt.  Breakup Dynamics of Capillary Bridges on Hydrophobic Stripes: Research data, 2021.  doi:10.48328/tudatalib-425.
  * [30] M. Hartmann and S. Hardt.  Stability of Evaporating Droplets on Chemically Patterned Surfaces.  Langmuir, 35(14):4868–4875, 2019.  doi:10.1021/acs.langmuir.9b00172.
  * [31] I. M. Hauner, A. Deblais, J. K. Beattie, H. Kellay, and D. Bonn.  The Dynamic Surface Tension of Water.  The Journal of Physical Chemistry Letters, 8(7):1599–1603, 2017\.  doi:10.1021/acs.jpclett.7b00267.
  * [32] L. He, X. Sui, W. Liang, Z. W., and A. Akbarzadeh.  Numerical analysis of anisotropic wetting of chemically striped surfaces.  RSC Advances, 8(55):31735–31744, 2018.  doi:10.1039/C8RA06626D.
  * [33] C. Huh and L. E. Scriven.  Hydrodynamic model of steady movement of a solid/liquid/fluid contact line.  Journal of Colloid and Interface Science, 35(1):85–101, 1971.  doi:10.1016/0021-9797(71)90188-3.
  * [34] H. P. Jansen, K. Sotthewes, J. van Swigchem, H. J. W. Zandvliet, and E. S. Kooij.  Lattice Boltzmann modeling of directional wetting: Comparing simulations to experiments.  Physical Review E, 88(1):013008, 2013.  doi:10.1103/PhysRevE.88.013008.
  * [35] H. P. Jansen, K. Sotthewes, H. J. W. Zandvliet, and E. S. Kooij.  Potential of lattice Boltzmann to model droplets on chemically stripe-patterned substrates.  Applied Surface Science, 361:122–132, 2016.  doi:10.1016/j.apsusc.2015.11.120.
  * [36] J. B. Keller and M. J. Miksis.  Surface Tension Driven Flows.  SIAM Journal on Applied Mathematics, 43(2):268–277, 1983.  doi:10.1137/0143018.
  * [37] T. A. Kowalewski.  On the separation of droplets from a liquid jet.  Fluid Dynamics Research, 17(3):121–145, feb 1996.  doi:10.1016/0169-5983(95)00028-3.
  * [38] J. Léopoldès, A. Dupuis, D. G. Bucknall, and J. M. Yeomans.  Jetting Micron-Scale Droplets onto Chemically Heterogeneous Surfaces.  Langmuir, 19(23):9818–9822, 2003.  doi:10.1021/la0353069.
  * [39] Y. Li, P. Lv, C. Diddens, H. Tan, H. Wijshoff, and D. Versluis, M.and Lohse.  Evaporation-Triggered Segregation of Sessile Binary Droplets.  Physical Review Letters, 120(22):224501, jun 2018.  doi:10.1103/PhysRevLett.120.224501.
  * [40] Y. Li and J. E. Sprittles.  Capillary breakup of a liquid bridge: identifying regimes and transitions.  Journal of Fluid Mechanics, 797:29–59, 2016.  doi:10.1017/jfm.2016.276.
  * [41] A. Lippert.  Direct numerical simulations of thermocapillary driven motions in two-phase flows.  PhD thesis, Darmstadt, 2016.  URL: https://tuprints.ulb.tu-darmstadt.de/5817/.
  * [42] H. Matsui, Y. Noda, and T. Hasegawa.  Hybrid Energy-Minimization Simulation of Equilibrium Droplet Shapes on Hydrophilic/Hydrophobic Patterned Surfaces.  Langmuir, 28(44):15450–15453, 2012.  doi:10.1021/la303717n.
  * [43] C. Neto, D. R. Evans, E. Bonaccurso, H.-J. Butt, and V. S. J. Craig.  Boundary slip in newtonian liquids: A review of experimental studies.  Reports on Progress in Physics, 68(12):2859–2897, 2005.  doi:10.1088/0034-4885/68/12/R05.
  * [44] J. K. Park, J. Ryu, B. C. Koo, S. Lee, and K. H. Kang.  How the change of contact angle occurs for an evaporating droplet: effect of impurity and attached water films.  Soft Matter, 8(47):11889, 2012.  doi:10.1039/c2sm26559a.
  * [45] A. R. Parker and C. R. Lawrence.  Water capture by a desert beetle.  Nature, 414(6859):33–34, 2001.  doi:10.1038/35102108.
  * [46] S. Popinet.  An accurate adaptive solver for surface-tension-driven interfacial flows.  Journal of Computational Physics, 228(16):5838–5866, 2009.  doi:10.1016/j.jcp.2009.04.042.
  * [47] M. Renardy, Y. Renardy, and J. Li.  Numerical Simulation of Moving Contact Line Problems Using a Volume-of-Fluid Method.  Journal of Computational Physics, 171(1):243–263, 2001.  doi:10.1006/jcph.2001.6785.
  * [48] W. J. Rider and D. B. Kothe.  Reconstructing Volume Tracking.  Journal of Computational Physics, 141(2):112–152, 1998.  doi:10.1006/jcph.1998.5906.
  * [49] M. Rieber.  Numerische Modellierung der Dynamik freier Grenzflächen in Zweiphasenströmungen, volume 459 of Fortschritt-Berichte VDI: Reihe 7, Strömungstechnik.  VDI-Verlag, Düsseldorf, 2004.
  * [50] M. Rieber and A. Frohn.  A numerical study on the mechanism of splashing.  International Journal of Heat and Fluid Flow, 20(5):455–461, oct 1999.  doi:10.1016/S0142-727X(99)00033-8.
  * [51] D. Song, B. Song, H. Hu, X. Du, and F. Zhou.  Selectively splitting a droplet using superhydrophobic stripes on hydrophilic surfaces.  Physical Chemistry Chemical Physics, 17(21):13800–13803, 2015.  doi:10.1039/C5CP01530H.
  * [52] R. L. Speth and E. Lauga.  Capillary instability on a hydrophilic stripe.  New Journal of Physics, 11(7):075024, 2009.  doi:10.1088/1367-2630/11/7/075024.
  * [53] J. E. Sprittles and Y. D. Shikhmurzaev.  Finite element framework for describing dynamic wetting phenomena.  International Journal for Numerical Methods in Fluids, 68(10):1257–1298, 2012.  doi:10.1002/fld.2603.
  * [54] G. Strang.  On the Construction and Comparison of Difference Schemes.  SIAM Journal on Numerical Analysis, 5(3):506–517, 1968.  doi:10.1137/0705041.
  * [55] T. Tolle, D. Gründing, D. Bothe, and T. Marić.  Computing volume fractions and signed distances from triangulated surfaces immersed in unstructured meshes, 2021.  arXiv:2101.08511.
  * [56] G. Tryggvason, B. Bunner, O. Ebrat, and W. Tauber.  Computations of Multiphase Flows by a Finite Difference/Front Tracking Method. I. Multi-Fluid Flows.  Lecture Series - von Karman Institute For Fluid Dynamics, 1998.
  * [57] G. Tryggvason, R. Scardovelli, and S. Zaleski.  Direct Numerical Simulations of Gas–Liquid Multiphase Flows.  Cambridge University Press, 2011.  doi:10.1017/CBO9780511975264.
  * [58] J. Z. Wang, Z. H. Zheng, H. W. Li, W. T. S. Huck, and H. Sirringhaus.  Dewetting of conducting polymer inkjet droplets on patterned surfaces.  Nature Materials, 3(3):171–176, 2004.  doi:10.1038/nmat1073.
  * [59] Y. Wu, F. Wang, M. Selzer, and B. Nestler.  Investigation of Equilibrium Droplet Shapes on Chemically Striped Patterned Surfaces Using Phase-Field Method.  Langmuir, 35:8500–8516, 2019.  doi:10.1021/acs.langmuir.9b01362.
  * [60] D. L. Youngs.  An Interface Tracking Method for a 3D Eulerian Hydrodynamics Code.  Technical Report AWRE/44/92/35, 1984.
  * [61] L. Zhang, J. Wu, M. N. Hedhili, X. Yang, and P. Wang.  Inkjet printing for direct micropatterning of a superhydrophobic surface: toward biomimetic fog harvesting surfaces.  Journal of Materials Chemistry A, 3(6):2844–2852, 2015.  doi:10.1039/C4TA05862C.
  * [62] L. Zou, H. Wang, X. Zhu, Y. Ding, R. Chen, and Q. Liao.  Droplet splitting on chemically striped surface.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 537(August 2017):139–148, 2018.  doi:10.1016/j.colsurfa.2017.09.052.