2.1. Thermodynamics and contact line conditions

Based on (2.4)–(2.7), we now derive appropriate contact line conditions following the principles of non-equilibrium thermodynamics (Ren, Hu & E Reference Ren, Hu and E.2010; Ren & E 2011). We start with the system total energy under the lubrication approximation:

(2.8) \begin{equation} \begin{aligned} \mathcal{E} & = \mathcal{E}_s + \mathcal{E}_b, \\ \mathcal{E}_s & = \frac {1}{{\textit{Ca}}} \int _0^{A(t)}\left (1 + \frac {\delta ^2}{2} \left |\partial _XH\right |^2\right) \mathrm {d} X \\ &\quad {}+ \tilde {\gamma }_1 \int _0^{A(t)}\left (1 + \frac {\delta ^2}{2} \left |\partial _XG\right |^2 \right)\mathrm {d} X + \tilde {\gamma }_2 \int _{A(t)}^{X_r} \left (1 + \frac {\delta ^2}{2} \left |\partial _XG\right |^2\right) \mathrm {d} X, \\ \mathcal{E}_b & = \frac {\tilde {C}_b}{2} \int _0^{X_r} \delta ^2 \left |\partial _{XX} G\right |^2 \mathrm {d} X, \end{aligned} \end{equation}

where $\mathcal{E}_s$ denotes the total interfacial energy, and $\mathcal{E}_b$ is the bending energy of the elastic sheet. To avoid an infinite interfacial energy, here we apply the far-field condition (2.6) at the fixed boundary $X=X_r$ (instead of $X=+\infty$ ) to suppress any energy dissipation at this point. (It is more rigorous to deal with a finite energy, although the derivation in this subsection remains valid, at least formally, with $X_r = \infty$ .)

Next, we compute the energy dissipation rate. For the bending energy, we have

(2.9) \begin{equation} \begin{aligned} \frac {\mathrm {d}\mathcal{E}_b}{\mathrm {d} T} ={} &\delta ^2 \tilde {C}_b \left ( - [\![ \partial _{XX}G ]\!] {\cdot }\partial _{XT}G + [\![ \partial _{X}^3G ]\!] {\cdot }\partial _TG \right)_{X=A(T)} \\ &{}+ \delta ^2 \frac {\tilde {C}_b}{2} [\![ \left |\partial _{XX}G\right |^2 ]\!]\, \dot {A}(T) \\ &{}+ \delta ^2 \tilde {C}_b \left \{\int _0^{A(T)} \partial _X^4G {\cdot }\partial _TG \ \mathrm {d} X + \int _{A(T)}^{X_r} \partial _X^4G {\cdot }\partial _TG \ \mathrm {d} X\right \}\!, \end{aligned} \end{equation}

where we have used the symmetry condition (2.5) at $X=0$ , and the far-field condition (2.6) at $X=X_r$ . The symbol $[\![ {\cdot }]\!]$ denotes the jump of a quantity across the contact line from $X=A(T)^-$ to $X=A(T)^+$ .

The dissipation rate of the interfacial energy can be decomposed into three contributions:

(2.10) \begin{align} \frac {\mathrm {d} \mathcal{E}_s}{\mathrm {d} T} ={}& \dot {\mathcal{E}}_{s,1} + \dot {\mathcal{E}}_{s,2} + \dot {\mathcal{E}}_{s,3}, \notag \\ \dot {\mathcal{E}}_{s,1} ={}& \delta ^2 \left \{ \frac {1}{{\textit{Ca}}}\, \partial _XH {\cdot }\partial _TH + (\tilde {\gamma }_1 - \tilde {\gamma }_2)\, \partial _XG {\cdot }\partial _TG \right \}_{X=A(T)}, \notag \\ \dot {\mathcal{E}}_{s,2} ={}& \dot {A}(T) \left \{ \frac {1}{{\textit{Ca}}} \left (1 + \frac {\delta ^2}{2}\left \vert \partial _XH\right \vert ^2\right) + \left ( \tilde {\gamma }_1 - \tilde {\gamma }_2 \right) \left ( 1 + \frac {\delta ^2}{2}\left \vert \partial _XG\right \vert ^2 \right) \right \}_{X=A(T)},\notag \\ \dot {\mathcal{E}}_{s,3} ={}& -\delta ^2 \left \{ \frac {1}{{\textit{Ca}}}\int _0^{A(T)} \partial _{XX}H {\cdot }\partial _TH \ \mathrm {d} X \right . \notag \\ &\left . {}+ \tilde {\gamma }_1 \int _0^{A(T)} \partial _{XX}G {\cdot }\partial _TG \ \mathrm {d} X + \tilde {\gamma }_2 \int _{A(T)}^{X_r} \partial _{XX}G {\cdot }\partial _TG \ \mathrm {d} X \right \}\!. \end{align}

We take the time derivative of the kinematic compatibility condition (2.7) and use the resulting relation to simplify $\dot {\mathcal{E}}_{s,1}$ and $\dot {\mathcal{E}}_{s,2}$ . After some manipulation, we obtain

(2.11) \begin{equation} \begin{aligned} \dot {\mathcal{E}}_{s,1} + \dot {\mathcal{E}}_{s,2} ={} & \delta ^2 \left \{ \left ( \frac {1}{{\textit{Ca}}}\,\partial _XH + (\tilde {\gamma }_1 - \tilde {\gamma }_2)\, \partial _XG \right)\partial _TG \right \}_{X=A(T)} \\ &{} + \frac {\delta ^2}{{\textit{Ca}}}\, \dot {A}(T) \left \{ \frac {1}{2}\varTheta _Y^2 - \frac {1}{2}\left \vert \partial _XH - \partial _XG\right \vert ^2 \right \}_{X=A(T)} + \operatorname {O}\big(\delta ^4\big), \end{aligned} \end{equation}

where $\varTheta _Y$ is the rescaled Young’s angle defined via Young’s relation (Young Reference Young1805)

(2.12) \begin{equation} {\textit{Ca}}\left (\tilde {\gamma }_2 - \tilde {\gamma }_1\right) = 1 - \frac {\delta ^2}{2} \varTheta _Y^2, \end{equation}

the right-hand side being the small-angle expansion of the cosine function. To handle $\dot {\mathcal{E}}_{s,3}$ , we use the thin film (2.4), which leads to the identity

(2.13) \begin{equation} \begin{aligned} &\int _0^{A(T)} \partial _{XX}H {\cdot }\partial _TH \ \mathrm {d} X \\ &= \int _0^{A(T)} \partial _{XX}H {\cdot }\partial _TG \ \mathrm {d} X + \int _0^{A(T)} \frac {1}{{\textit{Ca}}} \left (H-G\right)^2\left (H-G+\lambda \right)\big \vert \partial _{X}^3 H\big \vert ^2 \mathrm {d} X, \end{aligned} \end{equation}

where we have used integration by parts and the fact that $F=0$ at $X=0$ and $X = A(T)$ .

Finally, we combine (2.9), (2.10), (2.11) and (2.13), and apply the sheet (2.4c ) to obtain the dissipation rate for the total energy:

(2.14a) \begin{align} \frac {\mathrm {d} \mathcal{E}}{\mathrm {d} T} = & - \frac {\delta ^2}{{\textit{Ca}}^2} \int _0^{A(T)} \left (H-G\right)^2\left (H-G+\lambda \right)\big \vert \partial _{X}^3 H\big \vert ^2 \mathrm {d} X \end{align}

(2.14b) \begin{align} &{} - \delta ^2 \tilde {C}_b \bigl \{ [\![ \partial _{XX}G ]\!] {\cdot }\partial _{XT}G \bigr \}_{X=A(T)} + \delta ^2 \frac {\tilde {C}_b}{2}\, [\![ \left \vert \partial _{XX}G\right \vert ^2 ]\!]\, \dot {A}(T) \end{align}

(2.14c) \begin{align} & {}+ \delta ^2 \left \{ \left (\frac {1}{{\textit{Ca}}}\,\partial _XH + \left (\tilde {\gamma }_1 - \tilde {\gamma }_2\right)\partial _XG + \tilde {C}_b\, [\![ \partial _{X}^3 G ]\!] \right) \partial _TG \right \}_{X=A(T)} \end{align}

(2.14d) \begin{align} &{} + \frac {\delta ^2}{{\textit{Ca}}} \left \{\frac {1}{2}\varTheta _Y^2 - \frac {1}{2}\left \vert \partial _XH - \partial _XG\right \vert ^2\right \}_{X=A(T)} \dot {A}(T) \end{align}

(2.14e) \begin{align} & {}+ \operatorname {O}\big(\delta ^4\big). \end{align}

