2.1. Governing equations

Following previous studies (Ben Amar & Bonn Reference Ben Amar and Bonn2005; Lindner et al. Reference Lindner, Derks and Shelley2005; Sinha et al. Reference Sinha, Dutta and Tarafdar2008; Nase et al. Reference Nase, Derks and Lindner2011; Dias & Miranda Reference Dias and Miranda2013a), we assume that the motion of the fluid is governed by Darcy's law, which is a two-dimensional approximation of the Navier–Stokes system obtained by averaging the equations over the narrow gap between the plates. Following Lamb (Reference Lamb1932), He & Belmonte (Reference He and Belmonte2011) and Anjos, Dias & Miranda (Reference Anjos, Dias and Miranda2017), the two-dimensional approximation holds when $\mathcal {E}={\tilde {b}(\tilde {t})}/{\tilde {R}(\tilde {t})}\ll 1$ and Reynolds number $Re=({\rho \tilde {b}(\tilde {t}) \dot {\tilde {b}}(\tilde {t})})/{12 \mu }\ll 1$, where $\tilde {R}$ is the equivalent radius of the fluid drop, $\dot {\tilde {b}}(\tilde {t})={\textrm {d}\tilde {b}(\tilde {t})}/{\textrm {d}\tilde {t}}$ is the lifting speed of the upper plate, $\rho$ is the density of the fluid and $\mu$ is the viscosity of the fluid. The equations are

(2.1)\begin{equation} \tilde{\boldsymbol{u}}={-}\frac{\tilde{b}^2(\tilde{t})}{12\mu}\tilde{\boldsymbol{\nabla}} \tilde{P} \quad \text{in} \quad\tilde{\varOmega}(\tilde{t}),\end{equation}

where $\tilde {\boldsymbol {u}}$ is the velocity and $\tilde {P}$ is the pressure.

From volume conservation, we obtain the gap-averaged incompressibility condition as

(2.2)\begin{equation} \tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\tilde{\boldsymbol{u}}={-}\frac{\dot{\tilde{b}}(\tilde{t})}{\tilde{b}(\tilde{t})}\quad \text{in} \ \tilde{\varOmega}(\tilde{t}). \end{equation}

The pressure jump $[\tilde {P}]$ across the interface is given by the Laplace–Young condition, which is the product of surface tension $\sigma$ and the curvature $\mathcal {K}$ (Park & Homsy Reference Park and Homsy1984; Dias & Miranda Reference Dias and Miranda2013a; Anjos & Miranda Reference Anjos and Miranda2014), as follows:

(2.3)\begin{equation} [\tilde{P}]=\tilde{P}|_{\partial\varOmega^+}-\tilde{P}|_{\partial\varOmega^-}=\sigma\mathcal{K} \quad \text{on}\quad \partial\tilde{\varOmega}(\tilde{t}), \end{equation}

where the $+/-$ superscripts denote the limit from the exterior and interior of $\varOmega$, respectively. We approximate the curvature as $\mathcal {K}=\tilde {\kappa }+({2}/{\tilde {b}(\tilde {t})})\cos (\theta _c)$, where $\tilde {\kappa }$ denotes the curvature of the planar interfacial curve and $\theta _c$ represents the contact angle measured between the plates and the curved meniscus. The second term accounts for the curvature across the gap width, which is associated with the interface profile in the z-direction. When $\tilde {b}(\tilde {t})$ increases, this term decreases indicating that surface tension effects across the gap become weaker. However, $\tilde {\kappa }$ increases since the interface shrinks. Even though the second term varies with time, it is uniform in space if we assume $\theta _c$ is a constant. Thus it can be absorbed in the pressure as a redefined pressure. We use this redefined pressure in the rest of this paper. The scaled normal derivative $({1}/{\mu })({\partial \tilde {P}}/{\partial \boldsymbol {n}})$ is continuous across $\partial \varOmega$ and the normal velocity of the interface is thus

(2.4)\begin{equation} \tilde{V}={-}\frac{\tilde{b}^2(\tilde{t})}{12\mu}\frac{\partial \tilde{P}}{\partial \boldsymbol{n}}\quad \text{on}\ \partial\tilde{\varOmega}(\tilde{t}), \end{equation}

where $\boldsymbol {n}$ is the unit normal vector pointing into $\tilde {\varOmega }(\tilde {t})$.

We non-dimensionalize the system using a characteristic length $L_0$, time $T={\tilde {b}_0}/{\dot {\tilde {b}}_0}$ and pressure $P_0=({12\mu L_0^2})/({T\tilde {b}_0^2})$, where $L_0$ is a characteristic size of the initial fluid domain $\varOmega (0)$, and $\tilde {b}_0$ and $\dot {\tilde {b}}_0$ are the initial values of $\tilde {b}$ and $\dot {\tilde {b}}$. Further, we define the non-dimensional modified pressure as ${P}=\tilde {P}/P_0- ({\dot {b}(t)}/{4b^3(t)})|\boldsymbol {x}|^2$, where $b(t)=\tilde {b}(\tilde {t})/b_0$, $t=\tilde {t}/T$, ${\dot {b}(t)={\textrm {d}b}/{\textrm {d}t}}$ and $\boldsymbol {x}=\tilde {\boldsymbol {x}}/L_0$. Then, the non-dimensional version of (2.1)–(2.4) becomes

