# Decreasing contact angles at accelerating three-phase moving contact lines.

Carlos A. Galeano-Rios1<EMAIL_ADDRESS>1School of Mathematics,
University of Birmingham, Birmingham, B15 2TT, United Kingdom

###### Abstract

We consider a liquid drop placed on a smooth homogeneous solid substrate as it
spreads from rest to its eventual equilibrium state. The problem is studied
numerically in the framework of a model where the contact angle formed by the
drop’s free surface with the substrate is not prescribed as a function of the
contact-line speed and has to be found as part of the solution. It is shown
that in this spontaneous spreading, as the drop starts moving and the contact
line accelerates, the dynamic contact angle decreases with the contact-line
speed, which is in line with what experimental observations on spontaneous
wetting report, though it is contrary to what is assumed in most models aimed
at describing experiments on the “forced” spreading. Nontrivial aspects of the
implementation of the model in a finite-element-based algorithm are discussed
in detail.

###### keywords:

Dynamic wetting, moving contact lines, contact angle.

##  1 Introduction

The relationship between the dynamic contact angle and the contact-line speed
in dynamic wetting has been intensively studied both experimentally (e.g.
(Hoffman, 1975; Ngan & Dussan, 1982; Seebergh & Berg, 1992; Šikalo et al.,
2005)) and theoretically (e.g. Ranabothu et al. (2005); Ren & E (2007);
Legendre & Maglio (2013)). In theoretical studies, almost invariably the
contact angle is assumed to be a function of the contact-line speed, which
assumes the contact angle increases as the speed increases. This function,
assumed the same for all flows, is used as a fixed input into the models.

However, recent experiments (Mohammad Karim et al., 2016) have shown that the
contact angle-vs-speed relationship is quantitatively different for the so-
called “spontaneous” and “forced” contact-line motion, and the discrepancy
goes beyond the experimental error. Furthermore, experiments by Bayer &
Megaridis (2006) also show that the dynamic contact angle can decrease as the
contact-line speed increases with the authors coming to the same conclusion
that no one-to-one map from the contact-line velocity to the dynamic contact
angle can account for all types of flows and even for different regimes of
qualitatively the same flow. The situation has been qualitatively summarised
in a recent review (Shikhmurzaev, 2020).

The direct measurements of the dependence of the dynamic contact angle on
contact line velocity (Mohammad Karim et al., 2016; Bayer & Megaridis, 2006),
as well as earlier experiments on the flow dependence of the contact angle (De
Ruijter et al., 1999) and their analysis, (Davies et al., 2006) indicate that
the dynamic contact angle should be an outcome, part of the solution, rather
than an input in the form of a given function into the model. Such sensitivity
of the contact angle and hence the entire free surface shape to the flow make
any realistic mathematical model much more complex than the simpler ‘slip
models’ explored so far (see Shikhmurzaev (2020) for a review) and its
adequate numerical implementation a formidable problem in its own right.

The only known model to date that treats the contact angle as an outcome
rather than an input is the model introduced in Shikhmurzaev (1993) and
further developed in Shikhmurzaev (2007); Sprittles & Shikhmurzaev (2012b),
where dynamic wetting is described as a process by which a fresh liquid-solid
interface is continuously created. Near a moving contact line this interface
is yet out of equilibrium and the out-of-equilibrium surface tensions acting
on the contact line “negotiate”, via the Young equation, the value of the
contact angle. This contact angle thus depends on the stage of the process,
and on the flow, and is, in principle, different from the static (equilibrium)
one. Furthermore, it has been shown that for the capillary rise flow starting
from rest the dynamic contact angle can go down from $90$ as the contact line
accelerates (Sprittles & Shikhmurzaev, 2012b).

The main challenge in the present work consists in fully implementing the
interface formation model for cases in which the contact angle can vary in a
wide range, in particular, for flows that involve obtuse contact angles. The
issue at hand is then that, as shown in Sprittles & Shikhmurzaev (2011a), in
all models generalising the Navier slip condition, for obtuse contact angles,
a parasitic eigensolution appears, and this eigensolution gives rise to
numerical artefacts in the pressure distribution, namely to sharp mesh-
dependent spikes with opposite signs near the contact line. An additional
difficulty is that the physics of the problem introduces disparate length
scales so that accurate numerical implementation of the model, even without
addressing the artefact, poses a difficult problem.

In the present work, we address the problem of the implementation of the full
interface formation model in the entire range of variation of the dynamic
contact angle, in particular addressing the numerical artefact mentioned
above, and we show that, for the spreading of a two-dimensional drop from rest
over a horizontal solid substrate, there exists a stage of the process where
the dynamic contact angle decreases as the contact line accelerates. The
structure of the paper is as follows. In section 2, we summarise the interface
formation equations and discuss the split- domain formulation introduced for
its numerical implementation, used to control the parasitic eigensolution in
the vicinity of the contact line. Section 3 discusses the way the split-domain
formulation is integrated into the weak form of the equations, and covers
important technical aspects of the numerical implementation, including the
improvements to previous implementations of the interface formation model. In
section 4, the results of our modelling are discussed, and in section 5 we
give a brief review of the implications of these findings and analyse some
possible directions of further research.

##  2 Problem formulation

We consider the spreading of a 2D drop of an incompressible Newtonian fluid of
density $\rho$ and kinematic viscosity $\nu$ over a flat homogeneous solid
substrate. If a drop, initially placed on a substrate, is gently pressed
against the latter, e.g. by gravity, so as to form a small area of contact, it
will begin to spread until the contact angle formed between its free surface
and the solid boundary reaches its equilibrium value.

The initial configuration of the droplet $\Omega(t=0)$ is that of a circular
arc with a contact angle of $160^{o}$ (see the schematic representation in
figure 1). The solid occupies the bottom half of the plane and the remaining
space is filled by dynamically passive gas with constant pressure $p_{g}$. The
motion of the droplet starts form rest and fluid flow in the droplet is
subject to gravitational acceleration $\boldsymbol{g}$ and capillary effects.

Figure 1:  Schematics of initial configuration of the droplet.

We introduce Cartesian coordinates $(x,z)$ with the $x$ axis coinciding with
the solid surface and the $z$ axis given by the symmetry line in the 2D
droplet, pointing away from the solid. We consider the problem in
dimensionless form, using the droplet radius $L\coloneqq R$ as unit length,
$U\coloneqq\sigma_{1,e}/\left(\rho\nu\right)$ as unit speed (where
$\sigma_{1,e}$ is the equilibrium surface tension of the free surface),
$P\coloneqq\rho\nu U/R$ as the unit stress.

The flow is governed by the Navier-Stokes equation for incompressible flow,
i.e.

$\Rey\left[\frac{\partial\boldsymbol{u}}{\partial
t}+\boldsymbol{u}\cdot\nabla\boldsymbol{u}\right]=\nabla\cdot\mathsfbi{P}+S\kern-0.39993ptt\,\boldsymbol{g}$
(1)

and

$\nabla\cdot\boldsymbol{u}=0,$ (2)

where $\boldsymbol{u}$ is the fluid velocity,
$\boldsymbol{g}=-g\hat{\boldsymbol{z}}$, $\Rey=\nu U/L^{2}$ and
$S\kern-0.39993ptt\,=gL^{2}/(\nu U)$, with $U$ and $L$ being the velocity and
length scales. Moreover, $\mathsfbi{P}$ is the stress tensor, defined as

$\mathsfbi{P}=-p\mathsfbi{I}+[\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T}],$
(3)

where $p$ is the pressure and $\mathsfbi{I}$ is the metric tensor.

###  2.1 Boundary conditions

The fluid domain $\Omega(t)$ is bounded by the free surface $S_{1}$, the
liquid-solid surface $S_{2}$ and the axis of symmetry $S_{3}$ (figure 1). As
required by the interface formation model, we introduce surface velocity and
surface density variables $\boldsymbol{v}^{s}_{i}$ and $\rho^{s}_{i}$,
respectively; which correspond to boundary $S_{i}$ with $i=1,2$.

####  2.1.1 The free surface

The tangential and normal dynamic boundary conditions (DBC) and the kinematic
boundary condition (KBC) are respectively given by

$C\kern-0.39993pta\,\boldsymbol{n}_{1}\cdot\left[\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T}\right]\cdot(\mathsfbi{I}-\boldsymbol{n}_{1}\boldsymbol{n}_{1})+\nabla\sigma_{1}=0,$
(4)
$C\kern-0.39993pta\,\boldsymbol{n}_{1}\cdot\left\\{p_{g}-p+\boldsymbol{n}_{1}\cdot\left[\nabla\boldsymbol{u}+(\nabla\boldsymbol{u})^{T}\right]\right\\}=\sigma_{1}\nabla\cdot\boldsymbol{n}_{1}$
(5)

and

$\left(\boldsymbol{c}-\boldsymbol{v}^{s}_{1}\right)\cdot\boldsymbol{n}_{1}=0;$
(6)

where $\boldsymbol{c}=(u^{c},w^{c})$ is the velocity of the surface
coordinates, $\boldsymbol{n}_{1}$ is the normal to the free-surface that
points into $\Omega$ and $C\kern-0.39993pta\,=\rho\nu U/\sigma_{1,e}$, with
$\sigma_{1,e}$ being the equilibrium surface tension coefficient of the free
surface, i.e. the surface tension coefficient that is measured when the
droplet is static. Following Sprittles & Shikhmurzaev (2012b), the two dynamic
boundary conditions can be combined into a single vector boundary condition,
given by

$\left(p_{g}+\mathsfbi{P}\right)\cdot\boldsymbol{n}_{1}=-\frac{\nabla^{s}\cdot\left[\sigma_{1}(\mathsfbi{I}-\boldsymbol{n}_{1}\boldsymbol{n}_{1})\right]}{C\kern-0.39993pta\,}.$
(7)

The free-surface and the liquid-solid surface are treated as 2-dimensional
phases that exchange mass with the 3-dimensional bulk phase. A surface density
is introduced for each interface, which reflects the fact that interfaces are
transitional regions with average properties that are different from those of
the bulk. An equilibrium surface density exists for each surface, and the
deviation of surface tension from its nominal equilibrium value depends
(linearly) on how much the surface density deviates from its equilibrium
condition. The model also incorporates mass-transport effects along the
interfaces and relates surface tension gradients to the slip between bulk and
interface phases. Consequently, five further equations are required for the
free surface; namely:

$\left(\mathsfbi{I}-\boldsymbol{n}_{1}\boldsymbol{n}_{1}\right)\cdot\left(\boldsymbol{v}^{s}_{1}-\boldsymbol{u}\right)=\frac{1+4\bar{\alpha}_{g}\bar{\beta_{g}}}{4\bar{\beta_{g}}}\nabla\sigma_{1},$
(8)
$\left(\boldsymbol{u}-\boldsymbol{v}^{s}_{1}\right)\cdot\boldsymbol{n}_{1}=Q_{g}\left(\rho^{s}_{1}-\rho^{s}_{1,e}\right),$
(9) $\sigma_{1}=\lambda_{g}\left(1-\rho^{s}_{1}\right),$ (10)
$\epsilon_{g}\left[\frac{\partial\rho^{s}_{1}}{\partial
t}+\nabla\cdot\left(\rho^{s}_{1}\boldsymbol{v}^{s}_{1}\right)\right]=\rho^{s}_{1,e}-\rho^{s}_{1},$
(11)

and

$\rho^{s}_{1}=\rho^{s}_{1,e},\qquad\text{at }t=0;$ (12)

where $\bar{\alpha}_{g}=\alpha_{g}\sigma_{1,e}/(UL)$,
$\bar{\beta}_{g}=\beta_{g}UL/\sigma_{1,e}$, $Q_{g}=\rho^{s}_{(0)}/(\rho
U\tau_{g})$, $\epsilon_{g}=\tau_{g}U/L$ and
$\lambda_{g}=\gamma_{g}\rho^{s}_{(0)}/\sigma_{1,e}$; with $\alpha_{g}$,
$\sigma_{1,e}$, $\beta_{g}$, $\rho^{s}_{(0)}$, $\rho^{s}_{1,e}$, $\tau_{g}$
and $\gamma_{g}$ being physical properties of the surface which are specified
in section 3.5.

Equations (8-10) respectively correspond to the slip condition between the
1-dimensional (1D) free-surface phase and the bulk phase, the equation for
mass exchange between bulk and surface phase, the equation of state that
relates surface tension and surface density, and the equation for mass
transport along the interface.

