# The effect of contact angle hysteresis on a droplet in a viscoelastic two-
phase system

Kazem Bazesefidpar Corresponding author<EMAIL_ADDRESS>Outi Tammisola  SeRC
(Swedish e-Science Research Centre) and FLOW, Dept. of Engineering Mechanics,
KTH Royal Institute of Technology, SE-10044 Stockholm, Sweden

## abstract

We investigate the dynamic behaviour of a two-dimensional (2D) droplet
adhering to a wall in Poiseuille flow at low Reynolds numbers, in a system
where either the droplet is viscoelastic (V/N) or the surrounding medium
(N/V). The fluid viscoelasticity has been modeled by the Giesekus constitutive
equation, and the Cahn–Hilliard Phase-Field method is used to capture the
interface between the two phases. The contact angle hysteresis is represented
by an advancing contact angle $\theta_{A}$ and a receding contact angle
$\theta_{R}$. The results reveal that the deformation of the viscoelastic drop
over time is changed due to the presence of polymeric molecules, and it can be
categorized in two stages prior to depinning of the contact lines. In the
first stage, the viscoelastic droplet speeds up and deforms faster, while in
the second stage, the Newtonian counterpart accelerates and its deformation
outpaces the viscoelastic droplet. The deformation of viscoelastic drop is
retarded significantly in the second stage with increasing Deborah number
$De$. In the V/N case, the viscous bending is enhanced on the receding side
for small $De$, but it is weakened by further increase in $De$, and this non-
monotonic behavior brings about an increase in the receding contact line
velocity at small $De$ and a decrease at large $De$. On the advancing side,
the viscous bending is decreased monotonically for $Ca<0.25$, and hence the
advancing contact line velocity is decreased with increasing $De$; however,
the viscous bending presents a non-monotonic behavior for $Ca=0.25$ at the
advancing side. The non-monotonic behavior on the receding side is attributed
to the emergence of outward pulling stresses in the vicinity of the receding
contact line and the inception of strain-hardening at higher $De$, while the
reduction in the viscous bending at the advancing side is the result of just
strain-hardening due to the presence of dominant extensional flow on the
advancing side. Finally, in the N/V system, the viscoelasticity of the medium
suppresses the droplet deformation on both receding and advancing sides, and
this effect is more pronounced with increasing De; the weakening effect of
viscous bending is enhanced significantly at the advancing side by increasing
the Giesekus mobility parameter. These results give a thorough understanding
of viscoelastic effect on both drop deformation and depinning at both contact
lines over a surface with contact angle hysteresis.

Keywords: Wetting, viscoelasticity, contact angle hysteresis, superhydrophobic
surface, Diffuse-interface method

##  1 Introduction

The wetting of a solid substrate by a polymeric drop is an important process
in many applications such as inkjet printing, coating, and self-cleaning
surfaces. Among these applications, superhydrophobic surfaces have recently
drawn the attention of many researchers due to their application in drag
reduction and self-cleaning [Cheng and Rodak, 2005, Neinhuis and Barthlott,
1997]. These surfaces are inspired by natural surfaces, such as gecko skin
[Watson et al., 2015] and lotus leaves that give rise to the so-called lotus
effect [Shirtcliffe et al., 2010]. Superhydrophobic surfaces are usually
characterized (on a macroscopic length scale) by the advancing contact angle
$\theta_{A}$ and the receding contact angle $\theta_{R}$, and the difference
between these two contact angles called Contact Angle Hysteresis (CAH). The
properties of a surface such as microscopic roughness and inhomogeneities may
lead to CAH [Joanny and Robbins, 1990], and the movement of a droplet on the
surface can be significantly affected by CAH. Furthermore, the non-Newtonian
properties of a droplet affect its motion on a solid surface. For example, the
elasticity of a droplet moderately speeds up the contact-line motion on a
hydrophilic flat and chemically homogeneous surface at late stages of wetting
and for a high-viscosity liquid [Yue and Feng, 2012, Wang et al., 2015, 2017],
but a recent study [Yada et al., 2023] revealed that increasing the polymer
concentration does not modify the contact-line friction for aqueous polymer
(shear-thinning) solutions on the smooth surface in the initial wetting
regime. Thus, the elasticity of a drop has a negligible effect on the early
stages of the wetting process, and the contact line motion at early times can
be controlled by changing the viscosity of the solvent or adding
roughness/inhomogeneities to the surface.

In contrast to the smooth hydrophilic and hydrophobic surfaces, the elasticity
of a droplet plays an important role on the superhydrophobic surfaces.
Experiment results suggest that the superhydrophobic surface may lose its
self-cleaning property when they interact with highly elastic fluids [Xu et
al., 2018]. The motion of a Boger fluid drop, which shows viscoelasticity
without shear thinning, on the superhydrophobic surface is strikingly
different from a Newtonian drop. The velocity of the viscoelastic drop is
notably reduced in comparison with the Newtonian drop, and complex branch-like
patterns are left on the tilted superhydrophobic surface. Previous studies
also showed that the polyacrylamide drops at higher polymer concentrations
were noticeably more elongated when sliding down a hydrophilic surface with
small CAH [Varagnolo et al., 2017]. The wetting area of a viscoelastic drop in
a Newtonian surrounding fluid has been studied by numerical simulations using
the Oldroyd-B model [Wang et al., 2023]. They found that the wetting area
decreased when the elasticity of the fluid increased, and the same was
observed when the polymeric viscosity ratio decreased. Also droplet motion on
lubricant surfaces has been analysed. A lubricant surface is constructed by
trapping a suitable low surface tension lubricating liquid inside a texture
[Lafuma and Quéré, 2011], and these surfaces hardly pin drops. An oscillatory
motion has been reported for the viscoelastic drop having both shear-thinning
and elasticity on lubricant surfaces [Sartori et al., 2022].

Our aim is to model the effect of CAH on viscoelastic fluids on a continuum
level, neglecting the microscopic details of the roughness, using a model that
can handle a dynamic contact angle, but still allows receding and advancing
contact angles to have different values. The most popular approach has been
developed by Spelt [2005] for a level-set method. An intermediate contact
angle is computed, and the contact line pins if it is within the hysteresis
window $[\theta_{R},\theta_{A}]$, otherwise it moves. This approach has then
been extended to the phase-field method [Ding and Spelt, 2008], volume-of-
fluid method [Dupont and Legendre, 2010], lattice Boltzmann method [Ba et al.,
2013], and front-tracking method [Shin et al., 2018]. It should be noted that
the above methods need ghost cells outside the boundary to impose the contact
angle condition and a slip model to circumvent the well-known singularity in
the contact line [Batchelor et al., 2003]; however, the Cahn–Hilliard equation
automatically regularizes the singularity at the contact line [Jacqmin, 2000].
A CAH model capable of capturing pinning, advancing and receding automatically
without needing any explicit knowledge of contact line velocity was proposed
by Yue [2020], and it has been further extended to level-set method [Zhang and
Yue, 2020]. Recently, an angle-dependent line friction model for the
Cahn–Hilliard equation has been proposed by Amberg [2022], and the numerical
results of the proposed model were in a very good agreement with the
experimental data.

In the present work, we investigate the effect of a fluid’s viscoelastic
properties on the deformation and movement of a droplet with contact line
hysteresis in a 2D Poiseuille flow at low Reynolds number ($Re$) by direct
numerical simulation. The deformation and yield condition (depinning of
contact lines) of both 2D and three-dimensional (3D) Newtonian drops were the
subject of many previous studies at low and high $Re$, which are summarized
below.

The yield condition of a 2D droplet in a shear flow and stokes regime was
examined for a broad range of capillary numbers (Ca), viscosity ratios
($\lambda_{\mu}$), and CAH ($\Delta\theta=\theta_{A}-\theta_{R}$)
[Dimitrakopoulos and Higdon, 1997]. It is found that increasing $\Delta\theta$
results in an increase in the critical capillary number ($Ca_{cr}$), above
which the droplet starts to move. Their study also demonstrated that a droplet
with a higher viscosity (higher $\lambda_{\mu}$) on a surface with a specific
CAH starts to move at a smaller $Ca_{cr}$ than a droplet with a lower
viscosity.

The effect of $Ca$, $\lambda_{\mu}$, slip velocity, and surfactants on the 2D
droplet deformation and yield condition in Stokes flow regime was studied by
Schleizer and Bonnecaze [1999] in both shear-driven and pressure-driven flows.
They demonstrated that surfactants reduce the deformation of a drop due to the
generation of marangoni stresses along the droplet surface; in addition, a
decrease in the deformation of the drop when the surface is hydrophilic and an
increase in the deformation when the surface is hydrophobic have been observed
in the presence of a slip velocity without considering CAH.

Later, 2D droplet deformation in shear flow with pinned and moving contact
lines was studied at moderate $Re$ using a level-set method by Spelt [2006].
The author showed that the critical Weber number ($We_{cr}$) approaches a
constant value when $Re$ increases. The droplet becomes unstable for
$We>We_{cr}$, and the transient behaviour is very different at different
Reynolds numbers. It was demonstrated that the slip at the wall increases the
$We_{cr}$, and the slip length has a significant effect on the contact line
speed. It should be mentioned that the behavior of the droplet with slip
velocity was reported to be dependent on the slip length.

Afterward, the critical conditions for the motion of a 3D droplet at both low
and moderate $Re$ in the presence of CAH were studied [Ding and Spelt, 2008].
The inertial effects led to a complicated structure around the deformed drop
with a strong 3D effect for the moderate Re, whereas their results were in
agreement with previous 2D studies at low Re.

In the present study, we conduct direct numerical simulations to study the
effect of fluid viscoelasticity on the deformation and depinning of a 2D drop
either in V/N or N/V systems over a surface in the presence of CAH at low
$Re$. In particular, we focus on the physical effects of polymers on the
deformation and depinning processes. We apply the Cahn–Hilliard equation with
a dynamic contact angle boundary condition for capturing the interface between
the two phases, and the Giesekus constitutive equation to model the
viscoelasticity of the drop. First, the qualitative effect of the elasticity
is investigated over a wide range of the $Ca$ at a low constant $De=1$. We
compare the time evolution of both dynamic contact angle and velocity of the
receding and advancing contact lines for a viscoelastic droplet (V/N) and a
Newtonian droplet in a Newtonian medium over surfaces with a varying but
significant degree of contact line hysteresis
($\Delta\theta=[100^{\circ},130^{\circ},140^{\circ}]$). Next, the effect of
$De$ is studied in the V/N system. Finally, the effect of $De$ and $\alpha$
has been investigated on the drop’s deformation and the depping of contact
lines in the N/V system.

##  2 Equations

###  2.1 Governing equations

Consider an incompressible two-phase system composed of a Newtonian fluid with
viscosity $\mu_{n}=\mu_{s}$ and a viscoelastic (Giesekus) fluid with solvent
viscosity $\mu_{s}$, and polymeric viscosity $\mu_{p}$. A phase-field variable
$\phi$ varies smoothly from phase one ($\phi=1$) to another phase ($\phi=-1$)
with $\phi=0$ denoting the fluid/fluid interface. This two-phase system is
modeled by the following coupled system of equations [Yue et al., 2004, Abels
et al., 2012]:
the linear momentum and continuity equations

