4.1 Inner-region dynamics

It turns out that, contrary to the outer-region analysis, which relies on assuming a nearly circular contact line, the inner-region asymptotics of $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h$ as the outer region is approached can be determined analytically up to $O(\unicode[STIX]{x1D706})$ , provided that $a(\unicode[STIX]{x1D719},t)$ varies sufficiently smoothly that terms like $\unicode[STIX]{x1D706}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}a$ and $\unicode[STIX]{x1D706}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}^{2}a$ can be neglected. To study the dynamics in the inner region, we consider the stretching transformation (see also Lacey Reference Lacey1982),

(4.2a,b ) $$\begin{eqnarray}\boldsymbol{x}=\boldsymbol{c}+\frac{\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D717}_{\ast }}(\unicode[STIX]{x1D701}\unicode[STIX]{x1D749}-\unicode[STIX]{x1D709}\unicode[STIX]{x1D742}),\quad h=\unicode[STIX]{x1D706}\unicode[STIX]{x1D6F6},\end{eqnarray}$$

where $\unicode[STIX]{x1D701}$ and $\unicode[STIX]{x1D709}$ are measures of the distances along the unit tangent, $\unicode[STIX]{x1D749}$ , and unit outward normal vectors, $\unicode[STIX]{x1D742}$ , at the contact line $\boldsymbol{x}=\boldsymbol{c}$ , respectively. Under this change of variables, $\unicode[STIX]{x1D6F6}=\unicode[STIX]{x1D709}$ satisfies the local contact angle condition (3.10), so that a quasi-static expansion of the form $\unicode[STIX]{x1D6F6}=\unicode[STIX]{x1D709}+\tilde{\unicode[STIX]{x1D6F6}}+\cdots \,$ with $\unicode[STIX]{x1D709}\gg \tilde{\unicode[STIX]{x1D6F6}}$ yields the following equation for $\tilde{\unicode[STIX]{x1D6F6}}$ :

(4.3) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}^{3}\tilde{\unicode[STIX]{x1D6F6}}=-\frac{\unicode[STIX]{x2202}_{t}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D742}}{\unicode[STIX]{x1D717}_{\ast }^{3}(\unicode[STIX]{x1D709}^{2}+1)},\end{eqnarray}$$

where we have neglected $O(\unicode[STIX]{x1D706})$ terms in accordance with the above-mentioned assumptions about the contact line. Equation (4.3) is of the form $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}^{3}\tilde{\unicode[STIX]{x1D6F6}}=A/(\unicode[STIX]{x1D709}^{2}+1)$ , where $A$ is independent of $\unicode[STIX]{x1D709}$ and has been previously encountered in the context of spreading in two dimensions (see e.g. Eggers Reference Eggers2004, Savva & Kalliadasis Reference Savva and Kalliadasis2011). Its asymptotic structure is investigated with conditions $\tilde{\unicode[STIX]{x1D6F6}}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}\tilde{\unicode[STIX]{x1D6F6}}=0$ at $\unicode[STIX]{x1D709}=0$ including $\tilde{\unicode[STIX]{x1D6F6}}/\unicode[STIX]{x1D709}^{2}\rightarrow 0$ as $\unicode[STIX]{x1D709}\rightarrow \infty$ , which ensures compatibility with the behaviour of the outer solution as the contact line is approached. Although the solution to this problem is easily obtained analytically, we are interested in the behaviour of $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}\tilde{\unicode[STIX]{x1D6F6}}$ as we approach the outer region. In the end, we find

(4.4) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}\tilde{\unicode[STIX]{x1D6F6}}\sim \frac{\unicode[STIX]{x2202}_{t}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D742}}{\unicode[STIX]{x1D717}_{\ast }^{3}}\ln (\text{e}\unicode[STIX]{x1D709})\quad \text{as }\unicode[STIX]{x1D709}\rightarrow \infty ,\end{eqnarray}$$

where we have retained the two leading-order terms in the asymptotic expansion. Clearly, $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h=-\unicode[STIX]{x1D717}_{\ast }\unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}h$ , which allows us to write down the behaviour of the normal derivative of the inner solution as the outer region is approached as

(4.5) $$\begin{eqnarray}-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{in}\sim \unicode[STIX]{x1D717}_{\ast }+\frac{\unicode[STIX]{x2202}_{t}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D742}}{\unicode[STIX]{x1D717}_{\ast }^{2}}\ln \frac{\text{e}\unicode[STIX]{x1D717}_{\ast }(\boldsymbol{c}-\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D742}}{\unicode[STIX]{x1D706}}\quad \text{as }\frac{(\boldsymbol{c}-\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D742}}{\unicode[STIX]{x1D706}}\rightarrow \infty ,\end{eqnarray}$$

where $h_{in}$ denotes the inner solution. Equation (4.5) is a generalisation of the inner-region behaviour obtained in two-dimensional geometries, noting that its derivation relies only on assuming quasi-static dynamics and weakly varying contact lines. Casting this equation in terms of the $r$ and $\unicode[STIX]{x1D719}$ variables, we obtain

(4.6) $$\begin{eqnarray}-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{in}\sim \unicode[STIX]{x1D717}_{\ast }+\frac{\unicode[STIX]{x2202}_{t}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D742}}{\unicode[STIX]{x1D717}_{\ast }^{2}}\ln \frac{\text{e}a^{2}\unicode[STIX]{x1D717}_{\ast }(1-r)}{\unicode[STIX]{x1D706}\sqrt{(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}a)^{2}+a^{2}}}\quad \text{as }\frac{a(1-r)}{\unicode[STIX]{x1D706}}\rightarrow \infty .\end{eqnarray}$$

If we further assume that the contact line is nearly circular, so that only linear terms in $b_{m}$ are retained, we use (3.2) and (3.6) to expand $\unicode[STIX]{x2202}_{t}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D742}$ according to