####  2.1.2 The liquid-solid surface

On the liquid-solid surface, the impermeability condition and the generalised
Navier-slip condition are respectively given by

$\left(\boldsymbol{v}^{s}_{2}-\boldsymbol{u}_{s}\right)\cdot\boldsymbol{n}_{2}=0,$
(13)

and

${\boldsymbol{n}_{2}\cdot\mathsfbi{P}\cdot(\mathsfbi{I}-\boldsymbol{n}_{2}\boldsymbol{n}_{2})+\frac{\nabla\sigma^{2}}{2C\kern-0.39993pta\,}}\\\
={\bar{\beta}_{s}(\boldsymbol{u}-\boldsymbol{u}_{s})\cdot(\mathsfbi{I}-\boldsymbol{n}_{2}\boldsymbol{n}_{2})},$
(14)

where $\boldsymbol{n}_{2}$ is the normal to the liquid-solid interface that
points into $\Omega$ and $\boldsymbol{u}_{s}$ is the velocity of the solid
surface. Here, we include a generic $\boldsymbol{u}_{s}$ for added generality,
despite all cases here considered satisfying $\boldsymbol{u}_{s}=0$. Moreover,
$\bar{\beta}_{s}=\beta_{s}L/(\rho\nu U)$, with $\beta_{s}$ being a surface
property specified in section 3.5.

As was the case in the free surface, here we also need further equations,
these are

$\left[\boldsymbol{v}^{s}_{2}-\frac{1}{2}\left(\boldsymbol{u}+\boldsymbol{u}_{s}\right)\right]\cdot\left(\mathsfbi{I}-\boldsymbol{n}_{2}\boldsymbol{n}_{2}\right)=\bar{\alpha}_{s}\nabla\sigma_{2},$
(15)
$\left(\boldsymbol{u}-\boldsymbol{v}^{s}_{2}\right)\cdot\boldsymbol{n}_{2}=Q_{s}\left(\rho^{s}_{2}-\rho^{s}_{2,e}\right),$
(16) $\sigma_{2}=\lambda_{s}\left(1-\rho^{s}_{2}\right),$ (17)
$\epsilon_{s}\left[\frac{\partial\rho^{s}_{2}}{\partial
t}+\nabla\cdot\left(\rho^{s}_{2}\boldsymbol{v}^{s}_{2}\right)\right]=\rho^{s}_{2,e}-\rho^{s}_{2},$
(18)

and

$\rho^{s}_{2}=1-\frac{\sigma_{g-s}-\sigma_{1,e}\cos(\theta_{c})}{\gamma_{s}},\qquad\text{at
}t=0;$ (19)

where $\bar{\alpha}_{s}=\alpha_{s}\sigma_{1,e}/(UL)$,
$Q_{s}=\rho^{s}_{(0)}/(\rho U\tau_{s})$, $\epsilon_{s}=\tau_{s}U/L$ and
$\lambda_{s}=\gamma_{s}\rho^{s}_{(0)}/\sigma_{1,e}$; with $\alpha_{s}$,
$\sigma_{1,e}$, $\rho^{s}_{(0)}$, $\rho^{s}_{2,e}$, $\tau_{s}$ and
$\gamma_{s}$ being surface properties specified in section 3.5. The initial
condition for equation (18), given by equation (19), is obtained combining
equations (12), (10), (17) and the Young equation (see equation 22 below).

####  2.1.3 The symmetry axis

On the axis of symmetry the flow is subject to the usual impermeability and
no-tangential-stress conditions, given respectively by

$\boldsymbol{u}\cdot\boldsymbol{n}_{3}=0$ (20)

and

$\boldsymbol{n}_{3}\cdot\mathsfbi{P}\cdot\left(\mathsfbi{I}-\boldsymbol{n}_{3}\boldsymbol{n}_{3}\right)=0,$
(21)

where $\boldsymbol{n}_{3}$ is the normal to the axis of symmetry that points
into $\Omega$.

####  2.1.4 Contact line conditions

On the free-surface and liquid-solid surface, we have boundary conditions
which themselves are partial differential equations (PDEs) that require
boundary conditions of their own. In particular, at the contact line, where
these two surfaces meet, they must satisfy Young’s equation, namely

$\sigma_{1}\cos(\theta_{c})+\sigma_{2}=\sigma_{g-s},$ (22)

where $\sigma_{g-s}$ is the solid-gas surface tension.

Furthermore, mass flux at the contact line must satisfy

$\rho^{s}_{1}\left(\boldsymbol{v}^{s}_{1}-\boldsymbol{U}_{c}\right)\cdot\boldsymbol{m}_{1}+\rho^{s}_{2}\left(\boldsymbol{v}^{s}_{2}-\boldsymbol{U}_{c}\right)\cdot\boldsymbol{m}_{2}=0,$
(23)

where $\boldsymbol{m}_{i}$, defined only on the boundary of $S_{i}$, is the
unit tangent to boundary $S_{i}$ that points inward; and $\boldsymbol{U}_{c}$
is the velocity of the contact line.

####  2.1.5 Conditions at the ends of $S_{3}$

At the apex of the droplet, i.e. where the free surface meets the axis of
symmetry, we have

$\theta_{a}=\frac{\pi}{2},$ (24)

where $\theta_{a}$ is the angle between the axis of symmetry and the free
surface. Moreover, we have

$\boldsymbol{v}^{s}_{1}\cdot\boldsymbol{n}_{3}=0,\qquad\text{at }x=0.$ (25)

We note that the combination of equations (8), (20), (24) and (25) implies
that, at the apex, we have

$\nabla\sigma_{1}=0,\qquad\text{at }x=0,$ (26)

since, at this point, $\boldsymbol{n}_{3}$ is parallel to
$\boldsymbol{t}_{1}$, the unit tangent to the free surface.

Similarly, where the liquid-solid surface meets the axis of symmetry, we also
have

$\boldsymbol{v}^{s}_{2}\cdot\boldsymbol{n}_{3}=0,\qquad\text{at }x=0,$ (27)

which, by an argument entirely analogue to the one used above, implies

$\nabla\sigma_{2}=0,\qquad\text{at }x=0.$ (28)

Equations (1)-(3), subject to boundary conditions (6)-(28), form the closed
system we need to solve.

##  3 Split domain, weak formulation and numerical implementation

We proceed to implement an approximate solution to the system of equations
presented above. We first state the equations in an arbitrary Lagrangian-
Eulerian reference system, we then introduce a split-domain formulation, which
is necessary to solve flows in which the contact angle can be obtuse. Finally,
we introduce the weak formulation of the problem and we proceed to apply the
finite element method using Galerkin’s method.

###  3.1 Arbitrary Lagrangian-Eulerian reference system

The free-boundary problem at hand is naturally suited for the use of an
arbitrary Lagragian-Eulerian (ALE) reference system, which requires that those
equations that involve temporal derivatives, i.e. (1), (11) and (18), be re-
written as

$\Rey\left[\frac{\partial\boldsymbol{u}}{\partial
t}+\left(\boldsymbol{u}-\boldsymbol{c}\right)\cdot\nabla\boldsymbol{u}\right]-\nabla\cdot\mathsfbi{P}+S\kern-0.39993ptt\,\boldsymbol{e}_{z}=0,$
(29) $\epsilon_{g}\left\\{\frac{\partial\rho^{s}_{1}}{\partial
t}+\rho^{s}_{1}\nabla\cdot\boldsymbol{c}+\nabla\cdot\left[\rho^{s}_{1}\left(\boldsymbol{v}^{s}_{1}-\boldsymbol{c}\right)\right]\right\\}=\rho^{s}_{1,e}-\rho^{s}_{1},$
(30)

and

$\epsilon_{s}\left\\{\frac{\rho^{s}_{2}}{\partial
t}+\rho^{s}_{2}\nabla\cdot\boldsymbol{c}+\nabla\cdot\left[\rho^{s}_{2}\left(\boldsymbol{v}^{s}_{2}-\boldsymbol{c}\right)\right]\right\\}=\rho^{s}_{2,e}-\rho^{s}_{2},$
(31)

where $\boldsymbol{c}$ is the velocity of the ALE coordinates. We note that,
in the latter two equations above, we have used the identity
$\nabla^{s}(f\boldsymbol{h})=\boldsymbol{h}\cdot\nabla^{s}f+f\nabla^{s}\cdot\boldsymbol{h}$,
to express the equations in a more convenient format.

It is also worth noting that, in this formulation, we can re-write (23) as

$\rho^{s}_{1}\left(\boldsymbol{v}^{s}_{1}-\boldsymbol{c}\right)\cdot\boldsymbol{m}_{1}+\rho^{s}_{2}\left(\boldsymbol{v}^{s}_{2}-\boldsymbol{c}\right)\cdot\boldsymbol{m}_{2}=0.$
(32)

###  3.2 The split-domain formulation

When the contact angle is acute, the equations listed above can be solved
numerically, for instance, using the finite element method (FEM), as done in
Sprittles & Shikhmurzaev (2013). However, as was shown by Sprittles &
Shikhmurzaev (2011b), when the contact angle is obtuse, a straightforward
computational implementation of this set of equations leads to nonphysical
solutions that dominate the flow around the contact line. These effects are
indeed observed with any of the standard FEM implementations for dynamic
wetting, regardless of the inclusion of the interface formation equations.

In Sprittles & Shikhmurzaev (2011a), the origin of such nonphysical behaviour
is found, and a solution to the problem is provided. The undesired behaviour
is caused by the numerical method picking up a member of a one-dimensional
space of solutions to the problem of Stokes’ flow in a wedge (subject to no-
tangential-stress conditions on its boundaries). The removal of these
nonphysical features is achieved by expressing the velocity field, in the
vicinity of the contact line, as a sum of one of this eigensolution (whose
amplitude is to be determined), and a second unknown component of the velocity
field. Effectively, this means that we remove the eigensolution contribution
when calculating the remaining components of the flow. The eigensolution
contribution is later put back, completing the correct physical picture.

We can thus write

$\boldsymbol{u}=\bar{\boldsymbol{u}}+A(t)\check{\boldsymbol{u}},$ (33)

where $\check{\boldsymbol{u}}$ is the solution to the problem

$\nabla\cdot\mathsfbi{P}=0,$ (34) $\nabla\cdot\check{\boldsymbol{u}}=0$ (35)

on the plane, which is subject to

$\boldsymbol{n}\cdot\mathsfbi{P}\cdot\left(\mathsfbi{I}-\boldsymbol{n}\boldsymbol{n}\right)=0,$
(36)

along the lines than bound a wedge of angle $\theta_{c}$.

The conditions above determine the eigensolution to be given by the stream
function

$\psi=B\left(\check{x}^{2}+\check{z}^{2}\right)^{\lambda/2}\sin\left(\lambda\left[\theta+\theta_{c}-\pi\right]\right),$
(37)

where $\lambda=\pi/\theta_{c}$, $B$ is an arbitrary factor (to be set as
$B=1$), $\check{x}=x-x_{c}$, and $\check{z}=z-z_{c}$, with $(x_{c},z_{c})$
being the coordinates of the contact line. We note that this is a solution on
the entire plane, not just the wedge, though boundary conditions are satisfied
at the boundary of the wedge. This fact is important, as it allows us to
evaluate the solution slightly outside the wedge, which can be necessary in
this implementation.

A factor $A(t)$, the “amplitude” of the eigensolution, multiplied by $\psi$,
becomes one more unknown in the system, and an additional equation is required
to determine it, namely

$\lim_{r\to 0}\frac{\partial p}{\partial\theta}=0,$ (38)

where $r=\sqrt{\check{x}^{2}+\check{z}^{2}}$,
$\theta=\tan^{-1}\left(-\check{z}/\check{x}\right)$, with $0<\theta<\pi$.

Figure 2:  Schematic representation of the split domain. Compatibility
conditions for the bulk equations are applied along the $S_{4}$ line, and
continuity is assumed for the interface formation variables on points $J$ and
$K$.

The introduction of the eigensolution is only necessary in the vicinity of the
contact line, and it is therefore convenient to split the $\Omega$ domain into
two sub-domains, $\Omega_{\text{n}}$ and $\Omega_{\text{f}}$, respectively
spanning the near-field of the contact line and the far-field (see figure 2).

In this split-domain formulation we solve two sets of PDEs; which, at the
separatrix of the respective domains ($S_{4}$), are subject to the following
compatibility conditions