$\displaystyle\rho(\frac{\partial{\mathbf{u}}}{\partial{t}}+({\mathbf{u}}\cdot\nabla){\mathbf{u}})+{\mathbf{J}}\cdot\nabla{\mathbf{u}}$
$\displaystyle=$
$\displaystyle-\nabla{p}+\nabla\cdot\mu(\nabla{{\mathbf{u}}}+\nabla{{\mathbf{u}}^{T}})+G\nabla\phi$
$\displaystyle+\nabla\cdot({\frac{(1\pm\phi)}{2}\boldsymbol{\tau}_{p}}),$ (1)
$\displaystyle\nabla\cdot{\mathbf{u}}=0,$ (2)

the Giesekus constitutive model for the viscoelastic fluid

$\displaystyle\boldsymbol{\tau}_{p}+\lambda_{H}(\frac{\partial{\boldsymbol{\tau}_{p}}}{\partial{t}}+{\mathbf{u}}\cdot\nabla{\mathbf{\boldsymbol{\tau}_{p}}}-\boldsymbol{\tau}_{p}\nabla{{\mathbf{u}}}-\nabla{{\mathbf{u}}^{T}}\boldsymbol{\tau}_{p})+\frac{\alpha\lambda_{H}}{\mu_{p}}(\boldsymbol{\tau}_{p}\cdot\boldsymbol{\tau}_{p})$
$\displaystyle=$
$\displaystyle\mu_{p}(\nabla{{\mathbf{u}}}+\nabla{{\mathbf{u}}^{T}}),\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$
(3)

and the Cahn-Hilliard equation for the movement of the interface

$\displaystyle\frac{\partial{\phi}}{\partial{t}}+\nabla\cdot({{\mathbf{u}}\phi})=M\nabla^{2}G,$
(4) $\displaystyle G=-\lambda\nabla^{2}{\phi}+h(\phi).$ (5)

In the above equations, $\mathbf{u}(\mathbf{x},t)$ is the velocity vector,
$p(\mathbf{x},t)$ is the pressure, and $\boldsymbol{\tau}_{p}$ is the extra
stress due to the polymers. In the Giesekus equation, $\lambda_{H}$ is the
polymer relaxation time, and $\alpha$ is the Giesekus mobility parameter. In
the Cahn-Hilliard equation, $G$ is the chemical potential, $M$ is the mobility
parameter, $\lambda$ is the mixing energy density which is related to the
surface tension in the sharp-interface limit [Yue et al., 2004] as

$\displaystyle\sigma=\frac{2\sqrt{2}}{3}\frac{\lambda}{\eta}.$ (6)

The function $h(\phi)$ in Eq. (5) is defined as

$\displaystyle
h(\phi)=F^{\prime}(\phi)=\frac{\lambda}{\eta^{2}}\phi(\phi^{2}-1),$ (7)

with $F(\phi)=\frac{\lambda}{4\eta^{2}}(\phi^{2}-1)^{2}$ being the double-well
potential, and where $\eta$ is the capillary width of the interface. The
density $\rho$ and the dynamic viscosity $\mu$ of the mixture are a linear
function of the phase-field variable as:

$\displaystyle\rho=\frac{(1+\phi)}{2}{\rho_{1}}+\frac{(1-\phi)}{2}{\rho_{2}},$
(8)
$\displaystyle\mu=\frac{(1+\phi)}{2}{\mu_{s1}}+\frac{(1-\phi)}{2}{\mu_{s2}}.$
(9)

The total viscosity of the non-Newtonian phase is $\mu=\mu_{s}+\mu_{p}$. The
flux of density due to the diffusion of components $\mathbf{J}$ [Abels et al.,
2012] in Eq. (1) is given by

$\displaystyle\mathbf{J}=-\frac{(\rho_{1}-\rho_{2})}{2}M\mathbf{\nabla}{G}.$
(10)

The equations (1), (3), and (4) should be supplemented with appropriate
boundary and initial conditions. The following boundary conditions are
imposed:

$\displaystyle\boldsymbol{u}=\mathbf{0},$ (11)
$\displaystyle{\mathbf{n}}\cdot\nabla{G}=0,$ (12)
$\displaystyle-\eta\mu_{f}(\frac{\partial{\phi}}{\partial{t}}+\mathbf{u}\cdot\nabla\phi)=\frac{3}{2\sqrt{2}}\sigma\eta\mathbf{n}\cdot\nabla{\phi}+f^{\prime}_{w}(\phi)\hskip
28.45274pton\hskip 14.22636pt\partial\Omega_{w},$ (13)
$\displaystyle\mathbf{n}\cdot\nabla{\phi}=0\hskip 28.45274pton\hskip
14.22636pt\partial\Omega\backslash\partial\Omega_{w}$ (14)

Equations (11) and (12), respectively, impose the no-slip and zero diffusive
flux at the wall, where $\mathbf{n}$ is the outward pointing normal vector to
the boundary. Equations (13) and (14) correspond to the wall energy relaxation
[Jacqmin, 2000, Carlson et al., 2009] and neutral wettability condition,
respectively. $\partial\Omega\backslash\partial\Omega_{w}$ denotes the
boundaries other than the solid walls intersected with contact line, and
$\mu_{f}$ is the contact line friction. The surface chemical potential
$L({\phi,\nabla\phi})$ is defined as

$\displaystyle
L({\phi,\nabla\phi})=\lambda\mathbf{n}\cdot\nabla{\phi}+f^{\prime}_{w}(\phi),$
(15)

where $f_{w}=-\sigma\frac{\phi(3-\phi^{2})}{4}cos\theta$ is the wall energy.

##  3 Numerical method

The governing equations are discretized by a second-order finite difference
method on a dual-resolution grid. The numerical method used in the present
study has been described in detail and validated in [Bazesefidpar et al.,
2022], except for the Cahn-Hilliard equation which is described in the
following section.

###  3.1 Cahn-Hilliard equation

Different approaches have been proposed to deal with the non-linear term in
the Cahn-Hilliard equation (Eq. 4) embedded in the chemical potential $G$ [Gao
and Wang, 2014, Shen and Yang, 2010, Guillén-González and Tierra, 2013, Shen
et al., 2018], but scalar auxiliary variable (SAV) proposed by [Shen et al.,
2018] has the advantages that the resulting discretized equations are linear,
unconditionally energy stable, and relatively straight-forward to implement.
Challenges still remain when using this method, and the accuracy of SAV scheme
declines if too large time steps are chosen so that the results are generally
no longer accurate. One strategy to overcome this issue is to use adaptive
time stepping strategy [Shen et al., 2019] or the relaxed-SAV (RSAV) method
[Jiang et al., 2022, Zhang and Shen, 2022], and the RSAV method proposed by
Jiang et al. [2022] has been used in the present study. The Equations (4-5)
with the boundary conditions (12-14) can be reformulated and discretized in
time by using RSAV-BDF2 scheme [Jiang et al., 2022] as

$\displaystyle\frac{\gamma_{0}\phi^{n+1}-\hat{\phi}}{\Delta
t}+\nabla\cdot({{\mathbf{u^{*,n+1}}}\phi^{*,n+1}})=M\nabla^{2}G^{n+1},$ (16)
$\displaystyle
G^{n+1}=-\lambda\nabla^{2}{\phi^{n+1}}+\frac{\tilde{r}^{n+1}}{\sqrt{E(\phi^{*,n+1})}}h(\phi^{*,n+1}),$
(17) $\displaystyle\frac{\gamma_{0}\tilde{r}^{n+1}-\hat{r}}{\Delta
t}=\int_{\Omega}\frac{h(\phi^{*,n+1})}{2\sqrt{E(\phi^{*,n+1})}}\frac{\gamma_{0}\phi^{n+1}-\hat{\phi}}{\Delta
t},$ (18) $\displaystyle{\mathbf{n}}\cdot\nabla{G}^{n+1}=0,$ (19)
$\displaystyle\eta\mu_{f}(\frac{\gamma_{0}\phi^{n+1}-\hat{\phi}}{\Delta
t})=-L^{n+1}_{A,R}\hskip 28.45274pton\hskip 14.22636pt\partial\Omega_{w},$
(20) $\displaystyle\mathbf{n}\cdot\nabla{\phi^{n+1}}=0\hskip
28.45274pton\hskip 14.22636pt\partial\Omega\backslash\partial\Omega_{w},$ (21)
$\displaystyle
r^{n+1}=\zeta_{0}\tilde{r}^{n+1}+(1-\zeta_{0})\sqrt{E(\phi^{n+1})},$ (22)

where $E(t)=\int_{\Omega}F(\phi)$, $g^{*,n+1}=2g^{n}-g^{n-1}$ represents the
second-order explicit approximation of $g^{n+1}$,
$\hat{g}=2g^{n}-\frac{1}{2}g^{n-1}$, and $\gamma_{0}=\frac{3}{2}$. The
$L^{n+1}_{A,R}$ on the right-hand side of Eq. 20 is defined as [Yue, 2020]

$L^{n+1}_{A,R}=\begin{cases}L_{A}({\phi^{*,n+1},\nabla\phi^{n+1})}\qquad
if\qquad 0>L^{*,n+1}_{A}>L^{*,n+1}_{R},\\\
L_{R}({\phi^{*,n+1},\nabla\phi^{n+1})}\qquad if\qquad
L^{*,n+1}_{A}>L^{*,n+1}_{R}>0,\\\ 0\qquad otherwise,\end{cases}$ (23)

where the explicit second order approximation $\phi^{*,n+1}$ of the phase
field function has been used in computing $L_{A}$ and $L_{R}$ with A and R
corresponding to advancing and receding contact angles respectively. If the
contact angle is between the hysteresis window $[\theta_{R},\theta_{A}]$, the
contact line is pinned; otherwise it moves. The speed of the receding contact
line (when it is not pinned) depends on the $\theta_{d,R}$, and the advancing
contact line speed on the $\theta_{d,A}$, where $\theta_{d}$ is the dynamic
contact angle defined on the wall ($y=0$). The above system of equations
(16-21) can be transformed into two coupled equations [Yang et al., 2019] by
introducing new field functions as

$\displaystyle\phi^{n+1}=\phi_{1}^{n+1}+z^{n+1}\phi_{2}^{n+1},$ (24)
$\displaystyle G^{n+1}=G_{1}^{n+1}+z^{n+1}G_{2}^{n+1},$ (25)

resulting in the following final set of equations.

The first coupled system of equations for $\phi_{1}$ and $G_{1}$ is

$\displaystyle\nabla^{2}G_{1}^{n+1}-\frac{\gamma_{0}}{M\Delta
t}\phi_{1}^{n+1}=\frac{-\hat{\phi}}{M\Delta
t}+\frac{1}{M}\nabla\cdot({{\mathbf{u^{*,n+1}}}\phi^{*,n+1}}),$ (26)
$\displaystyle\nabla^{2}\phi_{1}^{n+1}+\frac{1}{\lambda}G_{1}^{n+1}=\frac{1}{\lambda\gamma_{0}}(\hat{r}-\frac{1}{2}\int_{\Omega}\frac{h(\phi^{*,n+1})}{\sqrt{E(\phi^{*,n+1})}}\hat{\phi})\frac{h(\phi^{*,n+1})}{\sqrt{E(\phi^{*,n+1})}},$
(27)

with the following boundary conditions

