2.1. The squeezing problem $(F>0)$

Consider what happens when a drop of viscous liquid of volume $V$ and dynamic viscosity $\mu$ is placed between two horizontal plates $z=0$ and $z=h(t)>0$, the upper plate being impelled towards the fixed lower plate by a constant force $F$. It is observed in this situation that, when viewed from above, the drop tends to become of circular form of expanding radius $r=a(t)$, where, with $h(0)=h_{0}, a(0)=a_{0}$, by conservation of volume,

(2.1)\begin{equation} {\rm \pi}\, a^{2}(t)\,h(t) ={\rm \pi}\, a_{0}^{2}\,h_{0}=V= \text{const.} \end{equation}

We shall describe this evolving circular situation, shown in figure 1, as the basic state. In this state, provided $h(t) \ll a(t)$, standard lubrication theory implies that (i) the pressure $p(r,t)$ is effectively independent of $z$ and (ii) the radial velocity component $u(r,z,t)$ satisfies $|\partial u/\partial z| \gg |\partial u/\partial r|$ and is driven by the pressure gradient according to the equation

(2.2)\begin{equation} \frac{ \partial p}{\partial r }= \mu \frac{\partial^{2}u}{\partial \,z\, ^{2}}, \end{equation}

with boundary conditions $u(r,0,t)=u$[$r, h(t), t$]$=0$. Hence

(2.3)\begin{equation} u(r,z,t)=\frac{1}{2\mu}\frac{ \partial p}{\partial r }\,z\,(z-h), \end{equation}

and, averaging over $z$,

(2.4)\begin{equation} \bar{u}(r,t) \equiv \frac{1}{h}\int_{0}^{h}u(r,z,t)\, \text{d}z={-}\frac{h^{2}}{12\mu}\, \frac{ \partial p}{\partial r }. \end{equation}

Mass conservation and incompressibility indicate that the radial flux $Q(r,t)=2{\rm \pi} r\, h\, \bar {u}(r,t)$ at each radius $r$ is equal to the flux $-{\rm \pi} r^2\, \text {d}h/\text {d}t$ driven downwards by the upper plate. It follows that

(2.5)\begin{equation} \frac{ \partial p}{\partial r}= \frac{6\,\mu\, r}{h^3}\, \frac{\text{d}h}{\text{d}t} \end{equation}

and we note that

(2.6)\begin{equation} \nabla^{2}p=\frac{1}{r}\frac{\partial}{\partial r} r\, \frac{ \partial p}{\partial r} = \frac{12\mu}{h^3}\, \frac{\text{d}h}{\text{d}t}. \end{equation}

This is the particular form of ‘Reynolds’ equation’ applicable to this problem (see for example § 3.6 of Moffatt Reference Moffatt1977).

Let $\gamma$ denote the surface tension at the liquid/air boundary. Assuming a circular meniscus with zero contact angle, the radius of curvature is $h(t)/2$, so that there is a jump of pressure $2\gamma /h(t)$ with corresponding suction across this boundary. (If the contact angle is greater than zero, then the radius of curvature is increased, and the surface tension effect correspondingly reduced. Note also that, for the moment, we neglect the relatively small curvature $a(t)^{-1}$ of the boundary of the drop in the $\{x,y\}$ plane; this is, however, retained in the stability calculation of § 5 below.) Under these conditions, the boundary condition on $r=a(t)$ is

(2.7)\begin{equation} p=p_{a}-2\gamma/h(t) ,\quad \text{at}\ r=a(t), \end{equation}

where $p_{a}$ is the atmospheric pressure. With this boundary condition the pressure in the drop is then given by

(2.8)\begin{equation} p(r,t)-p_{a}= \frac{3\mu}{h^{3}(t)}\frac{\text{d}h}{\text{d}t}\left[r^2 -a^{2}(t)\right] - \frac{2\gamma}{h(t)}= p_{0}(r,t),\quad \text{say}. \end{equation}

The total vertical force exerted by the fluid on the upper plate is then

(2.9)\begin{equation} F_{p} =2{\rm \pi}\int_{0}^{a} p_{0}(r,t)\,r \,\text{d}r ={-}\frac{3{\rm \pi}\,\mu}{2h^3}\frac{\text{d}h}{\text{d}t}\, a^{4}-\frac{2{\rm \pi}\gamma}{h(t)}a^2. \end{equation}

Using (2.1) and neglecting plate inertia, the force balance $F=F_{p}$ gives

(2.10)\begin{equation} F=\frac{3\mu V^2}{8{\rm \pi}}\frac{\text{d}}{\text{d}t}\left(\frac{1}{h^4}\right)-\frac{2\gamma V}{h^2}. \end{equation}

With dimensionless variables

(2.11)\begin{equation} X=h_{0}^{2}/h^{2}(t) =a^{4}(t)/a_{0}^{4}\quad \text{and}\quad \tau=t/t_{0}, \quad \text{where} \ t_0=\frac{3}{8{\rm \pi}} \frac{\mu}{F}\,\frac{V^2}{h_0^4} = \frac{3{\rm \pi}^3}{8}\frac{\mu}{F}\frac{a_{0}^{8}}{V^{2}}, \end{equation}