$\boldsymbol{u}_{\text{f}}=\bar{\boldsymbol{u}}+A\check{\boldsymbol{u}},$ (39)

where $\boldsymbol{u}_{\text{f}}$ is the velocity as calculated in the far
field; and

$\mathsfbi{P}_{\text{f}}\cdot\boldsymbol{n}_{4}=\mathsfbi{P}_{\text{n}}\cdot\boldsymbol{n}_{4},$
(40)

where $\boldsymbol{n}_{4}$ points into the far field domain.

Equations for flow in the vicinity of an obtuse contact angle are obtained by
substituting (33), into the formulation given by (1)-(23). The equations on
$\Omega_{\text{f}}$ and $\Omega_{\text{n}}$ involve functions that are defined
in both sub-domains, such as $p$, $\boldsymbol{v}^{s}_{i}$, $\rho^{s}_{i}$ and
$\sigma_{i}$, for $i=1,2$. Naturally, these must be continuous across $S_{4}$.

It should be emphasised that the split-domain formulation is only necessary
while the contact angle is obtuse. However, in this problem (as well as in
many other dynamic wetting ones), the contact angle is free to fluctuate
between acute and obtuse; consequently, it is more convenient to always act as
if we are solving two sets of PDEs in their respective domains; the far field
problem being fixed and the near field problem being switched from the obtuse-
contact-angle formulation to the acute one, simply by forcing $A$ to be zero
and removing equation (38) when needed.

###  3.3 Weak formulation

Once the split-domain problem is clearly defined, we introduce its weak
formulation, in preparation for the use of the finite element method.

The split domain formulation re-shapes our original bulk PDEs in the following
form

$\displaystyle\chi_{\Omega_{\text{f}}}\mathcal{P}_{\text{f}}\left(\boldsymbol{u},p,\boldsymbol{u}^{s}_{1},\rho^{s}_{1},\sigma_{1},\boldsymbol{u}^{s}_{2},\rho^{s}_{2},\sigma_{2}\right)\
\ \ \ \ \ \ \ \ \ $
$\displaystyle+\chi_{\Omega_{\text{n}}}\mathcal{P}_{\text{n}}\left(\bar{\boldsymbol{u}},p,\boldsymbol{u}^{s}_{1},\rho^{s}_{1},\sigma_{1},\boldsymbol{u}^{s}_{2},\rho^{s}_{2},\sigma_{2}\right)=0,$
(41)

where $\mathcal{P}_{\text{f}}$ and $\mathcal{P}_{\text{n}}$ represent the
differential operators in the far-field and near-field bulk, respectively.
Consequently, the weak form can be written as

$\int\limits_{\Omega_{\text{f}}}{\phi_{i}\mathcal{P}_{\text{f}}}+\int\limits_{\Omega_{\text{n}}}{\phi_{i}\mathcal{P}_{\text{n}}}=0,$
(42)

where $\phi_{i}$ is an arbitrary element of our set of test functions and
$\Omega_{\text{f}}$ and $\Omega_{\text{n}}$ are the far-field and near-field
sub-domains, respectively.

It is important to highlight that the PDEs being solved on both sub-domains
are the same. That is, we are still solving a smooth PDE in the entire domain.
It is only the expression of the bulk velocity that changes, and the form in
which we write the equations changes as a consequence of this.

Similarly, for each of the surface PDEs in the interface formation model we
have

$\displaystyle\chi_{S_{j,\text{f}}}\mathcal{P}^{j}_{\text{f}}\left(\boldsymbol{u},p,\boldsymbol{u}^{s}_{j},\rho^{s}_{j},\sigma_{j}\right)\
\ \ \ \ \ \ \ \ \ $
$\displaystyle+\chi_{S_{j,\text{n}}}\mathcal{P}^{j}_{\text{n}}\left(\bar{\boldsymbol{u}},p,\boldsymbol{u}^{s}_{j},\rho^{s}_{j},\sigma_{j}\right)=0,$
(43)

with $j=1,2$. For a given element $\phi^{j}_{i}$ of our set of test functions
on the boundary, we have

$\int\limits_{S_{j,\text{f}}}{\phi^{j}_{i}\mathcal{P}^{j}_{\text{f}}}+\int\limits_{S_{j,\text{n}}}{\phi^{j}_{i}\mathcal{P}^{j}_{\text{n}}}=0.$
(44)

The form above implies that residuals associated to test functions that are
supported on the union of the near-field and far-field sub-domains will, in
general, include a contribution from each. The expressions for each residual
contribution in the far and near fields are given in Appendix A.

###  3.4 Numerical Implementation

We use Galerkin’s method with the FEM on the weak form of the split-domain
formulation of the problem in an Arbitrary Lagrangian-Eulerian system of
reference. We follow closely the methodology of Sprittles & Shikhmurzaev
(2012b), including their choice to use the spine method, the V6P3 taylor-Hood
triangular elements, and their meshing algorithm based on the bi-polar system
of coordinates, where the separatrix between near-field and far-field sub-
domains happens along one of the circular spines. However, the exact features
of our mesh are somewhat different, so as to optimise it for our applications.

Close to the contact line, the spines are separated using the criteria given
by Sprittles & Shikhmurzaev (2012b), i.e. the distance between consecutive
spines along the horizontal liquid-solid line increases in a geometric
sequence that uses the ratio recommended in Sprittles & Shikhmurzaev (2012b)
(ratio = 1.07). Next to this region, we use a very similar geometric criteria,
this time with a ratio of 1.1. We also introduce a geometric progression for
the separation between the spines along the horizontal axis in the vicinity of
the axis of symmetry, that is, we define an initial separation for the spines
closest to the vertical axis, and we place neighbouring spines with increasing
separation in a geometric sequence of ratio $1.2$. The remaining middle
section of the liquid-solid line is filled with equally spaced spines, the
spacing in this middle section is chosen as close as possible to the average
of the two spacings to its left and right. As the flow progresses, the length
of the liquid-solid line (which is considered a spine in the problem) changes,
and the separation between points of intersection of the circular spines with
the horizontal axis change in proportion.

Given the location of the contact line, the separation between the axis of
symmetry and its neighbouring spine, and the separation between the contact
line and its neighbouring circular spine; all we need to place every spine is
the number of spines in each of the sections where the geometric progression
for the spacing applies. The number of spines in the intermediate section will
be determined by the size of the spacing which is chosen as described above.
See figure 3 for an schematic representation.

Figure 3: Schematic representation of the spine spacing along the horizontal
axis. Bullets correspond to intersection between the spines and the liquid-
solid line, point $O$ corresponds to the axis of symmetry and point $C$ to the
contact line.

Another modification to the meshing algorithm presented on Sprittles &
Shikhmurzaev (2012b) is necessary in order to verify equation (38). The
numerical implementation of this limit is done by ensuring that the two
pressure nodes on each interface (one on $S_{1}$ and one on $S_{2}$) that are
at the closest (strictly positive) distance to the contact line, are located
on a circumference centred at the contact line. This is solved assigning the
centre of the circular spine with the two smallest radii to be given the
contact line itself. When the mesh is chosen in this way, the numerical
verification of equation (38) is simply done by equating the pressure at these
two nodes, as is recommended in Sprittles & Shikhmurzaev (2011a).

The intersection of each circular spine with the horizontal axis is enough to
define the radius of each spine (level set of the bipolar coordinate system),
as explained in Sprittles & Shikhmurzaev (2012b); since this is determined by
the initial distribution and the length of the liquid-solid line, the length
of all spines (including the liquid-solid line) is all that is needed to fully
determine the shape of the domain.

Figure 4: Detail of the mesh in the vicinity of the contact line.

Spines are numbered from $1$ to $n_{s}$, with the first spine being the
liquid-solid line and the subsequent ones listed in increasing order of
distance to the contact line. Spine $n_{s}$ is on the axis of symmetry.

The mesh in Sprittles & Shikhmurzaev (2012b) had a single element to have a
vertex on the contact line. This was perfectly justified in the cases there
considered, as they only ever dealt contact angles that are less than or equal
to $90^{o}$. However, in our case, when the angle can be close to $180^{o}$,
this convention would deform this element almost to the point of degeneracy.
We thus allow a generic number of elements at the contact line, three being
sufficient for all purposes here considered.

As in Sprittles & Shikhmurzaev (2012b), our elements are placed along circular
bands around the contact line,as can be seen in figure 4. Each band has nodes
along three spines. As shown in figure 4, our mesh has three triangular
elements in touch with the contact line the contact line. The number of
elements placed along the following band of elements is of seven and each
subsequent band includes two elements more than the number of elements in the
prior band. This pattern is maintained over the first portion of the region
where we use the factor of $1.07$ to determine the progression of the spacing
between spines (see the $\times$ markers next to the contact line in figure
3). The first band formed by the remaining spines in this first section
($\triangle$ markers) includes a single element more that the prior band, and
all subsequent bands in the entire mesh keep the exact same number of elements
as the prior band, as shown in dashed blue lines in figure 4.

The separatrix between near and far field is given by the last solid black
spine in figure 4, and in practice we use $38$ circular spines with the $1.07$
spacing ratio, with $14$ of those in the region next to the contact line,
$115$ spines in the region with the factor of $1.1$, and 21 in the region with
the factor of $1.2$. Moreover, for the configuration shown in the sections
below, we use a minimum length that $2^{12}$ times smaller than the minimum
recommended by Sprittles & Shikhmurzaev (2013).

Gaussian integration is used to numerically find the residuals and their
derivatives. We use an 8th order rule with 16 points based on permutation
stars on a triangle, following Zhang et al. (2009). The exact eigensolution is
sampled at the Gaussian integration points, as opposed to it being sampled at
the nodes and interpolated by the hat functions, following the advise given in
Sprittles & Shikhmurzaev (2011a). To capture the singular nature of the
pressure at the contact line, we follow Sprittles & Shikhmurzaev (2011b) and
use an interpolating function for pressure that is singular at that point.

Temporal derivatives are discretised using the implicit, second-order,
variable-step-size, backward differentiation formula (VS-BDF2) following
Celaya et al. (2014) and Denner (2014). The initial time step is obtained
using the implicit Euler method. The first time step size is taken to be the
characteristic time divided by a power of $2$, which is chosen to guarantee
that the change in angle is less than one degree in this step. Subsequent
time-steps are set fixing a maximum and minimum desired change in the contact
angle. If the prior time step violated either condition, the next time step is
chosen as double or half the prior one, accordingly. In the case discussed in
the sections below we used a maximum change in dynamic contact angle of
$1^{o}$ and a minimum change in contact angle of $0.1^{o}$, over the first
portion of the simulation; and when the Newton method started taking more than
five iterations to converge, we changed this numbers to $0.2^{o}$ and
$0.02^{o}$, respectively. To verify whether our choice of time-step is
adequate, we run the same simulation with half the chosen time-step until it
reaches maximum velocity. We then compare the resulting predictions for the
dynamic contact angle, and we accept the simulation if they differ by less
than $0.1^{o}$.

Five improvements, with respect to what was done in Sprittles & Shikhmurzaev
(2013) and Sprittles & Shikhmurzaev (2012a), to the weak form of the interface
formation residuals are introduced here. First, equations are re-cast so that
the surface tension functions do not strictly need to be differentiable;
secondly, certain terms in the momentum equations are exchanged for boundary
terms using Gauss’ theorem and the incompressibility condition (reducing the
computations needed); thirdly, residuals for equations (30) and (31) are
expressed in way that does not require the free-surface normal to be
differentiated (thus removing all need for the calculation of second
derivatives to find the residuals). More importantly, the fourth improvement
relates to the implementation of equation (32), which is done in a way that is
symmetric and more robust. Finally, the fifth improvement consists in
equations (9) and (16) being combined with (11) and (18) to produce more
numerically tractable residuals. The residual equations used in this work are
given in full detail in appendix A, where each of these improvements is shown.

Out of these improvements, we highlight the symmetric implementation of
equation (32). In prior implementations, this equation was used exclusively to
feed information from the contact line into the liquid-solid surface
(Sprittles & Shikhmurzaev, 2013, see equation 43); whereas, in the present
work, information can flow from one boundary to the other as dictated by the
flow.

The eigensolution is switched off (i.e. $A\equiv 0$) for contact angles lower
than $\pi/2$, which can happen halfway through a Newton method iteration, and
therefore the standard Newton algothim has to be adapted to select the correct
system of equations to solve.