This equation shows that the energy dissipation consists of dominant terms of $\operatorname {O} (\delta ^2)$ and some higher-order contributions. Among the dominant terms, the term in (2.14a ) arises from the dynamics of the bulk fluid, while the remaining terms (2.14b )–(2.14d ) are associated with the contact line dynamics. Under the physical constraint $H-G \geqslant 0$ for $0 \lt X \lt A(T)$ , the term in (2.14a ) is always non-positive, and the thin-film motion always dissipates energy. To ensure that the total system energy decays, at least to leading order in $\operatorname {O} (\delta ^2)$ , we will derive and impose appropriate conditions at the moving contact line. For this purpose, we have written the contact-line-related terms in the form of a product of a generalised force and a generalised flux. Next, we examine these terms one by one.

In the first term of (2.14b ), the generalised force is the jump in the curvature across the contact line, while the generalised flux is the rate of sheet rotation about the contact line. We assume that the sheet rotation here does not induce any energy change, so we set the generalised force to zero, and obtain a continuity condition for the curvature:

(2.15) \begin{equation} [\![ \partial _{XX}G ]\!] = 0 \quad \text{at } X=A(T). \end{equation}

The second term of (2.14b ) vanishes given the continuity of the sheet curvature established above. Therefore, this term does not contribute to the energy dissipation.

In (2.14c ), the generalised force is the resultant vertical force acting at the contact line, while the generalised flux is the vertical velocity of the sheet at the contact line. (Note that this is not a material velocity, however.) We assume that the sheet motion at the contact line alone does not induce any energy change. This yields the following balance condition at the contact line:

(2.16) \begin{equation} \frac {1}{{\textit{Ca}}}\,\partial _XH + \left (\tilde {\gamma }_1 - \tilde {\gamma }_2\right) \partial _XG = -\tilde {C}_b\, [\![ \partial _{X}^3 G ]\!] \quad \text{at } X=A(T). \end{equation}

In (2.14d ), the generalised force is an approximation of the unbalanced Young stress in the horizontal direction, while the generalised flux is the contact line speed. Assuming a linear constitutive relation between the two, we obtain

(2.17) \begin{equation} \frac {1}{2} \left \vert \partial _XH - \partial _XG\right \vert ^2 - \frac {1}{2}\varTheta _Y^2 = {\textit{Ca}}\, \tilde {\mu }_{\varLambda }\, \dot {A}(T). \end{equation}

Here, $\tilde {\mu }_{\varLambda } \gt 0$ is a dimensionless friction coefficient, representing the resistance to contact line motion. It is a physical parameter whose value in general depends on fluid properties, substrate roughness and microscopic characteristics.

To conclude this section, we have derived a model for the motion of a thin film on an elastic sheet based on the lubrication approximation. The primary unknowns of this model are the liquid–air interface profile $z=h(x,t)$ and the elastic sheet profile $z=g(x,t)$ , with the contact line located at $x=a(t)$ . The dimensionless governing equations are summarised below. For better readability, we revert to lowercase variables and omit overhead tildes:

(2.18a) \begin{align} \partial _t\left (h-g\right) + {\textit{Ca}}^{-1}\, \partial _x \left ( \left (h-g\right)^2 \left (h-g+\lambda \right) \partial _{x}^3 h\right) = 0, &\quad 0 \lt x \lt a(t), \end{align}

(2.18b) \begin{align} C_b\, \partial _x^4 g - \gamma _1\, \partial _{xx}g = {\textit{Ca}}^{-1}\, \partial _{xx} h, &\quad 0 \lt x \lt a(t), \end{align}

(2.18c) \begin{align} C_b\, \partial _x^4 g - \gamma _2\, \partial _{xx}g = 0, &\quad a(t) \lt x. \end{align}

The boundary conditions are

(2.19a) \begin{align} \partial _x h = \partial _{x}^3 h = 0 \quad &\text{at } x=0, \end{align}

(2.19b) \begin{align} \partial _x g = \partial _{x}^3 g = 0 \quad &\text{at } x=0, \end{align}

(2.19c) \begin{align} g = \partial _{xx} g = 0 \quad &\text{at } x = +\infty . \end{align}

The contact line conditions at $x=a(t)$ are

(2.20a) \begin{gather} [\![ g ]\!] = [\![ \partial _{x} g ]\!] = [\![ \partial _{xx} g ]\!] = 0, \end{gather}

(2.20b) \begin{gather} {\textit{Ca}}^{-1}\,\partial _xh + \left (\gamma _1-\gamma _2\right)\partial _xg = -C_b\, \big[\!\big[ \partial _{x}^3 g \big]\!\big], \end{gather}

(2.20c) \begin{gather} \frac {1}{2} \left \vert \partial _xh - \partial _xg\right \vert ^2 - \frac {1}{2}\theta _Y^2 = {\textit{Ca}}\, \mu _{\varLambda \,} \dot {a}(t), \end{gather}

(2.20d) \begin{gather} h(a(t), t) = g(a(t), t). \end{gather}

In this model, the total energy of the system changes as

(2.21) \begin{align} \frac {\mathrm {d} \mathcal{E}}{\mathrm {d} t} = - \frac {\delta ^2}{{\textit{Ca}}^2} \int _0^{a(t)} \left (h-g\right)^2\left (h-g+\lambda \right)\big \vert \partial _{x}^3 h\big \vert ^2 \mathrm {d} x - \delta ^2 \mu _{\varLambda } \left \vert \dot {a}(t)\right \vert ^2 + \operatorname {O}\big(\delta ^4\big), \end{align}

where the first and second terms on the right-hand side represent the dissipation rates due to the thin-film motion and the contact line motion, respectively.

3.1. Outer region

We expand the contact line speed in powers of $\epsilon$ :

(3.4) \begin{align} \dot {a}(t) = \dot {a}_0(t) + \epsilon\, \dot {a}_1(t) + \cdots , \end{align}

where $\epsilon \dot {a}_1$ is $\operatorname {o} (1)$ for the expansion to be well ordered. Similarly, we expand the profiles

(3.5) \begin{align} \begin{aligned} &h(x,t) = h_0(x,t) + \epsilon\, h_1(x,t) + \cdots , \\ &g(x,t) = g_0(x,t) + \epsilon\, g_1(x,t) + \cdots . \end{aligned} \end{align}

Substituting these expansions into the governing equation (3.2), we obtain the leading-order equations for $h_0$ and $g_0$ :

(3.6a) \begin{align} \partial _x\left ( \left (h_0-g_0\right)^3 \partial _{x}^3 h_0 \right) = 0, \quad & 0 \lt x \lt a(t), \end{align}

(3.6b) \begin{align} -\gamma _1\, \partial _{xx}g_0 = \partial _{xx}h_0, \quad & 0 \lt x \lt a(t), \end{align}

(3.6c) \begin{align} \partial _{xx} g_0 = 0, \quad & a(t) \lt x. \end{align}

The leading-order boundary conditions are given by

(3.7) \begin{equation} \begin{aligned} &\partial _x h_0 = \partial _{x}^3 h_0 = 0 \quad \text{at } x=0, \\ &\partial _x g_0 = 0 \quad \text{at } x=0, \\ &g_0 = 0 \quad \text{at } x = +\infty , \end{aligned} \end{equation}

and the leading-order contact line conditions are

(3.8) \begin{equation} \begin{aligned} &{[\![} g_0 ]\!] = 0, \\ &h_0(a(t), t) = g_0(a(t), t). \end{aligned} \end{equation}

To solve this system, we first deduce from (3.6c ) and the far-field condition that $g \equiv 0$ when $a(t) \leqslant x$ . We then solve (3.6a ) and (3.6b ) with the boundary and contact line conditions. To fully determine the solution, we need to impose the conservation of the liquid volume $V$ . The leading-order solution is then given by

(3.9) \begin{equation} \begin{gathered} h_0(x,t) = \frac {p_0}{2} \big(a^2 - x^2\big), \quad 0 \leqslant x \leqslant a(t), \\ g_0(x,t) = \begin{cases} \dfrac {p_0}{2\gamma _1}\left (x^2 - a^2\right)\!, & 0 \leqslant x \leqslant a(t), \\ 0, & a(t) \lt x, \end{cases} \end{gathered} \end{equation}

where we have introduced the pressure $p_0$ for notational convenience:

(3.10) \begin{equation} p_0 = \frac {3V}{a^3\big(1 + \gamma _1^{-1}\big)}. \end{equation}

We note that this solution automatically satisfies the higher-order boundary conditions $\partial _x^3g_0 = 0$ at $x=0$ , and $\partial _x^2g_0 = 0$ at $x=+\infty$ . For later use, we compute the leading-order film thickness:

(3.11) \begin{equation} h_0 - g_0 = \frac {3V}{2a^3} \big(a^2 - x^2\big). \end{equation}

The next-order equations from (3.2) are

(3.12a) \begin{align} \partial _t\left (h_0 - g_0\right) + \partial _x \big((h_0-g_0)^3 \partial _{x}^3 h_1\big) = 0, \quad & 0 \lt x \lt a(t), \end{align}

(3.12b) \begin{align} -\gamma _1\, \partial _{xx}g_1 = \partial _{xx}h_1, \quad & 0 \lt x \lt a(t), \end{align}

(3.12c) \begin{align} \partial _{xx}g_1 = 0, \quad & a(t) \lt x. \end{align}