(2.5)\begin{gather} \nabla^2 {P}=0\quad \text{in}\ \varOmega, \end{gather}

(2.6)\begin{gather}{}[{P}]=\tau \kappa-\frac{\dot{b}(t)}{4b^3(t)}|\boldsymbol{x}|^2\quad \text{on}\ \partial\varOmega, \end{gather}

(2.7)\begin{gather} V={-}b^2(t)\frac{\partial{P}}{\partial \boldsymbol{n}}+\frac{\dot{b}(t)}{2b(t)}\boldsymbol{x}\boldsymbol{\cdot}\boldsymbol{n}\quad\textrm{on}\ \partial\varOmega, \end{gather}

where $\tau ={\sigma \tilde {b}_0^3}/({12\mu \dot {\tilde {b}}_0L_0^3})$ is a non-dimensional surface tension. Taking $L_0$ to be the equivalent radius of $\varOmega (0)$, e.g. the radius of a circle with the same enclosed area, the non-dimensional volume of the fluid is ${\rm \pi}$ since the initial non-dimensional gap is $b(0)=1$.

2.2. Linear theory

In this section, we briefly review the linear stability analysis in Shelley et al. (Reference Shelley, Tian and Wlodarski1997) and Zhao et al. (Reference Zhao, Li, Yin, Belmonte, Lowengrub and Li2018). We consider the interface to be a slightly perturbed circle,

(2.8)\begin{equation} r(\alpha,t)=R(t)+\epsilon \delta(t)\cos(k\alpha), \end{equation}

where $\epsilon$ $\ll 1$, the perturbation mode $k\geq 2$ is an integer, $\alpha \in [0,2{\rm \pi} ]$ is the polar angle and $\delta (t)$ is the amplitude of the perturbation. The shape factor $({\delta }/{R})(t)={\delta (t)}/{R(t)}$ can be used to characterize the size of the perturbation relative to the underlying circle (Mullins & Sekerka Reference Mullins and Sekerka1963). Then, it can be shown that the shape factor evolves according to

(2.9)\begin{equation} \left(\frac{\delta}{R}\right)^{{-}1}\frac{\textrm{d}}{\textrm{d}t}\left(\frac{\delta}{R}\right)=\frac{\dot{b}k}{2b}- \tau\frac{b^2(k^3-k)}{R^3}. \end{equation}

Using the relationship $R(t)={1}/{\sqrt {b(t)}}$, which arises from volume conservation at the level of linear theory, the shape factor grows only when

(2.10)\begin{equation} \dot{b}(t)>2\tau (k^2-1)b^{9/2}(t). \end{equation}

This indicates that the band of modes satisfying ${|k|<\sqrt {{\dot {b}b^{-9/2}}/{2\tau }+1}}$ are unstable and the interface can develop fingering patterns. Note that in the special cases that $\dot {b}=1$, or even $\dot {b}=b$, all modes $|k|>1$ become stable in the long time limit. This occurs because ${\dot {b}b^{-9/2}}\to 0$ as $t\to \infty$. On the other hand, if we consider much faster rates of gap increase,

(2.11)\begin{equation} b(t)=\left(1-\tfrac{7}{2}\tau \mathcal{C}t\right)^{-{2}/{7}} , \end{equation}

where $\mathcal {C}$ is a constant, then the band of unstable modes is fixed in time and depends on $\mathcal {C}$ where ${|k|<\sqrt {\mathcal {C}/2+1}}$. Note this gap tends to infinity at the finite time $T_{\mathcal {C}}=2/(7\tau \mathcal {C})$.

We find the fastest growing mode $k_{max}$ by setting

(2.12)\begin{equation} \frac{\textrm{d}}{\textrm{d}k}\left[\left(\frac{\delta }{R}\right)^{{-}1}\frac{\textrm{d}}{\textrm{d}t}\left(\frac{\delta}{R}\right) \right]=0 \end{equation}

and solving for $k_{max}$. The maximum perturbation mode $k^*$ has the largest magnitude relative to its initial value. We get $k^*$ by solving

(2.13)\begin{equation} \frac{\textrm{d}}{\textrm{d}k}\left(\left(\frac{\delta}{R}\right)_0^{{-}1}\boldsymbol{\cdot}\left(\frac{\delta}{R}\right)\right)=0, \end{equation}

where

(2.14)\begin{equation} \frac{\delta}{R}(t)=\left(\frac{\delta}{R}\right)_0b^{{k}/{2}}(t)\exp\left({\tau (k-k^3)\int_0^tb^{7/2}(s)\,\textrm{d}s}\right) \end{equation}

and $({\delta }/{R})_0$ is the initial perturbation. Thus, mode $k^*$ is given by

(2.15)\begin{equation} k^*=\sqrt{\frac{\ln{b(t)}}{\displaystyle 6\tau \int_0^t b^{7/2}(s)\,\textrm{d}s}+\frac{1}{3}}. \end{equation}

When $({\delta }/{R})_0$ is independent of $k$, then $k^*$ corresponds to the mode with the largest amplitude. When $({\delta }/{R})_0$ depends on $k$, such as the exponential decay in (3.1), the mode with largest amplitude can be determined by solving $({\textrm {d}}/{\textrm {d}k})({\delta }/{R})=0$. Referring to this mode as $\tilde {k}^*$, we obtain