It is extremely important to highlight that when using the eigensolution in
the near field, the angle that one should use to calculate variable $\lambda$
is the angle measured on the mesh $\theta_{c,m}$, i.e. the one given by the
solid surface and the tangent to the free-surface line element at the contact
line. This angle converges in the Newton method iterations to the value of
$\theta_{c}$, i.e. the variable that appears in the momentum and other
equations, but if $\theta_{c}$ is used to calculate $\lambda$ at the start of
the Newton method iterations, the method shows a tendency to escape
convergence. Convergence of these two angles is absolutely vital for a
consistent simulation. In fact, it is the criterion used in Sprittles &
Shikhmurzaev (2012b) to identify a sufficiently refined mesh, which was also
used here.

Care should be taken when calculating this angle and any other quantity that
depends on the tangents to the free surface. In particular, in the vicinity of
the contact line, elements are extremely small, so it is often convenient to
translate their coordinates to a point in the element before operating on
them. Moreover, when finding lengths in these elements, small quantities are
squared, then added, and the result is raised to the $1/2$ power; hence, it is
often convenient to normalise the numbers by a typical size of them before
squaring them. One can then find the lengths and, finally, re-scale them.

We note that the advice for the minimum length that needs to be resolved,
given in Sprittles & Shikhmurzaev (2013), appears to be sufficient only for
acute angles. The larger the angle, the more this length needs to be reduced.
Here we reduce the length by trial and error, until the result of the
simulations shows that the contact angle, as measured in from the mesh, and
the dynamic contact angle variable $\theta_{c}$ are within $0.2^{o}$ of each
other.

All entries in the Jacobian of the system of equations are calculated
analytically, with the only exception of the entries in the column that
correspond to the derivatives with respect to the length of the liquid-solid
line; for these entries, derivatives are calculated using the chain rule, with
the derivative of the coordinates of the nodes found numerically, while every
other derivative in the chain is found analytically.

All code was written in Matlab and is made available at
https://github.com/grcarlosa/IFM_repo. This includes the meshing algorithms,
construction of the residuals and the Jacobian of the system, as well as the
implementation of the modified Newton method. Extensive notes on these
derivations are available at the public repository in
https://github.com/grcarlosa/IFM_notes_repo.

###  3.5 Physical constants

Bulk fluid constants are taken from standard values for water at $25^{o}$C,
i.e.

$\rho=1\times 10^{3}\,\rm{kg/m}^{3}$ (45) $\nu=1\times
10^{-6}\,\rm{m}^{2}/\rm{s}$ (46)

and an equilibrium surface tension for the free surface of

$\sigma_{1,e}=7.286\times 10^{-2}\,N/m.$ (47)

Gravitational effects are not considered, i.e. we take $g=0\,\rm{m/s}^{2}$,
and we assume a constant atmospheric pressure of the gas phase, which we take
as $p_{g}=0\,\rm{N}/\rm{m}^{2}$.

The characteristic units are given by $L=R=0.5\,\mu$m ,
$U=\sigma_{1,e}/\left(\rho\nu\right)=72.86$ m$/$s and $P=\rho\nu
U/R=145,720\,$ N$/$m2. Interface formation constants are taken following the
recommendations of Sprittles & Shikhmurzaev (2013), Sprittles & Shikhmurzaev
(2012a). In particular, we assumed a width of the interfacial layer of

$\upsilon=1\times 10^{-9}\rm{m},$ (48)

also following (Bazhlekov & Chesters, 1996), and the relaxation times for the
interfaces are given by

$\tau_{s}=7\times 10^{-6}\frac{\rm{m}\,\rm{s}^{2}}{\rm{kg}}=7\times
10^{-12}\,\rm{s},$ (49)

with $\tau_{g}=\tau_{s}/100$.

From the assumptions above, it follows that

$\rho^{s}_{(0)}=\rho\upsilon=1\times 10^{-6}\,\rm{kg}/\rm{m}^{2},$ (50)
$\alpha_{g}=\alpha_{s}=\frac{\upsilon}{\rho\nu}=1\times
10^{-6}\,\rm{m}^{2}\rm{s}/kg,$ (51)
$\beta_{g}=\beta_{s}=\frac{\rho\nu}{\upsilon}=1\times
10^{6}\,\rm{kg}/(\rm{m}^{2}\rm{s}),$ (52)

and we take

$\rho^{s}_{1,e}=0.6\rho^{s}_{(0)}=6\times 10^{-7}\,\rm{kg}/\rm{m}^{2},$ (53)

which, together with equations (10), (47) and (50), implies

$\gamma_{g}=\frac{\sigma_{1,e}}{\rho^{s}_{(0)}-\rho^{s}_{1,e}}=1.82\times
10^{5}\,\rm{m}^{2}/\rm{s}^{2},$ (54)

and we take $\gamma_{s}=\gamma_{g}$.

We assume an equilibrium contact angle of

$\theta_{c,e}=\pi/3,$ (55)

and a gas-solid surface tension $\sigma_{g-s}=0$, which from the Young
equation (22) and (47) imply

$\sigma_{2,e}=-3.64\times 10^{-2}\,N/m.$ (56)

In turn, from equation (17), it follows that

$\rho^{s}_{2,e}=1.2\times 10^{-6}\,\rm{kg}/\rm{m}^{2}.$ (57)

We these physical parameters, we obtain the following dimensionless numbers
$\Rey=36.43$, $S\kern-0.39993ptt\,=0$, $C\kern-0.39993pta\,=1$,
$\bar{\alpha}_{g}=0.002$, $\bar{\beta}_{g}=5\times 10^{5}$,
$\epsilon_{g}=0.0102$, $L_{g}=Q_{g}\epsilon_{g}=0.002$ (see equation 76),
$\bar{\lambda}_{g}=2.5$, $\bar{\beta}_{s}=500$, $\bar{\alpha}_{s}=0.002$,
$\epsilon_{s}=1.02$, $L_{s}=Q_{s}\epsilon_{s}=0.002$, $\bar{\gamma}_{s}=2.5$.

##  4 Results

Simulations show an initial acceleration of the contact line, followed a
comparatively slow descend of the contact-line velocity. The contact angle is
monotonically decreasing from start to end of the contact-line motion.
Consequently, we observe a decreasing contact angle at an accelerating contact
line over the initial portion of the simulation, followed by a longer regime
in which the contact line velocity drops to zero as the contact angle
continues to relax to its equilibrium value.

|
---|---
|
Figure 5: Large-scale view of the spreading process. Arrows show the flow
velocity in the frame of reference of the stationary solid and their lengths
are re-scaled at each time to better illustrate relative distribution of the
flow at each time.

Motion of the droplet is initiated as the out-of-equilibrium surface density
along the liquid-solid line drives mass exchange with the bulk and mass
transport along the surface, driving the liquid-solid surface density towards
its equilibrium. The free surface is sucked into the liquid-solid surface,
thus perturbing the free-surface density in the vicinity of the contact line,
away from its initial equilibrium. These changes in surface densities imply
changes in surface tension, which control the contact angle and induce
Marangoni flow in the bulk, resulting in a fully capillary-driven spontaneous
flow of the droplet towards its equilibrium configuration.

Flow in the bulk is initially concentrated near the contact line (see figure
5), where we do not observe the noisy behaviour associated with the artefact
described in Sprittles & Shikhmurzaev (2011a), indicating that the removal of
such artefacts, by means of the use of our split-domain formulation and
eigensolution, was successful (see figure 6). A video animation of the flow
field is provided as supplementary material.

|
---|---
|
Figure 6: Close-up of the velocity fields in the frame of reference of the moving contact line. Arrow lengths are re-scaled at each time to better illustrate spatial distribution of velocity. |
---|---
|
Figure 7: Contact line velocity and contact angle as functions of time. Panels
correspond to the microsecond-scale process for velocity (a) and contact angle
(c), and the nanosecond-scale regime of contact line acceleration (b) and
initial drop in the contact angle (d).

Figure 7 shows the evolution of the contact-line velocity and contact angle
over the full process (panel a for velocity and c for contact angle), as well
as a detail of the early stages (panel b for velocity and d for contact
angle). Panels a and b also show the vertical velocity of the free surface at
the apex of the droplet. Similarly, panels c and d also compare the dynamic
contact angle (i.e. the result of the solution of the system of equations for
variable $\theta_{c}$) with the contact angle measured from the tangent to the
free surface at the contact line $\theta_{c,m}$. The latter angle is
ultimately determined by the length of the two circular spines that are
closest to the contact line. Moreover, panels c and d, also show the amplitude
of the eigensolution contribution needed to satisfy the limit condition of
axial symmetry of pressure as one approaches the contact line.

Oscillations in the contact line velocity and contact angle curve in figures
7a and 7c, correspond to the apex of the droplet suddenly descending faster
and slowing down soon after that. This reflects the ability of the model to
capture effects that the global configuration of the flow has on the contact
line behaviour.

As can be seen in figures 7c and 7d, the two measurements of the contact angle
($\theta_{c}$ and $\theta_{c,m}$) are virtually indistinguishable. The
convergence of these two ways of measuring the contact angle was used as the
criterion to determine that our mesh was sufficiently refined, as recommended
by Sprittles & Shikhmurzaev (2013).

Figure 7c illustrates how massively important the contribution of the
eigensolution becomes for angles slightly above $90^{o}$, in agreement with
what was reported in Sprittles & Shikhmurzaev (2011a). The eigensolution
contribution is present from the very start of the simulation, though not with
amplitudes comparable to those found when the contact angle is slightly beyond
$90^{o}$.

The resulting relation between contact line velocity and contact angle is
multi-valued, as predicted in Shikhmurzaev (2020), and it shows the influence
of global flow patterns in the oscillations in the curve (which correspond to
the oscillations seen in figure 7a and 7c). Figure 8 summarises the process
trajectory at the contact line, where dashed lines indicate the equilibrium
values of the surface densities and contact angle, and arrows show the
direction of advancing time. It can be clearly seen that the liquid-solid
surface density at the contact line increases monotonically towards its
equilibrium; while the free-surface density starts from equilibrium, falls in
the vicinity of the contact line (as a consequence of the flux of mass into
the liquid-solid surface). As the contact line reaches its maximum velocity,
the free-surface density continues to fall for some time, then gradually
returns to its equilibrium value.

Figure 8 best demonstrates the impossibility to capture the complete picture
of the contact angle behaviour with a single function relating contact line
velocity to contact angle, as we have two contact angles for all but one
velocity. These two “branches” of the contact angle versus contact line
velocity curve are traversed in very different time-scales, as can be seen in
figure 7. Consequently, a function approximating well the relation between
velocity and angle for the lower branch would still capture the greatest part
of the spreading behaviour; however, it would still be lacking the physical
picture of how the droplet lands on that branch.

Figure 8: Process trajectories for contact angle and surface densities at the
contact line.

Away from the contact line, the free-surface density remains essentially
undisturbed. In the far field, the liquid-solid surface density moves, in an
approximately uniform manner, towards its equilibrium, reaching equilibrium in
the far field region before it does in the vicinity of the contact line.

We further consider the degenerate case of $\epsilon_{g}=0$, which can be
thought of as the case when the free-surface relaxation is instantaneous. This
implies that the surface tension on the free surface will be constant
everywhere and throughout the simulation, and from this fact it follows that
there will be no slip between free surface and bulk (see equation 8).

This case is of particular interest to understand the consequences for the
global flow of the perturbation of the free-surface equilibrium in the
vicinity of the contact line. In fact, figure 9 shows that the qualitative
picture is rather similar to the one shown in figure 8, with the obvious
exception of the free-surface density at the contact line being kept constant.
However, the maximum velocity reached by the contact line is only of about
$40\%$ of that reached in the first case, which indicates that the
perturbation of the free surface from its equilibrium near the contact line,
contributes substantially to the mechanics of spontaneous spreading.

It should be noted that, even when the case $\epsilon_{g}=0$ is essentially
equivalent to the relaxation of the free surface being instantaneous, if we
simply took the relaxation time of the free surface $\tau_{g}=0$, coefficient
$Q_{g}$ in equation (9) would become infinite. This issue disappears naturally
with one of the improvements introduced in this work, as the dimensionless
number $Q_{g}$ is effectively removed from the residual equations (see the
discussion around equation 76 in the appendix for the full details). Instead,
the constant $L_{g}=Q_{g}\epsilon_{g}=0$, which is independent of $\tau_{g}$,
is introduced. In this improved implementation, the right-hand side of the
equation controlling mass exchange between the free surface and the bulk (76)
is identically null when the relaxation time is zero, as
$\rho^{s}_{1}=\rho^{s}$ everywhere and at all times. Nevertheless, when
considering the case $\epsilon_{g}=0$, it is still important to impose
$L_{g}=0$, to avoid undesired numerical issues. The imposition of $L_{g}=0$
has no negative consequences in the implementation for this case, as it does
not enter the residuals elsewhere.