(4.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D742} & = & \displaystyle {\dot{b}}_{0}+\left({\dot{b}}_{1}-\frac{b_{2}}{b_{0}}{\dot{b}}_{1}^{\ast }\right)\text{e}^{\text{i}\unicode[STIX]{x1D719}}+\left({\dot{b}}_{2}-\frac{3b_{3}}{2b_{0}}{\dot{b}}_{1}^{\ast }\right)\text{e}^{2\text{i}\unicode[STIX]{x1D719}}\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{m=3}^{M}\left({\dot{b}}_{m}+\frac{(m-1)b_{m-1}}{2b_{0}}{\dot{b}}_{1}-\frac{(m+1)b_{m+1}}{2b_{0}}{\dot{b}}_{1}^{\ast }\right)\text{e}^{\text{i}m\unicode[STIX]{x1D719}}+\cdots \,,\end{eqnarray}$$

where the dots indicate higher-order corrections which we omit. Here we defined

(4.8) $$\begin{eqnarray}b_{m}(t)=\left\{\begin{array}{@{}ll@{}}x_{c}(t)-\text{i}y_{c}(t), & m=1,\\ a_{m}(t), & m\neq 1,\end{array}\right.\end{eqnarray}$$

$b_{M+1}=0$ and ${\dot{b}}_{1}^{\ast }$ denotes the complex conjugate of ${\dot{b}}_{1}$ , which appears in (4.7) due to the presence of ${\dot{x}}_{c}\sin \unicode[STIX]{x1D719}-{\dot{y}}_{c}\cos \unicode[STIX]{x1D719}=({\dot{b}}_{1}\text{e}^{\text{i}\unicode[STIX]{x1D719}}-{\dot{b}}_{1}^{\ast }\text{e}^{-\text{i}\unicode[STIX]{x1D719}})/(2\text{i})$ in $\unicode[STIX]{x2202}_{t}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D742}$ ; see also (3.7). Note that $b_{1}$ is defined as above, because the term ${\dot{x}}_{c}\cos \unicode[STIX]{x1D719}+{\dot{y}}_{c}\sin \unicode[STIX]{x1D719}=\text{Re}({\dot{b}}_{1}\text{e}^{\text{i}\unicode[STIX]{x1D719}})$ in (3.3) appears in the orders at which matching is performed.

Based on the expansion above, we see that the logarithmically singular terms in (4.6) do not contain $O({\dot{b}}_{0}b_{m})$ terms, thus justifying their omission in the outer-region analysis that follows; had these terms been included, they would have merely contributed a higher-order correction to the evolution of the harmonics of the contact line. Note also that (4.7) contains $O({\dot{b}}_{1}b_{m})$ and $O({\dot{b}}_{1}^{\ast }b_{m})$ terms which are nevertheless generally smaller compared to $O({\dot{b}}_{m})$ terms since we assumed that $b_{m}\ll 1$ . Hence, the corrections that result by including such terms in a typical situation are expected to be of lesser importance, compared to, say, the more restrictive assumption that the shape of the contact line is nearly circular. However, we shall see later on that there can also exist situations where the contributions of these terms need to be taken into account, due to their tendency to strengthen the coupling between $b_{m}$ and the centroid motion. Although in what follows we will omit $O({\dot{b}}_{1}b_{m})$ and $O({\dot{b}}_{1}^{\ast }b_{m})$ terms, it is worthwhile to stress that we will briefly revisit our analysis in § 6.4 to deduce a more involved system of equations which includes these terms, mostly to highlight their effect when we have a droplet undergoing a long excursion from its original position. Hence, neglecting $O({\dot{b}}_{1}b_{m})$ and $O({\dot{b}}_{1}^{\ast }b_{m})$ terms for now, (4.7) simplifies further to

(4.9) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{c}\boldsymbol{\cdot }\unicode[STIX]{x1D742}=\mathop{\sum }_{m=0}^{M}{\dot{b}}_{m}\text{e}^{\text{i}m\unicode[STIX]{x1D719}}+\cdots\end{eqnarray}$$

and (4.6) reduces to

(4.10) $$\begin{eqnarray}-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{in}\sim \unicode[STIX]{x1D717}_{\ast }+\frac{1}{\unicode[STIX]{x1D717}_{\ast }^{2}}\ln \left(\text{e}b_{0}\bar{\unicode[STIX]{x1D717}}_{\ast }\frac{1-r}{\unicode[STIX]{x1D706}}\right)\mathop{\sum }_{m=1}^{M}{\dot{b}}_{m}\text{e}^{\text{i}m\unicode[STIX]{x1D719}}\quad \text{as }\frac{b_{0}(1-r)}{\unicode[STIX]{x1D706}}\rightarrow \infty ,\end{eqnarray}$$

where

(4.11) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D717}}_{\ast }=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D717}(x_{c}+a\cos \unicode[STIX]{x1D719},y_{c}+a\sin \unicode[STIX]{x1D719})\,\text{d}\unicode[STIX]{x1D719}\end{eqnarray}$$

denotes the average value of $\unicode[STIX]{x1D717}_{\ast }$ along the contact line. Equation (4.10) needs to be matched with the corresponding outer solution as the contact line region is approached, which is determined next.

4.2 Outer-region dynamics

Since slip effects manifest themselves only in the vicinity of the contact line, we neglect slip in the bulk of the droplet. Because of this, the outer solution cannot possibly satisfy (3.10) and we only require that the solution vanishes along the contact line, (3.9), the constraint (3.14), as well as the appropriate pole conditions at the origin of the polar coordinate system. Denoting the solution in the outer region as $h_{out}$ , (3.3) simplifies to

(4.12) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x2202}_{t}h_{out}-\frac{1}{a}\left[{\dot{x}}_{c}\cos \unicode[STIX]{x1D719}+{\dot{y}}_{c}\sin \unicode[STIX]{x1D719}+r\unicode[STIX]{x2202}_{t}a+\frac{\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}a}{a}({\dot{x}}_{c}\sin \unicode[STIX]{x1D719}-{\dot{y}}_{c}\cos \unicode[STIX]{x1D719})\right]\unicode[STIX]{x2202}_{r}h_{out}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\frac{{\dot{x}}_{c}\sin \unicode[STIX]{x1D719}-{\dot{y}}_{c}\cos \unicode[STIX]{x1D719}}{ar}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}h_{out}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h_{out}^{3}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FB}^{2}h_{out})=0.\end{eqnarray}$$