(2.10) takes the convenient dimensionless form

(2.12)\begin{equation} \text{d}\kern0.06em X^{2}/\text{d}\tau \equiv 2X\,\text{d}\kern0.06em X/\text{d}\tau =1+\lambda\, X\quad\text{with}\ X(0)=1, \end{equation}

where

(2.13)\begin{equation} \lambda =(2\gamma/F) \left(V/h_{0}^2\right)=(2\gamma/F) \left({\rm \pi}^{2} a_{0}^{4}/V\right). \end{equation}

Here, $\lambda$ is the relevant dimensionless measure of surface tension. It is small $(\sim 0.006)$ in the experiments described in § 7, but surface tension effects can evidently still become important as and when $X$ increases to O$(\lambda ^{-1})$, i.e. when $a/a_{0}$ increases to O$(\lambda ^{-1/4})$.

Equation (2.12) integrates to give

(2.14)\begin{equation} \tau=\frac{2}{\lambda^2} \left((X-1)\lambda -\log{\frac{1+\lambda X}{1+\lambda}}\right), \end{equation}

thus determining $X$ implicitly as a function of $\tau$ and $\lambda$ (cf. Ward Reference Ward2006). For $\lambda \ll 1$, this gives

(2.15)\begin{equation} \tau= X^2-1 - {\tfrac{2}{3}} \left( X^{3} -1\right) \lambda +\text{O}\left(\lambda^2\right)\,. \end{equation}

Figure 6(a) shows $X^{2}(\tau )$ as determined by (2.14) for three small values of $\lambda$; as expected from (2.15), the curves all asymptote to $X^{2}\sim 1+\tau$ as $\lambda \rightarrow 0$.

Figure 6. (a) Plot of $X^{2}(\tau ; \lambda )=a^{8}(t)/a_{0}^{8}$, with $\tau \equiv t/t_{0}$, as determined by (2.11) and (2.14), for $\lambda =0.001,0.01,0.02$; note that $X^{2}\sim 1+\tau$ as $\lambda \rightarrow 0$; (b) plot of $X^{2}(\tau ; |\lambda |)$, with $\tau =t/|t_{0}|$, for the case of contraction $F<0$ (as given by (2.19)); here $X^{2}\sim 1-\tau$ as $\lambda \rightarrow 0$.

If surface tension is neglected (i.e. $\lambda =0$), then (2.15) gives $X(\tau )=(1+\tau )^{1/2}$ or equivalently

(2.16a,b)\begin{equation} h(\tau)= h_{0}(1+\tau)^{{-}1/4}\,,\quad a(\tau)=a_{0}(1+\tau)^{1/8}. \end{equation}

The result for $a(\tau )$ is the ‘one-eighth power law’ (as indicated in the introduction). Figures 1 and 7 show the results of a simple experiment (details in § 7); the one-eighth power law shows up quite accurately in this experiment, for which $t_{0}\approx 2.15$ s and $\lambda \approx 3\times 10^{-3}$. The time taken to approach this one-eighth power law is evidently very short in this experiment. This ‘adjustment time’ has been recently investigated by Ball & Huppert (Reference Ball and Huppert2019) and Webber & Huppert (Reference Webber and Huppert2020) for the case of gravitational spread of a viscous fluid (Huppert Reference Huppert1982). They found that this time, somewhat independent of initial conditions, can be just fractions of a second in laboratory experiments, but many days in geological situations, where the extremely viscous lava domes from volcanic eruptions can eventually spread horizontally to hundreds of kilometres.

Figure 7. (a) Plot of $\hat {a}(\tau )\equiv a(t)/a_{0}$ for a squeezed drop; the red curve shows the function defined by (2.16a,b); the dotted lines correspond to the three images shown in figure 1; (b) plot of $\hat {a}^8(\tau )$, showing clearly the one-eighth power law over the full range of $\tau$.

2.2. The contraction problem $F<0$

If $F<0$ (i.e. if the upper plate is pulled away from the lower plate with a constant traction force $|F|$) leading to contraction of the area of the drop, then $t_{0}<0$, and it is natural therefore to redefine $\tau$ as $\tau =t/|t_{0}|$; with this modification, (2.12) becomes

(2.17)\begin{equation} \text{d}\kern0.06em X^{2}/\text{d}\tau \equiv 2X\,\text{d}\kern0.06em X/\text{d}\tau ={-}1+|\lambda|\, X\quad\text{with again}\ X(0)=1. \end{equation}

Now, a steady state with $X=1$ is possible if $|\lambda |=1$; but, if $|\lambda |\ne 1$,

(2.18)\begin{equation} \left|\text{d}\kern0.06em X^{2}/\text{d}\tau\right|_{\tau=0} ={-}1+|\lambda|, \end{equation}

so that $X^{2}$ increases or decreases from 1 according as $|\lambda | >$ or $<1$. If $|\lambda | >1$, the reduction in pressure in the drop due to surface tension is sufficient to overcome the traction force, thus decreasing $h(\tau )$ and increasing $a(\tau )$, ensuring good adhesion.

Equation (2.18) now integrates to give