Figure 9: Process trajectories for contact angle and surface densities at the
contact line in the limit case $\epsilon_{g}=0$.

##  5 Discussion

We implement the full set of methods, needed to model spontaneous spread of a
droplet from a small initial contact area. These include, the removal of
numerical artefacts known to appear when dealing with dynamic wetting with
obtuse contact angles, the capture of the singular nature of pressure at the
origin by the inclusion of singular elements at the contact line, and the
resolution of dynamics at the scale of the interfacial layer in a continuum
mechanics framework. Moreover, we use a method that produces the contact angle
as an outcome of the application of first principles; which has the added
advantage of allowing the contact angle to be influenced by global flow
patterns. We are thus able to capture the initial regime of contact line
acceleration with a decreasing contact angle,as well as the second, longer
regime of the commonly modelled behaviour.

The proper implementation of these approaches required careful design of a
mesh that is suited to tackle the varied scales of the problem (ranging from
the width of the interfacial layer to the radius of the droplet, and from the
relaxation times of the surfaces to the duration of the spreading), while also
remaining computationally tractable.

In comparison to simpler approaches, e.g. those based on the imposition a
priori of a one-to-one map relating contact-line speed and contact angle, our
model does require a substantially greater amount of work to be implemented;
however, the result is well worth the effort, as cases that can not conform to
the assumptions of other models are now at hand. Approaches that do not
resolve the interfacial scale remain useful to model regimes where the flow
results in the type of contact angle behaviour that they assume.

One difficulty with the present method is found in the implementation of the
Young equation. As its derivative with respect to the contact angle is very
close to zero when $\theta_{c}$ approaches $180^{o}$, thus yielding Jacobians
that turn ill-posed for very large contact angles. It is possible, however,
that this type of numerical issue can be made less pressing by the use of some
form of preconditioning.

The main contribution of the present work is to bring together in a single
implementation all the ingredients that were needed for the simulation of
these spontaneous spreading flows, as well as making the necessary adaptations
in each one. Each of these ingredients had been implemented independently in
separate prior works (Sprittles & Shikhmurzaev, 2011a, b, 2012b, 2013), but
they needed to be used simultaneously to deal with more general cases.
Moreover, in the process of combining these methods, we were also able to
improve some aspects of the implementation, as discussed in section 3.4. We
highlight, in particular, our symmetric implementation of the mass balance
condition for the contact line, which is introduced here and allows for the
natural implementation of this boundary condition for transport along the
surface, without resorting to additional assumptions on the flow of
information, as was done in prior works. Further contributions include
extensive notes on the derivations, which should facilitate future works to
implement this model.

In view of the experimental evidence (De Ruijter et al., 1999; Mohammad Karim
et al., 2016; Bayer & Megaridis, 2006) showing the need for methods that can
produce a contact angle prediction that depends on the global configuration of
the flow, as well as the evidence from lattice Boltzman simulations pointing
in the same direction (Davies et al., 2006); it is clear that there was a
missing _continuum mechanics_ piece of this puzzle. The present work
constitutes a first fully coherent step towards providing that lacking piece
in the modelling picture. Combining, adapting and improving the relevant
existing strategies in a single code proved sufficient to solve the model
problem at hand.

The methods used in this work are naturally relevant to most, if not all,
dynamic wetting problems; thus opening countless possibilities for
ramifications of the present work. More importantly, the present formulation
can address problems that are known to be outside the scope of any other known
method, as discussed in Shikhmurzaev (2020), while also being capable of
addressing the traditional “forced” spreading scenarios.

Ongoing work includes the study of the dependence of the behaviour of the
evolution of contact angle on physical parameters of the fluids for spreading
droplets, the consideration of axisymmetric configurations, spreading over
curved surfaces and spreading from a non-zero initial velocity (such as the
following a droplet impact). Moreover, as the problem here considered involves
a variety of scales, it is likely that its implementation can benefit from the
use of different non-dimensionalisations in different parts of the domain.
This and other strategies to optimise the implementation are being considered.

We highlight the relevance of the modelling here presented for drop-on-demand
additive manufacturing; as this manufacturing technique involves placing
individual on a substrate, where they are allow to spread and then solidify.
Moreover, the framework used for this work has been shown to be capable of
naturally integrating heat transfer and phase transition effects
(Shikhmurzaev, 2021; Belozerov & Shikhmurzaev, 2022), making the addition of
such effects a natural extension of the present work. The methods here
presented are most directly applicable to the cases when drop-on-demand
3D-printing is done with Newtonian fluids, such as in the case of printing
with liquid metal (Simonelli et al., 2019; Gilani et al., 2021, 2022).

The author is very thankful to Prof. Yulii Shikhmurzaev for the support and
feedback, which were both vital for this project, and to Prof. James Sprittles
for helpful discussions. The author gratefully acknowledges the support of
EPSRC project EP/P031684/1.

##  Appendix A Residual equations

We introduce test functions $\phi_{i}$, which yield the residual equations
detailed below.

The $i$-th residual equation for momentum in the $\xi$-direction (where
$\xi=x,y$) is given by

$\mathcal{M}^{\xi}_{i}=M^{\xi}_{i}+\bar{M}^{\xi}_{i},$ (58)

where $M^{\xi}_{i}$ and $\bar{M}^{\xi}_{i}$ are the contribution of the far
and near field, respectively.

For the far field, we have

$M^{\xi}_{i}=M^{\xi,0}_{i}+M^{\xi,1}_{i}+M^{\xi,2}_{i}+M^{\xi,3}_{i}+M^{\xi,4}_{i}$
(59)

where

$\displaystyle M^{x,0}_{i}=$
$\displaystyle-S\kern-0.39993ptt\,\int\limits_{\Omega_{\text{f}}}{\phi_{i}\hat{g}_{x}}+\Rey\int\limits_{\Omega_{\text{f}}}{\phi_{i}\partial_{t}u}+\Rey\int\limits_{\Omega_{\text{f}}}{\left[\phi_{i}u\partial_{x}u+w\partial_{z}u\right]}-\Rey\int\limits_{\Omega_{\text{f}}}{\phi_{i}\left[u^{c}\partial_{x}u+w^{c}\partial_{z}u\right]}$
(60)
$\displaystyle+\int\limits_{\Omega_{\text{f}}}{\left[\partial_{x}u\partial_{x}\phi_{i}+\partial_{z}u\partial_{z}\phi_{i}\right]}-\int\limits_{\Omega_{\text{f}}}{p\partial_{x}\phi_{i}},$
$\displaystyle M^{z,0}_{i}=$
$\displaystyle-S\kern-0.39993ptt\,\int\limits_{\Omega_{\text{f}}}{\phi_{i}\hat{g}_{z}}+\Rey\int\limits_{\Omega_{\text{f}}}{\phi_{i}\partial_{t}w}+\Rey\int\limits_{\Omega_{\text{f}}}{\left[\phi_{i}u\partial_{x}w+w\partial_{z}w\right]}-\Rey\int\limits_{\Omega_{\text{f}}}{\phi_{i}\left[u^{c}\partial_{x}w+w^{c}\partial_{z}w\right]}$
(61)
$\displaystyle+\int\limits_{\Omega_{\text{f}}}{\left[\partial_{x}\phi_{i}\partial_{x}w+\partial_{z}\phi_{i}\partial_{z}w\right]}-\int\limits_{\Omega_{\text{f}}}{p\partial_{z}\phi_{i}},$
$\displaystyle M^{\xi,1}_{i}=$
$\displaystyle-\frac{\sigma^{1}(x_{J},z_{J})\phi_{i}(x_{J},z_{J})}{C\kern-0.39993pta\,}m^{1,\text{f}}_{\xi}(x_{J},z_{J})+\frac{\sigma^{1}(x_{a},z_{a})\phi_{i}(x_{a},z_{a})}{C\kern-0.39993pta\,}m^{1}_{\xi}(x_{a},z_{a})$
(62)
$\displaystyle-\int\limits_{S_{1,\text{f}}}{\phi_{i}\left[n^{1}_{x}\partial_{\xi}u+n^{1}_{z}\partial_{\xi}w\right]}-\int\limits_{S_{1,\text{f}}}{\phi_{i}p^{g}n^{1}_{\xi}}+\frac{1}{C\kern-0.39993pta\,}\int\limits_{S_{1,\text{f}}}{t^{1}_{\xi}\sigma^{1}\partial_{s}\phi_{i}},$
$\displaystyle M^{\xi,2}_{i}=$
$\displaystyle-\int\limits_{S_{2,\text{f}}}{\phi_{i}\left[n^{2}_{x}\partial_{\xi}u+n^{2}_{z}\partial_{\xi}w\right]}+\int\limits_{S_{2,\text{f}}}{\lambda^{2}\phi_{i}n^{2}_{\xi}}+\left(\frac{1}{4C\kern-0.39993pta\,\bar{\alpha}_{s}}+\bar{\beta}_{s}\right)\int\limits_{S_{2,\text{f}}}{\phi_{i}\left[ut^{2}_{x}+wt^{2}_{z}\right]t^{2}_{\xi}}$
(63)
$\displaystyle+\left(\frac{1}{4C\kern-0.39993pta\,\bar{\alpha}_{s}}-\bar{\beta}_{s}\right)\int\limits_{S_{2,\text{f}}}{\phi_{i}\left[u^{s}t^{2}_{x}+w^{s}t^{2}_{z}\right]t^{2}_{\xi}}-\frac{1}{2C\kern-0.39993pta\,\bar{\alpha}_{s}}\int\limits_{S_{2,\text{f}}}{\phi_{i}\left[u^{s}_{2}t^{2}_{x}+w^{s}_{2}t^{2}_{z}\right]t^{2}_{\xi}},$
$\displaystyle M^{\xi,3}_{i}=$
$\displaystyle-\int\limits_{S_{3}}{\phi_{i}\left[n^{3}_{x}\partial_{\xi}u+n^{3}_{z}\partial_{\xi}w\right]}+\int\limits_{S_{3}}{\phi_{i}\left[\gamma^{3}t^{3}_{\xi}+\lambda^{3}n^{3}_{\xi}\right]},$
(64) $\displaystyle M^{\xi,4}_{i}=$
$\displaystyle-\int\limits_{S_{4}}{\phi_{i}\left[n^{4}_{x}\partial_{\xi}u+n^{4}_{z}\partial_{\xi}w\right]}+\int\limits_{S_{4}}{\phi_{i}\left[\gamma^{4}t^{4}_{\xi}+\lambda^{4}n^{4}_{\xi}\right]},$
(65)

where $\hat{\boldsymbol{g}}=(\hat{g}_{x},\hat{g}_{z})$,
$\boldsymbol{u}=(u,w)$, $\boldsymbol{c}=(u^{c},w^{c})$,
$\boldsymbol{n}_{j}=(n^{j}_{x},n^{j}_{z})$,
$\boldsymbol{m}_{j}=(m^{j}_{x},m^{j}_{z})$ and
$\boldsymbol{m}_{1}^{\text{f}}=(m^{1,\text{f}}_{x},m^{1,\text{f}}_{z})$ (with
the f super-index referring to the far-field half of the free surface),
$p_{g}$ is the gas pressure on the free-surface, the sub-index $a$ indicates
the apex of the droplet (intersection of $S_{1}$ and $S_{3}$) and the sub-
index $J$ refers to the intersection of $S_{1}$ and $S_{4}$ (see figure 2).
Moreover, $\boldsymbol{t}_{j}=(t^{j}_{x},t^{j}_{z})$ is the unit tangent to
boundary $S_{j}$ and $\partial_{s}$ stands for the derivative with respect to
the arc-length variable, which increases along $S_{j}$ in the direction of
tangent $\boldsymbol{t}_{j}$, and $\lambda_{j}\boldsymbol{n}_{j}$ and
$\gamma_{j}\boldsymbol{t}_{j}$ are the normal and tangential stress,
respectively. Moreover;
$\mathbb{P}\cdot\boldsymbol{n}_{i}=\lambda^{i}\boldsymbol{n}^{i}+\gamma^{i}\boldsymbol{t}_{i}$,
$\boldsymbol{u_{s}}=(u_{s},w_{s})$, and
$\boldsymbol{v}^{s}_{i}=(u^{s}_{i},w^{s}_{i})$.