The assumption of slow dynamics allows us to introduce the quasi-static expansion

(4.13) $$\begin{eqnarray}h_{out}(r,\unicode[STIX]{x1D719},t)=h_{0}(r,\unicode[STIX]{x1D719},t)+h_{1}(r,\unicode[STIX]{x1D719},t)+\cdots \,,\end{eqnarray}$$

where $h_{0}\gg h_{1}$ . This means that the capillary forces dominate the viscous ones in the bulk and the leading-order term, $h_{0}$ , is a quasi-equilibrium solution cast in the form

(4.14) $$\begin{eqnarray}h_{0}(r,\unicode[STIX]{x1D719},t)=h_{0,0}(r,b_{0})+\mathop{\sum }_{m=2}^{M}b_{m}(t)h_{0,m}(r,b_{0})\text{e}^{\text{i}m\unicode[STIX]{x1D719}}+\cdots\end{eqnarray}$$

and $h_{1}$ is a correction associated with the motion of the contact line, written as

(4.15) $$\begin{eqnarray}h_{1}(r,\unicode[STIX]{x1D719},t)=\mathop{\sum }_{m=0}^{M}{\dot{b}}_{m}(t)h_{1,m}(r,b_{0})\text{e}^{\text{i}m\unicode[STIX]{x1D719}}+\cdots \,,\end{eqnarray}$$

noting that in (4.14) the first harmonic is omitted since $a_{1}(t)\equiv 0$ . Moreover, we omitted $O({\dot{b}}_{0}b_{m})$ and $O({\dot{b}}_{1}b_{m\pm 1})$ terms in (4.15), consistently with (4.10). We also emphasise the quasi-static treatment of the dynamics, given that (4.14) and (4.15) do not depend explicitly on $t$ , but through the functions $b_{m}$ and their time derivatives.

At leading order, we have zero fluid velocity within the droplet, so that the first term in (4.13), $h_{0}$ , satisfies the Poisson equation with a homogeneous Dirichlet condition, namely

(4.16a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}h_{0}(r,\unicode[STIX]{x1D719},t)=q(t),\quad 0<r<1,0\leqslant \unicode[STIX]{x1D719}<2\unicode[STIX]{x03C0}, & \displaystyle\end{eqnarray}$$

(4.16b ) $$\begin{eqnarray}\displaystyle & \displaystyle h_{0}(1,\unicode[STIX]{x1D719},t)=0, & \displaystyle\end{eqnarray}$$

$q(t)$

$h_{0}$

(4.17) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}(\cdot )=\frac{1}{b_{0}^{2}r}\unicode[STIX]{x2202}_{r}[r\unicode[STIX]{x2202}_{r}(\cdot )]+\frac{1}{b_{0}^{2}r^{2}}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}^{2}(\cdot )-\frac{1}{b_{0}^{3}r}\mathop{\sum }_{m=2}^{M}b_{m}\{2\unicode[STIX]{x2202}_{r}[r\unicode[STIX]{x2202}_{r}(\cdot )]-m^{2}\unicode[STIX]{x2202}_{r}(\cdot )\}\text{e}^{\text{i}m\unicode[STIX]{x1D719}}+\cdots \,.\end{eqnarray}$$

Plugging (4.14) into (4.16a ) and matching the $\text{e}^{\text{i}m\unicode[STIX]{x1D719}}$ terms in the series expansion yields

(4.18) $$\begin{eqnarray}{\mathcal{L}}_{0}h_{0,0}=qb_{0}^{2}\end{eqnarray}$$

for the axisymmetric component ( $m=0$ ), whereas the higher harmonics ( $m>1$ ) satisfy

(4.19) $$\begin{eqnarray}{\mathcal{L}}_{m}h_{0,m}=\frac{1}{b_{0}r}\unicode[STIX]{x2202}_{r}(2r\unicode[STIX]{x2202}_{r}h_{0,0}-m^{2}h_{0,0}),\end{eqnarray}$$

where

(4.20) $$\begin{eqnarray}{\mathcal{L}}_{m}(\cdot )=\frac{1}{r}\unicode[STIX]{x2202}_{r}[r\unicode[STIX]{x2202}_{r}(\cdot )]-\frac{m^{2}}{r^{2}}(\cdot ),\quad m\geqslant 0.\end{eqnarray}$$

Solving (4.18) and (4.19) with the appropriate conditions is rather straightforward. Equation (4.18) is solved with $h_{0,0}(1,t)=0$ and (3.14), giving $q=-16/b_{0}^{4}$ and

(4.21) $$\begin{eqnarray}h_{0,0}=\frac{4}{b_{0}^{2}}(1-r^{2}).\end{eqnarray}$$

Equation (4.19) is clearly an equi-dimensional equation and has the solution

(4.22) $$\begin{eqnarray}h_{0,m}=\frac{8}{b_{0}^{3}}(r^{m}-r^{2})\quad \text{for }m\geqslant 2,\end{eqnarray}$$

found by requiring $h_{0,m}(1,t)=0$ and $h_{0,m}$ to be finite as $r\rightarrow 0$ . Thus, for given $a(\unicode[STIX]{x1D719},t)$ , (4.14) with (4.21) and (4.22) form an approximate solution to (4.16).

To determine $h_{1}$ using the series expansion (4.15), we need to retain terms that are linear in ${\dot{b}}_{m}$ in (4.12), with $m=0,1,2,\ldots ,M$ . We find that $h_{1}$ must satisfy

(4.23) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}h_{0}-\frac{1}{b_{0}}[\text{Re}({\dot{b}}_{1}\text{e}^{\text{i}\unicode[STIX]{x1D719}})+r\unicode[STIX]{x2202}_{t}a]\unicode[STIX]{x2202}_{r}h_{0,0}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h_{0,0}^{3}\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FB}^{2}h_{1})+\cdots =0.\end{eqnarray}$$