$\begin{cases}\mathbf{n}\cdot\nabla{\phi_{1}}^{n+1}+\frac{\gamma_{0}\eta\mu_{f}}{\lambda\Delta
t}\phi_{1}^{n+1}=\frac{\eta\mu_{f}}{\lambda\Delta t}\hat{\phi}-\\\ \hskip
85.35826pt\frac{1}{\lambda}f^{\prime}_{w}(\phi^{*,n+1})_{A,R}\qquad if\qquad
L^{*,n+1}_{A}L^{*,n+1}_{R}>0,\\\
\phi_{1}^{n+1}=\frac{1}{\gamma_{0}}\hat{\phi}\qquad otherwise\qquad
on\qquad\partial\Omega_{w}.\end{cases}$ (28)

$\displaystyle{\mathbf{n}}\cdot\nabla{G}_{1}^{n+1}=0\qquad
on\qquad\partial\Omega_{w},$ (29)

and the second coupled system of equations for $\phi_{2}$ and $G_{2}$ is

$\displaystyle\nabla^{2}G_{2}^{n+1}-\frac{\gamma_{0}}{M\Delta
t}\phi_{2}^{n+1}=0,$ (30)
$\displaystyle\nabla^{2}\phi_{2}^{n+1}+\frac{1}{\lambda}G_{2}^{n+1}=\frac{1}{2\lambda}\frac{h(\phi^{*,n+1})}{\sqrt{E(\phi^{*,n+1})}},$
(31)

with the following boundary conditions

$\begin{cases}\mathbf{n}\cdot\nabla{\phi_{2}}^{n+1}+\frac{\gamma_{0}\eta\mu_{f}}{\lambda\Delta
t}\phi_{2}^{n+1}=0\qquad if\qquad L^{*,n+1}_{A}L^{*,n+1}_{R}>0,\\\
\phi_{2}^{n+1}=0\qquad otherwise\qquad on\qquad\partial\Omega_{w},\end{cases}$
(32)

$\displaystyle{\mathbf{n}}\cdot\nabla{G}_{2}^{n+1}=0\qquad
on\qquad\partial\Omega_{w},$ (33)

First, the coupled equations (26-27) with the boundary conditions (28-29) are
solved, and then the coupled equations (30-31) with the boundary conditions
(32-33) are solved. The resultant coupled equations are non-symmetric linear
system so that GMRES iterative solver with geometric multigrid preconditioner
with the Message-Passing Interface (MPI) for distributed-memory
parallelization provided by PETSc library [Balay et al., 2019] has been used
to solve them. Then, $\phi^{n+1}$ and $G^{n+1}$ is computed by using Eq.
(24-25). At the end, the relaxation technique Eq. (22) is used to improve the
accuracy and consistency of the SAV scheme. The interested readers are
referred to [Jiang et al., 2022] for the definition of $\zeta_{0}$.

Figure 1: The Initial configuration for the droplet in the 2D with
$\theta_{0}$. $\theta_{D}$ has been extracted at the height $y=2\hskip
1.42271ptCn$ above the wall.

(a)

(b)

(c)

Figure 2: The steady-state deformed droplet profile $\phi=0$ at $Re=0.01$,
$\theta_{R}=10^{\circ}$, and $\theta_{A}=170^{\circ}$ (a) $Ca=0.05$ (b)
$Ca=0.1$ (c) $Ca=0.15$. The simulation data are compared with the results in
[Schleizer and Bonnecaze, 1999] and [Yue, 2020]

##  4 Validation

###  4.1 Droplet deformation with pinned contact lines

The implementation of the Cahn-Hilliard equation with CAH has been validated
against a test case from the literature [Schleizer and Bonnecaze, 1999, Yue,
2020]: a 2D Newtonian droplet in the Newtonian medium (N/N) with a large
hysteresis window $[\theta_{R},\theta_{A}]$. The computational domain is a 2D
rectangle, $\Omega=[0,8H]\times[0,2H]$, with the droplet placed on the bottom
wall, see Fig. 1. The no-slip and no-penetration conditions are imposed on the
both walls, and a fully developed parabolic velocity
$u=6\overline{U}((\frac{y}{L_{y}})-(\frac{y}{L_{y}})^{2})$ is imposed at the
inlet and outlet of the channel. The droplet initially is semicircle with
contact angle $\theta_{0}=60^{\circ}$, radius $R=0.9023H$, and height
$h=0.4511H$. This test case is defined by the length scale $L_{ref}=H$, which
is the half channel width; $u_{ref}=(\frac{3\overline{U}}{H})h$, which is the
velocity scale. This results in the following dimensionless numbers [Schleizer
and Bonnecaze, 1999, Dupont and Legendre, 2010, Liu et al., 2015, Yue, 2020].
From the physical parameters of the problem, we obtain the Reynolds number,
$Re=\frac{\rho_{1}u_{ref}H}{\mu_{1}}$, the ratio between inertial and viscous
forces; the Capillary number, $Ca=\frac{\mu_{1}u_{ref}}{\sigma}$, the ratio
between viscous and surface tension forces, density ratio,
$\lambda_{\rho}=\frac{\rho_{2}}{\rho_{1}}$, the ratio between ambient density
to droplet density, and the viscosity ratio,
$\lambda_{\mu}=\frac{\mu_{2}}{\mu_{1}}$, which is the ratio between ambient
viscosity to droplet viscosity. Furthermore, from the phase field model, we
obtain the Peclet number $Pe=\frac{2\sqrt{2}}{3}\frac{u_{ref}H\eta}{M\sigma}$,
the ratio between the advection and diffusion; the Cahn number,
$Cn=\frac{\eta}{H}$, which represents the ratio between the interface width
and the characteristic length scale; $\Lambda=\frac{\mu_{f}Cn}{\mu_{1}}$,
which characterizes the wall energy relaxation rate.

Following Yue [2020], we consider the following values of the relevant
dimensionless numbers:

$Re=0.01,\quad\lambda_{\mu}=1,\quad\lambda_{\rho}=1,\quad Cn=0.01,$
$Pe=\frac{2\sqrt{2}}{3}\frac{Ca}{Cn},\quad\frac{\mu_{f}}{\mu_{1}}=\frac{0.1}{Cn},\quad\theta_{R}=10^{\circ},\quad\theta_{A}=170^{\circ}.$

The comparison of the droplet interface $\phi=0$ with the results by Schleizer
and Bonnecaze [1999], Yue [2020] at different Capillary numbers when they
reach steady state represents a good agreement except at $Ca=0.15$ at the
receding contact angle, see Fig. 2; the same agreement and error at $Ca=0.15$
is indeed reported by [Yue, 2020], and it is due to a very small receding
contact angle.

##  5 Results

We scrutinize the effect of the viscoelasticity on both deformation and
depinning of a drop over a surface with different degrees of the hysteresis
window, $\Delta\theta=140^{\circ}$ with the choice of $\theta_{R}=30^{\circ}$
and $\theta_{A}=170^{\circ}$, $\Delta\theta=130^{\circ}$ with the choice of
$\theta_{R}=30^{\circ}$ and $\theta_{A}=160^{\circ}$, and
$\Delta\theta=100^{\circ}$ with the choice of $\theta_{R}=40^{\circ}$ and
$\theta_{A}=140^{\circ}$. These hysteresis windows are larger than encountered
in physical experiments, and are chosen to prolong the pinning of the contact
lines, following previous numerical works of [Dimitrakopoulos and Higdon,
1997, Schleizer and Bonnecaze, 1999, Dimitrakopoulos, 2007, Ding et al., 2010,
Dupont and Legendre, 2010]. A rectangular domain, $\Omega=[0,2H]\times[0,8H]$,
is used for the V/N system with the droplet placed on the middle of the bottom
wall, and a larger domain, $\Omega=[0,2H]\times[0,16H]$, for the N/V system in
order to allow a fully developed polymeric stresses in the inlet and outlet
boundaries. It should be noted that the velocity inside the droplet has been
set to zero in the beginning of the simulation, but a parabolic velocity
profile is imposed elsewhere.

In addition to the Newtonian dimensionless numbers introduced in section 4.1,
the following dimensionless numbers have been used to study the effect of the
viscoelasticity : the dimensionless cross-sectional area
$A=({4A_{d}}/{L^{2}_{y}})$ representing the ratio between the droplet area
($A_{d}$) and the square of channel half height, the Deborah number
$De=(\lambda_{H}u_{ref}/H)$ representing the ratio between elastic forces to
the viscous forces, and the polymeric viscosity ratio $\beta=(\mu_{s}/\mu)$
which is the ratio between the solvent viscosity and the total viscosity. The
dimensionless cross-sectional area is kept constant at $A=0.5$ in the
following sections (section 5.1 and section 5.2), and the semi-circular
droplet’s initial radius
$r_{0}=\sqrt{A_{d}H/(\theta_{0}-\sin{\theta_{0}}\cos{\theta_{0}})}$ and the
droplet’s initial height $h=r_{0}(1-\cos{\theta_{0}})$ are computed for a
specific initial contact angle ($\theta_{0}=90^{\circ}$). Furthermore, the
following dimensionless numbers are kept constant in all simulations: $Re=0.5$
(so that the inertial effects are negligible), $\lambda_{\mu}=1$,
$\lambda_{\rho}=1$, and $\beta=0.1$. It should be noted that $\lambda_{\mu}$
is the ratio between total viscosities of the two fluids; Since polymeric
viscosity increases with polymer concentration, the viscosity of the Newtonian
fluid also needs to increase if the polymer concentration changes (for example
by adding more viscous fluids).

Regarding numerical parameters related to the phase field model, we choose an
affordable but small capillary length ($Cn=\eta/H=0.00375$), and the Peclet
number is set to $Pe=(2\sqrt{2}/3)(Ca/Cn)$ according to the guidelines [Yue et
al., 2010, Yue, 2020] in order to reach the so-called sharp-interface limit;
the contact line friction is set to $\mu_{f}=(\Lambda/Cn)\mu_{1}$ [Carlson et
al., 2011, Yue, 2020] with $\Lambda=0.1$. The contact line is pinned when
$\frac{\partial{\phi}}{\partial{t}}|^{t^{n+1}}=0$, so the droplet’s
deformation and the contact line’s pinning/depinning are independent of the
choice of $\mu_{f}$ prior to depinning.

To analyse the results further, we have extracted different quantities in the
simulation, e.g., the locations of the receding contact line $X_{cl,R}$, and
the advancing contact line $X_{cl,A}$, and the values of the dynamic contact
angle at the receding side $\theta_{D,R}$, and at the advancing side
$\theta_{D,A}$. Here, the dynamic contact angle has been computed using the
geometric expression
$\cos{\theta_{D}}=(\mathbf{n}\cdot\nabla{\phi})/(|\nabla{\phi}|)$ at the
interface ($\phi=0$) at the height $y=2\hskip 1.42271ptCn$ above the wall
[Wang et al., 2015]. The velocity of the contact line ($U_{cl}=\frac{\sigma
U^{*}_{cl}}{\mu_{f}}$) is computed by using second-order central difference
with extracted $X_{cl}$ in the simulation, and the average contact line
velocity ($\overline{U}_{cl}$) is computed by taking average over a nominal
interval $[0,t_{end}]$.