We highlight that the first integrals on the right-hand side (RHS) of
equations (62)-(65), though they are line integrals over the boundary, involve
derivatives with respect to $x$ and $z$, and therefore the algorithm used to
calculate these terms must resort to data from triangular elements in the
bulk. These integrals originate from boundary terms that appear through the
use of the Gauss’ divergence theorem on terms that result from the divergence
of the stress tensor. Bulk terms that result from this transformation are
eliminated using the incompressibility condition. This results in an exact
implementation of the incompressibility condition in this these bulk terms,
and it has proven to eliminate a tendency to have noise in the velocities at
the boundary.

For the near field, we have

$\bar{M}^{\xi}_{i}=\bar{M}^{\xi,0}_{i}+\bar{M}^{\xi,1}_{i}+\bar{M}^{\xi,2}_{i}+\bar{M}^{\xi,5}_{i},$
(66)

where once again we have $\xi=x,z$ and where

$\displaystyle\bar{M}^{x,0}_{i}=$
$\displaystyle\Rey\int\limits_{\Omega^{n}}{\phi_{i}\partial_{t}u}-S\kern-0.39993ptt\,\int\limits_{\Omega^{n}}{\phi_{i}\hat{g}_{x}}+A\int\limits_{\Omega^{n}}{\left[\partial_{x}\phi_{i}\partial_{x}\check{u}+\partial_{z}\phi_{i}\partial_{z}\check{u}\right]}$
(67)
$\displaystyle+\Rey\int\limits_{\Omega^{n}}{\phi_{i}\left[\bar{u}\partial_{x}\bar{u}+\bar{w}\partial_{z}\bar{u}-u^{c}\partial_{x}\bar{u}-w^{c}\partial_{z}\bar{u}\right]}$
$\displaystyle+\Rey
A\int\limits_{\Omega^{n}}{\phi_{i}\left[\check{u}\partial_{x}\bar{u}+\check{w}\partial_{z}\bar{u}+\bar{u}\partial_{x}\check{u}+\bar{w}\partial_{z}\check{u}\right]}-\Rey
A\int\limits_{\Omega^{n}}{\phi_{i}\left[u^{c}\partial_{x}\check{u}+w^{c}\partial_{z}\check{u}\right]}$
$\displaystyle+\int\limits_{\Omega^{n}}{\left[\partial_{x}u\partial_{x}\phi_{i}+\partial_{z}u\partial_{z}\phi_{i}\right]}+\Rey\left(A\right)^{2}\int\limits_{\Omega^{n}}{\phi_{i}\left[\check{u}\partial_{x}\check{u}+\check{w}\partial_{z}\check{u}\right]}-\int\limits_{\Omega^{n}}{p\partial_{x}\phi_{i}},\
\ \ \ $ $\displaystyle\bar{M}^{z,0}_{i}=$
$\displaystyle\Rey\int\limits_{\Omega_{\text{n}}}{\phi_{i}\partial_{t}w}-S\kern-0.39993ptt\,\int\limits_{\Omega_{\text{n}}}{\phi_{i}\hat{g}_{z}}+A\int\limits_{\Omega_{\text{n}}}{\left[\partial_{x}\phi_{i}\partial_{x}\check{w}+\partial_{z}\phi_{i}\partial_{z}\check{w}\right]}$
(68)
$\displaystyle+\Rey\int\limits_{\Omega_{\text{n}}}{\phi_{i}\left[\bar{u}\partial_{x}\bar{w}+\bar{w}\partial_{z}\bar{w}-u^{c}\partial_{x}\bar{w}-w^{c}\partial_{z}\bar{w}\right]}$
$\displaystyle+\Rey
A\int\limits_{\Omega_{\text{n}}}{\phi_{i}\left[\check{u}\partial_{x}\bar{w}+\check{w}\partial_{z}\bar{w}+\bar{u}\partial_{x}\check{w}+\bar{w}\partial_{z}\check{w}\right]}$
$\displaystyle-\Rey
A\int\limits_{\Omega_{\text{n}}}{\phi_{i}\left[u^{c}\partial_{x}\check{w}+w^{c}\partial_{z}\check{w}\right]}+\int\limits_{\Omega_{\text{n}}}{\left[\partial_{x}\phi_{i}\partial_{x}w+\partial_{z}\phi_{i}\partial_{z}w\right]}$
$\displaystyle+\Rey\left(A\right)^{2}\int\limits_{\Omega_{\text{n}}}{\phi_{i}\left[\check{u}\partial_{x}\check{w}+\check{w}\partial_{z}\check{w}\right]}-\int\limits_{\Omega_{\text{n}}}{p\partial_{z}\phi_{i}},\
\ \ \ $ $\displaystyle\bar{M}^{\xi,1}_{i}=$
$\displaystyle\frac{\sigma^{1}(x_{J},z_{J})\phi_{i}(x_{J},z_{J})}{C\kern-0.39993pta\,}m^{1,n}_{\xi}(x_{J},z_{J})+\frac{\sigma^{1}(x_{c},z_{c})\phi_{i}(x_{c},z_{c})}{C\kern-0.39993pta\,}m^{1}_{\xi}(x_{c},z_{c})$
(69)
$\displaystyle-A\int\limits_{S_{1,\text{n}}}{\phi_{i}\left[n^{1}_{x}\partial_{\xi}\check{u}+n^{1}_{z}\partial_{\xi}\check{w}\right]}-\int\limits_{S_{1,\text{n}}}{\phi_{i}\left[n^{1}_{x}\partial_{\xi}\bar{u}+n^{1}_{z}\partial_{\xi}\bar{w}\right]}$
$\displaystyle-\int\limits_{S_{1,\text{n}}}{\phi_{i}p^{g}n^{1}_{\xi}}+\frac{1}{C\kern-0.39993pta\,}\int\limits_{S_{1,\text{n}}}{t^{1}_{\xi}\sigma^{1}\partial_{s}\phi_{i}},$
$\displaystyle\bar{M}^{\xi,2}_{i}=$
$\displaystyle-\int\limits_{S_{2,\text{n}}}{\phi_{i}\left[A\left(n^{2}_{r}\partial_{\xi}\check{u}+n^{2}_{z}\partial_{\xi}\check{w}\right)+n^{2}_{r}\partial_{\xi}\bar{u}+n^{2}_{z}\partial_{\xi}\bar{w}\right]}$
(70)
$\displaystyle{+{\left(\frac{1}{4C\kern-0.39993pta\,\bar{\alpha}_{s}}+\bar{\beta}_{s}\right)A\int\limits_{S_{2,\text{n}}}\phi_{i}\left[\check{u}t^{2}_{r}+\check{w}t^{2}_{z}\right]t^{2}_{\xi}}}$
$\displaystyle+\left(\frac{1}{4C\kern-0.39993pta\,\bar{\alpha}_{s}}+\bar{\beta}_{s}\right)\int\limits_{S_{2,\text{n}}}{\phi_{i}\left[\bar{u}t^{2}_{r}+\bar{w}t^{2}_{z}\right]t^{2}_{\xi}}$
$\displaystyle+\left(\frac{1}{4C\kern-0.39993pta\,\bar{\alpha}_{s}}-\bar{\beta}_{s}\right)\int\limits_{S_{2,\text{n}}}{\phi_{i}\left[u^{s}t^{2}_{r}+w^{s}t^{2}_{z}\right]t^{2}_{\xi}}$
$\displaystyle-\frac{1}{2C\kern-0.39993pta\,\bar{\alpha}_{s}}\int\limits_{S_{2,\text{n}}}{\phi_{i}\left[u^{s}_{2}t^{2}_{r}+w^{s}_{2}t^{2}_{z}\right]t^{2}_{\xi}}+\int\limits_{S_{2,\text{n}}}{\lambda^{2}\phi_{i}n^{2}_{\xi}},$

and

$\displaystyle\bar{M}^{\xi,4}_{i}=$ $\displaystyle
A\int\limits_{S_{4}}{\phi_{i}\left[n^{4}_{r}\partial_{\xi}\check{u}+n^{4}_{z}\partial_{\xi}\check{w}\right]}-\int\limits_{S_{4}}{\phi_{i}\left[n^{4}_{r}\partial_{\xi}\bar{u}+n^{4}_{z}\partial_{\xi}\bar{w}\right]}$
(71)
$\displaystyle-\int\limits_{S_{4}}{\phi_{i}\left[\lambda^{4}n^{4}_{\xi}+\gamma^{4}t^{4}_{\xi}\right]};$

where $\boldsymbol{\bar{u}}=(\bar{u},\bar{w})$,
$\check{\boldsymbol{u}}=(\check{u},\check{w})$ and
$\boldsymbol{m}^{\text{n}}_{1}=(m^{1,\text{n}}_{x},m^{1,\text{n}}_{z})$ (where
the n super-index referring to the near half of the free surface). We also
highlight that $\boldsymbol{n}^{1}$ and $\boldsymbol{n}^{2}$ point into the
near-field domain, while $\boldsymbol{n}^{4}$ points out of it.

It is important to stress that the integrand in the first term on the RHS of
equation (67) involves the full velocity of the fluid in the $x$-direction
(i.e. $u=\bar{u}+A\check{u}$). In this term, it is more convenient not to
separate these two components of the velocity field at this point. In this
way, when the finite difference approximation for $\partial_{t}u$ is done,
only the velocity at the time-step being calculated needs to be split into its
two components, making the algorithm independent of the contribution of the
eigensolution at previous time steps.

We highlight that the third integral on the RHS of equation (67) comes from
applying Guass’ divergence theorem on boundary integrals of the stress caused
by the eigensolution. We note that the tangential stress caused by the
eigensolution is only null along the boundaries of the wedge of angle
$\theta_{c}$, not along the curved free surface that bounds the fluid domain.

We also note that the first term on the RHS of equation (62) and the last
integral on the RHS of equation (65) will be cancelled by the analogue terms
for the near field given in (69) and (71); which means that, in the
computational implementation, there is no need to introduce variables
$\lambda^{4}$, $\gamma^{4}$ and $\boldsymbol{m}_{j}(x_{J},z_{J})$, as they are
cancelled in the residual equations once near-field and far-field
contributions are included.

Following Sprittles & Shikhmurzaev (2012b), we introduce linear test functions
for the continuity equation, which we call $\psi_{i}$. The continuity
residuals are thus given by

$C_{i}=\int\limits_{\Omega_{\text{f}}}{\psi_{i}\partial_{x}u+\int\limits_{\Omega_{\text{f}}}\psi_{i}\partial_{z}w}.$
(72)

We note that functions $\psi_{i}$, which are linear in the master element
coordinates, will also be used in the linear combination to expand the
pressure in the bulk.

###  A.1 Free surface

We introduce surface test functions $\phi^{j}_{i}$, which coincide with the
bulk test functions $\phi_{i}$ when sampled on boundary $S_{j}$, but are re-
numbered along the boundary, so as to skip those which are identically zero on
the interfaces.

The KBC residuals are given by $\mathcal{K}_{i}=K_{i}+\bar{K}_{i}$, where

$K_{i}=\int\limits_{S_{1,\text{f}}}{\phi^{1}_{i}\left[n^{1}_{x}\left(u^{s}_{1}-u^{c}\right)+n^{1}_{z}\left(w^{s}_{1}-w^{c}\right)\right]},$
(73)

and the expression for $\bar{K}_{i}$ is identical, except for the domain of
integration.

The condition for slip between the bulk and the surface phase yield the
residuals $\mathcal{S}^{1}_{i}=S^{1}_{i}+\bar{S}^{1}_{i}$, where

$\displaystyle S^{1}_{i}=$
$\displaystyle\int\limits_{S_{1,\text{f}}}{\phi^{1}_{i}t^{1}_{x}\left[u^{s}_{1}-u\right]}+\int\limits_{S_{1,\text{f}}}{\phi^{1}_{i}t^{1}_{z}\left[w^{s}_{1}-w\right]}+\frac{1+4\bar{\alpha}_{g}\bar{\beta}_{g}}{4\bar{\beta}_{g}}\int\limits_{S_{1,\text{f}}}{\sigma^{1}\partial_{s}\phi^{1}_{i}}$
(74)
$\displaystyle-\frac{1+4\bar{\alpha}_{g}\bar{\beta}_{g}}{4\bar{\beta}_{g}}\phi^{1}_{i}(x_{a},z_{a})\sigma^{1}(x_{a},z_{a})+\frac{1+4\bar{\alpha}_{g}\bar{\beta}_{g}}{4\bar{\beta}_{g}}\phi^{1}_{i}(x_{J},z_{J})\sigma^{1}(x_{J},z_{J})$