Using the approximation

(4.24) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}h_{0}(r,\unicode[STIX]{x1D719},t)\approx {\dot{b}}_{0}\unicode[STIX]{x2202}_{b_{0}}h_{0,0}(r,b_{0})+\mathop{\sum }_{m=2}^{M}{\dot{b}}_{m}(t)h_{0,m}(r,b_{0})\text{e}^{\text{i}m\unicode[STIX]{x1D719}}+\cdots\end{eqnarray}$$

in (4.23) and collecting $\text{e}^{\text{i}m\unicode[STIX]{x1D719}}$ terms with $m\geqslant 0$ gives the following set of decoupled linear differential equations,

(4.25a ) $$\begin{eqnarray}\displaystyle & \displaystyle 8r(2r^{2}-1)+\unicode[STIX]{x2202}_{r}[r(1-r^{2})^{3}\unicode[STIX]{x2202}_{r}P_{0}]=0,\quad m=0, & \displaystyle\end{eqnarray}$$

(4.25b ) $$\begin{eqnarray}\displaystyle & \displaystyle 8r^{m+2}+r\unicode[STIX]{x2202}_{r}[r(1-r^{2})^{3}\unicode[STIX]{x2202}_{r}P_{m}]-m^{2}(1-r^{2})^{3}P_{m}=0,\quad m\geqslant 1, & \displaystyle\end{eqnarray}$$

(4.26) $$\begin{eqnarray}P_{m}=\frac{64}{b_{0}^{7}}{\mathcal{L}}_{m}h_{1,m}.\end{eqnarray}$$

Equations (4.25) need to be considered with the appropriate pole conditions, together with a vanishing mass flux along the contact line and $h_{1,m}(1,t)=0$ . As in the inner region, we only need to obtain the asymptotic behaviour of $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{out}$ in the overlap region so that matching with (4.10) can be carried out. It is straightforward to see that compatibility with (4.10) requires that

(4.27) $$\begin{eqnarray}P_{m}\sim \frac{1}{1-r}\quad \text{as }r\rightarrow 1^{-},\end{eqnarray}$$

so that we can write

(4.28) $$\begin{eqnarray}-\unicode[STIX]{x2202}_{r}h_{1,m}\sim \frac{b_{0}^{7}}{64}[\ln (1-r)+\unicode[STIX]{x1D6FD}(m)]+\cdots \quad \text{as }r\rightarrow 1^{-},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}(m)$ is a function of $m$ that needs to be determined in order to complete the analysis of the outer region.

The asymptotic behaviour of the case when $m=0$ can be found by integrating (4.25a ) once, which gives

(4.29) $$\begin{eqnarray}\unicode[STIX]{x2202}_{r}P_{0}=\frac{4r}{(1-r^{2})^{2}}.\end{eqnarray}$$

Then, by multiplying each side of (4.29) by $r^{2}(1-r^{2})/2$ , making use of (4.26) and integrating from 0 to $1-\unicode[STIX]{x1D716}$ with $0<\unicode[STIX]{x1D716}\ll 1$ , integration by parts gives

(4.30) $$\begin{eqnarray}1-\ln (\unicode[STIX]{x1D716})-\unicode[STIX]{x1D6FD}(0)+O(\unicode[STIX]{x1D716}\ln \unicode[STIX]{x1D716})=\int _{0}^{1-\unicode[STIX]{x1D716}}\frac{2r^{3}}{1-r^{2}}\,\text{d}r,\end{eqnarray}$$

after applying the aforementioned conditions on $h_{1,m}$ , its asymptotics, (4.28), and the condition $\int _{0}^{1}rh_{1,0}\,\text{d}r=0$ . Using

(4.31) $$\begin{eqnarray}\ln \unicode[STIX]{x1D716}=-\int _{0}^{1-\unicode[STIX]{x1D716}}\frac{1}{1-r}\,\text{d}r\end{eqnarray}$$

and rearranging terms in (4.30) gives the value of $\unicode[STIX]{x1D6FD}(0)$ ,

(4.32) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}(0)=1+\int _{0}^{1}\frac{1}{1-r}\left(1-2\frac{r^{3}}{1+r}\right)\,\text{d}r=2+\ln 2,\end{eqnarray}$$

as we take the limit $\unicode[STIX]{x1D716}\rightarrow 0^{+}$ . This result is identical to Hocking’s work on axisymmetric droplets in the special case when gravity effects become negligible (Hocking Reference Hocking1983).

Contrary to the $m=0$ case, the values of $\unicode[STIX]{x1D6FD}(m)$ for $m\geqslant 1$ can only be specified in terms of integrals of transcendental functions. First, we determine the solution to (4.25b ) as

(4.33) $$\begin{eqnarray}P_{m}=\frac{4r^{m}}{(m+4)(1-r^{2})^{2}}\left[\frac{f_{m}(r^{2})}{f_{m}(1)}-1\right],\end{eqnarray}$$

where $f_{m}(r)$ denotes the Gauss hypergeometric function

(4.34) $$\begin{eqnarray}f_{m}(r)=_{2}F_{1}\left(\frac{m-1+\sqrt{m^{2}+9}}{2},\frac{m-1-\sqrt{m^{2}+9}}{2};m+1;r\right),\end{eqnarray}$$

noting the agreement of the asymptotics of (4.33) as $r\rightarrow 1^{-}$ with (4.27), as expected. The function $f_{m}$ is bounded for $0\leqslant r\leqslant 1$ ; as $m\rightarrow \infty$ , we have that $f_{\infty }\rightarrow \sqrt{1-r}$ , so that (4.27) can no longer be satisfied.

To proceed, we substitute (4.33) in (4.26), multiply both sides by $r^{m+1}$ and integrate, as before, from 0 to $1-\unicode[STIX]{x1D716}$ . After successive integrations by parts, we find