The boundary conditions and contact line conditions take the same form as in (3.7) and (3.8), with $h_0$ , $g_0$ replaced by $h_1$ , $g_1$ , respectively. As before, we have $g_1 \equiv 0$ when $a(t) \leqslant x$ . To solve for $h_1$ , we note that the dependence of the leading-order profiles $h_0$ and $g_0$ on time is only through $a(t)$ . We integrate (3.12a ) once with respect to $x$ to obtain

(3.13) \begin{align} \dot {a}_0 \frac {3V}{2a^4} x \big(x^2 - a^2\big) + \partial _{x}^3 h_1 \left [\frac {3V}{2a^3}\big(a^2-x^2\big)\right ]^3 = \text{const}. \end{align}

We set the integration constant to zero, as it will yield a solution that matches the solutions in the other regions and is consistent with the numerical results. We proceed by further integrations. Using the boundary and contact line conditions as well as volume conservation, we can determine the next-order solutions:

(3.14a) \begin{gather} h_1 = \dot {a}_0 \frac {a^4}{9V} \biggl \{ \left (a+x\right) \log \left (a+x\right) + \left (a-x\right) \log \left (a-x\right) - 2a \log {2a} + \frac {3}{2a} \big(a^2-x^2\big) \biggr \}, \end{gather}

(3.14b) \begin{gather} g_1 = \begin{cases}-\dfrac {1}{\gamma _1}h_1, & 0 \leqslant x \leqslant a(t), \\ 0, & a(t) \lt x. \end{cases} \end{gather}

Before moving on, we introduce a moving coordinate system centred at the contact line. Let $x = a(t) - y$ , where $y \leqslant a(t)$ . In this coordinate system, the contact line is located at $y=0$ , and the governing equation for $h = h(y,t)$ becomes

(3.15) \begin{equation} \epsilon \left [ \partial _t \left (h-g\right) + \dot {a}\, \partial _y \left (h-g\right) \right ] + \partial _y \left ( \left (h-g\right)^2\left (h-g+\lambda \right) \partial _{y}^3 h \right) = 0. \end{equation}

The governing equations for $g = g(y,t)$ remain the same as in (3.2b )–(3.2c ), but with $\partial _x$ replaced by $\partial _y$ . The symmetry conditions are now imposed at $y=a(t)$ , the far-field conditions at $y=-\infty$ , and the contact line conditions at $y=0$ . Their forms remain unchanged except for $\partial _x$ being replaced by $\partial _y$ . Near the contact line as $y\rightarrow 0^+$ , the outer solutions in terms of $y$ are

(3.16) \begin{equation} \begin{aligned} h_0 \sim & \alpha _{\textit{app}} y - \frac {\alpha _{\textit{app}}}{2a}y^2, \\ h_1 \sim & \frac {\dot {a}_0}{\theta _{\textit{app}}^2} \left \{y\log {y} + \left (2 - \log {2a} \right) y + \cdots \right \}\!, \end{aligned} \end{equation}

and

(3.17) \begin{equation} \begin{aligned} g_0 \sim & \beta _{\textit{app}} y - \frac {\beta _{\textit{app}}}{2a} y^2, \\ g_1 \sim & -\frac {1}{\gamma _1} h_1. \end{aligned} \end{equation}

Here, $\alpha _{\textit{app}}$ , $\beta _{\textit{app}}$ and $\theta _{\textit{app}}$ are the apparent contact angles for the interface, for the sheet, and between them, respectively. They are defined as

(3.18) \begin{equation} \begin{aligned} \alpha _{\textit{app}} =& \frac {\partial h_0(0,{\cdot })}{\partial y} = \frac {3V}{a^2 \big(1 + \gamma _1^{-1}\big)}, \\ \beta _{\textit{app}} =& \frac {\partial g_0(0,{\cdot })}{\partial y} = -\frac {1}{\gamma _1} \frac {3V}{a^2\big(1 + \gamma _1^{-1}\big)}, \\ \theta _{\textit{app}} =& \alpha _{\textit{app}} - \beta _{\textit{app}} = \frac {3V}{a^2}. \end{aligned} \end{equation}

Next, we analyse the problem in the bending region.

3.2. Bending region

We rescale the variables in the bending region as $y = \epsilon \tilde {y}$ , $h = \epsilon \tilde {h}$ and $g = \epsilon \tilde {g}$ . The governing equations, written in terms of these variables, are

(3.19a) \begin{gather} \epsilon ^2\,\partial _t\big(\tilde {h}-\tilde {g}\big) + \epsilon \dot {a}\,\partial _{\tilde {y}}\big(\tilde {h}-\tilde {g}\big) + \partial _{\tilde {y}}\big[\big(\tilde {h}-\tilde {g}\big)^2\big(\tilde {h}-\tilde {g}+\lambda \big) \partial _{\tilde {y}}^3\tilde {h} \big] = 0, \quad 0 \lt \tilde {y}, \end{gather}

(3.19b) \begin{gather} C_b\,\partial _{\tilde {y}}^4\tilde {g} - \gamma _1\, \partial _{\tilde {y}\tilde {y}}\tilde {g} = \partial _{\tilde {y}\tilde {y}} \tilde {h}, \quad 0 \lt \tilde {y}, \end{gather}

(3.19c) \begin{gather} C_b\,\partial _{\tilde {y}}^4\tilde {g} - \gamma _2\, \partial _{\tilde {y}\tilde {y}}\tilde {g} = 0, \quad \tilde {y} \lt 0. \end{gather}

To determine the appropriate boundary conditions as $\tilde {y} \rightarrow \pm \infty$ , we rewrite the outer solutions (3.16) and (3.17) in terms of $\tilde {y}$ , and arrange terms according to powers of $\epsilon$ . Accordingly, we impose the far-field behaviour

(3.20a) \begin{gather} \partial _{\tilde {y}}\tilde {h} \sim \alpha _{\textit{app}} - \epsilon \frac {\alpha _{\textit{app}}}{a} \tilde {y} + \epsilon \frac {\dot {a}_0}{\theta _{\textit{app}}^2} \left ( \log \tilde {y} + 3 + \log \epsilon - \log {2a} + \cdots \right)\!, \end{gather}

(3.20b) \begin{gather} \partial _{\tilde {y}} \tilde {g} \sim \beta _{\textit{app}} - \epsilon \frac {\beta _{\textit{app}}}{a} \tilde {y} - \epsilon \frac {\dot {a}_0}{\gamma _1\theta _{\textit{app}}^2} \left (\log \tilde {y} + 3 + \log \epsilon - \log {2a} + \cdots \right)\!, \end{gather}

as $\tilde {y} \rightarrow + \infty$ . In addition, from the outer solutions we have

(3.21) \begin{equation} \tilde {g} = \partial _{\tilde {y}} \tilde {g} = \partial _{\tilde {y}\tilde {y}} \tilde {g} = 0 \end{equation}

as $\tilde {y} \rightarrow -\infty$ . At the contact line, we impose

(3.22) \begin{equation} \begin{gathered} {[\![} \tilde {g} ]\!] = [\![ \partial _{\tilde {y}} \tilde {g} ]\!] = [\![ \partial _{\tilde {y}\tilde {y}} \tilde {g} ]\!] = 0, \\ \tilde {h}(0, {\cdot }) = \tilde {g}(0, {\cdot }). \end{gathered} \end{equation}

We omit the contact line conditions (3.3b ) and (3.3c ) at this stage, as they are only relevant in the immediate vicinity of the contact line, i.e. the inner region.

We expand $\tilde {h}$ and $\tilde {g}$ in powers of $\epsilon$ :

(3.23) \begin{align} \begin{aligned} \tilde {h}(\tilde {y},t) \sim & \tilde {h}_0(\tilde {y},t) + \epsilon \tilde {h}_1(\tilde {y},t) + {\cdots } , \\ \tilde {g}(\tilde {y},t) \sim & \tilde {g}_0(\tilde {y},t) + \epsilon \tilde {g}_1(\tilde {y},t) + {\cdots } . \end{aligned} \end{align}

The leading-order equations for $\tilde {h}_0$ and $\tilde {g}_0$ are

(3.24a) \begin{align} \partial _{\tilde {y}}\big [\big(\tilde {h}_0 - \tilde {g}_0\big)^3 \partial _{\tilde {y}}^3 \tilde {h}_0 \big] = 0, \quad & 0 \lt \tilde {y}, \end{align}

(3.24b) \begin{align} C_b\, \partial _{\tilde {y}}^4 \tilde {g}_0 - \gamma _1\, \partial _{\tilde {y}\tilde {y}} \tilde {g}_0 = \partial _{\tilde {y}\tilde {y}} \tilde {h}_0, \quad & 0 \lt \tilde {y}, \end{align}

(3.24c) \begin{align} C_b\, \partial _{\tilde {y}}^4 \tilde {g}_0 - \gamma _2\, \partial _{\tilde {y}\tilde {y}} \tilde {g}_0 = 0, \quad & \tilde {y} \lt 0. \end{align}