and

$\displaystyle\bar{S}^{1}_{i}=$
$\displaystyle\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\left[u^{s}_{1}-\bar{u}\right]t^{1}_{r}}+\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\left[w^{s}_{1}-\bar{w}\right]t^{1}_{z}}-A\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\check{u}t^{1}_{r}}-A\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\check{w}t^{1}_{z}}$
(75)
$\displaystyle+\frac{1+4\bar{\alpha}_{g}\bar{\beta}_{g}}{4\bar{\beta}_{g}}\int\limits_{S_{1,\text{n}}}{\sigma^{1}\partial_{s}\phi^{1}_{i}}$
$\displaystyle+\frac{1+4\bar{\alpha}_{g}\bar{\beta}_{g}}{4\bar{\beta}_{g}}\phi^{1}_{i}(r_{c},z_{c})\sigma^{1}(r_{c},z_{c})-\frac{1+4\bar{\alpha}_{g}\bar{\beta}_{g}}{4\bar{\beta}_{g}}\phi^{1}_{i}(r_{J},z_{J})\sigma^{1}(r_{J},z_{J}).$

We note that the last term of the RHS of the two equations above will cancel
each other, sa $\phi^{1}_{i}$ is continuous.

To derive the mass exchange residuals, we combine equation (9) with (11),
obtaining

$\left(\boldsymbol{u}-\boldsymbol{v}^{s}_{1}\right)\cdot\boldsymbol{n}_{1}=-L_{g}\left[\frac{\partial\rho^{s}_{1}}{\partial
t}+\nabla\cdot\left(\rho^{s}_{1}\boldsymbol{v}^{s}_{1}\right)\right],$ (76)

where $L_{g}=\epsilon_{g}Q_{g}=\rho^{s}_{(0)}/(\rho L)$. This substitution is
required as the original form of (9) yields residuals that lead to numerical
issues (unnecessary operations between very small and very large numbers).
These problems are naturally avoided once the equation is reformulated as
shown above. Thus, using equation (30), we have
$\mathcal{E}^{1}_{i}=E^{1}_{i}+\bar{E}^{1}_{i}$, where

$\displaystyle E^{1}_{i}=$
$\displaystyle\left\\{L_{g}\phi_{i}\rho^{s}_{1}\left[m^{1,\text{f}}_{r}\left(u_{c}-u^{s}_{1}\right)+m^{1,\text{f}}_{z}\left(w_{c}-w^{s}_{1}\right)\right]\right\\}|_{(x_{J},z_{J})}$
(77)
$\displaystyle+L_{g}\int\limits_{S_{1,\text{f}}}{\phi^{1}_{i}\partial_{t}\rho^{s}_{1}}-L_{g}\int\limits_{S_{1,\text{f}}}{\rho^{s}_{1}\left(\partial_{s}\phi^{1}_{i}\right)\left(u^{s}_{1}t^{1}_{r}+w^{s}_{1}t^{1}_{z}\right)}+L_{g}\int\limits_{S_{1,\text{f}}}{\rho^{s}_{1}\left(\partial_{s}\phi^{1}_{i}\right)\left(t^{1}_{r}u_{c}+t^{1}_{z}w_{c}\right)}$
$\displaystyle+L_{g}\int\limits_{S_{1,\text{f}}}{\phi^{1}_{i}\rho^{s}_{1}\left(t^{1}_{r}\partial_{s}u_{c}+t^{1}_{z}\partial_{s}w_{c}\right)}+\int\limits_{S_{1,\text{f}}}{\phi^{1}_{i}n^{1}_{x}\left[u-u^{s}_{1}\right]}+\int\limits_{S_{1,\text{f}}}{\phi^{1}_{i}n^{1}_{z}\left[w-w^{s}_{1}\right]}\
\ $

and

$\displaystyle\bar{E}^{1}_{i}=$
$\displaystyle\left\\{L_{g}\phi_{i}\rho^{s}_{1}\left[m^{1,\text{f}}_{r}\left(u_{c}-u^{s}_{1}\right)+m^{1,\text{f}}_{z}\left(w_{c}-w^{s}_{1}\right)\right]\right\\}|_{(x_{J},z_{J})}$
(78)
$\displaystyle-\left\\{L_{g}\phi_{i}\rho^{s}_{1}\left[m^{2}_{r}\left(u_{c}-u^{s}_{2}\right)+m^{2}_{z}\left(w_{c}-w^{s}_{2}\right)\right]\right\\}|_{(x_{c},z_{c})}$
$\displaystyle+L_{g}\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\partial_{t}\rho^{s}_{1}}-L_{g}\int\limits_{S_{1,\text{n}}}{\rho^{s}_{1}\left(\partial_{s}\phi^{1}_{i}\right)\left(u^{s}_{1}t^{1}_{r}+w^{s}_{1}t^{1}_{z}\right)}+L_{g}\int\limits_{S_{1,\text{n}}}{\rho^{s}_{1}\left(\partial_{s}\phi^{1}_{i}\right)\left(t^{1}_{r}u_{c}+t^{1}_{z}w_{c}\right)}$
$\displaystyle+L_{g}\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\rho^{s}_{1}\left(t^{1}_{r}\partial_{s}u_{c}+t^{1}_{z}\partial_{s}w_{c}\right)}+A\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\check{u}n^{1}_{x}}+A\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\check{w}n^{1}_{z}}$
$\displaystyle\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\left[\bar{u}-u^{s}_{1}\right]n^{1}_{x}}+\int\limits_{S_{1,\text{n}}}{\phi^{1}_{i}\left[\bar{w}-w^{s}_{1}\right]n^{1}_{z}},$

where we have also used the mass balance at the contact line, given by
equation (32), to express flows at the contact line in terms of the flux
entering boundary $2$. This substitution, and others very similar to it, are,
in fact, the only way in which equation (32) is imposed.

We note that the first term of each of the latter two equations (boundary
contributions assessed at point $J$) will cancel each other out when the far-
field and near-field contributions are added together. We also highlight that
the boundary terms for flux at the apex of the droplet is omitted as it is
equals to zero, which follows from equations (25) and the fact that $u^{c}=0$
on boundary $3$.

The mass transport along the free surface yields

$\displaystyle\mathcal{D}^{1}_{i}=$
$\displaystyle-\left\\{\epsilon_{g}\phi_{i}\rho^{s}_{1}\left[m^{2}_{r}\left(u_{c}-u^{s}_{2}\right)+m^{2}_{z}\left(w_{c}-w^{s_{2}}_{c}\right)\right]\right\\}|_{(x_{c},z_{c})}$
(79)
$\displaystyle-\rho^{s}_{1,e}\int\limits_{S_{1}}{\phi^{1}_{i}}+\int\limits_{S_{1}}{\phi^{1}_{i}\rho^{s}_{1}}+\epsilon_{g}\int\limits_{S_{1}}{\phi^{1}_{i}\partial_{t}\rho^{s}_{1}}-\epsilon_{g}\int\limits_{S_{1}}{\rho^{s}_{1}\left(\partial_{s}\phi^{1}_{i}\right)\left(u^{s}_{1}t^{1}_{r}+w^{s}_{1}t^{1}_{z}\right)}$
$\displaystyle+\epsilon_{g}\int\limits_{S_{1}}{\rho^{s}_{1}\left(\partial_{s}\phi^{1}_{i}\right)\left(t^{1}_{r}u_{c}+t^{1}_{z}w_{c}\right)}+\epsilon_{g}\int\limits_{S_{1}}{\phi^{1}_{i}\rho^{s}_{1}\left(t^{1}_{r}\partial_{s}u_{c}+t^{1}_{z}\partial_{s}w_{c}\right)}.$

For simplicity, the equation above has not been divided into near-field and
far-field contributions, since the expression for both are essentially
identical. Moreover, had we separated the contributions, we would have had
boundary terms at the separatrix, which once again would be present in both
contributions, but would always cancel each other out. Moreover, boundary
terms for the apex are identically null for the same reasons as mentioned for
the mass exchange residuals, and the boundary terms for the contact line in
(79) have been written in terms of the liquid-solid surface variables, using
equation (32). We note that this substitution was not performed in Sprittles &
Shikhmurzaev (2013) nor Sprittles & Shikhmurzaev (2012a). However, we find
that this additional substitution is necessary and sufficient to satisfy
equation (32) in all cases considered here.

Finally, the state equation relating surface density and surface tension along
boundary $1$ is satisfied in its strong form, yielding

$T^{1}_{i}=\sigma^{1}_{i}+\lambda_{g}\left(\rho^{s}_{1}-1\right).$ (80)

###  A.2 The liquid-solid surface

The impermeability equation leads to the following residual

$I_{i}=\int\limits_{S_{2}}{\phi^{2}_{i}n^{2}_{x}\left(u^{s}_{2}-u_{s}\right)}+\int\limits_{S_{2}}{\phi^{2}_{i}n^{2}_{z}\left(w^{s}_{2}-w_{s}\right)},$
(81)

where, once again there is no need to separate the integration domains, as the
expression for this particular equation are identical in both domains.

The equation for slip between the bulk and the surface leads to
$\mathcal{S}_{i}=S_{i}+\bar{S}_{i}$, where

$\displaystyle S^{2}_{i}=$
$\displaystyle\bar{\alpha}_{s}\phi^{2}_{i}(x_{L},z_{L})\sigma^{2}(x_{L},z_{L})-\bar{\alpha}_{s}\phi^{2}_{i}(x_{o},z_{o})\sigma^{2}(x_{o},z_{o})$
(82)
$\displaystyle+\int\limits_{S_{2,\text{f}}}{\phi^{2}_{i}\left(u^{s}_{2}t^{2}_{x}-\frac{1}{2}u-\frac{1}{2}u_{s}\right)}+\int\limits_{S_{2,\text{f}}}{\phi^{2}_{i}\left(w^{s}_{2}-\frac{1}{2}w-w_{s}\frac{1}{2}\right)t^{2}_{z}}$
$\displaystyle+\bar{\alpha}_{s}\int\limits_{S_{2,\text{f}}}{\sigma^{2}\partial_{s}\phi^{2}_{i}},$

and

$\displaystyle\bar{\mathcal{S}}^{2}_{i}=$
$\displaystyle\bar{\alpha}_{s}\phi^{2}_{i}(x_{c},z_{c})\sigma^{2}(x_{c},z_{c})-\bar{\alpha}_{s}\phi^{2}_{i}(x_{L},z_{L})\sigma^{2}(x_{L},z_{L})$
(83)
$\displaystyle+\int\limits_{S_{2,\text{n}}}{\phi^{2}_{i}\left(u^{s}_{2}-\frac{1}{2}\bar{u}-\frac{1}{2}u_{s}\right)t^{2}_{r}}+\int\limits_{S_{2,\text{n}}}{\phi^{2}_{i}\left(w^{s}_{2}-\frac{1}{2}\bar{w}-\frac{1}{2}w_{s}\right)t^{2}_{z}}$
$\displaystyle-\frac{1}{2}A\int\limits_{S_{2,\text{n}}}\check{u}t^{2}_{r}-\frac{1}{2}A\int\limits_{S_{2,\text{n}}}\check{w}t^{2}_{z}+\bar{\alpha}_{s}\int\limits_{S_{2,\text{n}}}{\sigma^{2}\partial_{s}\phi^{2}_{i}},$

where $(x_{L},z_{L})$ is the intersection between boundary $2$ and the
separatrix. We highlight that once again, the boundary term corresponding to
the separatrix will be cancelled out whenever the test functions (in this case
$\phi^{2}_{i}$) are continuous.

Following transformations that are entirely analogous to the ones performed
for equation (77), the condition for mass exchange between surface and bulk
yields $\mathcal{E}^{2}_{i}=E^{2}_{i}+\bar{E}^{2}_{i}$, where