(4.35) $$\begin{eqnarray}\int _{0}^{1-\unicode[STIX]{x1D716}}\frac{4r^{2m+1}}{(m+4)(1-r^{2})^{2}}\left[\frac{f_{m}(r^{2})}{f_{m}(1)}-1\right]\,\text{d}r=-\text{ln}\,\unicode[STIX]{x1D716}-\unicode[STIX]{x1D6FD}(m)+O(\unicode[STIX]{x1D716}\ln \unicode[STIX]{x1D716}),\end{eqnarray}$$

where we have made use of the conditions on $h_{1,m}$ and its asymptotics, (4.28). Clearly, the integrand on the left-hand side has a logarithmic singularity as $\unicode[STIX]{x1D716}\rightarrow 0$ which is balanced by the logarithmic singularity of the right-hand side. Using (4.31) and rearranging allows us to remove the singular terms thus obtaining the following integral expression for $\unicode[STIX]{x1D6FD}(m)$ :

(4.36) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}(m)=\int _{0}^{1}\left\{\frac{1}{1-r}-\frac{4r^{2m+1}}{(m+4)(1-r^{2})^{2}}\left[\frac{f_{m}(r^{2})}{f_{m}(1)}-1\right]\right\}\,\text{d}r,\quad m\geqslant 1.\end{eqnarray}$$

Figure 2.

Plot of $\unicode[STIX]{x1D6FD}$ as a function of $m\geqslant 1$ . (a) The values of $\unicode[STIX]{x1D6FD}$ for the first 40 harmonics. (b) Plot on a logarithmic scale for $m$ (solid curve), showing also the empirically deduced large- $m$ asymptotics of $\unicode[STIX]{x1D6FD}$ (dashed curve). In both (a) and (b) the discrete values of $\unicode[STIX]{x1D6FD}$ are joined for clarity.

This expression cannot be simplified further and needs to be evaluated numerically. We performed this computation once and the values of $\unicode[STIX]{x1D6FD}(m)$ are stored and retrieved in subsequent simulations as necessary. Some implementation details on how $\unicode[STIX]{x1D6FD}(m)$ and the associated Gauss hypergeometric function are computed can be found in appendix A. Figure 2 shows the result of this computation, where we see a gradual divergence of $\unicode[STIX]{x1D6FD}$ as $m\rightarrow \infty$ , which is a manifestation of the breakdown of the theory when $m$ becomes sufficiently large. The divergence of $\unicode[STIX]{x1D6FD}(m)$ is logarithmic and we see that $\unicode[STIX]{x1D6FD}(m)\sim \ln (2\text{e}m)$ fits the numerical data excellently as $m\rightarrow \infty$ , noting that this result was obtained empirically. We will return to the behaviour of $\unicode[STIX]{x1D6FD}(m)$ for large values of $m$ after matching is performed, as this will allow us to set the order of magnitude of $M$ , the number of harmonics retained in the truncated Fourier series, (4.1).

Knowing the asymptotics of $\unicode[STIX]{x2202}_{r}h_{1,m}$ allows us to obtain the normal derivative of the free surface, $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{out}$ as the contact line is approached. We have

(4.37) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{out}=\frac{1}{b_{0}}\left(1-\mathop{\sum }_{m=2}^{M}\frac{b_{m}}{b_{0}}\text{e}^{\text{i}m\unicode[STIX]{x1D719}}\right)\unicode[STIX]{x2202}_{r}(h_{0}+h_{1})+\cdots \,,\end{eqnarray}$$

which needs to be evaluated as $r\rightarrow 1^{-}$ . We have

(4.38) $$\begin{eqnarray}-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{out}\sim \unicode[STIX]{x1D703}+\frac{b_{0}^{6}}{64}\mathop{\sum }_{m=0}^{M}{\dot{b}}_{m}[\ln (1-r)+\unicode[STIX]{x1D6FD}(m)]\text{e}^{\text{i}m\unicode[STIX]{x1D719}}\quad \text{as }r\rightarrow 1^{-},\end{eqnarray}$$

where $\unicode[STIX]{x1D703}$ is the (approximate) apparent contact angle

(4.39) $$\begin{eqnarray}\unicode[STIX]{x1D703}=\frac{8}{b_{0}^{3}}\left[1+\mathop{\sum }_{m=2}^{M}(1-m)\frac{b_{m}}{b_{0}}\text{e}^{\text{i}m\unicode[STIX]{x1D719}}\right]+\cdots \,.\end{eqnarray}$$

Having obtained the asymptotics of the inner, (4.10), and outer regions, (4.38), the final step of the analysis is to match the two solutions.

4.3 Matching

Just as in the corresponding analysis in two dimensions by Vellingiri et al. (Reference Vellingiri, Savva and Kalliadasis2011), matching is performed with the cubes of the asymptotics of the inner and outer normal derivatives, referring the interested reader to the work of Sibley, Nold & Kalliadasis (Reference Sibley, Nold and Kalliadasis2015a ) which offers a generalisation of how matching can be formally performed for cases when the cubes of the slopes do not match like they do here.

Hence, by considering the third powers of (4.10) and (4.38) we can eliminate the $r$ -dependent logarithmic terms, and by requiring that the remaining terms match we find

(4.40) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D703}^{3}-\unicode[STIX]{x1D717}_{\ast }^{3}}{3}=\mathop{\sum }_{m=0}^{M}{\dot{b}}_{m}\left[\ln \frac{\text{e}b_{0}\bar{\unicode[STIX]{x1D717}}_{\ast }}{\unicode[STIX]{x1D706}}-\unicode[STIX]{x1D6FD}(m)\right]\text{e}^{\text{i}m\unicode[STIX]{x1D719}},\end{eqnarray}$$

which is reminiscent of the Voinov–Cox–Hocking law and offers an extension of the previous results in two dimensions (Vellingiri et al. Reference Vellingiri, Savva and Kalliadasis2011) to the three-dimensional setting. It constitutes a simplified problem for the slow motion of a drop, since, as we shall see, the effort required to solve (4.40) is considerably less compared to solving (2.12). Using (4.40), we can write down the evolution equations for $b_{m}$ , by collecting the corresponding Fourier coefficients of the left-hand side. Hence, we can write