To solve these equations, we first consider (3.24a ) for $\tilde {h}_0$ . We integrate the equation once, set the constant of integration to zero, then integrate further and match the far-field condition in (3.20a ), to obtain

(3.25) \begin{equation} \partial _{\tilde {y}} \tilde {h}_0 \equiv \alpha _{\textit{app}}. \end{equation}

To solve (3.24b ) and (3.24c ), it is convenient to introduce

(3.26) \begin{align} \begin{gathered} l_i = \sqrt {C_b / \gamma _i}, \quad i=1, 2, \\ \tilde {y}_1 = \tilde {y} / l_1, \quad \tilde {y} \geqslant 0, \\ \tilde {y}_2 = \tilde {y} / l_2, \quad \tilde {y} \leqslant 0. \end{gathered} \end{align}

The general solutions to the fourth-order (3.24b ) and (3.24c ) are given by

(3.27) \begin{align} \tilde {g}_0 = C^+_i\, e^{\tilde {y}_i} + C^-_i\, e^{-\tilde {y}_i} + B_i \tilde {y} + A_i, \quad i=1, 2, \end{align}

where $i=1, 2$ corresponds to the solution when $\tilde {y}\gt 0$ and $\tilde {y}\lt 0$ , respectively. The eight constants are to be determined using the far-field conditions (3.20b ) and (3.21), and the contact line conditions in (3.22). From the far-field conditions, we find that $C_1^+ = 0$ , $C_2^- = 0$ , $A_2 = B_2 = 0$ and $B_1 = \beta _{\textit{app}}$ . We then apply the continuity conditions in (3.22) to obtain

(3.28) \begin{equation} \tilde {g}_0 = \begin{cases} C_1^-\, e^{-\tilde {y}_1} + \beta _{\textit{app}} \tilde {y} + A_1, & \tilde {y} \geqslant 0, \\ C_2^+\, e^{\tilde {y}_2}, & \tilde {y} \lt 0, \end{cases} \end{equation}

where

(3.29) \begin{equation} C_1^- = \frac {\sqrt {C_b\gamma _2}}{\sqrt {\gamma _1}\left (\sqrt {\gamma _1} + \sqrt {\gamma _2}\right)} \beta _{\textit{app}}, \quad C_2^+ = \frac {\sqrt {C_b\gamma _1}}{\sqrt {\gamma _2}\left (\sqrt {\gamma _1} + \sqrt {\gamma _2}\right)} \beta _{\textit{app}}, \quad A_1 = C_2^+ - C_1^-. \end{equation}

For future reference, we compute the leading-order slope of the sheet at the contact line:

(3.30) \begin{equation} \tilde {\beta } = \frac {\partial \tilde {g}_0(0,{\cdot })}{\partial \tilde {y}} = \frac {\sqrt {\gamma _1}}{\sqrt {\gamma _1} + \sqrt {\gamma _2}} \beta _{\textit{app}}. \end{equation}

This slope is independent of $C_b$ .

We then consider the next-order problem:

(3.31a) \begin{align} \dot {a}_0\, \partial _{\tilde {y}}\big(\tilde {h}_0 - \tilde {g}_0\big) + \partial _{\tilde {y}}\big[ \big(\tilde {h}_0 - \tilde {g}_0\big)^3 \partial _{\tilde {y}}^3 \tilde {h}_1 \big] = 0, & \quad 0 \lt \tilde {y}, \end{align}

(3.31b) \begin{align} C_b\, \partial _{\tilde {y}}^4 \tilde {g}_1 - \gamma _1\, \partial _{\tilde {y}\tilde {y}} \tilde {g}_1 = \partial _{\tilde {y}\tilde {y}}\tilde {h}_1, & \quad 0 \lt \tilde {y}, \end{align}

(3.31c) \begin{align} C_b\, \partial _{\tilde {y}}^4 \tilde {g}_1 - \gamma _2\, \partial _{\tilde {y}\tilde {y}} \tilde {g}_1 = 0, & \quad \tilde {y} \lt 0. \end{align}

We integrate (3.31a ) for $\tilde {h}_1$ once, and set the constant of integration to zero, to obtain

(3.32) \begin{equation} \partial _{\tilde {y}}^3 \tilde {h}_1 = -\frac {\dot {a}_0}{\big(\tilde {h}_0 - \tilde {g}_0\big)^2} = -\frac {\dot {a}_0}{\left (\theta _{\textit{app}} \tilde {y} + C_1^- - C_1^-\, e^{-\tilde {y}_1}\right)^2}. \end{equation}

This equation does not admit a closed-form solution; however, it is possible to find approximate solutions in the far field as $\tilde {y} \rightarrow +\infty$ , and the near field as $\tilde {y}\rightarrow 0$ . To find the far-field behaviour, we neglect the exponentially small term on the right-hand side, and consider the far-field approximation $\tilde {h}_1^\infty$ satisfying

(3.33) \begin{equation} \partial _{\tilde {y}}^3 \tilde {h}^{\infty }_1 = -\frac {\dot {a}_0}{\theta ^2_{\textit{app}}\left (\tilde {y} + \tilde {c}\right)^2}, \end{equation}

where the constant is $\tilde {c} = C_1^-/\theta _{\textit{app}} \lt 0$ . We integrate this equation and use the far-field condition (3.20a ) to determine the constant of integration. We obtain

(3.34) \begin{equation} \begin{aligned} \partial _{\tilde {y}} \tilde {h}_1^{\infty } = - \frac {\alpha _{\textit{app}}}{a} \tilde {y} + \frac {\dot {a}_0}{\theta _{\textit{app}}^2} \big \{\log (\tilde {y}+\tilde {c}) + E_1 \big \}, \end{aligned} \end{equation}

where $E_1 = 3 + \log (\epsilon /2a)$ .

Similarly, we consider the far-field approximation $\tilde {g}_1^{\infty }$ for $\tilde {g}_1$ , satisfying

(3.35) \begin{equation} C_b\, \partial _{\tilde {y}}^4 \tilde {g}_1^{\infty } - \gamma _1\, \partial _{\tilde {y}\tilde {y}} \tilde {g}_1^{\infty } = \partial _{\tilde {y}\tilde {y}} \tilde {h}_1^{\infty } = \frac {\dot {a}_0}{\theta _{\textit{app}}^2} \frac {1}{\tilde {y} + \tilde {c}} - \frac {\alpha _{\textit{app}}}{a}. \end{equation}

To solve this equation, we integrate it twice and use the far-field condition (3.20b ) to determine the constant of integration. We obtain

(3.36) \begin{align} \partial _{\tilde {y}\tilde {y}} \tilde {g}_1^{\infty } = D_1\, e^{-\tilde {y}_1} - \frac {\beta _{\textit{app}}}{a} + \frac {1}{\sqrt {\gamma _1 C_b}} \frac {\dot {a}_0}{\theta _{\textit{app}}^2}\, \varphi \left (\frac {\tilde {y}+\tilde {c}}{l_1}\right)\!, \end{align}

where $D_1$ is a constant to be determined, and $\varphi$ is a specific solution to a related second-order ordinary differential equation; see (B3) in Appendix B. We integrate once more to obtain

(3.37) \begin{equation} \partial _{\tilde {y}}\tilde {g}_1^{\infty } = -D_1l_1\, e^{-\tilde {y}_1} - \frac {\beta _{\textit{app}}}{a} \tilde {y} + \frac {\dot {a}_0}{\gamma _1\theta _{\textit{app}}^2} \left \{ \varphi ^{\prime}\left (\frac {\tilde {y}+\tilde {c}}{l_1}\right) - \log \left (\tilde {y}+\tilde {c}\right) - E_1 \right \}\!, \end{equation}

where the constant of integration is determined by matching to the far-field condition (3.20b ). We leave the constant $D_1$ undetermined, as it does not contribute to the far-field behaviour.

So far, we have obtained the far-field behaviour of $\tilde {h}_1$ and $\tilde {g}_1$ . We next turn to their near-field behaviour as $\tilde {y} \rightarrow 0^+$ . To solve (3.32) for $\tilde {h}_1$ in the near field, we first notice from the leading-order solutions (3.25)–(3.28) that $\tilde {h}_0 - \tilde {g}_0$ is analytic in a neighbourhood of $\tilde {y} = 0$ . Using a power series to represent the right-hand side of (3.32), we have

(3.38) \begin{align} \partial _{\tilde {y}}^3 \tilde {h}_1 = -\frac {\dot {a}_0}{\big(\alpha _{\textit{app}} - \tilde {\beta }\big)^2 \tilde {y}^2} + \operatorname {O}\left (\frac {1}{\tilde {y}}\right)\!, \quad \tilde {y} \rightarrow 0^+, \end{align}

where the remainder is a power series consisting of terms $\tilde {y}^k\ (k\geqslant -1)$ . We perform a term-by-term integration to obtain the near-field approximation