(2.19)\begin{equation} \tau=\frac{2}{|\lambda|^2} \left((X-1)|\lambda|+\log{\frac{1-|\lambda| X}{1-|\lambda|}}\right), \end{equation}

again determining $X$ implicitly as a function of $\tau$. Figure 6(b) shows $X^{2}(\tau )$ as determined by (2.19) for several values of $|\lambda |$, both less than and greater than 1. For $|\lambda |<1$, $X^2$ goes to zero at a time $\tau ^{*}(|\lambda |)$, indicating a finite-time singularity [$h(\tau )=\infty$] at this time. If $|\lambda |=0$, the solution is

(2.20)\begin{equation} X=(1-\tau)^{1/2}\quad \text{so}\quad h(\tau)=h_{0}\left(1-\tau\right)^{{-}1/4}\,,\quad a(\tau)=a_{0} (1- \tau)^{1/8}, \end{equation}

indicating the finite-time singularity at $\tau =1$ anticipated in the introduction. These results hold only for so long as $h(\tau )\ll a(\tau )$, and cease to apply as $\tau =t/|t_{0}|$ approaches the singularity time, specifically when $1-\tau \sim (h_{0}/a_{0})^{8/3}$. At this stage, it is to be expected that the two plates will rapidly separate; however, as will be shown in § 4 below, the evolution for $\tau >0$ is subject to a fingering instability.

More generally, if $F=F(t)$ is time dependent, then (2.10) still holds, giving