In the following sections, we present results for non-Newtonian cases, where
either the droplet is laden with polymers and therefore viscoelastic (V/N), or
the doplet is Newtonian but the ambient fluid is viscoelastic (N/V). In both
cases, polymers affect the droplet deformation and dynamics in two ways: the
polymeric stresses affect the interfacial deformation, and they also induce
changes in the flow field which can affect the droplet dynamics [Ramaswamy and
Leal, 1999, Yue et al., 2005b, Aggarwal and Sarkar, 2007]. The stress tensor
excluding the term due to mixing energy can be written as [Yue et al., 2005a]

$\displaystyle\boldsymbol{\tau}=-p\boldsymbol{I}+\mu(\nabla{{\mathbf{u}}}+\nabla{{\mathbf{u}}^{T}})+\frac{(1\pm\phi)}{2}\boldsymbol{\tau}_{p}.$
(34)

Components of the stress tensor normal (${\mathbf{n}}_{s}$) and tangential
(${\mathbf{t}}$) to the interface ($\phi=0$), see Fig. 1, for the N/N, V/N,
and N/V systems are shown along the interface in order to perceive the
dominant forces there.

Furthermore, our numerical experiments showed that the polymeric stresses
inside the droplet were growing without bound for $De\geq 5$ when using
Oldroyd-B model ($\alpha=0$). It is known that results from Oldroyd-B model
may not be accurate in extensionally dominated flows, and it will be shown
later that there are such regions in these flows. Therefore, we employ the
Giesekus model by setting $\alpha=0.05$ for the V/N and $\alpha=[0.05,0.1]$
for N/V systems in order to prevent this issue while producing sufficiently
large polymeric stresses to be able to observe viscoelastic effects.

(a)

(b)

(c)

(d)

Figure 3: The time evolution of $\theta_{D}$ for the Newtonian and Giesekus
drops with $De=1$. N ( ) and VE ( ) refer to Newtonian and Giesekus drops
respectively. (a) $\theta_{D,R}$ over the surface with
$\Delta\theta=100^{\circ}$ (b) $\theta_{D,A}$ over the surface with
$\Delta\theta=100^{\circ}$ (c) $\theta_{D,R}$ over the surface with
$\Delta\theta=140^{\circ}$ (d) $\theta_{D,A}$ over the surface with
$\Delta\theta=140^{\circ}$.

(a)

(b)

(c)

(d)

(e)

Figure 4: Different components of the stresses along $\phi=0$ at time
$t=0.1\hskip 1.42271pt$s for drops with $Ca=0.2$ and $De=1$ over the surface
with $\Delta\theta=140^{\circ}$. (a) $\tau^{n}_{vi}$ (b) $\tau^{t}_{vi}$ (c)
$\tau^{n}_{p}$ (d) $\tau^{t}_{p}$ (e) $\tau^{n}_{pre}$.

(a)

(b)

(c)

(d)

Figure 5: The time evolution of $U_{cl,R}$ ($X^{t=0}_{cl,R}=3.436$) and
$U_{cl,A}$ ($X^{t=0}_{cl,A}=4.564$) for the Newtonian and Giesekus drops with
$De=1$. N and VE refer to Newtonian and Giesekus drops respectively. (a)
$U_{cl,R}$ for the surface with $\Delta\theta=100^{\circ}$ (b) $U_{cl,A}$ for
the surface with $\Delta\theta=100^{\circ}$ (c) $\overline{U}^{*}_{cl,R}$ vs
$Ca$ (d) $\overline{U}^{*}_{cl,A}$ vs $Ca$.

###  5.1 Elastic effects in the V/N system

In this section, we investigate the case where a viscoelastic droplet is
surrounded by a Newtonian fluid (V/N system), while the capillary number,
contact line hysteresis window and Deborah number are varied.

####  5.1.1 Effect of the capillary number

In this section, we investigate the effect of the capillary number on the
deformation of the viscoelastic droplet, at a fixed Deborah number ($De=1$)
and at different values of receding and advancing contact angles ($\theta_{R}$
and $\theta_{A}$). Figure 3 presents the time evolution of the instantaneous
contact angle $\theta_{D}$ at both the receding and advancing sides over two
different surfaces with $\Delta\theta=100^{\circ}$ and
$\Delta\theta=140^{\circ}$, and we observe that viscoelasticity of the droplet
affects them differently. Figures 3(a) and 3(c) show that the elasticity of
the droplet causes $\theta_{D,R}$ to attain a smaller steady-state value in
comparison to the Newtonian counterpart for small $Ca$ when the contact line
is pinned, and $\theta_{D,R}$ gains negligibly larger value when contact line
is moving for ($De\leq 1$), see Fig. 5(a). The advancing contact angle
$\theta_{D,A}$ reaches a smaller equilibrium value for $Ca\leq 0.2$ and a
larger value for $Ca>0.2$, see Figs. 3(b) and 3(c); It can be observed that
$\theta^{VE}_{D,A}$ is approaching $\theta^{N}_{D,A}$ and surpassing it as the
advancing contact line moves with considerable velocity which depends on both
$Ca$ and $\Delta\theta$, see Fig. 5(b). It is worth mentioning that for larger
values of $De$, the behaviour of $\theta_{D,R}$ is different, as will be
discussed in section 5.1.2.

The transient behavior of the $\theta_{D}$ reveals that the difference between
dynamic contact angles ($\theta^{VE}_{D}-\theta^{N}_{D}$) for a specific $Ca$
is not monotonic in time, and this trend is more pronounced for the
($\theta_{D,R}$) in the early times, see Figs. 3(a) and 3(c) for $t\leq
2\hskip 1.42271pt$s. The reason is that the Newtonian and viscoelastic
droplets deform at different rates. The viscoelastic drop speeds up and
deforms faster than the Newtonian counterpart at early times, and this
behavior has been observed in all of the cases in our simulation regardless of
the $Ca$ and $\Delta\theta$. Since the steady contact angle value is lower for
viscoelastic droplets, this causes the contact angle curves in fig. 3 to
intersect.

To find the reason for the faster deformation of viscoelastic droplets at
early times, we show different components of the stresses in Eq. (34) normal
and tangential to the interface ($\phi=0$) of the Newtonian and viscoelastic
droplets at $Ca=0.2$ over the surface with $\Delta\theta=140^{\circ}$ in
figure 4. We show the instantaneous stresses along the interface at
$t=0.1\hskip 1.42271pt$s, where the angle $\varphi$ along the interface is
measured clockwise from the negative $x$-direction. Let us first consider the
action of viscous stress in the Newtonian case (black markers in figure Fig. 4
(a)). The normal component of the viscous stress $\tau^{n}_{vi}$ at
$t=0.1\hskip 1.42271pt$s pushes the droplet inward at the back side for
$\varphi\lesssim 100^{\circ}$ and outward in the front side
$\varphi\approx[100^{\circ},172^{\circ}]$, both of which are in favor of the
droplet deformation. Near the advancing contact line however, the viscous
normal stress pushes the droplet inward at $\varphi\gtrsim 172^{\circ}$, which
resists the increase of $\theta_{D,A}$ near the wall. Let us now consider the
viscoelastic case, green markers in figure 4 a. Viscous normal stress has the
same effects as in the Newtonian case, but weaker. Also, the tangential
component of viscous stress (Fig. 4 (b)) which is related to the acceleration
of the fluid elements near the interface is weaker in the viscoelastic flow.
Polymeric stresses act similarly to viscous stresses, but have a significantly
smaller amplitude (Figs. 4 (c,d)).

Thus, neither the viscous stress nor the polymeric stress is influential in
the faster deformation of the viscoelastic droplet at the early times. The
reason is instead found from the pressure distribution at the interface (Fig.
4 (e)). The pressure pulls the droplet inward for all values of the $\varphi$,
and its magnitude is larger in comparison to other components due to the
presence of the surface tension in the interface (Young–Laplace equation). The
magnitude of the pressure is larger for $\varphi\lesssim 20^{\circ}$ and
smaller for $\varphi\gtrsim 158^{\circ}$ for the viscoelastic droplet in
comparison to the Newtonian counterpart; thus, the pressure is responsible for
the faster decrease of the dynamic contact angle in the receding side, and
faster increase on the advancing side, both of them increasing the deformation
of the viscoelastic droplet at $t=0.1\hskip 1.42271pt$s. The change in the
pressure is the consequence of the modification of the overall flow field by
polymeric stresses.

Next, we will look at the effect of contact line hysteresis window
$\Delta\theta=\theta_{A}-\theta_{R}$ on pinning and velocity of the contact
lines. Firstly, we observe that as hysteresis increases, both Newtonian and
viscoelastic droplet are pinned at higher capillary numbers. Secondly, for
droplets that are not pinned, the steady state contact line speed is higher
for viscoelastic droplets when $Ca$ is larger ($Ca\geq 0.25$). Figs. 5(a,b)
display the time evolution of the contact line velocity for different $Ca$
over surfaces with $\Delta\theta=100^{\circ}$, and Figs. 5(c,d) present the
variation of averaged dimensionless contact line velocity
($\overline{U}^{*}_{cl}$) vs. $Ca$ over surfaces with different
$\Delta\theta$. To characterise the depinning process, the time for the
inception of motion ($U_{cl}(t_{in})\approx 10^{-4}\hskip
1.42271pt\frac{m}{s}$) in both contact lines have been marked on their
corresponding velocity curve for the hysteresis window
$\Delta\theta=100^{\circ}$, Figs. 5(a,b). Both of the contact lines are pinned
when $Ca<0.2$. When $Ca\geq 0.2$, the advancing contact line starts to move
first (when the $\theta_{D,A}|_{y=0}\leq\theta_{A}$ at $t\leq t_{in}$) for the
Newtonian drop; the receding contact line of the Giesekus drop moves first in
this cases. At this moderate $De$, the both advancing and receding contact
lines of the viscoelastic droplet move faster over all surfaces (different
$\Delta\theta$) with $Ca=0.25$, but the advancing contact line of the Giesekus
drop moves faster when $Ca=0.25$ and slower when $Ca=0.2$ in comparison with
its Newtonian counterpart, see Figs. 5(c,d).

The contact line speed can be approximated by the difference between the
dynamic contact angle $\theta$ and the equilibrium contact angle $\theta_{e}$
($\theta_{A}$ and $\theta_{R}$ in the present of CAH) for minor variation of
$\theta$ from $\theta_{e}$ [Yue and Feng, 2011, Amberg, 2022, Yada et al.,
2023]: $\frac{\mu_{f}}{\sigma}U_{cl}\propto\theta-\theta_{e}$, where the
contact line friction $\mu_{f}$ is a constant determined by surface and fluid
properties. Hence, the contact line speed is enhanced if the droplet
experiences more bending causing the dynamic contact angle to differ from the
equilibrium contact angle. At the receding contact line, more bending implies
a decrease in $\theta_{D,R}$.The elasticity ($De=1$) of the droplet enhances
the viscous bending at the receding contact line
($\theta^{VE}_{D,R}\leq\theta^{N}_{D,R}$) after its depinning over a large
span of time, see Figs. 3(a,c), for the simulations in this section and
therefore $\theta^{VE}_{D,R}\leq\theta^{N}_{D,R}$ implies that
$U^{Ve}_{cl}\geq U^{N}_{cl}$. However, the viscous bending at the advancing
contact line is initially weakened by the elasticity of the droplet
($\theta^{VE}_{D,A}\leq\theta^{N}_{D,A}$). Thus, the elasticity of the droplet
causes $U^{Ve}_{cl}\leq U^{N}_{cl}$ when $\theta^{VE}_{D,A}<\theta^{N}_{D,A}$
and $U^{Ve}_{cl}\geq U^{N}_{cl}$ when $\theta^{VE}_{D,A}>\theta^{N}_{D,A}$.