(4.41) $$\begin{eqnarray}{\dot{b}}_{m}=\left(\ln \frac{\text{e}b_{0}\bar{\unicode[STIX]{x1D717}}_{\ast }}{\unicode[STIX]{x1D706}}-\unicode[STIX]{x1D6FD}(m)\right)^{-1}\int _{0}^{2\unicode[STIX]{x03C0}}\frac{\unicode[STIX]{x1D703}^{3}-\unicode[STIX]{x1D717}_{\ast }^{3}}{3\unicode[STIX]{x03C0}}\text{e}^{-\text{i}m\unicode[STIX]{x1D719}}\,\text{d}\unicode[STIX]{x1D719},\end{eqnarray}$$

where the real and imaginary parts of (4.41) give, respectively, ordinary differential equations (ODEs) for the cosine and sine coefficients of the truncated Fourier series in (4.1), except the case when $m=1$ which gives ODEs for $x_{c}(t)$ and $y_{c}(t)$ .

Note that in (4.40) and (4.41) we have not written $\unicode[STIX]{x1D703}^{3}$ explicitly in its linearised form for illustration purposes, because the analogy with the two-dimensional geometry is more readily seen this way and also because we will refer to (4.41) later on. We also chose not to assume anything about the strength of the heterogeneities, because, under certain conditions, it is possible to have a nearly circular contact line even when heterogeneities are strong. If we do use (4.39) and retain only the linear terms in $b_{m}$ we find

(4.42a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{b}}_{0}=\frac{1}{3}\left(\ln \frac{b_{0}\bar{\unicode[STIX]{x1D717}}_{\ast }}{2\unicode[STIX]{x1D706}\text{e}}\right)^{-1}\left[\frac{512}{b_{0}^{9}}-\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D717}_{\ast }^{3}\,\text{d}\unicode[STIX]{x1D719}\right], & \displaystyle\end{eqnarray}$$

(4.42b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{x}}_{c}=-\frac{1}{3\unicode[STIX]{x03C0}}\left(\ln \frac{\text{e}b_{0}\bar{\unicode[STIX]{x1D717}}_{\ast }}{\unicode[STIX]{x1D706}}-\unicode[STIX]{x1D6FD}(1)\right)^{-1}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D717}_{\ast }^{3}\cos \unicode[STIX]{x1D719}\,\text{d}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$

(4.42c ) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{y}}_{c}=-\frac{1}{3\unicode[STIX]{x03C0}}\left(\ln \frac{\text{e}b_{0}\bar{\unicode[STIX]{x1D717}}_{\ast }}{\unicode[STIX]{x1D706}}-\unicode[STIX]{x1D6FD}(1)\right)^{-1}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D717}_{\ast }^{3}\sin \unicode[STIX]{x1D719}\,\text{d}\unicode[STIX]{x1D719}, & \displaystyle\end{eqnarray}$$

(4.42d ) $$\begin{eqnarray}\displaystyle & \displaystyle {\dot{b}}_{m}=\left(\ln \frac{\text{e}b_{0}\bar{\unicode[STIX]{x1D717}}_{\ast }}{\unicode[STIX]{x1D706}}-\unicode[STIX]{x1D6FD}(m)\right)^{-1}\left[\frac{512(1-m)b_{m}}{b_{0}^{10}}-\frac{1}{3\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D717}_{\ast }^{3}\text{e}^{-\text{i}m\unicode[STIX]{x1D719}}\,\text{d}\unicode[STIX]{x1D719}\right],\quad m\geqslant 2.\quad \qquad & \displaystyle\end{eqnarray}$$

$\unicode[STIX]{x1D717}_{\ast }\equiv 1$

${\dot{b}}_{m}$

$b_{m}$

The last statement above is valid provided that the first term on the right-hand side of (4.42d ) is always positive for $m$ sufficiently large, namely

(4.43) $$\begin{eqnarray}\ln \frac{\text{e}b_{0}\bar{\unicode[STIX]{x1D717}}_{\ast }}{\unicode[STIX]{x1D706}}-\unicode[STIX]{x1D6FD}(m)>0.\end{eqnarray}$$

This inequality, however, cannot be true for all values of $m$ . The plots in figure 2 show that $\unicode[STIX]{x1D6FD}(m)$ grows logarithmically with $m$ , which justifies why we considered a truncated Fourier series for $a(\unicode[STIX]{x1D719},t)$ ; see (4.1). To specify approximately how many terms we must retain, take $b_{0}=2$ (the equilibrium angle when $\unicode[STIX]{x1D717}_{\ast }\equiv 1$ ) and use the empirically obtained asymptotics of $\unicode[STIX]{x1D6FD}\sim \ln (2\text{e}m)$ for $m\gg 1$ , which gives the criterion $M\ll \unicode[STIX]{x1D706}^{-1}$ for the analysis to hold for all terms in (4.1). This, however, is not particularly restrictive for our purposes, because $\unicode[STIX]{x1D706}$ is generally a small parameter.

It is important at this stage to contrast our analysis with the work of Greenspan (Reference Greenspan1978) and its extension by Greenspan & McCay (Reference Greenspan and McCay1981). Greenspan (Reference Greenspan1978) used (1.1) and the volume constraint in place of (2.12b ) and (2.12d ). However, enforcing the volume constraint at all times is a result of having zero mass flux along $C$ , which, in turn, is also equivalent to the kinematic condition (2.12b ) along $C$ (Young Reference Young, Brown and Davis1994). Hence, the condition (1.1) utilised by Greenspan (Reference Greenspan1978) can be viewed as a replacement for our contact angle condition (2.12d ) which become equivalent to each other as $k\rightarrow \infty$ . Yet, Greenspan (Reference Greenspan1978) treated (1.1) as a kinematic condition rather than a contact angle condition. This different treatment of (1.1) allows the extraction of evolution equations for the contact line just from the leading-order outer solution, but this type of analysis completely ignores the influence of slip, which is precisely the effect responsible for the contact line motion (see also Hocking Reference Hocking1992, for further comments on this point).