(2.16)\begin{equation} \tilde{k}^*=\sqrt{\frac{\ln{b(t)-2\beta}}{\displaystyle 6\tau \int_0^t b^{7/2}(s)\,\textrm{d}s}+\frac{1}{3}}, \end{equation}

where $\beta$ is the decay rate of initial perturbation with respect to $k$. Comparing the predictions from $k_{max}$, $k^*$ and $\tilde {k}^*$ with experimental observations for the number of fingers, we find that at early times the number of fingers predicted from $k^*$ is more consistent with the experimental results. This suggests that the largest initial perturbation modes have nearly uniform magnitudes. At later times, both $k^*$ and $\tilde {k}^*$ behave similarly, decay in time and are more consistent with experimental results than $k_{max}$.

Note that because the growth rate of the shape factor depends on time, $k^*$ need not be equal to $k_{max}$. Further, these have been used to estimate the number of fingers (Nase et al. Reference Nase, Derks and Lindner2011; Dias & Miranda Reference Dias and Miranda2013a).

When the special gap dynamics in (2.11) is used, the fastest growing mode $k_{max}$ can be prescribed by taking $\mathcal {C}=2(3k_{max}^2-1)$. Under this gap dynamics, the mode $k_{max}$ will remain the fastest growing at all times so that $k_{max}=k^*$ and the shape factor is

(2.17)\begin{equation} \frac{\delta(t)}{R(t)}= \frac{\delta(0)}{R(0)}b(t)^{{k_{max}^3}/({3k^2_{max}-1})} =\frac{\delta(0)}{R(0)}R(t)^{-({2k_{max}^3}/({3k^2_{max}-1}))}, \end{equation}

which diverges as $t\to T_C$. Alternatively, if $\mathcal {C}=2(\bar {k}^2-1)$, then the shape factor of mode $\bar {k}$ does not change in time. In other words, under linear theory, a $\bar {k}$-mode perturbation of the interface would evolve self-similarly. Both of these conditions suggest that in the special gap regime, there may be mode selection and that perturbations may persist (and even grow) as the fluid domain and interface shrinks.

2.3. Numerical method

Because of the strong nonlinearity and non-locality of the lifting plate equations, numerical methods are needed to characterize the nonlinear dynamics. However, the simulations are very challenging because of severe time and space step restrictions introduced by the rapid evolution and shrinking of the interface. To overcome these numerical issues, we have developed a rescaled boundary integral scheme (Zhao et al. Reference Zhao, Li, Yin, Belmonte, Lowengrub and Li2018); the method is briefly described here. The idea is to map the original time and space $({\boldsymbol {x}},t)$ into new coordinates $(\bar {\boldsymbol {x}},\bar {t})$ such that the interface can evolve at an arbitrary speed in the new rescaled frame (see also Li et al. Reference Li, Lowengrub and Leo2007; Zhao et al. Reference Zhao, Belmonte, Li, Li and Lowengrub2016, Reference Zhao, Yin, Lowengrub and Li2017). Introduce a new frame $(\bar {\boldsymbol {x}},\bar {t})$ such that

(2.18)\begin{gather} \boldsymbol{x}=\bar{R}(\bar{t})\bar{\boldsymbol{x}}(\bar{t},\alpha), \end{gather}