In this subsection, we have studied the effects of moderate viscoelasticity on
the deformation of the droplet and movement of both contact lines by comparing
a Newtonian droplet and a viscoelastic droplet with $De=1$. Some of the
conclusions will change when viscoelasticity is further increased in the next
subsection.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6: The effect of $De$ on the time evolution of $\theta_{D}$ at constant
capillary numbers in the V/N system. N and VE refer to Newtonian and Giesekus
drops respectively. (a) $Ca=0.15$ and $\Delta\theta=100^{\circ}$ (b) $Ca=0.15$
and $\Delta\theta=140^{\circ}$ (c) $Ca=0.2$ and $\Delta\theta=100^{\circ}$ (d)
$Ca=0.2$ and $\Delta\theta=140^{\circ}$ (e) $Ca=0.25$ and
$\Delta\theta=100^{\circ}$ (f) $Ca=0.25$ and $\Delta\theta=140^{\circ}$.

####  5.1.2 Effect of increasing Deborah number

The effect of increasing the droplet’s elasticity at a fixed capillary number
has been investigated by changing $De$ over the surfaces with
$\Delta\theta=100^{\circ}$ and $\Delta\theta=140^{\circ}$. The main effect
observed is that bending of both contact lines is decreased significantly with
increasing Deborah number to $De=5-10$. Thus, we have a non-monotonic behavior
with respect to $De$ in the receding contact line; the elasticity of the
droplet enhances the viscous bending for droplets with small elasticity, but
the viscous bending is decreased for droplets with larger $De$. This non-
monotonic behaviour in the receding contact angle also changes the receding
contact line’s velocity trend, and the velocity of the receding contact line
decreases for $De=[5,10]$ in comparison to the Newtonian counterpart in all
cases, see Fig. 7. Further increasing the elasticity of the droplet depending
on the hysteresis window of the surface and capillary number brings the
contact lines to a halt, see Figs. 7(a) and 7(b) and for $De=10$, $Ca=0.2$,
and $\Delta\theta=100^{\circ}$.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 7: The effect of $De$ on the time evolution of $U_{cl,R}$
($X^{t=0}_{cl,R}=3.436$) and $U_{cl,A}$ ($X^{t=0}_{cl,A}=4.564$) at constant
capillary numbers in the V/N system. N and VE refer to Newtonian and Giesekus
drops respectively. (a) $U_{cl,R}$ for $Ca=0.2$ with
$\Delta\theta=100^{\circ}$ (b) $U_{cl,A}$ for $Ca=0.2$ with
$\Delta\theta=100^{\circ}$ (c) $U_{cl,R}$ with $Ca=0.25$ and
$\Delta\theta=100^{\circ}$ (d) $U_{cl,A}$ with $Ca=0.25$ and
$\Delta\theta=100^{\circ}$ (e) $\overline{U}^{*}_{cl,R}$ vs $Ca$ (f)
$\overline{U}^{*}_{cl,A}$ vs $Ca$.

The advancing dynamic contact angle of the viscoelastic droplet generally
attains a smaller steady-state value with further increasing the elasticity to
$De=10$ (Fig. 6), and the velocity of the advancing contact line decreases
with increasing the elasticity of the droplet, see Figs. 7(b), 7(d) and 7(f).
The decrease in the velocity of the advancing contact line is in line with the
decreasing $\theta_{D,A}$. The same trend for $U_{cl,R}$ exists in the
receding contact line for $De=[5-10]$ but $\theta_{cl,R}$ gains a larger
value, see Figs. 7(a), 7(c) and 7(d).

(a)

(b)

(c)

(d)

Figure 8:  Time evolution of viscoelastic drops with $Ca=0.2$ over the surface
with $\Delta\theta=100^{\circ}$. The interface $\phi=0$ of the Newtonian and
Giesekus drops are depicted with a white and black, respectively. The trace of
the conformation tensor has been visualized on the background, and the
principal eigenvector of the conformation tensor is represented by purple line
segments.(a) $De=10$ at $t\approx 0.15\hskip 1.42271pt$s, (b) $De=10$ at
$t\approx 0.6\hskip 1.42271pt$s, (c) $De=10$ at $t\approx 1.8\hskip
1.42271pt$s, (d) $De=10$ at $t\approx 6\hskip 1.42271pt$s.

To gain a better understanding of the elastic effects, we have visualized the
trace of the conformation tensor $\mathbf{c}$ defined as
$\mathbf{\boldsymbol{c}}=(\lambda_{H}/\mu_{p})\mathbf{\boldsymbol{\tau}_{p}}+\mathbf{\boldsymbol{I}}$,
which is a measure of stretching of the polymer molecules, in figure 8 for the
highest Deborah number $De=10$ at $Ca=0.2$ and $\Delta\theta=100^{\circ}$
(corresponding to drops in Figs. 6(c) and 7(a,b)). The results for $De=5$ are
very similar to $De=10$. The interface (defined by $\phi=0$) of the
viscoelastic drop is shown with a black contour, and the corresponding
Newtonian drop by a white contour, with the trace of the conformation tensor
in the background, see Fig.8. Furthermore, we show by arrows the principal
direction of the conformation tensor, representing the orientation of the
dumbbells. The pull of the polymers affects the interface position most in
regions where the arrows are perpendicular to the interface ([Bird et al.,
1987b], pp. 64-69).

At early times (Figs. 8a), the viscoelastic drop speeds up and deforms faster
than the Newtonian counterpart. The polymer molecules are most stretched in
the vicinity of the receding and advancing contact lines at $t\approx
0.15\hskip 1.42271pt$s, but the dumbbells have acute angle with the interface
there, and hence do not influence the interface. At $t\approx 0.6\hskip
1.42271pt$s (Fig. 8(b)), the receding and advancing contact lines of the
Newtonian drop start moving already. However, the polymeric stresses have now
built up (which takes a finite time proportional to $\lambda_{H}$). Hence,
deformation and time evolution of the viscoelastic droplet have been retarded,
and both contact lines of the viscoelastic droplets remain pinned for
$De=5-10$. Furthermore, the dumbbells are becoming perpendicular to the
interface in $\varphi\approx[140^{\circ},160^{\circ}]$ and because polymers
are also very elongated here, the stresses should have a large influence on
the interface.

(a)

(b)

(c)

(d)

(e)

Figure 9: Different components of the stresses along $\phi=0$ at time
$t=0.6\hskip 1.42271pt$s for drops with $Ca=0.2$ over the surface with
$\Delta\theta=100^{\circ}$. (a) $\tau^{n}_{vi}$ (b) $\tau^{t}_{vi}$ (c)
$\tau^{n}_{p}$ (d) $\tau^{t}_{p}$ (e) $\tau^{n}_{pre}$

To compare the relative importance of polymeric and flow-induced stresses, we
again show them (Eq. 34) along the interface $\phi=0$, separated into normal
and tangential components, in Fig. 9. Increasing viscoelasticity to $De=5-10$
suppresses the tangential component (Fig. 9(b)) of the viscous stress along
the interface (
$\left[\frac{\tau^{t}_{vi}|_{VE}}{\tau^{t}_{vi}|_{N}}\right]_{max}\approx
0.3$), and this decrease is mainly responsible for the delay in the
deformation of the viscoelastic droplets. The change in viscous stresses due
to viscoelasticity has a larger magnitude than the previously analysed
polymeric stresses (Fig. 9(c,d)) which are only active in the segment
$\varphi\approx[110^{\circ},170^{\circ}]$. Further, the polymeric stresses are
largest for $De=5$ and decrease at $De=10$, while the contact line stays
pinned longest for $De=10$. Hence, the retardation and pinning of the
viscoelastic droplet is mainly caused by velocity modification in viscoelastic
flows around the interface, rather than by the polymeric stresses themselves.

Moreover, another important rheological property of the viscoelastic fluid is
the strain-hardening which usually occurs in extensional flow, and can be
attributed to the elongational viscosity (Bird et al. [1987a], p.p 132)
defined by $\mu_{e}\sim\mu_{e}(N_{1},\dot{\epsilon})$ where
$N_{1}=\tau_{xx}-\tau_{yy}$ is the first normal stress and $\dot{\epsilon}$ is
the elongational strain rate. Since the modification of the flow field is
considerable in particular for $De\geq 5$ when the polymeric stresses are
built up, we use the flow topology parameter [Lee et al., 2007] to find
extensional flow regions inside the droplets:

$\displaystyle\xi=\frac{|\boldsymbol{D}|-|\boldsymbol{\Omega}|}{|\boldsymbol{D}|+|\boldsymbol{\Omega}|}$
(35)

where
$\boldsymbol{D}=\frac{\nabla{{\mathbf{u}}}+\nabla{{\mathbf{u}}^{T}}}{2}=\frac{\dot{\boldsymbol{\gamma}}}{2}$
and
$\boldsymbol{\Omega}=\frac{\nabla{{\mathbf{u}}}^{T}-\nabla{{\mathbf{u}}}}{2}$
are the deformation and vorticity tensors respectively, and the magnitude of a
tensor $\boldsymbol{Q}$ can be computed by
$|\boldsymbol{Q}|=\sqrt{\frac{1}{2}\sum_{i}\sum_{j}Q^{2}_{ij}}$. The flow
parameter can vary between $\xi=-1$ corresponding to pure solid-like rotation
to $\xi=1$ corresponding to pure extensional flow, and $\xi=0$ presenting pure
shear flow.

To further contrast the effects at small and large elasticity, we compare the
flow topology of the viscoelastic droplets at $De=0$ (Newtonian case, Fig. 10
(a)), $De=1$ (Fig. 10 (b)) and $De=5$ (Fig. 10(c)). The flow type inside the
viscoelastic droplet with $De=5$ is more extensional (red color) at the front
where viscoelastic stresses pull the droplet interface, in comparison to the
other droplets. This is the extensional flow region, where extensional
viscosity may be higher. The extensional viscosity as a function of the
extensional strain rate at steady-state for the Giesekus model ($\alpha\neq
0$) can be expressed as: ([Bird et al., 1987a], p.p 368)

$\displaystyle\mu^{+}_{e}=\frac{\mu_{e}-\mu_{s}}{3\mu_{p}}=\frac{1}{6\alpha}\Big{[}3+\frac{1}{\lambda_{H}\dot{\epsilon}}\Big{[}\sqrt{1-4(1-2\alpha)\lambda_{H}\dot{\epsilon}+4(\lambda_{H}\dot{\epsilon})^{2}}-$
$\displaystyle\sqrt{1+2(1-2\alpha)\lambda_{H}\dot{\epsilon}+(\lambda_{H}\dot{\epsilon})^{2}}\Big{]}\Big{]}$
(36)