(3.39) \begin{equation} \tilde {h}_1^0 = \frac {\dot {a}_0}{\big(\alpha _{\textit{app}} - \tilde {\beta }\big)^2} \big( \tilde {y} \log \tilde {y} + \tilde {C}_1 \tilde {y}\big) + \operatorname {O}\big(\tilde {y}^2 \log \tilde {y}, \tilde {y}^2\big), \quad \tilde {y}\rightarrow 0^+, \end{equation}

where $\tilde {C}_1$ is a constant of integration, and the remainder series consists of the terms $\tilde {y}^k$ and $\tilde {y}^k \log \tilde {y} \ (k\geqslant 2)$ . We follow a similar procedure and get the equation for $\tilde {g}_1$ near $\tilde {y} = 0$ :

(3.40) \begin{align} C_b\, \partial _{\tilde {y}}^4 \tilde {g}_1 - \gamma _1\, \partial _{\tilde {y}\tilde {y}} \tilde {g}_1 = \frac {\dot {a}_0}{\big(\alpha _{\textit{app}} - \tilde {\beta }\big)^2} \frac {1}{\tilde {y}} + \operatorname {O}\big(\!\log \tilde {y}\big), \quad \tilde {y} \rightarrow 0^+ . \end{align}

We perform a term-by-term integration to obtain the near-field approximation

(3.41) \begin{equation} \tilde {g}_1^0 = \frac {\dot {a}_0}{C_b \big(\alpha _{\textit{app}} - \tilde {\beta }\big)^2} \tilde {y}^3 \log \tilde {y} + \operatorname {O}\big(\tilde {y}^3, \tilde {y}^4 \log \tilde {y}\big), \quad \tilde {y}\rightarrow 0^+, \end{equation}

where the remainder series consists of terms $\tilde {y}^k\ (k\geqslant 3)$ and $\tilde {y}^k\log \tilde {y}\ (k\geqslant 4)$ . In deriving this solution, we have imposed the continuity conditions at the contact line (see the next paragraph) and required matching with the intermediate solution in (3.59) (see § 3.4). Comparing (3.39) and (3.41), we see that $\tilde {h}_1^0$ exhibits a stronger singularity than $\tilde {g}_1^0$ .

For $\tilde {y} \lt 0$ , the next-order sheet solution is simply $\tilde {g}_1 \equiv 0$ , by noting the far-field condition (3.21) and the continuity conditions (3.22) up to the second derivative at the contact line. To summarise, we have obtained the interface and the sheet solutions up to $\operatorname {O} (\epsilon)$ in the bending region. The leading-order solutions are given in (3.25) and (3.28). The far-field behaviour of the next-order solutions is given in (3.34) and (3.37), while the near-field behaviour is given in (3.39) and (3.41).

3.3. Inner region

We now consider the inner region by introducing the rescaled variables $\bar {y} = y/\lambda$ , $\bar {h} =h/\lambda$ and $\bar {g} = g/\lambda$ . We rewrite the governing equation (3.2) in terms of these inner variables:

(3.42a) \begin{align} \lambda \epsilon\, \partial _t \left (\bar {h} - \bar {g}\right) + \epsilon \dot {a}\, \partial _{\bar {y}} \left (\bar {h} - \bar {g}\right) + \partial _{\bar {y}} \left [ \left (\bar {h}-\bar {g}\right)^2 \left (\bar {h}-\bar {g}+1\right) \partial _{\bar {y}}^3 \bar {h} \right ] = 0, \quad & 0 \lt \bar {y}, \end{align}

(3.42b) \begin{align} C_b \epsilon ^2\, \partial _{\bar {y}}^4 \bar {g} - \gamma _1 \lambda ^2\, \partial _{\bar {y}\bar {y}} \bar {g} = \lambda ^2\, \partial _{\bar {y}\bar {y}} \bar {h}, \quad & 0 \lt \bar {y}, \end{align}

(3.42c) \begin{align} C_b \epsilon ^2\, \partial _{\bar {y}}^4 \bar {g} - \gamma _2 \lambda ^2\, \partial _{\bar {y}\bar {y}} \bar {g} = 0, \quad & \bar {y} \lt 0. \end{align}

At the contact line $\bar {y} = 0$ , both the compatibility condition (3.3d ) and the continuity conditions in (3.3a ) apply. In addition, we impose the contact line conditions (3.3b ) and (3.3c ). They now become

(3.43a) \begin{gather} \partial _{\bar {y}} \bar {h} + \left (\gamma _1 - \gamma _2\right) \partial _{\bar {y}} \bar {g} = -C_b \lambda ^{-2} \epsilon ^2\, \big[\!\big[ \partial _{\bar {y}}^3 \bar {g} \big]\!\big], \end{gather}

(3.43b) \begin{gather} \frac {1}{2} \left \vert \partial _{\bar {y}} \bar {h} - \partial _{\bar {y}} \bar {g} \right \vert ^2 - \frac {1}{2}\theta _Y^2 = \epsilon \mu _{\varLambda } \dot {a}. \end{gather}

We note that since $\lambda$ is exponentially small compared to $\epsilon$ , the equations for $\bar {g}$ reduce to $\partial _{\bar {y}}^4 \bar {g} = 0$ , and the contact line condition (3.43a ) reduces to $[\![ \partial _{\bar {y}}^3 \bar {g} ]\!]=0$ . For the purpose of matching, we require the sheet curvature to vanish at infinity, and we obtain the solution

(3.44) \begin{equation} \partial _{\bar {y}} \bar {g} \equiv \beta _{{m}} = \text{const}, \end{equation}

which is accurate up to all orders of $\epsilon$ .

We only need to expand the interface profile

(3.45) \begin{align} \bar {h}(\bar {y},t) = \bar {h}_0(\bar {y},t) + \epsilon\, \bar {h}_1(\bar {y},t) + {\cdots } . \end{align}

From (3.42a ), we see that the leading-order equation for $\bar {h}$ is

(3.46) \begin{align} \partial _{\bar {y}} \left [ \left (\bar {h}_0 - \bar {g}_0\right)^2\left (\bar {h}_0 - \bar {g}_0 + 1\right) \partial _{\bar {y}}^3\bar {h}_0 \right ] = 0 \quad \text{for } 0 \lt \bar {y}, \end{align}

and from (3.43b ), the leading-order contact line condition is

(3.47) \begin{align} \partial _{\bar {y}} \bar {h}_0 - \partial _{\bar {y}} \bar {g} \equiv \theta _Y \quad \text{at } \bar {y} = 0. \end{align}

We further require the interface curvature to vanish at infinity, and hence obtain

(3.48) \begin{equation} \partial _{\bar {y}} \bar {h}_0 \equiv \alpha _{{m}} = \theta _Y + \beta _{{m}}, \quad 0 \lt \bar {y}. \end{equation}

For the next-order problem, the governing equation is

(3.49) \begin{equation} \dot {a}_0\, \partial _{\bar {y}} \left (\bar {h}_0 - \bar {g}_0\right) + \partial _{\bar {y}} \left [ \left (\bar {h}_0 - \bar {g}_0\right)^2 \left (\bar {h}_0 - \bar {g}_0 + 1\right) \partial _{\bar {y}}^3\bar {h}_1 \right ] = 0, \quad 0 \lt \bar {y}, \end{equation}

and the contact line condition is

(3.50) \begin{align} \partial _{\bar {y}} \bar {h}_1 = \frac {\mu _{\varLambda } \dot {a}_0}{\theta _Y} \quad \text{at } \bar {y} = 0. \end{align}

We integrate (3.49) and set the constant of integration to zero to obtain

(3.51) \begin{align} \partial _{\bar {y}}^3 \bar {h}_1 = -\frac {\dot {a}_0}{\theta _Y \bar {y} \left (\theta _Y \bar {y} + 1\right)}. \end{align}

We integrate further, require the curvature to vanish at infinity, and impose the contact line condition to obtain

(3.52) \begin{align} \partial _{\bar {y}} \bar {h}_1 = \frac {\dot {a}_0}{\theta _Y} \left \{ \frac {1}{\theta _Y}\log \left (\theta _Y \bar {y} + 1\right) + \bar {y} \log \left (1 + \frac {1}{\theta _Y \bar {y}}\right) + \mu _{\varLambda } \right \}\!. \end{align}

From here, we integrate one more time to obtain

(3.53) \begin{equation} \begin{aligned} \bar {h}_1 ={} & \frac {\dot {a}_0}{\theta _Y} \biggl \{ \left (\mu _{\varLambda } - \frac {1}{2\theta _Y}\right) \bar {y} + \frac {1}{2\theta _Y^2} \log \left (\theta _Y \bar {y} + 1\right) + \frac {1}{\theta _Y} \bar {y} \log \left (\theta _Y \bar {y} + 1\right) \\ & {}+ \frac {1}{2} \bar {y}^2 \left [\log \left (\theta _Y \bar {y} + 1\right) - \log \left (\theta _Y \bar {y}\right)\right ] \biggr \} \\ \sim{} & \frac {\dot {a}_0}{\theta _Y} \left \{ \frac {1}{\theta _Y}\bar {y}\log \bar {y} + \left (\mu _{\varLambda } + \frac {\log \theta _Y}{\theta _Y}\right) \bar {y} + {\cdots } \right \} \quad \text{as } \bar {y} \rightarrow +\infty . \end{aligned} \end{equation}