(2.19)\begin{gather} \bar{t}=\int_{0}^t \frac{1}{\rho(t')}\,\textrm{d}t', \end{gather}

where the space scaling $\bar {R}(\bar {t})$ captures the size of the interface, $\bar {\boldsymbol {x}}$ is the position vector of the scaled interface and $\alpha$ parameterizes the interface. The space scaling function maps the interface to its original size and the time scaling function $\rho (t)=\bar {\rho }(\bar {t})$ maps the original time $t$ to the new time $\bar {t}$, which has to be positive and continuous. If $\rho (t)<1$, then the evolution in the rescaled frame is slower than that in the original frame. Using mass conversation and requiring the area enclosed by the interface to be fixed in the new frame, we obtain the normal velocity in the new frame

(2.20)\begin{equation} \bar{V}= \frac{\bar{\rho}}{\bar{R}}{V}(t). \end{equation}

Representing the pressure $P$ as a double layer potential, with dipole density $\bar {\gamma }$, (2.5) and (2.6) can be written as the following Fredholm integral equation of the second kind:

(2.21)\begin{equation} \bar{\gamma}(\bar{\boldsymbol{x}})+\frac{1}{\rm \pi}\int_{\partial\bar{\varOmega}(\bar{t})}\bar{\gamma}(\bar{\boldsymbol{x}}')\left[\frac{\partial \ln|\bar{\boldsymbol{x}}-\bar{\boldsymbol{x}}'|}{\partial \boldsymbol{n}(\bar{\boldsymbol{x}}')}+\bar{R}(\bar{t})\right]d\bar{s}(\bar{\boldsymbol{x}}')=2\tau \bar{\kappa}-\frac{\dot{b}(t(\bar{t}))}{2b^3(t(\bar{t}))}\bar{R}^{3}|\bar{\boldsymbol{x}}|^2. \end{equation}

Once $\bar {\gamma }$ is determined, the normal velocity in the new frame $\bar {V}$ can be computed as

(2.22)\begin{equation} \bar{V}(\bar{\boldsymbol{x}})={-}\frac{b^2(t(\bar{t}))\bar{\rho}}{2{\rm \pi} \bar{R}^{3}}\int_{\partial\bar{\varOmega}}\bar{\gamma}_{\bar{s}}\frac{(\bar{\boldsymbol{x}}'-\bar{\boldsymbol{x}})^{{\perp}}\boldsymbol{\cdot}{\boldsymbol{n}}(\bar{s})}{|\bar{\boldsymbol{x}}'-\bar{\boldsymbol{x}}|^2}\,\textrm{d}\bar{s}', \end{equation}

where $\bar {\boldsymbol {x}}^{\perp }=(\bar {x}_2,-\bar {x}_1)$ and the interface evolution in the scaled frame is

(2.23)\begin{equation} \frac{\textrm{d}\bar{\boldsymbol{x}}(\bar{t},\alpha)}{\textrm{d}\bar{t}}\boldsymbol{\cdot} \boldsymbol{n}=\bar{V}(\bar{t},\alpha). \end{equation}

Here, we take the time scaling $\rho (t)\propto 1/\dot {b}$ so that in the rescaled frame, the gap increases linearly in time, e.g. $b(\bar {t})=1+c\bar {t}$, where $c$ is a constant. Using a spectrally accurate discretization method in space, a second-order accurate semi-implicit scheme in time and the generalized minimum residual scheme (known as GMRES) (Saad & Schultz Reference Saad and Schultz1986) to solve the Fredholm integral equation of the second kind, the method enables us to accurately compute the dynamics to far longer times than could previously be accomplished. Further details of the numerical method, including convergence studies, can be found in Zhao et al. (Reference Zhao, Li, Yin, Belmonte, Lowengrub and Li2018).

3.1. Comparison between linear theory, nonlinear simulations and experiments

We first compare the experimental results obtained in Nase et al. (Reference Nase, Derks and Lindner2011) with the results of linear theory (Dias & Miranda Reference Dias and Miranda2013a; Zhao et al. Reference Zhao, Li, Yin, Belmonte, Lowengrub and Li2018) and nonlinear simulations. In the experiments, a drop of viscous fluid, surrounded by air, was placed in a Hele-Shaw cell and the upper plate was pulled upward at a constant rate of increase, e.g. the non-dimensional gap width $b(t)=1+t$, where $t$ is the non-dimensional time. High resolution images were used to precisely determine the number of fingers over time. By varying the initial drop radius, the initial plate spacing, the rate of gap increase, and the viscosity, they systematically investigated how the number of fingers depends on the non-dimensional surface tension and confinement number (see Nase et al. (Reference Nase, Derks and Lindner2011) for details).

The results from two sets of experiments are shown in figure 2. In figure 2(a,b) the number of fingers is shown as a function of the non-dimensional time for different confinement numbers $C_0$, as labelled. The non-dimensional surface tensions are $\tau =3\times 10^{-5}$ (figure 2a) and $\tau =9.6\times 10^{-6}$ (figure 2b). The solid and dotted curves denote the number of fingers estimated from $k^*$ and $k_{max}$, respectively, obtained from linear theory. The stars denote the results from nonlinear simulations. The number of fingers in the simulations are calculated in exactly the same way as in Nase et al. (Reference Nase, Derks and Lindner2011) (e.g. number of air fingers penetrating the viscous fluid; see appendix A, figure 8). Characteristic drop morphologies using $\tau =9.6\times 10^{-6}$ are shown (at the non-dimensional times $t=1$, $2$, $3$) in figure 2(c–h) from simulations (c–e) and experiments (f–h) with $C_0=120$. The agreement between the experimental and simulation morphologies is striking. In our simulation, the interface eventually evolves into a circle, see figure 3 for the morphologies of the interface. Although the final gap width is approximately 16 times the initial gap width, Darcy approximation is valid ($\mathcal {E}\ll 1$) throughout the evolution.

Figure 2.

Comparison between linear theory, numerical simulations and experiments from Nase et al. (Reference Nase, Derks and Lindner2011) where the gap width is increased linearly in time. The numbers of fingers over time are shown in panels (a,b) predicted by linear theory ($k_{max}$ dotted curves; $k^*$ solid curves), numerical simulations (magenta stars) and experiments with different confinement numbers $C_0$, ranging from $30$ to $120$ as labelled. The dashed lines are exponential fits to the numerical simulations. In panel (a) the non-dimensional surface tension is $\tau =3\times 10^{-5}$ while in panel (b) $\tau =9.6\times 10^{-6}$. In panels (c–h), the simulated interfacial morphologies (c–e) and those from experiments (f–h), with $\tau =9.6\times 10^{-6}$, are shown at the dimensionless times $t=1,2,$ and $3$. There is an excellent agreement between the numerical simulations and experiments. Experimental results reprinted with permission.

Figure 3.

The simulated interfacial morphologies with $\tau =9.6\times 10^{-6}$ are shown at the dimensionless times $t=0$ (blue) ,$1$ (red), $12$ (black) and $16$ (green).

Linear theory is not able to fully describe the dynamics of the interface. But linear theory should give an insight into the dynamics. At early stages, when the perturbation is small, linear theory agrees well with nonlinear simulations. At long times, linear theory predicts that all modes vanish. In the experiments shown in Nase et al. (Reference Nase, Derks and Lindner2011), where the gap width increases linearly in time, the experiments demonstrate that the droplets vanish as a circle (see figures 3 and 4 in Nase et al. (Reference Nase, Derks and Lindner2011)), which is consistent with linear theory. In these figures, the confinement number $C_0=60$ while we choose to present the experimental observations with confinement number $C_0=120$ (figure 6 in Nase et al. (Reference Nase, Derks and Lindner2011)) as figure 2(c–h). The droplet with confinement number $C_0=120$ also vanishes as a circle although dynamics at late times are not presented in Nase et al. (Reference Nase, Derks and Lindner2011).

To produce the nonlinear simulation results, we took the initial shape of the drop to be a slightly perturbed circle,

(3.1)\begin{equation} r(\alpha,0)=1+\epsilon\sum_{k=k_{min}}^{k_N}\textrm{e}^{-\beta k}(a_k\cos(k\alpha)+b_k\sin(k\alpha)), \end{equation}

where each $k$, $a_k$ and $b_k$ are chosen randomly from a uniform distribution in the interval $(-1,1)$ and $\alpha$ parameterizes the interface. Varying $\epsilon$, $\beta$ and $k_N$ allows us to vary the amplitude of the perturbation and its modal content. Here, we took $\epsilon =0.05$, $\beta =0.2$, $k_{min}=2$, and $k_N$ is varied between 40 and 100. Since $b(t)\sim t$, we do not need to rescale time to slow down the evolution and take $\bar {t}=t$. We use N= 4096 mesh points along the interface, the time step $\Delta \bar {t}=1\times 10^{-4}$ and the surface tension $\tau =9.6\times 10^{-6}$ and $3\times 10^{-5}$.

As can be seen in figure 2, the number of fingers decreases in time and the drop eventually shrinks as a circle. The experimental results depend on both the confinement number $C_0$ and the non-dimensional surface tension $\tau$. At larger surface tensions (figure 2a), the effect of $C_0$ is more pronounced at later times while at small surface tensions (figure 2b), $C_0$ more significantly influences the early stages of the evolution. Generally, the larger $C_0$ is the more fingers that are observed.

The number of fingers predicted using $k_{max}$ significantly underpredicts the experimental (and simulation) results while the $k^*$ predictions agree better with the experimental data at small $C_0$. The numerical simulations predict more fingers than either linear estimate, and agree best with experimental data at large $C_0$. This illustrates the importance of accounting for nonlinear interactions. Indeed, in Zhao et al. (Reference Zhao, Li, Yin, Belmonte, Lowengrub and Li2018), it was shown that linear theory underpredicts the amplitudes of growing modes in the lifting plate problem.

For both surface tensions, the numerical simulations, and experimental data with large $C_0$, predict a biphasic, exponential decay of the number of fingers over time; rapid decay at early times (e.g. $\approx \textrm {e}^{-0.75t}$ in figure 2b) and slower decay over late times (e.g. $\approx \textrm {e}^{-0.26t}$ in figure 2b). The experimental data with small $C_0$ seem to predict a single rate of exponential decay (see appendix A, figure 9), consistent with the experimental results obtained in Ben Amar & Bonn (Reference Ben Amar and Bonn2005).

It is still not yet well understood why the confinement number influences the number of fingers. Note that the thickness $h$ of the wetting layer on the plates scales as $h \sim O(\textit {Ca}^{2/3})$ where $\textit {Ca}=\mu \tilde {V}/\sigma$ is the Capillary number (Park & Homsy Reference Park and Homsy1984; Park et al. Reference Park, Gorell and Homsy1984; Jackson et al. Reference Jackson, Stevens, Giddings and Power2015). Estimating $\tilde {V}\sim (\dot {b}_0/b_0) R_0= \dot {b}_0 C_0$ from (2.2) and the definition of $C_0$, we obtain $Ca\sim \dot {b}_0 C_0$. In the experiments in Nase et al. (Reference Nase, Derks and Lindner2011), the confinement number and lifting rate are changed such that the non-dimensional surface tension $\tau =\sigma \tilde {b}_0^3/(12\mu \dot {\tilde {b}}_0 R_0^3)$ is fixed. This implies that $\dot {b}_0\sim C_0^{-3}$ and therefore $Ca\sim C_0^{-2}$. Thus, the wetting layer thickness $h\sim C_0^{-4/3}$ increases as the confinement number decreases. As a simple test of this, we modified the lifting speed to reflect the fact that more fluid may be left on the plates when the confinement number is decreased, e.g. the lifting speed of the upper plate is increased to reflect the more rapid loss of fluid to the wetting layer. The results are shown in figure 10 in appendix A and indicate that with this modification, the nonlinear simulations are better able to fit experiments with smaller confinement numbers. Although, our approach for incorporating wetting effects in the model is ad-hoc, the results correlate with observed experimental trends. To more accurately account for wetting effects, we may follow the analysis in Park & Homsy (Reference Park and Homsy1984), Dias & Miranda (Reference Dias and Miranda2013a) and Anjos & Miranda (Reference Anjos and Miranda2014), where the Laplace–Young condition can be extended to account for wetting effects. This is beyond the scope of this paper and is not the subject of current research.

When wetting effects, as well as viscous stresses, are incorporated in the modified Laplace–Young condition, the range of agreement between linear predictions, using the maximum perturbation wavenumber $k^*$, and experiments can be extended to somewhat larger $C_0$. However, there is still disagreement between the two at large $C_0$ with weakly nonlinear theory underpredicting the number of fingers, see Dias & Miranda (Reference Dias and Miranda2013a). It is very likely that nonlinear effects also play a key role and the development of a numerical method to accurately account for nonlinearity, viscous stresses, wetting effects and fluid motion in three dimensions is a subject for future work.

3.2. Limiting dynamics in the special gap regime

We next study the interfacial dynamics in the special gap regime using $b(t)=(1-({7}/{2})\tau \mathcal {C}t)^{-{2}/{7}}$, which diverges at the finite time $T_C=2/(7\tau \mathcal {C})$. Recall that in this regime, according to linear theory (e.g. see § 2.2), the fastest growing mode $k_{max}$ and $\mathcal {C}$ are related by $k_{max}=\sqrt {(1+\mathcal {C}/2)/3}$. Further, $k_{max}=k^*$, the maximum perturbation mode, and the shape factor $\delta /R\propto R^{-2k_{max}^3/(3*k_{max}^2-1)}$, which monotonically increases over time, and diverges at $t=T_C$. Because of the very rapid dynamics, an infinitesimally small time step would be required as $t \to T_C$ if the time were not rescaled. This would make it virtually impossible to simulate the drop dynamics in the original frame. Consequently, we rescale time to slow down the evolution and take $\rho (t)=c/\dot {b}(t)$ so that in the new time scale $\bar {t}$ we have $b(\bar {t})=1+c\bar {t}$. This makes it possible to study the limiting dynamics in the special gap regime.

3.2.1. Examples of drop dynamics

In figure 4, we present results using $\mathcal {C}=52$ and surface tension $\tau =1\times 10^{-4}$, which makes $k_{max}=k^*=3$. Note that we reverse the orientation of the horizontal axis to reflect the fact that the radius $R$ decreases with time $t$. That is $R$ decreases from 1 as $t$ increases from 0. We use three different initial drop shapes (see insets in figure 4a with $R=1$): $r(\alpha ,0)=1+0.02(\cos (3\alpha )+\cos (5\alpha )+\cos (6\alpha ))$ (blue); $r(\alpha ,0)=1+0.02(\cos (3\alpha )+\sin (7\alpha )+\cos (15\alpha )+\sin (25\alpha ))$ (magenta); and $r(\alpha ,0)=1+0.02(\sin (6\alpha )+\cos (15\alpha )+\sin (25\alpha ))$ (red). Note that in the latter case, mode $k=3$ is not present initially. Here, we take $c=1/2$, and use the time step $\Delta \bar {t}=1\times 10^{-4}$ and $N=8192$ points along the interface. Using the temporal rescaling (Zhao et al. Reference Zhao, Li, Yin, Belmonte, Lowengrub and Li2018), we are able to dramatically reduce the equivalent original time step $\Delta t$ while the rescaled time step $\Delta \bar {t}$ is fixed. As a consequence, the equivalent original time step $\Delta t$ satisfies $\Delta t=({R^9}/{2\tau \mathcal {C}})\Delta \bar {t}$. When $R$ decreases from $1$ to $2\times 10^{-4}$ as in figure 4, $\Delta t$ decreases from $9.6\times 10^{-3}$ to $4.9\times 10^{-36}$. Such wide ranges of time steps are handled seamlessly using our rescaling approach.

Figure 4.

Simulations of drop dynamics in the special gap regime with $b_{\mathcal {C}}(t)=(1-({7}/{2})\tau \mathcal {C} t)^{-{2}/{7}}$ with ${\mathcal {C}=52}$ so that the fastest growing linear mode is $k_{max}=3$. Here $\tau =1\times 10^{-4}$. We reverse the orientation of the horizontal axis to reflect the fact that the radius $R$ decreases with time $t$. (a) The nonlinear shape perturbations for different initial drop shapes are plotted as functions of the equivalent drop radii, as labelled (see text for details). The drop shapes near the vanishing time $T_C=2/(7\tau \mathcal {C})\approx 54.9$ are shown as insets. (b) The dynamics of the drop with modes 3, 5 and 6 near the vanishing time. (c) The drop perimeters as a function of the effective drop radii as they vanish, together with linear fits near the vanishing time (slopes as labelled). The perimeters tend to a finite number as the drops vanish. (d) The fractal dimensions of the drops are shown as a function of the effective drop radii, together with drop morphologies (insets) at various stages of the evolution. These suggest that the limiting shapes have one-dimensional web-like morphologies.

The nonlinear shape factor $\delta /R$ is plotted in figure 4(a) as a function of the effective radius $R(t)$ in the original frame. Here, the nonlinear shape factor is calculated by ${\delta }/{R}=\max _{\alpha } |{|\bar {\boldsymbol {x}}(\alpha ,t)|}/{{\bar {R}}}-1 |$, where $\bar {\boldsymbol {x}}$ is the position vector measured from the centroid of the shape to the interface, ${\bar {R}}=\sqrt {{\bar {A}}/{{\rm \pi} }}$ is the effective radius of the drop in the rescaled frame and $\bar {A}$ is the constant area enclosed by the interface. Unlike the predictions of linear theory, the dynamics of the nonlinear shape factors are non-monotone due to nonlinear interactions among the modes, and the shape factors even decrease at early times (larger $R$). In particular, when mode $3$ is not present initially (red curve) the shape factor decreases until the drop becomes quite small ($R\approx 0.1$). However, as $R$ continues to decrease, eventually mode 3 starts to dominate and the shape factor grows rapidly. The shape perturbations grow throughout the dynamics (see also figure 5a below), which suggests that unlike the case for an expanding bubble (Li et al. Reference Li, Lowengrub and Leo2007, Reference Li, Lowengrub, Fontana and Palffy-Muhoray2009; Zhao et al. Reference Zhao, Belmonte, Li, Li and Lowengrub2016), the evolution does not become self-similar as the drop vanishes.

Figure 5.

Asymptotic behaviours of the drops from figure 4. (a) The nonlinear shape perturbations exhibit a biphasic power-law dependence on the effective drop radii and eventually diverge as $R^{-1}$ for $R\to 0$. (b) The maximum curvature of the drop also diverges as $R^{-1}$ for $R\to 0$. (c) The minimum neck width of the filaments tends to zero as $R^2$ for $R\to 0$.

The insets in figure 4(a) show the initial drop morphologies ($R=1$, top) and the final drop morphologies ($R$ as labelled, bottom) in the rescaled frame (the full dynamics can be found in appendix B figure 11) . Generally, the final morphologies are seen to have a 3-fold symmetric, one-dimensional, web-like network structure although the shapes are somewhat different as the drops vanish. The late time (small $R$) evolution in the original frame is plotted in figure 4(b), which shows the morphologies of the drop with initial condition containing modes $3$, $5$ and $6$ (blue drop in figure 4a). The drop dynamics in the rescaled frame can be found in appendix B (figure 12a). As the drop shrinks, we observe that the tips of the three fingers retract while the long filaments connecting the tips become thinner and tend to a finite length, as seen in the inset. The drop dynamics and morphologies when mode 3 is present initially (blue, magenta) are quite similar while in the third case (red) the shape has a less well-developed network structure because it takes some time for nonlinear interactions to generate mode 3 and then for mode 3 to dominate the shape. This is why the dominance of mode 3 emerges at much smaller $R$ than in the other cases (e.g. when the drop is approximately ${1}/{500}$ of its initial size).

As seen in figure 4(c), the interface perimeter $P\approx P_0+a R$ as $R$ tends to 0, where $a$ is a constant and $P_0$ is a finite number. The slope $a$ depends on the symmetry of the limiting shape ($a$ is a decreasing function of $k_{max}$) and the limiting perimeter $P_0$ depends on the initial shape; see figure 12(b) in appendix B for fits of the interface perimeter for other values of $k_{max}$. To test whether the limiting shapes are truly one-dimensional, we calculate the fractal dimension of the shapes. The fractal dimension $D_0$ can be approximated by a box counting algorithm: cover the pattern with a grid of square boxes of size $\zeta$ and define $N(\zeta )$ to be the total number of boxes of size $\zeta$ to cover the whole pattern (Praud & Swinney Reference Praud and Swinney2005), such that

(3.2)\begin{equation} D_0=\lim_{\zeta\rightarrow 0}-\frac{\log N(\zeta)}{\log \zeta}. \end{equation}

Figure 4(d) shows the fractal dimensions of the shapes as a function of effective drop radius. At early times, the drops remain compact especially for the case with initial modes 6, 15 and 25, which needs a long time for nonlinear interactions to create mode 3. Later there exists a transition as the fractal dimension decreases from approximately 2 down to approximately 1 as the drop vanishes. All together these results strongly suggest that the limiting shape is not a circle and but instead has a web-like network structure. Although the drop morphologies look similar to patterns of random fractals generated using a large unscreened angle threshold (Kaufman, Melroy & Dimino Reference Kaufman, Melroy and Dimino1989), the mechanism is different. In Kaufman et al. (Reference Kaufman, Melroy and Dimino1989), only the tip region is active, but in our case the whole interface is dynamic.

In figure 5 we analyse the properties of the limiting shapes in the original frame as the drop vanishes. In figure 5(a), the nonlinear shape factor is seen to diverge as $R$ tends to $0$. Interestingly, when $R$ is not so small, $\delta /R\sim R^{-2}$, which is consistent with the predictions of linear theory (e.g. see (2.17) with $k_{max}=3$). However, as $R$ decreases, nonlinear interactions increasingly dominate the evolution and the shape factor diverges more slowly, $\delta /R\sim R^{-1}$. This is due to the curvature of the drop tips, which diverges as $\kappa _{*}\sim R^{-1}$ as seen in figure 5(b). The curvature in the scaled frame $\bar \kappa _{*}$, on the other hand, is bounded and tends to a finite limit as $R\to 0$, see figure 12(c) in appendix B. The time dependence of the width $w$ of the neck region, shown within the boxed region in figure 4(b), is plotted in figure 5(c), which suggests the scaling $w\sim R^{2}$. This can be explained as follows. Let $L$ be the length of the neck region. Then, approximating the filament (neck region and drop tip) as a rectangle with a semicircular tip with radius $\kappa _{*}^{-1}$, the total area of the drop . Since $L$ tends to a finite constant as $R\to 0$, this suggests $w\sim R^2$. Further, taking the same approximation of the filament, the perimeter as suggested in figure 4(c).

Another interesting point is the possibility of drop topological changes. In Shelley et al. (Reference Shelley, Tian and Wlodarski1997), the authors have found theoretically that there exists droplet fission under an exponential gap when surface tension is not present. Using the same condition, we have investigated drop fission in a previous paper (Zhao et al. Reference Zhao, Li, Yin, Belmonte, Lowengrub and Li2018). High resolution simulations indicated the droplet with surface tension $\tau =10^{-4}$ investigated in Shelley et al. (Reference Shelley, Tian and Wlodarski1997) does not fission but rather shrinks as a circle. We found the same phenomenon (shrinkage to a circle) using an even smaller surface tension ($2\times 10^{-5}$). It is possible, however, other factors (even smaller surface tension, faster speed, initial set-up and three-dimensional effects) may contribute to fission. Using our special lifting strategy we do not find fission so far.

3.2.2. Mode selection and morphology diagram

Next, we investigate how the limiting shapes are selected by the parameter $\mathcal {C}$ using the special gap dynamics $b_{\mathcal {C}}(t)=(1-({7}/{2})\tau \mathcal {C}t)^{-{2}/{7}}$. We consider the dynamics using different initial shapes given by a perturbed circle with the initial radius $r(\alpha ,0)=1+2.5\times 10^{-3}\sum _{k=30}^{60}\exp (-0.2k)(a_k\cos (k\alpha )+b_k\sin (k\alpha ))$. The coefficients $a_k$ and $b_k$ are randomly selected using a uniform distribution in the interval $[-1,1]$. We generate two such initial shapes and we use the same initial shape for all $\mathcal {C}$, which we vary from 52 to 1000. According to linear theory, the fastest growing modes correspondingly ranges from $k_{max}=3$ to $12$. By considering initial shapes with these high modes, we guarantee that all the initial modes are decreasing (e.g. only modes $|k|<\sqrt {\mathcal {C}/2+1}$ are growing) and that the fastest growing mode is only generated by nonlinear interactions. The result is a morphology diagram given in figure 6, which shows that the dominant mode of the limiting shape (e.g. number of fingers) is an increasing, piecewise constant function of $\mathcal {C}$ where there are sharp transitions from $k$-fold to $(k+1)$-fold dominant limiting shapes. While the morphologies of the limiting shapes are not identical, the dominant mode is solely determined by the constant $\mathcal {C}$. For reference, we also plot the maximum growing mode $k_{max}$ (solid curve). Although linear theory provides a good approximation of the dominant mode of the limiting shape, nonlinear interactions are critical for determining where the transitions from $k$-fold to $(k+1)$-fold dominant limiting shapes occur. Further, the range of $\mathcal {C}$ for which the limiting shape is dominated by a particular mode $k$ is an increasing function of $\mathcal {C}$.

Figure 6.

A morphology diagram that relates the dominant mode of the limiting shapes (number of fingers) to the constant $\mathcal {C}$ using the special gap dynamics $b_{\mathcal {C}}(t)=(1-({7}/{2})\tau \mathcal {C}t)^{-{2}/{7}}$, with $\tau =1\times 10^{-4}$ and $\mathcal {C}$ is varied from 52 to 1000. According to linear theory, the fastest growing mode $k_{max}$ (solid curve) varies from 3 to 12. The dots represent nonlinear numerical simulations using the same initial data (see text for details). The results suggest that the dominant mode of the limiting shapes are selected by $\mathcal {C}$. Linear theory provides a good approximation of the dominant modes, but nonlinear interactions dictate where transitions from $k$-fold to $(k+1)$-fold dominant limiting shapes occur.

In appendix C, we present cases where the initial condition contains modes that grow. In these cases, the dominant mode of the limiting shape can be dependent on the initial condition as well as $\mathcal {C}$ (see figure 13). However, if the initial shape contains $k_{max}$ with magnitude comparable to the other initially growing modes, then the dominant mode of the limiting shape is still given by $k_{max}$, e.g. the limiting shape is still solely selected by $\mathcal {C}$ (see figure 14).

We have also investigated the angles between the thin fingers. In figure 7, we show the angles between pairs of fingers for drop morphologies from the morphology diagram (figure 6) as a function of symmetry mode $k$. Several characteristic drop morphologies are shown as insets, where the smallest and largest angles between pairs of fingers are marked. For each symmetry mode $k$, the average of angles is close to $360/k$ (green curve) although fingers do not seem to intersect at one common point when $k$ is large. Note there are some cases in which some angles deviate significantly from the average (blue morphologies), which reflects the importance of nonlinearities. These happen when the interface conditions are close to transitions from $k$ to $k+1$ mode symmetries. The red morphologies show cases when all angles are close to its average.

Figure 7.

The angles between pairs of fingers for drop morphologies from the morphology diagram (figure 6) are plotted as a function of the symmetry mode $k$. The green curve indicates the angle ${360}/{k}$. Several characteristic drop morphologies at $k=3$, $5$, $8$ and $10$ are shown as insets, where the smallest and largest angles between fingers are marked in the corresponding inset. The blue morphologies show cases when some angles deviate most from $360/k$ and the red morphologies show cases when all angles are close to $360/k$.