Fig. 10(d) presents $\mu^{+}_{e}$ as a function of strain rate111The magnitude
of strain rate in the extensional flow dominant regions in Figs. 10(b) and
10(c) is mostly less or equal to 2 ($|\dot{\boldsymbol{\gamma}}|\leq 2$)
$\dot{\epsilon}$ for three Deborah numbers. The magnitude of $\mu^{+}_{e}$ for
$De\geq 5$ increases steeply already at small elongational strain rates. The
strain rate for $De=5$ shown in figure 10 (f) is highest in shear-dominated
regions at the top, but also moderate in extensional flow areas at the front,
and therefore we should expect a large extensional viscosity resisting
deformation at the front. As observed, the shear strain rate is always higher
in the viscoelastic fluid than in the Newtonian fluid, and this is a
consequence of keeping the total viscosity ratio constant ($\lambda_{\mu}=1$);
the solvent viscosity of the viscoelastic droplet is smaller. The magnitude of
$N_{1}$ is largest in the top side where shear flow dominates and increases
with $De$,but those regions of the interface do not directly influence the
contact line movement.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 10: Different flow quantities of the Newtonian and Giesekus drops with
Ca = 0.2 over $\Delta\theta=140^{\circ}$ at $t\approx 6\hskip 1.42271pt$s.(a)
The flow parameter ($\xi$) on the background and the total stress vector
(${\mathbf{n}}_{s}\cdot\boldsymbol{\tau}$ in Eq. 34) along the interface
$\phi=0$ for the Newtonian droplet with $Ca=0.2$ and
$\Delta\theta=140^{\circ}$ at $t\approx 6\hskip 1.42271pt$s, (b) the same for
a Giesekus drop with $Ca=0.2$, $De=1$, (c) the same at $Ca=0.2$, $De=5$. (d)
The steady-state dimensionless elongational viscosity as a function of the
elongational strain rate for the Giesekus model with different $De$
corresponding to $Ca=0.2$ and $\Delta\theta=140^{\circ}$.(e) The magnitude of
strain rate ($|\dot{\boldsymbol{\gamma}}|$) for the droplet at $De=1$, and (f)
the same at $De=5$. (g) The first normal stress difference at $De=1$, and (h)
the same at $De=5$.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11: The time evolution of dynamic contact angle $\theta_{D}$ over a
surface with $\Delta\theta=140^{\circ}$ in the N/V system. N and VE refer to
Newtonian and Giesekus matrices respectively. (a) $\theta_{D,R}$ for $Ca=0.15$
(b) $\theta_{D,A}$ for $Ca=0.15$ (c) $\theta_{D,R}$ for $Ca=0.2$ (d)
$\theta_{D,A}$ for $Ca=0.2$ (e) $\theta_{D,R}$ for $Ca=0.25$ (f)
$\theta_{D,A}$ for $Ca=0.25$. The symbols ( ) and ($\star$) correspond to
$\alpha=0.05$ and $\alpha=0.1$ respectively.