We note that the inner solutions derived here are valid in an $\operatorname {O} (\lambda)$ neighbourhood of the contact line, or when $\bar {y} = \operatorname {O} (1)$ .

3.4. Intermediate region

We proceed to the intermediate region that bridges the inner and bending regions. We apply the transformation

(3.54) \begin{equation} y(s) = \exp \left (\frac {s-1}{\epsilon }\right)\!, \quad 0 \lt s \lt 1 + \epsilon \log \epsilon , \end{equation}

such that $y(0) = \lambda$ corresponds to the inner region, and $y(1 + \epsilon \log \epsilon) = \epsilon$ corresponds to the bending region. Following the usual practice, we also set

(3.55) \begin{equation} h(y(s),t) = H(s,t)\, y, \quad g(y(s),t) = G(s,t)\, y. \end{equation}

These intermediate variables are related to the inner variables via $\bar y = e^{s/\epsilon }$ , $\bar h = H(s, t)\,\bar y$ and $\bar g = G(s, t)\, \bar y$ .

To derive the equation for $H$ , we drop the exponentially small term in (3.42a ), integrate once, and set the constant of integration to zero. We then rewrite the equation in terms of the intermediate variables to obtain

(3.56) \begin{equation} \epsilon \dot {a} + \left (H-G\right)\left (H-G + e^{-s/\epsilon }\right) \big(-\epsilon\, \partial _sH + \epsilon ^3\,\partial _{s}^3H \big) = 0. \end{equation}

To derive the equation for $G$ , we rewrite (3.42b ) in terms of the intermediate variables:

(3.57) \begin{align}& C_b\epsilon ^2\,e^{2(1-s)/\epsilon } \big( 2\epsilon\, \partial _s G - \epsilon ^2\, \partial _{ss} G - 2\epsilon ^3\, \partial _{s}^3G + \epsilon ^4\, \partial _s^4G \big) - \gamma _1 \big(\epsilon\,\partial _sG + \epsilon ^2\, \partial _{ss}G\big) \notag\\&\quad = \epsilon\, \partial _sH + \epsilon ^2\, \partial _{ss} H. \end{align}

Here we note that when $0 \lt s \lt 1 + \epsilon \log \epsilon$ , the term $\epsilon ^2\, e^{2(1-s)/\epsilon }$ is positive and exponentially large. Therefore, the above equation reduces to

(3.58) \begin{equation} 2\, \partial _s G - \epsilon\, \partial _{ss}G - 2 \epsilon ^2\, \partial _{s}^3G + \epsilon ^3\, \partial _{s}^4 G = 0, \end{equation}

with some exponentially small corrections. This equation is linear, independent of $H$ , and can be solved exactly. We obtain $G \equiv \text{const}$ by discarding solutions that lead to an exponentially large curvature. By matching to the inner solution (3.44) and the bending solutions (3.28)–(3.30), we find

(3.59) \begin{equation} G \equiv \beta _{{m}} = \tilde {\beta }, \end{equation}

which is accurate to all orders of $\epsilon$ .

We then turn to solve (3.56). We drop the exponentially small term $e^{-s/\epsilon }$ . Since $\partial _sG \equiv 0$ , we can rewrite the equation in terms of $H-G$ :

(3.60) \begin{align} \dot {a}_0 + \epsilon \dot {a}_1 - \left (H-G\right)^2 \partial _{s}\left (H-G\right) = \operatorname {O}\big(\epsilon ^2, \epsilon ^2\log \epsilon \big). \end{align}

From here we solve for $H-G$ :

(3.61) \begin{align} \left (H-G\right)^3 = \left (3\dot {a}_0 s + c_0\right) + \epsilon \left (3\dot {a}_1 s + c_1\right) + \operatorname {O}\big(\epsilon ^2, \epsilon ^2\log \epsilon \big), \end{align}

where $c_0$ and $c_1$ are the integration constants to be determined. Further manipulation gives the two-term expansion

(3.62) \begin{equation} H-G = \left (3\dot {a}_0 s + c_0\right)^{1/3} + \epsilon \frac {3\dot {a}_1 s + c_1}{3\left (3\dot {a}_0 s + c_0\right)^{2/3}} + \operatorname {O}\big(\big(\epsilon \log \epsilon \big)^2\big). \end{equation}

3.5. Matching

The outer and bending solutions have been matched to each other. We now consider the matching between the intermediate and inner solutions. To this end, we rewrite the intermediate solutions (3.62) in terms of the inner variables and expand them in powers of $\epsilon$ :

(3.63) \begin{align} \begin{aligned} \bar {h} - \bar {g} &= \left (H - G\right) \bar {y} \\ &\sim \left (3\epsilon \dot {a}_0\log \bar {y} + c_0\right)^{1/3} \bar {y} + \epsilon \frac {3\epsilon \dot {a}_1 \log \bar {y} + c_1}{3 \left (3\epsilon \dot {a}_0 \log \bar {y} + c_0\right)^{2/3}} \bar {y} + {\cdots } \\ &\sim c_0^{1/3} \bar {y} + \frac {\epsilon }{c_0^{2/3}} \left (\dot {a}_0 \bar {y} \log \bar {y} + \frac {c_1}{3} \bar {y} \right) + {\cdots } . \end{aligned} \end{align}

This expansion is then matched to the inner solution

(3.64) \begin{align} \bar {h}-\bar {g} \sim \theta _Y \bar {y} + \frac {\epsilon \dot {a}_0}{\theta _Y} \left \{ \frac {1}{\theta _Y} \bar {y} \log \bar {y} + \left (\mu _{\varLambda } + \frac {\log \theta _Y}{\theta _Y}\right)\bar {y} \right \} + {\cdots } . \end{align}

Matching the two expansions term by term, we obtain

(3.65) \begin{equation} c_0 = \theta _Y ^ 3, \quad c_1 = 3 \dot {a}_0 \left ( \mu _{\varLambda } \theta _Y + \log \theta _Y \right)\!. \end{equation}

To match the intermediate solutions to the bending ones, we rewrite the intermediate solutions (3.62) in terms of the bending variables and re-expand the expression to get

(3.66) \begin{align} \begin{aligned} \tilde {h} - \tilde {g} &= \left ( H - G \right) \tilde {y} \\ &\sim \left (3\dot {a}_0\left (1 + \epsilon \log \epsilon \tilde {y}\right) + c_0\right)^{1/3} \tilde {y} + \epsilon \frac {3\dot {a}_1\left (1 + \epsilon \log \epsilon \tilde {y}\right) + c_1}{3 \left (3\dot {a}_0\left (1 + \epsilon \log \epsilon \tilde {y}\right) + c_0\right)^{2/3}} \tilde {y} + {\cdots } \\ &\sim \left (c_0 + 3\dot {a}_0\right)^{1/3} \tilde {y} + \frac {\epsilon }{\left (c_0 + 3\dot {a}_0\right)^{2/3}} \left \{ \dot {a}_0 \tilde {y} \log \tilde {y} + \left (\dot {a}_0 \log \epsilon + \dot {a}_1 + \frac {c_1}{3}\right) \tilde {y} \right \} + {\cdots } . \end{aligned} \end{align}

This is matched to the near-field bending solution (cf. (3.25)–(3.28), (3.39)–(3.41)):

(3.67) \begin{align} \tilde {h} - \tilde {g} \sim (\alpha _{\textit{app}} - \tilde {\beta }) \tilde {y} + \frac {\epsilon \dot {a}_0}{ (\alpha _{\textit{app}} - \tilde {\beta })^2} ( \tilde {y} \log \tilde {y} + \tilde {C}_1 \tilde {y} ) + {\cdots } . \end{align}

By matching the corresponding terms, we obtain

(3.68) \begin{equation} c_0 + 3 \dot {a}_0 = (\alpha _{\textit{app}} - \tilde {\beta })^3, \quad \frac {\dot {a}_0\log \epsilon + \dot {a}_1 + c_1/3}{\left (c_0 + 3\dot {a}_0\right)^{2/3}} = \frac {\dot {a}_0}{ (\alpha _{\textit{app}} - \tilde {\beta })^2}\tilde {C}_1. \end{equation}

Combining (3.65) and (3.68), we deduce

(3.69a) \begin{gather} 3 \dot {a}_0 = \left (\alpha _{\textit{app}} - \frac {\sqrt {\gamma _1}}{\sqrt {\gamma _1} + \sqrt {\gamma _2}} \beta _{\textit{app}}\right)^3 - \theta _Y^3, \end{gather}

(3.69b) \begin{gather} \dot {a}_1 = \big(\tilde {C}_1 - \mu _{\varLambda } \theta _Y - \log \theta _Y - \log \epsilon \big) \dot {a}_0, \end{gather}