$\displaystyle E^{2}_{i}=$
$\displaystyle\left\\{L_{s}\delta_{i,c}\rho^{s}_{2}\left[m^{2,\text{f}}_{x}\left(u_{c}-u^{s}_{2}\right)+m^{2,\text{f}}_{z}\left(w_{c}-w^{s}_{2}\right)\right]\right\\}_{x_{K},z_{K}}$
(84)
$\displaystyle+L_{s}\int\limits_{S_{2,\text{f}}}{\phi^{2}_{i}\partial_{t}\rho^{s}_{2}}-L_{s}\int\limits_{S_{2,\text{f}}}{\rho^{s_{2}}\left(\partial_{s}\phi^{2}_{i}\right)\left(u^{s_{2}}t^{2}_{r}+w^{s_{2}}t^{2}_{z}\right)}$
$\displaystyle+L_{s}\int\limits_{S_{2,\text{f}}}{\rho^{s_{2}}\left(\partial_{s}\phi^{2}_{i}\right)\left(t^{2}_{r}u_{c}-t^{2}_{z}w_{c}\right)}+L_{s}\int\limits_{S_{2,\text{f}}}{\phi^{2}_{i}\rho^{s}_{2}\left(t^{2}_{r}\partial_{s}u_{c}+t^{2}_{z}\partial_{s}w_{c}\right)}$
$\displaystyle\int\limits_{S_{2,\text{f}}}{\phi^{2}_{i}n^{2}_{x}\left[u-u^{s}_{2}\right]}+\int\limits_{S_{2,\text{f}}}{\phi^{2}_{i}n^{2}_{z}\left[w-w^{s}_{2}\right]}$

and

$\displaystyle\bar{E}^{2}_{i}=$
$\displaystyle\left\\{L_{s}\delta_{i,c}\rho^{s}_{2}\left[m^{2,\text{n}}_{x}\left(u_{c}-u^{s}_{2}\right)+m^{2,\text{n}}_{z}\left(w_{c}-w^{s}_{2}\right)\right]\right\\}_{x_{K},z_{K}}$
(85)
$\displaystyle-\left\\{L_{s}\delta_{i,c}\rho^{s}_{2}\left[m^{1}_{x}\left(u_{c}-u^{s}_{1}\right)+m^{1}_{z}\left(w_{c}-w^{s}_{1}\right)\right]\right\\}_{x_{c},z_{c}}$
$\displaystyle+L_{s}\int\limits_{S_{2,\text{n}}}{\phi^{2}_{i}\partial_{t}\rho^{s}_{2}}-L_{s}\int\limits_{S_{2,\text{n}}}{\rho^{s_{2}}\left(\partial_{s}\phi^{2}_{i}\right)\left(u^{s_{2}}t^{2}_{r}+w^{s_{2}}t^{2}_{z}\right)}$
$\displaystyle+L_{s}\int\limits_{S_{2,\text{n}}}{\rho^{s_{2}}\left(\partial_{s}\phi^{2}_{i}\right)\left(t^{2}_{r}u_{c}-t^{2}_{z}w_{c}\right)}+L_{s}\int\limits_{S_{2,\text{n}}}{\phi^{2}_{i}\rho^{s}_{2}\left(t^{2}_{r}\partial_{s}u_{c}+t^{2}_{z}\partial_{s}w_{c}\right)}$
$\displaystyle\int\limits_{S_{2,\text{n}}}{\phi^{2}_{i}\left(\bar{u}-u^{s}_{2}\right)n^{2}_{r}}+\int\limits_{S_{2,\text{n}}}{\phi^{2}_{i}\left(\bar{w}-w^{s}_{2}\right)n^{2}_{z}}$
$\displaystyle+A\int\limits_{S_{2,\text{n}}}{\phi^{2}_{i}\check{u}n^{2}_{r}}+A\int\limits_{S_{2,\text{n}}}{\phi^{2}_{i}\check{w}n^{2}_{z}},$

where $L_{s}=\epsilon_{g}Q_{g}=\rho^{s}_{(0)}/(\rho L)$. Here, once again the
boundary term that correspond to the separatrix will cancel each other out
when we add the contributions from both sub-domains.

The condition for mass transport along the surface results in

$\displaystyle\mathcal{D}^{2}_{i}=$
$\displaystyle-\left\\{\epsilon_{s}\delta_{i,c}\rho^{s}_{2}\left[m^{1}_{x}\left(u_{c}-u^{s}_{1}\right)+m^{1}_{z}\left(w_{c}-w^{s}_{1}\right)\right]\right\\}_{x_{c},z_{c}}$
(86)
$\displaystyle-\rho^{s}_{2,e}\int\limits_{S_{2}}{\phi^{2}_{i}}+\int\limits_{S_{2}}{\phi^{2}_{i}\rho^{s}_{2}}+\epsilon_{s}\int\limits_{S_{2}}{\phi^{2}_{i}\partial_{t}\rho^{s}_{2}}-\epsilon_{s}\int\limits_{S_{2}}{\rho^{s_{2}}\left(\partial_{s}\phi^{2}_{i}\right)\left(u^{s_{2}}t^{2}_{r}+w^{s_{2}}t^{2}_{z}\right)}$
$\displaystyle+\epsilon_{s}\int\limits_{S_{2}}{\rho^{s_{2}}\left(\partial_{s}\phi^{2}_{i}\right)\left(t^{2}_{r}u_{c}-t^{2}_{z}w_{c}\right)}+\epsilon_{s}\int\limits_{S_{2}}{\phi^{2}_{i}\rho^{s}_{2}\left(t^{2}_{r}\partial_{s}u_{c}+t^{2}_{z}\partial_{s}w_{c}\right)}.$

Here, in the same spirit as was done in equation (79), we have used equation
(32) to express the boundary terms as a function of the free-surface
variables.

As we did on boundary $1$, we can also satisfy the state equation in its
strong form, obtaining

$T^{2}_{i}=\sigma^{2}_{i}+\lambda_{s}\left(\rho^{s}_{2}-1\right).$ (87)

###  A.3 Pressure condition

When the contact angle is obtuse, we include $A$ as an unknown and the
equation

$p(x_{s},z_{s})-p(x_{g},z_{g})=0,$ (88)

where $(x_{s},z_{s})$ and $(x_{g},z_{g})$ are the two pressure nodes, one on
each interface, that are closest to the contact line, for which our mesh
ensures that they be located at the same distance from the contact line.

## References

  * Bayer & Megaridis (2006) Bayer, Ilker S & Megaridis, Constantine M 2006 Contact angle dynamics in droplets impacting on flat surfaces with different wetting characteristics. Journal of Fluid Mechanics 558, 415–449.
  * Bazhlekov & Chesters (1996) Bazhlekov, Ivan B & Chesters, Allan K 1996 Numerical investigation of the dynamic influence of the contact line region on the macroscopic meniscus shape. Journal of Fluid Mechanics 329, 137–146.
  * Belozerov & Shikhmurzaev (2022) Belozerov, A. A. & Shikhmurzaev, Y. D. 2022 The onset of solidification: from interface formation to the stefan regime. Journal of Chemical Physics Accepted (XXXX), XXXX.
  * Celaya et al. (2014) Celaya, E Alberdi, Aguirrezabala, JJ Anza & Chatzipantelidis, Panagiotis 2014 Implementation of an adaptive bdf2 formula and comparison with the matlab ode15s. Procedia Computer Science 29, 1014–1026.
  * Davies et al. (2006) Davies, Alan R, Summers, Jonathan L & Wilson, Mark CT 2006 On a dynamic wetting model for the finite-density multiphase lattice boltzmann method. International Journal of Computational Fluid Dynamics 20 (6), 415–425.
  * De Ruijter et al. (1999) De Ruijter, Michel J, Blake, TD & De Coninck, Joël 1999 Dynamic wetting studied by molecular modeling simulations of droplet spreading. Langmuir 15 (22), 7836–7847.
  * Denner (2014) Denner, Anja Katrin 2014 Experiments on temporal variable step bdf2 algorithms. Master’s thesis, University of Wisconsin Milwaukee.
  * Gilani et al. (2021) Gilani, N, Aboulkhair, NT, Simonelli, M, East, M, Ashcroft, I & Hague, RJM 2021 Insights into drop-on-demand metal additive manufacturing through an integrated experimental and computational study. Additive Manufacturing 48, 102402.
  * Gilani et al. (2022) Gilani, Negar, Aboulkhair, Nesma T, Simonelli, Marco, East, Mark, Ashcroft, Ian A & Hague, Richard JM 2022 From impact to solidification in drop-on-demand metal additive manufacturing using metaljet. Additive Manufacturing 55, 102827\.
  * Hoffman (1975) Hoffman, Richard L 1975 A study of the advancing interface. i. interface shape in liquid—gas systems. Journal of colloid and interface science 50 (2), 228–241.
  * Legendre & Maglio (2013) Legendre, Dominique & Maglio, Marco 2013 Numerical simulation of spreading drops. Colloids and Surfaces A: Physicochemical and Engineering Aspects 432, 29–37.
  * Mohammad Karim et al. (2016) Mohammad Karim, A, Davis, Stephen H & Kavehpour, HP 2016 Forced versus spontaneous spreading of liquids. Langmuir 32 (40), 10153–10158.
  * Ngan & Dussan (1982) Ngan, Clifton G & Dussan, Elizabeth 1982 On the nature of the dynamic contact angle: an experimental study. Journal of Fluid Mechanics 118, 27–40.
  * Ranabothu et al. (2005) Ranabothu, Srinivas R, Karnezis, Cassandra & Dai, Lenore L 2005 Dynamic wetting: hydrodynamic or molecular-kinetic? Journal of colloid and interface science 288 (1), 213–221.
  * Ren & E (2007) Ren, Weiqing & E, Weinan 2007 Boundary conditions for the moving contact line problem. Physics of fluids 19 (2), 022101\.
  * Seebergh & Berg (1992) Seebergh, Jill E & Berg, John C 1992 Dynamic wetting in the low capillary number regime. Chemical Engineering Science 47 (17-18), 4455–4464.
  * Shikhmurzaev (1993) Shikhmurzaev, Yulii D 1993 The moving contact line on a smooth solid surface. International Journal of Multiphase Flow 19 (4), 589–610.
  * Shikhmurzaev (2007) Shikhmurzaev, Yulii D 2007 Capillary flows with forming interfaces. CRC Press.
  * Shikhmurzaev (2020) Shikhmurzaev, Yulii D 2020 Moving contact lines and dynamic contact angles: a ‘litmus test’ for mathematical models, accomplishments and new challenges. The European Physical Journal Special Topics 229 (10), 1945–1977.
  * Shikhmurzaev (2021) Shikhmurzaev, Yulii D 2021 Solidification and dynamic wetting: A unified modeling framework. Physics of Fluids 33 (7), 072101.
  * Šikalo et al. (2005) Šikalo, Š, Tropea, C & Ganić, EN 2005 Dynamic wetting angle of a spreading droplet. Experimental Thermal and Fluid Science 29 (7), 795–802.
  * Simonelli et al. (2019) Simonelli, Marco, Aboulkhair, Nesma, Rasa, Mircea, East, Mark, Tuck, Chris, Wildman, Ricky, Salomons, Otto & Hague, Richard 2019 Towards digital metal additive manufacturing via high-temperature drop-on-demand jetting. Additive Manufacturing 30, 100930.
  * Sprittles & Shikhmurzaev (2011a) Sprittles, James E & Shikhmurzaev, Yulii D 2011a Viscous flow in domains with corners: Numerical artifacts, their origin and removal. Computer methods in applied mechanics and engineering 200 (9-12), 1087–1099.
  * Sprittles & Shikhmurzaev (2011b) Sprittles, James E & Shikhmurzaev, Yulii D 2011b Viscous flows in corner regions: Singularities and hidden eigensolutions. International Journal for Numerical Methods in Fluids 65 (4), 372–382.
  * Sprittles & Shikhmurzaev (2012a) Sprittles, James E & Shikhmurzaev, Yulii D 2012a The dynamics of liquid drops and their interaction with solids of varying wettabilities. Physics of Fluids 24 (8), 082001.
  * Sprittles & Shikhmurzaev (2012b) Sprittles, James E & Shikhmurzaev, Yulii D 2012b Finite element framework for describing dynamic wetting phenomena. International Journal for Numerical Methods in Fluids 68 (10), 1257–1298.
  * Sprittles & Shikhmurzaev (2013) Sprittles, James E & Shikhmurzaev, Yulii D 2013 Finite element simulation of dynamic wetting flows as an interface formation process. Journal of Computational Physics 233, 34–65.
  * Zhang et al. (2009) Zhang, Linbo, Cui, Tao & Liu, Hui 2009 A set of symmetric quadrature rules on triangles and tetrahedra. Journal of Computational Mathematics pp. 89–96.