As a result, the slip length does not appear anywhere in the set of equations obtained by Greenspan (Reference Greenspan1978) and Greenspan & McCay (Reference Greenspan and McCay1981). Moreover, it can be seen that the approach to equilibrium is algebraic rather than exponential. Had this condition been utilised properly (see e.g. Hocking Reference Hocking1992, Sibley et al. Reference Sibley, Nold, Savva and Kalliadasis2015b ) a treatment via matched asymptotics would have been required like here. In the end, however, as Sibley et al. (Reference Sibley, Nold, Savva and Kalliadasis2015b ) have shown, the equations obtained by matching can be made equivalent to each other with a suitable rescaling of $\unicode[STIX]{x1D706}$ . It is noteworthy that Hocking (Reference Hocking1992) criticised conditions of the form (1.1) arguing against their usefulness given that the dynamics can be formally mapped to the dynamics arising from using (2.12d ), which does not depend on additional material constants, such as $k$ in (1.1).

Apart from these key differences in the methodology, we should emphasise that our approach is general and can deal with droplets undergoing any type of motion. In contrast, Greenspan (Reference Greenspan1978) and Greenspan & McCay (Reference Greenspan and McCay1981) considered separately the cases where the local contact angle exhibits linear spatial variations and when the droplet does not move from its original position so that the contact line can be described in a fixed polar coordinate system at all times. In contrast, (4.40) can be used to treat both cases just as easily.

5.1 Reduced problem

The reduced problem, (4.40) and its linearised version (4.42) is generally a non-stiff system of ODEs, where we have explicit expressions for each harmonic, but coupling of each component can arise through the computation of $\bar{\unicode[STIX]{x1D717}}_{\ast }$ and the harmonics of $\unicode[STIX]{x1D717}_{\ast }^{3}$ . In general, we use $M=150$ , resulting in 301 unknowns (one axisymmetric component and 150 Fourier cosine and 150 sine modes). The input variables at each time step, namely $a_{m}(t)$ , $x_{c}(t)$ and $y_{c}(t)$ , are used to construct the contact line using the inverse fast Fourier transform, which allows us to obtain $\unicode[STIX]{x1D717}_{\ast }^{3}$ . Then $\unicode[STIX]{x1D717}_{\ast }^{3}$ is decomposed into its Fourier modes using the fast Fourier transform, which are then used in either (4.40) or (4.42) to extract the values for ${\dot{b}}_{m}$ and perform the time stepping using the routine ode15s in MATLAB. This is a rather straightforward implementation and a typical simulation takes about one minute or less on a standard personal computer.

5.2 Full problem

The full problem requires significantly more effort, both in terms of implementation and actual computing time, as it amounts to solving a set of coupled stiff PDEs for $h(r,\unicode[STIX]{x1D719},t)$ and $a(\unicode[STIX]{x1D719},t)$ as well as a set of ODEs for the coordinates of the centroid, $(x_{c}(t),y_{c}(t))$ . The numerical issues due to stiffness worsen as $\unicode[STIX]{x1D706}\rightarrow 0$ , because the boundary layers in $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h$ become narrower, thus requiring more mesh points near the contact line to be resolved.

Such difficulties need to be overcome in order to perform the computations accurately and efficiently, and to properly compare the results of our analysis with the predictions of the full problem. Here we opted for the pseudospectral collocation method, which has been successfully used in the past by the present authors to deal with contact line problems in two dimensions (see e.g. Savva & Kalliadasis Reference Savva and Kalliadasis2009, Reference Savva and Kalliadasis2011, Reference Savva and Kalliadasis2012, Reference Savva and Kalliadasis2013, Reference Savva and Kalliadasis2014). An attractive feature of this method is that it naturally clusters collocation points near the contact line, which allows us to capture boundary-layer dynamics with comparatively fewer collocation nodes as opposed to, say, a finite difference scheme based on a uniform grid. Although giving all the fine details of our numerical method is beyond the scope of this work, its basic principles are briefly outlined in appendix B.

5.3 Hybrid approach

The hybrid approach provides an attractive compromise between the speed of solving the low-order, reduced model, (4.40), and the accuracy of the full problem, (2.12). It is based on the observation that the apparent contact angle used in (4.40) is only approximately provided from low-order corrections about a circular contact line, obtained by solving Laplace’s equation by a truncated perturbation series; see (4.14).

In essence, the reduced model uses as input the contact line harmonics, from which the harmonics of the apparent contact angle are obtained by solving Laplace’s equation approximately. From the apparent contact angle, the reduced model then gives the rates of change of each harmonic of the contact line. Hence, as we shall see in the next section, the equilibrium obtained in the long-time limit by solving the reduced problem can often be different from that obtained by solving the full problem. This shortcoming of the reduced model can be remedied by solving numerically Laplace’s equation, thus obtaining improved approximations for the apparent contact angle at each time step.

The most appropriate method to solve Laplace’s equation in this setting is arguably the boundary integral method, as it outputs directly the apparent contact angle. This requires significantly less computational effort as opposed to traditional discretisation schemes, which determine the solution everywhere in the domain. This approach is similar to the work of Glasner (Reference Glasner2005), apart from the evolution of the contact line, which was computed by specifying the form of the normal velocity of the contact line as $\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D742}=\unicode[STIX]{x1D703}^{3}-\unicode[STIX]{x1D717}_{\ast }^{3}$ . A notable difference between our work and Glasner’s is that Glasner parametrised the contact line with respect to the arc length, unlike here where the parametrisation is done with respect to the polar angle. Glasner’s parametrisation allows the investigation of cases where the contact line cannot be described by a one-to-one function in polar coordinates. However, we do not require such flexibility since the presented asymptotic theory assumes weak contact line deformations.