where we have used (3.30) for the angle $\tilde {\beta }$ . Equation (3.69a ) provides an estimate for the leading-order contact line speed in terms of the apparent contact angles $\alpha _{\textit{app}}$ and $\beta _{\textit{app}}$ . We note that this equation was derived under the assumption ${\textit{Ca}} = \epsilon$ ; the general case, which can be obtained by a straightforward rescaling of the governing equations, is given in (1.3). Equation (3.69b ) gives an estimate for the next-order correction to the contact line speed; however, this correction is of limited practical use, as the constant $\tilde {C}_1$ (cf. (3.39)) remains undetermined. Nevertheless, $\tilde {C}_1$ can, in principle, be determined by solving the full next-order problem in the bending region either analytically or numerically, which we will not consider in this study.

Before concluding this section, we make three remarks concerning (3.69a ). First, the relation can be viewed as a generalisation of the classical Cox–Voinov relation. Taking the limit $\gamma _1 \rightarrow +\infty$ , we find that $\beta _{\textit{app}} \rightarrow 0$ (cf. (3.18)), and we formally recover the classical result

(3.70) \begin{equation} 3\dot {a}_0 = \theta _{\textit{app}}^3 - \theta _Y^3. \end{equation}

This limiting case corresponds to an infinitely pre-stretched sheet, which excludes any vertical displacement and behaves effectively as a rigid substrate. Second, as long as the friction coefficient $\mu _{\varLambda }$ is an $\operatorname {O} (1)$ parameter, it does not appear in the leading-order relation (3.69a ), but only affects the next-order relation (3.69b ) through the moving contact line condition applied in the inner region. Finally, (3.69a ) is independent of the bending modulus $C_b$ . This is because the slope of the sheet in the bending region near the contact line, $\tilde {\beta }$ (cf. (3.30)), does not depend on $C_b$ . This observation is further supported by the numerical results presented in the next section.

4.1. Interface and sheet profiles

We perform the simulations for the spreading and receding dynamics of a liquid film. The parameters in (2.18)–(2.20) are chosen as $\gamma _1 = 2 / \epsilon$ , $\gamma _2 = 2.95 / \epsilon$ and $C_b = \epsilon$ . The initial interface and sheet profiles are given by

(4.1) \begin{align} h(x,0) = c_V \left (1 + 0.2 \cos \frac {2\pi x}{a(0)}\right) \left \{1 - \left (\frac {x}{a(0)}\right)^2\right \}\!, \quad g(x,0) \equiv 0, \end{align}

where the constant $c_V$ is chosen so that the liquid volume is $V=1$ . We use the initial condition $a(0) = 1.0$ for a spreading film, and $a(0) = 2.0$ for a receding film. In figure 3, we present the simulation results of the interface and sheet profiles at several time instants.

Figure 3. Numerical solutions of the thin-film model (2.18)–(2.20), showing the interface and sheet profiles at different time instants during (a) spreading and (b) receding dynamics. Blue and orange solid lines denote the interface and sheet, respectively.

In the spreading case, the contact line advances to the right, the apex heights of both the interface and the sheet decrease, and the liquid film becomes flatter. By $t=32$ , the system has almost reached equilibrium, with the contact line speed dropping below $3 \times 10^{-3}$ . In the receding case, the contact line retreats to the left, and the liquid becomes more concentrated. By $t=24$ , the system has almost reached equilibrium, with the contact line speed falling below $5 \times 10^{-3}$ .

Next, we compare the computed interface and sheet profiles at $t=1$ and $t=4$ with the asymptotic solutions derived in § 3. The asymptotic solutions invlove the parameters $\gamma _i\ (i=1,2)$ , ${\textit{Ca}}$ , $C_b$ , $\theta _Y$ , $\mu _{\varLambda }$ , $\lambda$ , $\epsilon$ , and the liquid volume $V$ , together with the apparent contact angles $\alpha _{\textit{app}}$ and $\beta _{\textit{app}}$ , the contact line position $a(t)$ , and its velocity components $\dot {a}_0$ and $\dot {a}_1$ . All physical parameters ( $\gamma _i$ , ${\textit{Ca}}$ , $C_b$ , $\theta _Y$ , $\mu _{\varLambda }$ , $\lambda$ , $\epsilon$ , $V$ ) are prescribed in the simulations. The contact line position $a(t)$ at $ t=1, 4$ is obtained directly from the numerical solutions. Given $a(t)$ , the apparent contact angles $\alpha _{\textit{app}}$ and $\beta _{\textit{app}}$ are computed using (3.18). The leading-order prediction of the contact line speed $\dot {a}_0$ follows from (3.69a ). The next-order contact line speed $\dot {a}_1$ cannot be determined directly from (3.69b ) due to the undetermined constant $\tilde {C}_1$ , so we estimate it via $\dot {a}_1 = (\dot {a}_{{num}} - \dot {a}_0) / \epsilon$ , where $\dot {a}_{{num}}$ is the numerically obtained contact line speed. We note that no parameters or variables are adjusted by fitting in these comparisons.

We first focus on the spreading dynamics discussed above. In figure 4(a), the two-term outer solutions $h_0 + \epsilon h_1$ and $g_0 + \epsilon g_1$ from (3.9) and (3.14) are compared with the numerical solutions from the thin-film model at $t=1.0$ and $t=4.0$ . We observe good overall agreement between the asymptotic and numerical solutions, except that the outer solutions do not capture the smooth bending of the sheet near the contact line. We then examine the bending region. In figure 4(b), the leading-order bending solutions $\tilde {h}_0$ and $\tilde {g}_0$ , given by (3.25) and (3.28), respectively, are compared with the numerical results. These asymptotic solutions agree well with the numerical solutions near the contact line, even though only the leading-order terms are used.

Figure 4. Interface and sheet profiles for a spreading film at $t=1.0$ and $t=4.0$ , shown in (a) the outer region and (b) the bending region. Circles denote numerical solutions of the thin-film model (2.18)–(2.20), while solid lines show the asymptotic approximations: the two-term outer solutions from (3.9) and (3.14) in (a), and the leading-order bending solutions from (3.25) and (3.28) in (b).

We also compare the numerical and asymptotic solutions in all four regions; in particular, we check the agreement of their slopes. For the same spreading dynamics discussed above, we plot the interface and sheet slopes, $\partial _{y}h$ and $\partial _{y}g$ , against the intermediate variable $s = \epsilon \log y + 1$ in figures 5(a) and 5(b), respectively. In figure 5(a), we show the interface slopes in their respective regions of validity: (i) the near-field outer solution $\partial _y h_0 + \epsilon\, \partial _y h_1$ from (3.16); (ii) the far-field bending solution $\partial _{\tilde {y}}\tilde {h}_0 + \epsilon\, \partial _{\tilde {y}}\tilde {h}_1^{\infty }$ from (3.25) and (3.34); (iii) the near-field bending solution $\partial _{\tilde {y}}\tilde {h}_0 + \epsilon\, \partial _{\tilde {y}}\tilde {h}_1^0$ from (3.25) and (3.39); (iv) the intermediate solution $H + \epsilon\, H^{\prime}(s)$ from (3.59) and (3.62); and (v) the inner solution $\partial _{\bar {y}} \bar {h}_0 + \epsilon\, \partial _{\bar {y}} \bar {h}_1$ from (3.48) and (3.53). For the intermediate solution, the next-order contact line speed $\dot {a}_1$ cannot be directly obtained from (3.69b ), so we estimate it using the numerical contact line speed $\dot {a}_{{num}}$ via $\dot {a}_1 = (\dot {a}_{{num}} - \dot {a}_0) / \epsilon$ . The two-term bending solution $\tilde {h}_0 + \epsilon \tilde {h}_1^{\infty }$ is a rearrangement of the outer solution in the bending region (with some higher-order corrections). This explains the smooth connection observed between the outer and bending solutions. Overall, we observe good agreement between the numerical and asymptotic solutions at both $t=1.0$ and $t=4.0$ . In figure 5(b), we overlay the far-field solution $\partial _{\tilde {y}} \tilde {g}_0 + \epsilon\, \partial _{\tilde {y}} \tilde {g}_1^{\infty }$ (with $D_1=0$ ; cf. (3.37)) and the near-field solution $\partial _{\tilde {y}} \tilde {g}_0 + \epsilon\, \partial _{\tilde {y}} \tilde {g}_1^0$ from (3.28) and (3.41). These asymptotic solutions agree well with the numerical ones in their respective regions of validity. In particular, this validates our prediction of the sheet slope $\tilde {\beta }$ in the bending region, as given in (3.30). The discrepancy towards the inner region is possibly due to the fact that the bending length is not sufficiently small, which limits the accuracy of the asymptotic solution. We also note that the far-field and near-field solutions generally do not match each other near $\tilde {y} = \operatorname {O} (1)$ , since this lies outside their region of validity.