(2.21)\begin{equation} \frac{1}{h^4}=\left(\frac{8{\rm \pi}\,}{3\mu\,V^2}\right)\int_{0}^{t} F(t')\,\text{d}t' \,+\frac{1}{h_{0}^4}. \end{equation}

For prescribed $F(t)$, e.g. $F(t)\sim t^{\,\alpha }$, it is then a trivial matter to obtain results analogous to (2.16a,b); however, we shall not pursue these possibilities here.

3.1. Neglect of compressibility of air in the bubble

Note that, from (3.2), when $\text {d}h/\text {d}t<0$, $p_{b}(t)>p_{a}$, i.e. the pressure in the bubble is greater than atmospheric (this is essentially because $\bar {u}_{r} >0$, so that, according to (2.4), $p(r,t)$ is monotonic decreasing in $r$); the air in the bubble is therefore slightly compressed. This effect is, however, small in the experiments described in § 7 below; for, using (3.5) and (3.6), from (3.2) we obtain

(3.8)\begin{equation} p_{b}(\tau)=p_{a} +\frac{2F\,h_{0}}{\left(V+V_{b}\right)(1+\tau)^{9/4}}. \end{equation}

With the estimates $F\sim 10\ \textrm {kg}\ \textrm {m}\,\textrm {s}^{-2}$, $V+V_{b}\sim 10^{-5}\ \textrm {m}^{3}$ and $h_{0}\sim 2\times 10^{-3}\ \textrm {m}$, we find that $2Fh_{0}/(V+V_{b}) \sim 4\times 10^{3}\ \textrm {kg}\ \textrm {m}^{-1}\ \textrm {s}^{-2}$, small compared with atmospheric pressure $p_{a}\sim 10^{5}\ \textrm {kg}\ \textrm {m}^{-1}\ \textrm {s}^{-2}$. Thus at time $\tau$, $[p_{b}(\tau )-p_{a}]/p_{a}\sim 0.04 (1+\tau )^{-9/4}$ and Boyle's law implies a correspondingly small reduction in the volume of the bubble due to compression. The compressibility of the air in the bubble is therefore negligible in the conditions of our experiment, although other circumstances can be imagined (much larger force $F$ and correspondingly much shorter time scale $t_{0}$ as given by (2.11)) when the effect might become significant.

4.1. Elliptic drop

The situation is well illustrated by the case of an elliptic drop, for which $C(t)$ has the form

(4.4)\begin{equation} \frac{ x^{2}}{a^{2}(t)} + \frac{ y^{2}}{b^{2}(t)} =1,\quad\text{with}\ a(0)=a_{0},\,b(0)=b_{0}, \end{equation}

with

(4.5)\begin{equation} {\rm \pi}\,a(t)b(t)h(t)= V=\text{const.} ={\rm \pi} \,a_{0}b_{0}h_{0}. \end{equation}

We may assume $b_{0}\le a_{0}$. The pressure is given by

(4.6)\begin{equation} p(x,y,t)-p_{a}=\hat{p}(t)\left(\frac{ x^{2}}{a^{2}(t)} + \frac{ y^{2}}{b^{2}(t)} -1\right), \end{equation}

so that, using (4.3),

(4.7)\begin{equation} \nabla^{2}p= 2 \hat{p}(t)\left(\frac{ 1}{a^{2}(t)} + \frac{ 1}{b^{2}(t)} \right)=\frac{12\,\mu}{h^{3}}\frac{\text{d}h}{\text{d}t}. \end{equation}

Hence

(4.8)\begin{equation} \hat{p}(t)=\frac{6\,\mu}{h^{3}}\frac{\text{d}h}{\text{d}t} \frac{a^2 b^2}{a^2 +b^2 } = \frac{6\,\mu\,V^2}{{\rm \pi}^2 h^{5}}\frac{\text{d}h}{\text{d}t}\frac{1}{ a^2 +b^2 }. \end{equation}

The downward force $F$ applied to the upper plate is now

(4.9)\begin{align} F=\iint_{S(t)}\hat{p}(t)\left(\frac{ x^{2}}{a^{2}(t)} + \frac{ y^{2}}{b^{2}(t)} -1\right)\,\text{d}\kern0.06em x\,\text{d}y={-}(1/2) \hat{p}(t)\,{\rm \pi} a b ={-}(V/2)\hat{p}(t)/h(t). \end{align}

The pressure contours $p=$ const. are shown in figure 8. The pressure gradient at the boundary, as shown by the arrows, is greater where these contours are close together; the rate of expansion in the $y-$direction is therefore greater than in the $x$-direction. We can calculate these rates as follows. The $x$-component of velocity at the point $(a,0)$ is

(4.10)\begin{equation} \frac{\text{d}a}{\text{d}t}= \bar{u}(x,0)={-}\frac{h^2}{12\mu}\left.\frac{\partial p}{\partial x}\right|_{x=a,\,y=0}={-}\frac{h^2}{6\mu a}\hat{p}(t) ={-}\frac{1}{ a} \frac{V^2}{{\rm \pi}^2 h^{3}}\frac{\text{d}h}{\text{d}t}\frac{1}{ (a^2 +b^2 )}, \end{equation}

and hence

(4.11)\begin{equation} \frac{\text{d}a^2}{\text{d}t}= \frac{V^2}{{\rm \pi}^2}\frac{\text{d}}{\text{d}t}\left(\frac{1}{h^2}\right)\frac{1}{ \left(a^2 +b^2 \right)}. \end{equation}

Similarly,

(4.12)\begin{equation} \frac{\text{d}b^2}{\text{d}t}= \frac{V^2}{{\rm \pi}^2}\frac{\text{d}}{\text{d}t}\left(\frac{1}{h^2}\right)\frac{1}{ \left(a^2 +b^2 \right)}, \end{equation}

and it follows that $a^2 (t)-b^2 (t)=$ constant $=a_{0}^2-b_{0}^2$. We then have from (4.11)

(4.13)\begin{equation} \frac{\text{d}}{\text{d}t}\left[a^4 -\left(a_{0}^2-b_{0}^2\right)a^2\right]= \frac{V^2}{{\rm \pi}^2}\frac{\text{d}}{\text{d}t}\left(\frac{1}{h^2}\right), \end{equation}

so that, with the initial condition $a(0)=a_{0}$, $a^2$ satisfies the quadratic equation

(4.14)\begin{equation} a^4 -(a_{0}^2-b_{0}^2)a^2=V^2/{\rm \pi}^2h^2, \end{equation}

with relevant root

(4.15)\begin{equation} a^2(t)= {\tfrac{1}{2}} \left(a_{0}^2-b_{0}^2\right) + {\tfrac{1}{2}} \left[ \left(a_{0}^2 - b_{0}^2\right)^{2}+4V^2/{\rm \pi}^2 h^2\right]^{1/2}. \end{equation}

The corresponding expression for $b^2(t)$ is

(4.16)\begin{equation} b^2 (t)=a^2 (t)-\left(a_{0}^2-b_{0}^2\right)={-}{\tfrac{1}{2}} \left(a_{0}^2-b_{0}^2\right) + {\tfrac{1}{2}} \left[ \left(a_{0}^2-b_{0}^2\right)^{2}+4V^2/{\rm \pi}^2h^2\right]^{1/2}, \end{equation}

so that

(4.17)\begin{equation} a^2(t) + b^2 (t)= \left[ \left(a_{0}^2-b_{0}^2\right)^{2}+4V^2/{\rm \pi}^2h^2\right]^{1/2}. \end{equation}

From (4.8) we then have

(4.18)\begin{equation} \hat{p}(t)={-} \frac{6\,\mu\,V^2}{{\rm \pi} h^{4}}\frac{\text{d}h}{\text{d}t} \left[ \left(a_{0}^2-b_{0}^2\right)^{2}{\rm \pi}^2h^2+4V^2\right]^{{-}1/2}. \end{equation}

and from (4.19),

(4.19)\begin{equation} F=\frac{3\,\mu\,V^2}{8{\rm \pi}}\frac{\text{d}}{\text{d}t} \left(\frac{1}{h^4}\right) \left[ 1+\left\{\left(a_{0}^2-b_{0}^2\right){\rm \pi} h/2V\right\}^2\right]^{{-}1/2}, \end{equation}

in agreement with (2.10) when $a_{0}=b_{0}$.

Figure 8. Plan view of the elliptic drop in the $\{x,y\}$ plane: pressure contours $p=$ const. for the particular initial condition $a_{0}/b_{0}=3$, as given by (4.6); the magnitude of the pressure gradient at the boundary is indicated by the length of the arrows.

In dimensionless form, with $\tau =t/t_{0}$ [where $t_{0}$ is still defined by (2.11)] and $\hat {h}(\tau )=h(t)/h_{0}$, (4.19) becomes

(4.20)\begin{equation} \frac{\text{d}}{\text{d}\tau} \left(\frac{1}{\hat{h}^4}\right)= \left[ 1+\kappa^{2} \hat{h}^2\right]^{1/2}, \end{equation}

where

(4.21)\begin{equation} \kappa=\left(a_{0}^2-b_{0}^2\right) /2 a_{0}b_{0}. \end{equation}

Here, $\kappa \ge 0$ is a measure of the initial eccentricity of the ellipse. With this notation, (4.15) and (4.16) become

(4.22a,b)\begin{align} a^2(\tau)/a_{0}b_{0}= \kappa+ \left[\kappa^{2}+ \hat{h}^{{-}2}(\tau)\right]^{1/2},\quad b^2(\tau)/a_{0}b_{0}={-}\kappa + \left[\kappa^{2}+ \hat{h}^{{-}2}(\tau)\right]^{1/2}. \end{align}

The evolution is shown in figure 9 for initial conditions $a_{0}=1,b_{0}=0.5$ (so $\kappa =0.75$). In (a), $F>0$, and the drop expands, becoming gradually more circular in form. In (b), $F<0$, and the drop contracts with increasing ellipticity [as $\hat {h}(\tau )$ increases], becoming singular as $\tau \rightarrow \tau _c = 0.71672\ldots ,$ when $b(\tau )\rightarrow 0$ and $a(t)\rightarrow \surd {3}/2 =0.8660\ldots .$ However, this behaviour persists only for so long as $b(\tau )\gtrsim h_{0} \,\hat {h}(\tau )$, and it is moreover subject to a ‘fingering instability’. Surface tension must also become important where the curvature of the boundary becomes large. These effects will be considered in § 5 below.

Figure 9. Evolution of an elliptic drop, with initial condition $a_{0}=1,b_{0}=0.5,\kappa =0.75$ (plan views in the $\{x,y\}$ plane:); (a) expanding drop when $F > 0$, shown at times $\tau =0,1,10,10^2,10^3$ and $10^4$; the shape becomes gradually more circular in form; (b) contracting drop when $F < 0$; the minor axis of the ellipse $b(\tau )$ decreases to zero at the singularity time $\tau =0.71672356\ldots$; the contracting boundary, shown at times $\tau = 0,0.6,0.7,0.7166,0.71672356$, is subject to smaller-scale instabilities.

4.2. Linearised stability analysis

We here consider a circular drop with perturbed boundary

(4.23)\begin{equation} R(\theta, 0) =a_{0} +\epsilon \alpha_{0} \cos n\theta \end{equation}

at the initial instant, where $n$ is a positive integer, and $\epsilon \alpha _{0}$ an arbitrary initial amplitude with $0<\epsilon \ll 1$. When $n=1$, the circle is merely displaced in the $x$-direction without distortion; we may therefore suppose that $n \ge 2$, so that (4.2) is automatically satisfied. (When $n=2$, the circle is perturbed to an ellipse.) With the choice (4.23), $A(0)= {\rm \pi}a_{0}^{2} +{\frac {1}{2}}{\rm \pi} \, \epsilon ^{2} \alpha _{0}^2$.

With $0<\epsilon \ll 1$, the problem may be linearised. Then for $t>0$, $R(\theta ,t) =a(t) +\epsilon \, \alpha (t) \cos n\theta$, where $a(t)$ is given by (2.16a,b) and $A(t)={\rm \pi} a^{2}(t) +\text {O}(\epsilon ^{2})$. Thus, let $p(r,\theta ,t) =p_{0}(r,t) +\epsilon \, p_{1}(r,t) \cos n\theta$, where $p_{0}(r,t)$ is given by (2.8). At this level of approximation, $h(t)$ is still given by (2.16a,b), and since $\nabla ^{2}p_{0} =12\,\mu \, h^{-3}\text {d}h/\text {d}t$, (2.6) gives

(4.24)\begin{equation} \frac{1}{r}\frac{\partial}{\partial r} r \frac{\partial p_{1}}{\partial r} -\frac{n^2}{r^2}p_{1} =0. \end{equation}

The solution, finite at $r=0$, is $p_{1}(r,t) =k(t)r^{n}$, so that now $p(r,\theta ,t)= p_{0}(r,t) +\epsilon \, k(t)r^{n} \cos n\theta$.

The pressure boundary condition is now $p(r,\theta ,t)-p_{a}=0$ on $r=a(t)+\epsilon \,\alpha \cos n\theta$, and, noting that

(4.25)\begin{align} p_{0}(a(t)&+\epsilon\, \alpha(t)\cos n\theta,t)-p_{a} = p_{0}(a)-p_{a}\nonumber\\ &\quad +\epsilon\, \alpha(t)\cos n\theta \left. \frac{\partial p_{0}}{\partial r}\right|_{r=a} = \epsilon\, a(t)\,\alpha(t)\cos n\theta\frac{6\mu }{h^3}\frac{\text{d}h}{\text{d}t}, \end{align}

and that

(4.26)\begin{equation} \epsilon\, p_{1}(a(t)+\epsilon\, \alpha(t)\cos n\theta,t) =\epsilon\, k(t)a^{n} + \text{O}(\epsilon^2), \end{equation}

this boundary condition (at order $\epsilon$) determines $k(t)$, leading finally to

(4.27)\begin{equation} p(r, \theta,t)-p_{a}=\left[(r^2 -a^2) -2 \epsilon \,a \,(r/a)^{n}\,\alpha(t)\cos n\theta\right]\frac{3\mu}{h^3}\frac{\text{d}h}{\text{d}t} +\text{O}(\epsilon^2), \end{equation}

and so

(4.28)\begin{equation} \frac{\partial p}{\partial r} =\left [r- \epsilon\,n\, \alpha(t) (r/a)^{n-1}\cos n\theta\right]\frac{6\mu}{h^3}\frac{\text{d}h}{\text{d}t}+\text{O}(\epsilon^2). \end{equation}

On $r=R(\theta ,t)=a+\epsilon \, \alpha \cos n\theta$, this therefore gives

(4.29)\begin{equation} \frac{\partial p}{\partial r} = [a-\epsilon\, (n-1) \,\alpha \cos n\theta]\frac{6\mu}{h^3}\frac{\text{d}h}{\text{d}t} +\text{O}(\epsilon^2). \end{equation}

Now the curve $C(t)$ moves according to the equation

(4.30)\begin{equation} \frac{\partial R}{\partial t} = \frac{\text{d}a}{\text{d}t}+ \epsilon \, \frac{\text{d}\alpha}{\text{d}t}\cos n\theta = \bar u_{r} (R, t) +\text{O}(\epsilon^2), \end{equation}

and $\bar u_{r}$ still satisfies (2.4), i.e.

(4.31)\begin{equation} \bar u_{r} ={-}\frac{h^{2}}{12\mu}\,\left. \frac{ \partial p}{\partial r}\right|_{r=a(t)+\epsilon\, \alpha(t) \cos n\theta}\,. \end{equation}

At leading order, (4.30) merely confirms that $a^{2}(t)\,h(t)=$ const. At order $\epsilon$, with $h(t)$ given by (2.16a,b), we obtain the equation for $\alpha (t)$,

(4.32)\begin{equation} \frac{\text{d}\alpha}{\text{d}t}= \frac{(n-1)}{2h}\frac{\text{d}h}{\text{d}t}\,\alpha ={-}\frac{(n-1)\,\alpha}{8(t+t_{0})}, \end{equation}

which integrates to give

(4.33)\begin{equation} \alpha(t)=\alpha_{0}(1+t/t_{0})^{-(n-1)/8}. \end{equation}

Since $n>1$, this ensures that, when $F>0$ and so $t_{0}>0$, the perturbation decays to zero as $t\rightarrow \infty$; and in fact the higher the value of $n$, the more rapid is this decay.

The stability can be best understood from consideration of the sketch of figure 10(a) which shows the particular situation when $n=5$ (a fivefold hypotrochoid). Just as for the case of an elliptic drop (figure 8), the restoring force is greater where the isobars are closer, thus tending to restore a circular shape. Figure 10(c) shows the resulting evolution of the drop boundary in which the slow return to a circular shape is very evident. It is here that our conclusion differs from that of Savina & Nepomnyashchy (Reference Savina and Nepomnyashchy2015), as mentioned in the introduction.

Figure 10. (a) Sketch of isobars $p=$ const., indicating the stabilising mechanism for a squeezed drop; the drop expands more rapidly where the isobars are closer together, thus tending to restore a circular shape; (b) similar sketch indicating the destabilising mechanism at the boundary of an air bubble, in which the pressure $p_{b}$ is greater than atmospheric pressure $p_{a}$ (as given by (3.2)); the expansion is again more rapid where the isobars are closer together, leading now to growth of the perturbation amplitude; (c) boundary of squeezed drop, corresponding to the situation in (a); $n = 5,\, \epsilon = 0.3$; times $\tau = t/t_{0} = 0, 2, 4, 8, \ldots , 128$; (d) destabilising effect of air bubble, corresponding to the situation in (b); the outer boundary (blue) still spreads under squeezing, but the inner boundary (brown) is unstable; $n = 5, \,\epsilon = 0.0003,\, \eta =0.5,\, \beta (0)/\alpha (0) = 64;\, \tau = 0, 2, 4, 6, 8,10,12,14,16$.

We note that, for general $h(t)$, the first equality in (4.32) integrates to give

(4.34)\begin{equation} \frac{\alpha(t)}{\alpha(0)}=\left[\frac{h(t)}{h(0)}\right]^{(n-1)/2}=\left[\frac{a(t)}{a(0)}\right]^{-(n-1)/4}. \end{equation}

Thus the perturbation amplitude $\alpha (t)$ increases or decreases, implying instability or stability, according as $h(t)$ increases or decreases (or equivalently, according as $a(t)$ decreases or increases) respectively.

6.1. The situation when $F<0$

When $F<0$, with the basic state given by (3.7) there are still two perturbation modes, one stable and one unstable, but the unstable mode is now ‘driven’ by the fingering instability at the outer boundary $r=a(t)$, where the amplitude is much greater.

We note that, as already indicated in the introduction, the situation here is very different from that considered by Kanhurkar et al. (Reference Kanhurkar, Patankar, ul Islam and Gandhi2019), who conducted the plate-lifting experiment with a small hole in the upper plate, allowing for equalisation of pressure in the air bubble with atmospheric pressure. Air can flow through the hole, so the volume of the bubble is no longer constant as the lifting proceeds; this affects the pressure distribution throughout the viscous fluid annulus. Our situation is different in that the volume of the bubble is constant, and the pressure in the bubble is time dependent. This has an important effect on the stability problem in that, for either squeezing or lifting, there are two modes of perturbation of the basic time-dependent state, one stable and one unstable. In both cases, the unstable mode involves growth of the perturbation at both inner and outer boundaries of the viscous annulus.

7.1. Single expanding drop

As already mentioned, our ‘home experiments’ were carried out under the lockdown conditions of the current coronavirus epidemic. The first experiment was designed to verify the one-eighth power law (2.16a,b). The plates, of toughened glass, measured $420\ \textrm {mm} \times 280\ \textrm {mm}\times 4\ \textrm {mm}$ and weighed approximately 1.127 kg. They were certainly strong enough to resist any deformation resulting from pressure variations in the fluid. For fluid we used Lyle's black treacle (http://fiches.ranson.be/00006253nl-be.pdf) (density 1.41 $\textrm {g}\ \textrm {ml}^{-1}$, viscosity $\mu \approx 64.3\, \text {kg} \text { s}^{-1}\text {m}^{-1}$ at $22\,^{\circ }\textrm {C}$). A drop of volume $V= 5$ ml was placed on one horizontal plate and the second plate was placed on this and allowed to settle, so that the downward force was $F= Mg = 11.04 \,{\text {kg}\, \text {m}\,\text {s}^{-2}}$; thus $\mu /F = 5.82\ \text {m}^{-2}\ \text{s}$. The plates were illuminated from below, and a camera was set to take shots from above at 10 second intervals. the initial time $t=0$ was set when the drop radius was $a_{0}=30.72$ mm; thus $h_{0}= V/{\rm \pi} a_{0}^{2}=1.69$ mm and $h_{0}/a_{0}=0.055$. From (2.11), the time scale $t_{0}$ for this run was thus

(7.1)\begin{equation} t_0= \frac{3{\rm \pi}^3\,(5.82) (0.03072)^8}{8(5\times 10^{{-}6})^2} = 2.15\,\text{s}. \end{equation}

Figure 1 shows photos of the expanding drop at three times $\tau =t/t_{0}=$ 60, 172 and 321. The background grid allowed measurement of the diameter in each case. Figure 7(a) shows the growth of the drop radius for $0\le \tau \le 321$. The measurements agree well with the result (2.16a,b) as shown by the red curve; and figure 7(b) shows that the one-eighth power law is remarkably well satisfied over this full time range.

7.2. Contraction of drop by leverage

Figure 2 shows the effect of gently levering the plates open at the lower right-hand corner some time after the stage of figure 1(c) was reached. The fingering instability immediately develops and the fingers grow as the levering is gradually increased. The duration of course depends on the speed at which the angle between the plates is increased, and as the levering was applied by hand, exact repetition was not possible. This particular sequence was taken over a period of approximately one minute, and is quite typical. The characteristic tip-splitting and side-branching of fingers is evident, and cavitation bubbles can be detected in the very low-pressure region in panel (e). The reason for this cavitation may be understood with reference to the ‘hinged-plate’ problem briefly considered in § 2 of Moffatt (Reference Moffatt1964); if two plates $\theta =\pm \alpha (t)$ are hinged at $r=0$ with viscous fluid in the gap $-\alpha (t)<\theta <\alpha (t)$, and if $\alpha (t)$ is slowly increased from a very small value, then the pressure $p(r,t)$ behaves like $\mu \,\dot {\alpha }(t) \log r$ which becomes negative as $r$ decreases towards zero, so that cavitation is inevitable in this region. This is of course a two-dimensional idealisation, but the same physical mechanism is presumably responsible for the cavitation bubbles appearing in the present experiment.

Note that, with $F\sim 11\,{\text {kg}\, \text {m}\,\text {s}^{-2}}$ and taking $\gamma \sim 0.07$ kg s$^{-2}$ and $a_0 \sim 0.06\, \text {m}$ (as at the start of the leverage), the critical value $n_c$ defined in (5.7) evaluates to $n_{c} = (|F|/\gamma a_{0})^{1/3}\sim 14$. Thus we might expect of the order of 14 crests to emerge as the instability develops; given that this is simply an order of magnitude estimate and that the development soon becomes nonlinear, this is not inconsistent with what is evident in figure 2.

Note further that, as a consequence of the no-slip condition, a thin residual film is left on both plates in the regions invaded by the fingers. The thickness $s$ of these residual films is presumably given by $s/h \sim \mu U /\gamma$, the result obtained by Bretherton (Reference Bretherton1961) for a bubble rising in a capillary tube containing viscous liquid. In our present situation, with $\mu \sim 70\, \text {kg} \text { s}^{-1}\text {m}^{-1}$, $\gamma \sim 0.07$ kg s$^{-2}$ and the observed creeping velocity $U$ in the range $10^{-4}$–$10^{-5}~$m s$^{-1}$, this gives $s/h\sim 10^{-1}$–$10^{-2}$. When $h\sim 0.1$mm, the residual film thickness $s$ is therefore in the range 1–10 $\mu$m.

7.3. Evolution of an annular drop

Figure 3 shows an evolution sequence when the initial shape of the drop (panel a) is as near to a uniform circular annulus as could be achieved by hand. The curvature of the inner boundary of the annulus was inevitably not quite uniform, and where it was maximal a perturbation immediately developed, as anticipated by the stability analysis of § 6. This perturbation evolves into a secondary bubble (panel b) which is ejected, leaving the primary bubble (panel c) which deforms slowly to a distinctly non-circular shape (panel d); the outer boundary was only weakly perturbed from a circle, The secondary bubble detaches from the primary bubble, involving the breaking of a thin thread of fluid, and erupts through the outer boundary (panels d,e) while a new (tertiary) bubble grows very slowly from the residual cusp on the primary bubble; (thread breaking and cusp formation were possibilities anticipated by Shelley et al. Reference Shelley, Tiany and Wlodarski1997). After approximately 45 min, the situation (panel f) appears to be quasi-static, and very small wrinkles appear on the outer boundary, possibly associated with slight roughness of the two plates. After a further hour, leverage was introduced at the lower right-hand corner, and fingering rapidly developed (panels g–i) over about one minute. A growing finger punctures the tertiary bubble (panel g) at the point marked by the arrow, providing a connected path from the primary bubble to the exterior. Fingers also erupt from the bubble (panel j) and encounter the ‘invading’ fingers, forming a defensive ridge. The pattern continues to evolve round the outer boundary (panel k), and cavitation bubbles may be detected in the ‘north–western’ sector ahead of the advancing fronts. Panel (l) shows the ultimate fingering pattern just before rupture of the film when the plates separate completely; the ridge, which remains a prominent feature of this pattern, is breached only at the location of the puncture point.

7.4. Merging of four drops

We then carried out an experiment starting with four drops placed as near to the corners of a square as could be achieved, as shown in figure 4 (panel a). As the upper plate descends, these drops expand and ultimately merge, trapping an air bubble (panel e). The initial drop contact (panel c) involves an instantaneous geometrical singularity at a time $\tau =\tau _{c}$, which is subsequently (i.e. for $\tau >\tau _{c}$) resolved by surface tension. (We should distinguish such a singularity from the more familiar finite-time singularities, as treated for example by Eggers & Fontelos (Reference Eggers and Fontelos2015), for which a singularity develops as $\tau \rightarrow \tau _{c}$ and is resolved by surface tension or some other mechanism just before the singularity occurs, i.e. for $\tau <\tau _{c}$.) After the bubble is trapped, it shows no tendency to become circular; on the contrary, where the curvature on the air bubble is maximal, a secondary bubble begins to pinch off (panel g), again leaving the very thin residual films of treacle on both plates.

The formation of the secondary bubble is consistent with the qualitative description of the bubble instability in § 6 above: the pressure contours are compressed where the bubble curvature is large, leading to amplification of the perturbation in this region. The secondary bubble slowly propagates towards the outer boundary, separates from the primary bubble (a singular process analogous to the breaking of a viscous thread), and ultimately erupts through this boundary (panel i) – again a singular process. Up to this point, the total volume of trapped air is constant, but this falls very rapidly at the moment of eruption. At the same time, a tertiary bubble begins to form under the continuing influence of the pressure gradient; this suggests an iterative process, although much slower at each stage.

7.5. Contraction by levering

Figure 5 shows a continuation of the sequence of figure 4, showing first (panels a–c) the eruption of the secondary bubble and the very slow separation of the tertiary bubble from the primary bubble which is now much reduced in size. Note that the outer boundary remains relatively smooth but does not become circular; this is because the internal instability has a persistent weak disturbing effect on the outer boundary (see the comment following (6.17a,b)). In panel (d), levering has just commenced at the lower right-hand corner of the plates, and the first signs of the fingering instability are visible. In panel (e), fingers erupting from the outer boundary of the drop are evident, and tip-splitting and side-branching of the invading fingers are well developed; in panel (f) a ridge develops between the invading and defending cohorts of fingers. In panel (g) the defending cohort, now in retreat, is punctured by an invading finger; the pressure is now atmospheric ($p_{a}$) throughout the fingers and the bubbles, but less than atmospheric in the remaining bulk of treacle; in panel (h) the invading front of fingers advances round the boundary and cavitation bubbles appear ahead of the front. Finally, panel (i) shows that fingering has again extended throughout the drop, to the extent that the treacle is almost entirely contained in the narrow ‘tree branches and twigs’ that separate the fingers. The ultimate pattern retains an imprint of the original trapped bubbles, in the pronounced pentagonal ‘ridge’, still breached at only one point in the south–east sector.