Although the problem set-up is nearly identical to the work of Glasner (Reference Glasner2005), for the sake of completeness we reproduce its main results adapted to our notation and providing more details on the implementation in appendix C. Thus, for a given contact line shape, we wish to solve numerically the leading-order outer problem, (4.16). In this quasi-static setting, time is frozen, so $q(t)$ needs to be determined from the volume constraint, (3.14), cast in the form

(5.1) $$\begin{eqnarray}\int _{\unicode[STIX]{x1D6FA}(t)}h_{0}(\boldsymbol{x},t)\,\text{d}\boldsymbol{x}=V,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}(t)$ denotes the wet region of the substrate and $\boldsymbol{x}=(x-x_{c},y-y_{c})$ . Central to the boundary integral formulation is the free-space Green’s function,

(5.2) $$\begin{eqnarray}G(\boldsymbol{x},\boldsymbol{x}_{0})=-\frac{1}{2\unicode[STIX]{x03C0}}\ln |\boldsymbol{x}-\boldsymbol{x}_{0}|,\end{eqnarray}$$

which is the radially symmetric solution to

(5.3) $$\begin{eqnarray}-\unicode[STIX]{x1D6FB}^{2}G(\boldsymbol{x},\boldsymbol{x}_{0})=\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{0})\quad \text{for }\boldsymbol{x},\,\boldsymbol{x}_{0}\in \mathbb{R}^{2},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}$ denotes the Dirac delta function. To proceed, define

(5.4) $$\begin{eqnarray}w(\boldsymbol{x},t)=\frac{h_{0}(\boldsymbol{x},t)}{q(t)}+\frac{|\boldsymbol{x}|^{2}}{4},\end{eqnarray}$$

so that (4.16a ) and (4.16b ) prescribe the following problem for $w(\boldsymbol{x},t)$ :

(5.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FB}^{2}w(\boldsymbol{x},t)=0\quad \text{for }\boldsymbol{x}\in \unicode[STIX]{x1D6FA}(t), & \displaystyle\end{eqnarray}$$

(5.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle w(\boldsymbol{x},t)=\frac{|\boldsymbol{x}|^{2}}{4}\quad \text{for }\boldsymbol{x}\in C(t). & \displaystyle\end{eqnarray}$$

$h_{0}$

(5.6) $$\begin{eqnarray}\int _{C(t)}G(\boldsymbol{x},\boldsymbol{x}_{0})\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}w(\boldsymbol{x},t)\,\text{d}s(\boldsymbol{x})=\frac{|\boldsymbol{x}_{0}|^{2}}{8}+\frac{1}{4}\int _{C(t)}|\boldsymbol{x}|^{2}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}G(\boldsymbol{x},\boldsymbol{x}_{0})\,\text{d}s(\boldsymbol{x}),\end{eqnarray}$$

where $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}w$ and $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}G$ denote the normal derivatives of $w$ and $G$ , respectively, and $\boldsymbol{x}_{0}$ is a point on the contact line, $C(t)$ . Equation (5.6) can be solved for the unknown values of $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}w$ along $C(t)$ , from which we may determine $q(t)$ and then specify the apparent contact angle $\unicode[STIX]{x1D703}=-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{0}|_{C(t)}$ via (5.4), namely

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D703}(\boldsymbol{x},t)=-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}h_{0}(\boldsymbol{x},t)=q(t)\left[\frac{\boldsymbol{x}\boldsymbol{\cdot }\unicode[STIX]{x1D742}}{2}-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}w(\boldsymbol{x},t)\right],\quad \boldsymbol{x}\in C(t).\end{eqnarray}$$

Thus, to determine $q(t)$ , we follow Glasner (Reference Glasner2005) and consider (5.1) as

(5.8) $$\begin{eqnarray}V=\frac{1}{4}\int _{\unicode[STIX]{x1D6FA}(t)}h_{0}(\boldsymbol{x},t)\unicode[STIX]{x1D6FB}^{2}|\boldsymbol{x}|^{2}\,\text{d}\boldsymbol{x}.\end{eqnarray}$$

Then, utilising Green’s second integral identity, followed by the identity $\unicode[STIX]{x1D735}\boldsymbol{\cdot }(|\boldsymbol{x}|^{2}\boldsymbol{x})=4|\boldsymbol{x}|^{2}$ , (5.7) and the divergence theorem, (5.8) becomes

(5.9) $$\begin{eqnarray}V=\frac{q(t)}{4}\int _{C(t)}|\boldsymbol{x}|^{2}\left(\frac{\boldsymbol{x}\boldsymbol{\cdot }\unicode[STIX]{x1D742}}{4}-\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}w\right)\,\text{d}s(\boldsymbol{x}),\end{eqnarray}$$

which allows us to compute $q(t)$ once $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}w$ is known on $C(t)$ . Equation (5.9) is equivalent to equation (11) in Glasner (Reference Glasner2005), noting that it contains a typo in the first term inside the integral, which should read $|\boldsymbol{x}|^{2}\boldsymbol{x}\boldsymbol{\cdot }\unicode[STIX]{x1D742}/16$ instead of $|\boldsymbol{x}|^{2}\boldsymbol{x}/16$ . Summarising, the boundary integral formulation for obtaining the apparent contact angle involves solving an integral equation, (5.6), which results in a dense linear system for $\unicode[STIX]{x2202}_{\unicode[STIX]{x1D708}}w$ , computing $q(t)$ from (5.9) and $\unicode[STIX]{x1D703}(\boldsymbol{x})$ from (5.7). Further details about how these computations can be performed are given in appendix C. The hybrid analytical–numerical approach is closed by using the outcome of this calculation, $\unicode[STIX]{x1D703}(\boldsymbol{x})$ , in conjunction with the rates of change of the harmonics of the contact line, (4.40). The evolution of the contact line is then computed analogously to the reduced system. Some comparisons with Lacey’s leading-order asymptotic result are also offered in the present study, by solving (3.7), (3.12) and (3.13) with $W=(\unicode[STIX]{x1D717}_{\ast }^{3}-\unicode[STIX]{x1D703}^{3})\sqrt{a^{2}+(\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}a)^{2}}/(3\ln \unicode[STIX]{x1D706})$ .