Figure 5. (a) Interface slope $\partial _y h$ and (b) sheet slope $\partial _y g$ in a spreading film at $t=1.0$ and $t=4.0$ . Circles denote numerical solutions of the thin-film model (2.18)–(2.20), while lines indicate asymptotic solutions (see text for details). In (a), both far-field and near-field bending solutions are shown as black dashed lines, while in (b), only the bending region asymptotic solutions are displayed.

A similar comparison for the receding film is shown in figures 6(a) and 6(b). The slope profiles are visually similar at all time instants, so we only present the case $t=4.0$ . Again, we observe good agreement between the numerical results and the asymptotic predictions. We conclude that the asymptotic solutions in all four regions accurately capture the interface and sheet behaviours.

Figure 6. (a) Interface slope $\partial _y h$ and (b) sheet slope $\partial _y g$ for a receding film at $t=4.0$ . Circles denote numerical solutions of the thin-film model (2.18)–(2.20), while lines indicate asymptotic solutions (see text for details). In (a), both far-field and near-field bending solutions are shown as black dashed lines, while in (b), only the bending region asymptotic solutions are displayed.

It is interesting to note that the interface does not exhibit a clear bending region. In the numerical solutions, although the bending effect dominates the change in the sheet slope over the interval $s \in (0.5, 0.8)$ , the interface slope remains relatively unchanged in the same interval. This indicates that the interface is largely insensitive to small-scale variations in the substrate. We expect similar behaviour in liquid films spreading over substrates with variable height.

4.2. Contact line motion

We now compare the contact line speed predicted by the asymptotic relation (3.69a ) with that obtained from numerical solutions of the thin-film model (2.18)–(2.20). The initial conditions and parameters are the same as in the previous subsection, except that we now vary the tension: $\gamma _1 = 1/\epsilon , 2/\epsilon , 4/\epsilon , 8/\epsilon$ and $\gamma _2 = \gamma _1 + 0.95/\epsilon$ . Results for the spreading film and the receding film are shown in figures 7(a) and 7(b), respectively. Overall, the theoretical predictions agree well with the numerical results, especially for more rigid sheets, i.e. when $\gamma _1$ is large. The classical Cox–Voinov relation (3.70) is also shown in the same figures. As $\gamma _1$ increases, the results approach the classical one.

Figure 7. Contact line speed $\dot {a}$ versus film radius $a(t)$ for (a) a spreading film and (b) a receding film. Markers denote numerical solutions of the thin-film model (2.18)–(2.20). Solid lines show predictions from the asymptotic relation (3.69a ), while dashed lines are predictions using the classical Cox–Voinov relation (3.70) for rigid substrates.

Another important observation from these figures is that the spreading process becomes slower, and the receding process becomes faster, as the sheet becomes softer. The reason behind this behaviour can be understood from the dependence of the contact angle in the bending region, $\alpha _{\textit{app}} - \tilde {\beta }$ , on the sheet tension. Using (3.18) and (3.30), we can rewrite this angle as

(4.2a) \begin{gather} \alpha _{\textit{app}} - \tilde {\beta } = \rho (\gamma _1, \gamma _2)\, \theta _{\textit{app}}, \end{gather}

(4.2b) \begin{gather} \rho (\gamma _1, \gamma _2) = \frac {1}{1 + \gamma _1^{-1}} \left (1 + \frac {1}{\sqrt {\gamma _1} \left (\sqrt {\gamma _1} + \sqrt {\gamma _2}\right)}\right)\!, \end{gather}

where $\theta _{\textit{app}}$ depends only on the volume $V$ and the film radius $a(t)$ .

Recall that $\sigma _{\textit{pre}}$ is the stretching tension applied to the sheet (cf. (2.1)). When the stretching tension is large, we have

(4.3) \begin{align} \rho = 1-\frac {1}{2\sigma _{\textit{pre}}}+\operatorname {O}\left (\frac {1}{\sigma _{\textit{pre}}^2}\right)\!. \end{align}

In the limit of a rigid substrate, i.e. as $\sigma _{\textit{pre}}\rightarrow +\infty$ , we have $\rho \rightarrow 1$ , and recover the classical Cox–Voinov relation. For large but finite $\sigma _{\textit{pre}}$ , however, $\rho \lt 1$ and it increases monotonically with $\sigma _{\textit{pre}}$ . This reduction in the effective contact angle on softer substrates explains the observed phenomena of retarded spreading and enhanced receding in figures figures 7(a) and 7(b). It is worth noting that the retarded spreading of drops has also been observed on viscoelastic solids (Carré et al. Reference Carré, Gastel and Shanahan1996; Tamim & Bostwick Reference Tamim and Bostwick2023), where it is mainly attributed to viscous dissipation within the substrate. This mechanism, however, is not relevant here, since purely elastic solids do not exhibit such dissipation.

Next, we examine the effect of the bending modulus $C_b$ on contact line motion. In figure 8, we show numerical results from the thin-film model (2.18)–(2.20) for the contact line speed of a spreading film with $\gamma _1 = 1/\epsilon , 2/\epsilon , 4/\epsilon , 8/\epsilon$ and $C_b = 0.01\epsilon , 0.1\epsilon , \epsilon$ . We observe that the speed is barely affected by variations in $C_b$ . Similar behaviour is observed for the receding case (not shown). These results are consistent with the predictions from (3.69a ), where the contact line speed is independent of $C_b$ to leading order. We expect the relation to remain valid over a wide range of bending moduli, provided that the bending length satisfies $\lambda \ll l_1 \lesssim \epsilon$ .

Figure 8. Numerical results from the thin-film model (2.18)–(2.20) for the contact line speed $\dot {a}$ versus the film radius $a(t)$ in a spreading film. Circle, triangle and square markers correspond to $C_b = 0.01\epsilon$ , $0.1\epsilon$ and $\epsilon$ , respectively. Colours indicate sheet tension: blue for $\gamma _1 =1/\epsilon$ ; orange for $\gamma _1 = 2/\epsilon$ ; green for $\gamma _1=4/\epsilon$ ; and red for $\gamma _1 = 8/\epsilon$ .

Figure 9. Film radius $a(t)$ for (a) a spreading film with $a(0) = 1$ , and (b) a receding film with $a(0) = 2$ , with bending modulus $C_b = \epsilon$ . Markers denote numerical solutions of the thin-film model, while solid lines show the leading-order predictions by solving the generalised Cox–Voinov relation (3.69a ).

Finally, we plot the film radius against time in figure 9, with $C_b = \epsilon$ , varying tensions, and other parameters as prescribed at the beginning of the subsection. In this figure, the markers represent the numerical simulations of the thin-film model, while the solid lines correspond to the predictions from the asymptotic relation (3.69a ). Here, to get the radius $a(t)$ , we solve (3.69a ) as an ordinary differential equation for $a$ , with its right-hand side rewritten in terms of $a$ using (3.18). Overall, the thin-film dynamics is well captured by the asymptotic relation, particularly for large tensions $\gamma _1$ . The discrepancies are caused by the facts that the relation (3.69a ) is only a leading-order approximation, and that the bending length $l_1 = \sqrt {C_b/\gamma _1}$ is not sufficiently small for low values of $\gamma _1$ . Additionally, the assumption $\dot {a} \sim \operatorname {O} (1)$ is slightly violated during the initial transient dynamics. The errors accumulated in the transient stage carry over to the later stage of the evolution.

For a liquid film spreading on a rigid substrate, the film radius grows algebraically (Tanner’s law) and then approaches its equilibrium value exponentially (Thampi et al. Reference Thampi, Pagonabarraga, Adhikari and Govindarajan2016). A similar quantitative description can be obtained for the present problem. Let $a_{\infty }$ denote the equilibrium film radius predicted by (3.69a ). By substituting (3.18) and (4.2) into (3.69a ), we obtain

(4.4) \begin{equation} \dot {a}_0 = 9 \left (\rho (\gamma _1,\gamma _2)\,V\right)^3 \big(a_0^{-6} - a_{\infty }^{-6}\big). \end{equation}

At early times when $a_0$ is far below $a_\infty$ , (4.4) can be approximated by

(4.5) \begin{align} \dot {a}_0 \sim 9 \rho ^3V^3 a_0^{-6}\quad \Longrightarrow\quad a_0 \sim \big(63\rho ^3V^3t + a(0)\big)^{1/7}. \end{align}

This describes the algebraic regime characterised by the $t^{1/7}$ scaling. The exponent $1/7$ is the same as that in the case of a rigid substrate. At later times, when $\Delta a = a_0 - a_{\infty }$ is small, (4.4) can be approximated by

(4.6) \begin{align} \frac {{\rm d}}{{\rm d}t}\Delta a \sim -54\rho ^3V^3 a_{\infty }^{-7}\, \Delta a\quad \Longrightarrow\quad a_0(t) - a_{\infty } \sim \left (a(0) - a_\infty \right) \exp \big(-54\rho ^3V^3a_{\infty }^{-7}t\big). \end{align}

This is the exponential regime. Both regimes are quantitatively similar to the spreading on a rigid substrate, and differ in some constants only.