The magnitude and tendency (direction) of the total stress vector in Eq. 34
along the interface $\phi=0$ is depicted by arrows on top of figures 10
(a,b,c):(${\mathbf{n}}_{s}\cdot\boldsymbol{\tau}$, where ${\mathbf{n}}_{s}$ is
the unit vector normal to the interface. The stresses are mostly pulling the
interface inwards for both Newtonian (10a) and viscoelastic (10b,c) droplets.
However, in the vicinity of the receding contact line marked with circle, the
stresses are pulling the interface outward for the viscoelastic droplets. The
magnitude of these outward stresses increases with $De$ (elasticity of the
droplet) so that they are bringing about a larger $\theta_{D,R}$. The outward
stresses around the receding contact line were not observed when capillary
number was lowered to $Ca=0.15$, at least until $De=10$.

Concluding, the non-monotonic behavior in the receding dynamic contact angle
with increasing $De$ can be attributed to the combination of two effects: the
outward stresses in the vicinity of the receding contact line, and the
inception of the strain-hardening occurring for $De\geq 5$. Also the smaller
advancing dynamic contact angle at steady-state might be attributed to strain-
hardening.

(a)

(b)

Figure 12: The time evolution of $X_{cl,R}$ ($X^{t=0}_{cl,R}=7.436$) and
$X_{cl,A}$ ($X^{t=0}_{cl,A}=9.564$) at different $De$ in the N/V system. N and
VE refer to Newtonian and Giesekus matrices respectively. (a) $Ca=0.25$ over a
surface with $\Delta\theta=140^{\circ}$ (b) $Ca=0.25$ over a surface with
$\Delta\theta=100^{\circ}$. The symbols ( ) and ($\star$) correspond to
$\alpha=0.05$ and $\alpha=0.1$ respectively.

###  5.2 Effect of matrix viscoelasticity on a Newtonian droplet

In this section, we consider the effect of matrix viscoelasticity on the
deformation and depinning of a Newtonian droplet (N/V system) in the presence
of contact line hysteresis. The same boundary and initial conditions have been
applied as previously, with the polymeric stresses are set to zero in the
beginning of the simulation in the Giesekus medium. Figs. 11(a,c,e) present
the time evolution of $\theta_{D,R}$ of a Newtonian drop in the Giesekus
medium over a surface with $\Delta\theta=140^{\circ}$, and it can be seen that
the elasticity of the matrix has suppressed the drop’s deformation in the
receding side significantly for $De\geq 5$, in other words, the viscous
bending in the receding side has weakened drastically. The effect of changing
capillary number on the steady-state value of $\theta_{D,R}$ is limited and
small in the range $Ca=[0.15,0.25]$ used in our simulations when $De\geq 5$.
The same trend can be observed on the advancing side, where the
viscoelasticity of the medium decreases the viscous bending in particular for
$De=5-10$; Furthermore, the shear-thinning also plays an important role at
higher De by further reducing the viscous bending, see Figs. 11(b,d,f).

It is worth mentioning that the decreased viscous bending results in that both
contact angles are insensitive to capillary number for the Newtonian drop in
the viscoelastic matrix with high elasticity. For a Newtonian droplet in a
Newtonian matrix, when increasing capillary number from $Ca=0.15$ to
$Ca=0.25$, both receding and advancing contact angles vary by $\propto
30^{\circ}$. For a Newtonian droplet in a viscoelastic matrix (at $De=5$ and
$De=10$), the change with capillary number is less than $10^{\circ}$. Both
contact lines are pinned at $Ca=0.15-0.25$ over the surface with
$\Delta\theta=140^{\circ}$ in a viscoelastic matrix (Fig. 12(a)), while the
receding contact line of the Newtonian drop in the Newtonian medium moves with
a significant velocity at $Ca=0.25$. If the hysteresis window is decreased to
$\Delta\theta=100^{\circ}$ (Fig. 12(b)), the advancing contact line of the
least elastic medium ($De=1$) starts moving at $t\approx 0.85\hskip
1.42271pt$s and the receding contact line remains pinned, while both contact
lines start to move simultaneously in a Newtonian matrix with appreciable
velocities.

(a)

(b)

(c)

(d)

Figure 13:  Time evolution of the Newtonian drop in viscoelastic matrix with
$Ca=0.25$ over the surface with $\Delta\theta=100^{\circ}$. The interface
$\phi=0$ of the Newtonian and Giesekus drops are depicted with a white and
black lines, respectively. The trace of the conformation tensor has been
visualized on the background, and the principal eigenvector of the
conformation tensor is represented by purple line segments. (a) $De=1$ at
$t\approx 0.6\hskip 1.42271pt$s, (b) $De=5$ at $t\approx 0.6\hskip
1.42271pt$s, (c) $De=1$ at $t\approx 1.8\hskip 1.42271pt$s, (d) $De=5$ at
$t\approx 1.8\hskip 1.42271pt$s.

(a)

(b)

(c)

(d)

(e)

Figure 14: Different components of the stresses along $\phi=0$ at time
$t=0.6\hskip 1.42271pt$s for drops with $Ca=0.25$ over the surface with
$\Delta\theta=100^{\circ}$ in N/V system. (a) $\tau^{n}_{vi}$ (b)
$\tau^{t}_{vi}$ (c) $\tau^{n}_{p}$ (d) $\tau^{t}_{p}$ (e) $\tau^{n}_{pre}$.

The trace of conformation tensor along with principal stress orientation is
shown at $De=1$ 13(a,c) and $De=5$ 13(c,d), respectively, over a surface with
$\Delta\theta=100^{\circ}$. The stretching of polymer molecules is small in
the vicinity of the contact lines, while large stretching occurs around the
interface for $\varphi\approx[54^{\circ},138^{\circ}]$ since dumbbells are no
longer parallel to the interface in this portion (as they were in the V/N
case, fig. 8), this implies that large polymeric stresses is expected over a
large portion of the interface.

The normal component of the viscous stress $\tau^{n}_{vi}$ has been lessened
significantly due to the elasticity of the matrix at $t=0.6\hskip 1.42271pt$s
for $\varphi\approx[0,120^{\circ}]$ responsible for pulling the drop inward
and for $\varphi\approx[127^{\circ},164^{\circ}]$ responsible for pulling the
drop outward, see Fig. 14(a); in the same manner, the tangential component of
the viscous stress $\tau^{t}_{vi}$ has been weakened significantly for all
values of $\varphi$, see Fig. 14(b). Figs. 14(c) and 14(d) present that both
normal and tangential components of the polymer stresses corresponding to the
matrix with $De=5$ have been lessened in comparison to the medium with $De=1$.
It should be noted that the same trend exists between the the polymer stresses
at $t\approx 1.8\hskip 1.42271pt$s for mediums with $De=1$ and $De=5$. These
variations in the stresses along the interface ($\phi=0$) with increasing the
elasticity of the matrix reveals that the immobilization of the drop in the
matrix with $De\geq 5$ is due to the lessened viscous and polymeric stresses
along the interface.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 15: Different flow quantities of the Newtonian drop in viscoelastic
matrices with $Ca=0.25$ over the surface with $\Delta\theta=100^{\circ}$ at
$t\approx 1.8\hskip 1.42271pt$s. (a) flow parameter ($\xi$) with $De=1$, (b)
flow parameter ($\xi$) with $De=5$, (c) first normal stress ($N_{1}$) with
$De=1$, (d) first normal stress ($N_{1}$) with $De=5$, (e) the magnitude of
strain rate ($|\dot{\boldsymbol{\gamma}}|$) with $De=1$, (f) the magnitude of
strain rate ($|\dot{\boldsymbol{\gamma}}|$) with $De=5$.

The flow type becomes more extensionally dominated in the vicinity of the
receding and advancing contact lines with increasing the matrix elasticity
($De\geq 5$), in particular near the interface, see Figs. 15(a) and 15(b).
Large first normal stress $N_{1}$ is accompanied with dominating shear flow in
the matrix around the interface for $\varphi\approx[35^{\circ},145^{\circ}]$
so that the shear-thinning is the effective mechanism in this region of the
interface, and this shear-thinning effect is reflected in decrease in
$\theta_{A}$ with increasing $\alpha$, see Figs. 11(b), 11(d), 11(f). The
strain rate is very small and almost negligible within the matrix with $De\geq
5$ in the vicinity of the interface for $\varphi\approx[0^{\circ},22^{\circ}]$
and $\varphi\approx[170^{\circ},180^{\circ}]$, see Fig. 15(f) so that the
observed $tr(\boldsymbol{c})\approx 3$ in this region (Fig. 13(d)), is the
consequence of the small strain rate in these regions.

Summarizing, also viscoelasticity in the matrix phase in the N/V system
suppresses the deformation of the droplet and pins the contact lines, and does
so more efficiently than the viscoelasticity of the droplet in a V/N system.
In both cases, the effect of viscoelasticity mainly stems from changes in the
flow field and weakening the flow gradients around the droplet such as
tangential viscous shear, which results in decreasing viscous bending and
hence decreased movement of the contact lines.

##  6 Conclusions

In the present study, two-dimensional numerical simulations have been
conducted to investigate the effect of viscoelasticity either in the drop or
in the matrix on the depinning and deformation of a drop adhering to the wall
in Poiseuille flow at low Reynolds number, under the influence of contact
angle hysteresis. The results indicate that the viscoelastic properties have a
significant impact on the drop’s deformation and depping of the contact lines.

We first compare Newtonian and viscoelastic droplets at a low constant Deborah
number ($De=1$) and different capillary numbers. The elasticity of the drop
causes the viscoelastic drop to deform faster than the Newtonian counterpart
at the early initial times, but the Newtonian drop then catches up. The
viscoelastic drop continues to deform, and its receding dynamic contact angle
reaches a smaller final value than the Newtonian counterpart prior to
depinning; the advancing dynamic contact angle of a viscoelastic drop reaches
a smaller value if the receding contact angle is pinned, otherwise it attains
a larger final value. This final stage of the deformation depends on both $Ca$
and $\Delta\theta$ for a fixed $De$.

Secondly, we investigated the effect of increasing elasticity ($De=5$ and
$De=10$) over the surfaces with $\Delta\theta=140^{\circ}$ and
$\Delta\theta=100^{\circ}$. In the first stage of the deformation, the effect
of increasing $De$ is negligible, but in the second stage, the time evolution
of the dynamic contact angle at both receding and advancing sides is hindered
significantly with non-monotonic behavior in the receding dynamic contact
angle. This trend is attributed to the emergence of the pulling outward
stresses changing to inward along the interface depending on both $Ca$ and
$De$ and the inception of the strain-hardening inside the drop at both
receding and advancing sides.

Finally, we scrutinized the effect of matrix viscoelastic properties in N/V
system over surfaces with $\Delta\theta=140^{\circ}$ and
$\Delta\theta=100^{\circ}$. The viscoelasticity of the matrix lessens the
deformation of the Newtonian drop, in other words, the viscous bending has
been weakened in both receding and advancing sides due to the presence of
polymer molecules in the medium, and these reductions are more pronounced for
$De\geq 5$. These reductions are bring the drops in the N/V system to halt in
the matrix with large elasticity ($De\geq 5$). Furthermore, the shear-thinning
effect reduces the the viscous bending further in the advancing side by
increasing $\alpha$ from 0.05 to 0.1.

##  7 Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

##  8 Acknowledgments

This work was funded by European Research Council (ERC) through Starting grant
no. 852529 MUCUS, by the Swedish Research Council through grant VR 2017-0489,
and by Swedish e-Science Research Centre. We acknowledge the computing time on
the supercomputer Beskow at the PDC center, KTH provided by SNIC (Swedish
National Infrastructure for Computing),Sweden.

##  9 DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Appendix A:

The velocity of the contact line can be related to the static and dynamic
contact angles over a surface without hysteresis [Yue and Feng, 2011], and a
similar relation holds for the velocity of contact line over a surface with
$\theta_{A}$ and $\theta_{R}$. The phase-field variable in our simulation is a
function of $(x,y,t)$, and it is a function of two variable if we constrain it
to the wall $y=0$ so that it can be wrriten as

$\displaystyle\phi=\phi(x,t)\hskip 14.22636pton\hskip
14.22636pt\partial\Omega_{w}$ (A1)

it is initialized at $t=0$ ($\phi_{0}=\phi(x,0)$) on the wall and evolved over
time by Eq. 13. The level curves of $\phi(x,t)=C_{\phi}$, where $-1\lesssim
C_{\phi}\lesssim 1$, are curves defined on the $(x,t)$ plane for different
values of $C_{\phi}$. If we take the total differential of Eq. A1 and put it
equal zero
$d\phi=\frac{\partial{\phi}}{\partial{x}}dx+\frac{\partial{\phi}}{\partial{t}}dt=0$,
we can get a differential equation for the level curves on the $(x,t)$ plane

$\displaystyle\frac{dx}{dt}=-\frac{\frac{\partial{\phi}}{\partial{t}}}{\frac{\partial{\phi}}{\partial{x}}}$
(A2)

since we are interested on the time evolution of the contact line
($X_{cl}(t)$) on the $(x,t)$ plane and its velocity
$U_{cl}=\frac{dx}{dt}|_{\phi=0}$, we choose $C_{\phi}=0$ and evaluate the
$\frac{\partial{\phi}}{\partial{t}}|_{\phi=0}$ and
$\frac{\partial{\phi}}{\partial{x}}|_{\phi=0}$. We compute
$\frac{\partial{\phi}}{\partial{t}}|_{\phi=0}$ by using Eq. 13, and we assume
the interface intersecting the wall is a planar interface

$\displaystyle-\frac{\partial{\phi}}{\partial{t}}|_{\phi=0}=\frac{1}{\mu_{f}\eta}(\frac{3}{2\sqrt{2}}\sigma\eta\mathbf{n}\cdot\nabla{\phi}+f^{\prime}_{w}(\phi))|_{\phi=0}$
$\displaystyle=$
$\displaystyle\frac{3}{2\sqrt{2}}\frac{\sigma}{\mu_{f}}|\nabla{\phi}|^{A/R}_{\phi=0}(\cos{(\theta_{d})}|^{A/R}-\frac{1}{\sqrt{2}\eta|\nabla{\phi}|_{\phi=0}}\cos{(\theta_{A/R})})$
(A3)

where
$\mathbf{n}\cdot\nabla{\phi}|_{\phi=0}=|\nabla{\phi}|_{\phi=0}\cos{(\theta_{d}})$
is used. We then need to compute
$\frac{\partial{\phi}}{\partial{x}}|_{\phi=0}$

$\displaystyle\frac{\partial{\phi}}{\partial{x}}|_{\phi=0}=-|\nabla{\phi}|_{\phi=0}\cos{(\theta_{d}-\frac{\pi}{2})}=-|\nabla{\phi}|_{\phi=0}\sin{(\theta_{d}})$
(A4)

substituting Eqs. (A3-A4) in Eq. A2, we get the contact line velocity as

$\displaystyle\frac{\mu_{f}U_{cl}}{\sigma}=\frac{3}{2\sqrt{2}}\frac{1}{\sin{(\theta_{d}})}(\frac{1}{\sqrt{2}\eta|\nabla{\phi}|}\cos{(\theta_{A/R})}-\cos{(\theta_{d}}))$
(A5)

where $\theta_{R}$ is chosen when $\theta_{d}\leq\theta_{R}$ and $\theta_{A}$
if $\theta_{d}\geq\theta_{A}$ in the above equation. It should be noted that
$U_{cl}<0$ when the fluid with $\phi=-1$ replaces the fluid with $\phi=1$ if
not $U_{cl}>0$.

We take the derivative of Eq. A5 with respect to $\theta_{d}$ since we just
want to see the effect of varying $\theta_{d}$ on the $U_{cl}$ with having the
same $\theta_{A/R}$

$\displaystyle\frac{\mu_{f}}{\sigma}\frac{\partial{U_{cl}}}{\partial{\theta_{d}}}=(\frac{\partial{}}{\partial{\theta_{d}}}\frac{1}{\sin{(\theta_{d}})})(\frac{3}{2\sqrt{2}}(\frac{1}{\sqrt{2}\eta|\nabla{\phi}|}\cos{(\theta_{A/R})}-\cos{(\theta_{d}})))+$
$\displaystyle\frac{3}{2\sqrt{2}}\frac{1}{\sin{(\theta_{d})}}(\frac{\cos{(\theta_{A/R})}}{\sqrt{2}\eta}\frac{\partial{}}{\partial{\theta_{d}}}\frac{1}{|\nabla{\phi}|}-\frac{\partial{\cos{(\theta_{d}})}}{\partial{\theta_{d}}})$
$\displaystyle=$
$\displaystyle\frac{-\cos{(\theta_{d})}}{\sin^{2}{(\theta_{d})}}(\frac{3}{2\sqrt{2}}(\frac{1}{\sqrt{2}\eta|\nabla{\phi}|}\cos{(\theta_{A/R})}-\cos{(\theta_{d}})))+$
$\displaystyle\frac{3}{2\sqrt{2}}\frac{1}{\sin{(\theta_{d})}}(\frac{\cos{(\theta_{A/R})}}{\sqrt{2}\eta}\frac{\partial{}}{\partial{\theta_{d}}}\frac{1}{|\nabla{\phi}|})+\frac{3}{2\sqrt{2}}$
$\displaystyle=$
$\displaystyle-\frac{\mu_{f}U_{cl}}{\sigma}\cot{(\theta_{d})}+\frac{3}{2\sqrt{2}}(1+\frac{\cos({\theta_{A/R})}}{\sin{(\theta_{d})}}\frac{\partial{}}{\partial{\theta_{d}}}\frac{1}{\sqrt{2}\eta|\nabla{\phi}|})$
$\displaystyle=$
$\displaystyle-\frac{\mu_{f}U_{cl}}{\sigma}\cot{(\theta_{d})}+\frac{3}{2\sqrt{2}}(1+C_{\theta_{d}})\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad$
(A6)

The assumption that the $\phi=$constant profiles remain parallel (planar
interface) after infinitesimal variation in $\theta_{D}$ and reaching
equilibrium can simplify the Eq. A6. The gradient of phase-field variable
$|\nabla{\phi}|$ for a planar interface in equilibrium is given by

$\displaystyle|\nabla{\phi}|_{\phi=0}=\frac{1-\phi^{2}}{\sqrt{2}\eta}|_{\phi=0}=\frac{1}{\sqrt{2}\eta}$
(A7)

with this assumption
$\frac{\partial{}}{\partial{\theta_{d}}}\frac{1}{\sqrt{2}\eta|\nabla{\phi}|}$
should be small and negligible in comparison to other terms, and the Eq. A6
can be reduced to

$\displaystyle\frac{\mu_{f}}{\sigma}\frac{\partial{U_{cl}}}{\partial{\theta_{d}}}=-\frac{\mu_{f}U_{cl}}{\sigma}\cot{(\theta_{d})}+\frac{3}{2\sqrt{2}}$
(A8)

Eq. A8 anticipates a large variation in the velocity with a change in
$\theta_{d}$ due to the dependency on $\cot{}$ function over superhydrophobic
surfaces since these surfaces usually own large contact angles.

## References

  * Abels et al. [2012] H. Abels, H. Garcke, and G. Grün.  Thermodynamically consistent, frame indifferent diffuse interface models for incompressible two-phase flows with different densities.  _Math. Models Methods Appl. Sci._ , 22(03):1150013, 2012\.
  * Aggarwal and Sarkar [2007] N. Aggarwal and K. Sarkar.  Deformation and breakup of a viscoelastic drop in a Newtonian matrix under steady shear.  _J. Fluid Mech._ , 584:1–21, 2007.
  * Amberg [2022] G. Amberg.  Detailed modelling of contact line motion in oscillatory wetting.  _NPJ Microgr._ , 8(1):1, 2022.
  * Ba et al. [2013] Y. Ba, H. Liu, J. Sun, and R. Zheng.  Color-gradient lattice Boltzmann model for simulating droplet motion with contact-angle hysteresis.  _Phys. Rev. E_ , 88(4):043306, 2013.
  * Balay et al. [2019] S. Balay, S. Abhyankar, M. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, A. Dener, V. Eijkhout, W. Gropp, and D. Karpeyev.  PETSc users manual.  _Argonne National Laboratory_ , 2019.
  * Batchelor et al. [2003] G. Batchelor, H. Moffatt, and M. e. Worster.  _Perspectives in fluid dynamics: a collective introduction to current research_.  Cambridge University Press, 2003.  ISBN 978-0521780612.
  * Bazesefidpar et al. [2022] K. Bazesefidpar, L. Brandt, and O. Tammisola.  A dual resolution phase‐field solver for wetting of viscoelastic droplets.  _Internat. J. Numer. Methods Fluids_ , 94(9):1517–1541, 2022.
  * Bird et al. [1987a] R. Bird, R. Armstrong, and O. Hassager.  _Dynamics of polymeric liquids. Vol. 1: Fluid mechanics_.  Wiley, 1987a.
  * Bird et al. [1987b] R. Bird, C. Curtiss, R. Armstrong, and O. Hassager.  _Dynamics of polymeric liquids, volume 2: Kinetic theory_.  Wiley, 1987b.
  * Carlson et al. [2009] A. Carlson, M. Do-Quang, and G. Amberg.  Modeling of dynamic wetting far from equilibrium.  _Phys. Fluids_ , 21(12):121701, 2009.
  * Carlson et al. [2011] A. Carlson, M. Do-Quang, and G. Amberg.  Dissipation in rapid dynamic wetting.  _J. Fluid Mech._ , 682:213–240, 2011.
  * Cheng and Rodak [2005] Y. Cheng and D. Rodak.  Is the lotus leaf superhydrophobic?  _Appl. Phys. Lett._ , 86(14):144101, 2005.
  * Dimitrakopoulos [2007] P. Dimitrakopoulos.  Deformation of a droplet adhering to a solid surface in shear flow: onset of interfacial sliding.  _J. Fluid Mech._ , 580:451–466, 2007.
  * Dimitrakopoulos and Higdon [1997] P. Dimitrakopoulos and J. Higdon.  Displacement of fluid droplets from solid surfaces in low-Reynolds-number shear flows.  _J. Fluid Mech._ , 336:351–378, 1997.
  * Ding and Spelt [2008] H. Ding and P. Spelt.  Onset of motion of a three-dimensional droplet on a wall in shear flow at moderate Reynolds numbers.  _J. Fluid Mech._ , 599:341–362, 2008.
  * Ding et al. [2010] H. Ding, M. Gilani, and P. Spelt.  Sliding, pinch-off and detachment of a droplet on a wall in shear flow.  _J. Fluid Mech._ , 644:217–244, 2010.
  * Dupont and Legendre [2010] J. Dupont and D. Legendre.  Numerical simulation of static and sliding drop with contact angle hysteresis.  _J. Comput. Phys._ , 229(7):2453–2478, 2010.
  * Gao and Wang [2014] M. Gao and X. Wang.  An efficient scheme for a phase field model for the moving contact line problem with variable density and viscosity.  _J. Comput. Phys._ , 272:704–718, 2014.
  * Guillén-González and Tierra [2013] F. Guillén-González and G. Tierra.  On linear schemes for a Cahn–Hilliard diffuse interface model.  _J. Comput. Phys._ , 234:140–171, 2013.
  * Jacqmin [2000] D. Jacqmin.  Contact-line dynamics of a diffuse fluid interface.  _J. Fluid Mech._ , 402:57–88, 2000.
  * Jiang et al. [2022] M. Jiang, Z. Zhang, and J. Zhao.  Improving the accuracy and consistency of the scalar auxiliary variable (SAV) method with relaxation.  _J. Comput. Phys._ , 456:110954, 2022.
  * Joanny and Robbins [1990] J. Joanny and M. Robbins.  Motion of a contact line on a heterogeneous surface.  _J. Chem. Phys._ , 92(5):3206–3212, 1990.
  * Lafuma and Quéré [2011] A. Lafuma and D. Quéré.  Slippery pre-suffused surfaces.  _Eur. Phys._ , 96(5):56001, 2011.
  * Lee et al. [2007] J. Lee, R. Dylla-Spears, N. Teclemariam, and S. Muller.  Microfluidic four-roll mill for all flow types.  _Appl. Surf. Sci._ , 90(7), 2007.
  * Liu et al. [2015] H. Liu, Y. Ju, N. Wang, G. Xi, and Y. Zhang.  Lattice Boltzmann modeling of contact angle and its hysteresis in two-phase flow with large viscosity difference.  _Phys. Rev. E_ , 92(3):033306, 2015.
  * Neinhuis and Barthlott [1997] C. Neinhuis and W. Barthlott.  Characterization and distribution of water-repellent, self-cleaning plant surfaces.  _Ann. Bot._ , 79(6):667–677, 1997.
  * Ramaswamy and Leal [1999] S. Ramaswamy and L. Leal.  The deformation of a viscoelastic drop subjected to steady uniaxial extensional flow of a Newtonian fluid.  _J. Non-Newton. Fluid Mech._ , 85(2-3):127–163, 1999.
  * Sartori et al. [2022] P. Sartori, D. Ferraro, M. Dassie, A. Meggiolaro, D. Filippi, A. Zaltron, M. Pierno, and G. Mistura.  Oscillatory motion of viscoelastic drops on slippery lubricated surfaces.  _Commun. Phys._ , 5(1):81, 2022.
  * Schleizer and Bonnecaze [1999] A. Schleizer and R. Bonnecaze.  Displacement of a two-dimensional immiscible droplet adhering to a wall in shear and pressure-driven flows.  _J. Fluid Mech._ , 383:29–54, 1999.
  * Shen and Yang [2010] J. Shen and X. Yang.  Numerical approximations of Allen-Cahn and Cahn-Hilliard equations.  _Discrete Contin Dyn Syst_ , 28(4):1669, 2010\.
  * Shen et al. [2018] J. Shen, J. Xu, and J. Yang.  The scalar auxiliary variable (SAV) approach for gradient flows.  _J. Comput. Phys._ , 353:407–416, 2018.
  * Shen et al. [2019] J. Shen, J. Xu, and J. Yang.  A new class of efficient and robust energy stable schemes for gradient flows.  _SIAM Rev._ , 61(3):474–506, 2019.
  * Shin et al. [2018] S. Shin, J. Chergui, and D. Juric.  Direct simulation of multiphase flows with modeling of dynamic interface contact angle.  _Theor. Comput. Fluid Dyn._ , 32:655–687, 2018.
  * Shirtcliffe et al. [2010] N. Shirtcliffe, G. McHale, S. Atherton, and M. Newton.  An introduction to superhydrophobicity.  _Adv. Colloid Interface Sci._ , 161(1-2):124–138, 2010\.
  * Spelt [2005] P. Spelt.  A level-set approach for simulations of flows with multiple moving contact lines with hysteresis.  _J. Comput. Phys._ , 207(2):389–404, 2005.
  * Spelt [2006] P. Spelt.  Shear flow past two-dimensional droplets pinned or moving on an adhering channel wall at moderate Reynolds numbers: a numerical study.  _J. Fluid Mech._ , 561:439–463, 2006.
  * Varagnolo et al. [2017] S. Varagnolo, D. Filippi, G. Mistura, M. Pierno, and M. Sbragaglia.  Stretching of viscoelastic drops in steady sliding.  _Soft Matter_ , 13(17):3116–3124, 2017.
  * Wang et al. [2023] N. Wang, S. Li, L. Shi, X. Yuan, and H. Liu.  Viscoelastic effects on the deformation and breakup of a droplet on a solid wall in Couette flow.  _J. Fluid Mech._ , 963:A18, 2023.
  * Wang et al. [2015] Y. Wang, D. Minh, and G. Amberg.  Dynamic wetting of viscoelastic droplets.  _Phys. Rev. E_ , 92(4):043002, 2015.
  * Wang et al. [2017] Y. Wang, M. Do-Quang, and G. Amberg.  Impact of viscoelastic droplets.  _J. Non-Newton. Fluid Mech._ , 243:38–46, 2017.
  * Watson et al. [2015] G. Watson, D. Green, L. Schwarzkopf, X. Li, B. Cribb, S. Myhra, and J. Watson.  A gecko skin micro/nano structure–a low adhesion, superhydrophobic, anti-wetting, self-cleaning, biocompatible, antibacterial surface.  _Acta Biomater._ , 21:109–122, 2015.
  * Xu et al. [2018] H. Xu, A. Clarke, J. Rothstein, and R. Poole.  Viscoelastic drops moving on hydrophilic and superhydrophobic surfaces.  _J. Colloid Interface Sci._ , 513:53–61, 2018.
  * Yada et al. [2023] S. Yada, K. Bazesefidpar, O. Tammisola, G. Amberg, and S. Bagheri.  Rapid wetting of shear-thinning fluids.  _Phys. Rev. Fluids_ , 8(4):043302, 2023.
  * Yang et al. [2019] Z. Yang, L. Lin, and S. Dong.  A family of second-order energy-stable schemes for cahn–hilliard type equations.  _J. Comput. Phys._ , 383:24–54, 2019.
  * Yue [2020] P. Yue.  Thermodynamically consistent phase-field modelling of contact angle hysteresis.  _J. Fluid Mech._ , 899:A15, 2020.
  * Yue and Feng [2011] P. Yue and J. Feng.  Wall energy relaxation in the cahn–hilliard model for moving contact lines.  _Phys. Fluids_ , 23(1), 2011.
  * Yue and Feng [2012] P. Yue and J. Feng.  Phase-field simulations of dynamic wetting of viscoelastic fluids.  _J. Non-Newton. Fluid Mech._ , 189:8–13, 2012.
  * Yue et al. [2004] P. Yue, J. Feng, C. Liu, and J. Shen.  A diffuse-interface method for simulating two-phase flows of complex fluids.  _J. Fluid Mech._ , 515:293–317, 2004.
  * Yue et al. [2005a] P. Yue, J. Feng, C. Liu, and J. Shen.  Diffuse-interface simulations of drop coalescence and retraction in viscoelastic fluids.  _J. Non-Newton. Fluid Mech._ , 129(3):163–176, 2005a.
  * Yue et al. [2005b] P. Yue, J. Feng, C. Liu, and J. Shen.  Viscoelastic effects on drop deformation in steady shear.  _J. Fluid Mech._ , 540:427–437, 2005b.
  * Yue et al. [2010] P. Yue, C. Zhou, and Feng.  Sharp-interface limit of the Cahn–Hilliard model for moving contact lines.  _J. Fluid Mech._ , 645:279–294, 2010.
  * Zhang and Yue [2020] J. Zhang and P. Yue.  A level-set method for moving contact lines with contact angle hysteresis.  _J. Comput. Phys._ , 418:109636, 2020.
  * Zhang and Shen [2022] Y. Zhang and J. Shen.  A generalized SAV approach with relaxation for dissipative systems.  _J. Comput. Phys._ , 464:111311, 2022.