2.1 Experimental set-up

A schematic of the experimental set-up is shown in figure 1(a) (see also Pihler-Puzović et al. (Reference Pihler-Puzović, Illien, Heil and Juel2012) and Pihler-Puzović et al. (Reference Pihler-Puzović, Juel, Peng, Lister and Heil2015)). The experiments were performed in an elastic-walled Hele-Shaw cell with a rigid glass bottom and an elastic top boundary separated by a latex spacer with a circular cutout of radius $R_{cell}=180~\text{mm}$ . The top boundary was a thin square elastic sheet, made either of latex (Young’s modulus $E=2.1~\text{MPa}$ and Poisson’s ratio $\unicode[STIX]{x1D708}=0.5$ ) or polypropylene ( $E=3.7~\text{GPa}$ and $\unicode[STIX]{x1D708}=0.44$ ) with side length $2R_{lid}=400~\text{mm}$ and uniform thickness $d$ . The cell was initially filled with viscous liquid (silicone oil of viscosity $\unicode[STIX]{x1D707}=0.962~\text{Pa}~\text{s}$ , density $\unicode[STIX]{x1D70C}=961~\text{kg}~\text{m}^{-3}$ and surface tension $\unicode[STIX]{x1D6FE}=0.021~\text{N}~\text{m}^{-1}$ for the latex cell, or a glycerol–water mixture with $\unicode[STIX]{x1D707}=0.305~\text{Pa}~\text{s}$ , $\unicode[STIX]{x1D70C}=1235~\text{kg}~\text{m}^{-3}$ and $\unicode[STIX]{x1D6FE}=0.065~\text{N}~\text{m}^{-1}$ for the polypropylene cell) to a uniform height $h_{0}$ equal to the thickness of the spacer and radius $R_{fluid}\approx 150~\text{mm}$ .

Prior to the experiment, a small amount of gas was injected through a central hole in the base to form a bubble with radius $R_{init}\approx 5~\text{mm}$ . Nitrogen was then injected at a fixed flow rate $Q$ in the range $100$ – $1300~\text{ml}~\text{min}^{-1}$ , which caused the bubble to expand radially as well as inflate the sheet vertically. The evolution of the expanding, initially approximately circular gas–liquid interface and the development of the fingering instabilities were captured by a camera with a top-down view of the experiment. Retraction of the interface due to dewetting in the case of the glycerol/water solution on polypropylene occurred on time scales of the order of minutes, which were much longer than the experimental time scales ( ${\leqslant}15~\text{s}$ ). Hence, its influence on the fingering instability was negligible.

Table 1. Parameter values of $d$ , $h_{0}$ and $Q$ used in the experiments. The values $Q=Q_{35}$ yield zero growth rate at $\bar{R}=35~\text{mm}$ in the linear stability analysis (see § 4.1).

The values of $d$ , $h_{0}$ and $Q$ used in the experiments are shown in table 1. In the experiments with latex sheets, we either varied $d$ while $h_{0}$ was held constant at $0.56~\text{mm}$ (cells named D33–D97, respectively) or varied $h_{0}$ while $d$ was held constant at $0.56~\text{mm}$ (cells named H46–H79, respectively). The cell with $h_{0}=d=0.56~\text{mm}$ is listed as both D56 and H56 in table 1. In the experiments with polypropylene sheets, we investigated only one cell with $h_{0}=0.56~\text{mm}$ and $d=0.03~\text{mm}$ (named PP).

In order to quantify the variability inherent in the system, we performed at least six experiments for each value of $Q$ in each latex cell and three experiments for each value of $Q$ in the cell with the polypropylene sheet. In total, more than 650 experiments were performed, during which we observed the development of the instability. Results for lower values of $Q$ , for which the interface did not develop fingers, are not shown in table 1 but were reported in Pihler-Puzović et al. (Reference Pihler-Puzović, Juel, Peng, Lister and Heil2015). More details of the experimental protocol can be found in appendix A.

The typical time evolution of the patterns observed in the latex and polypropylene cells is shown in figure 1(b,c), using cells D56/H56 and PP with a volume flux $Q=500~\text{ml}~\text{min}^{-1}$ . At the start of the experiment, the perturbation amplitude is not uniform across modes and it is not a priori constant between experiments. Therefore, each pattern includes a wide range of finger widths – even at the very early stages of the system’s evolution. This is particularly visible in the innermost interface shape in figure 1(c), indicating that the system is very sensitive to initial perturbations. Tip splitting, most visible in the top right corner of the pattern in figure 1(c), occurs as the mean radius (circumference) of the pattern increases.

2.2 Image processing and extracted quantities

All images were processed in MATLAB 2013a by first converting each true colour picture to a binary image using an input value for the intensity threshold, resulting in a white image with black pixels representing the interface. Occasionally this procedure resulted in discontinuities in the interface of up to two pixels, particularly for the thicker (less translucent) latex sheets. Such discontinuities were filled in using a morphological closing algorithm (bwmorph in MATLAB). An active contour algorithm was then applied to extract the interface as a sequence of points $\boldsymbol{x}_{i}$ , with $1\leqslant i\leqslant M$ , going around the gas bubble. (The total number of points, $M$ , increased as the interface expanded.) In regions where the interface was thicker than a pixel, the active contour algorithm detected a meandering path, which we smoothed using a Savitzky–Golay filter with a 31-point window and a third-order polynomial. A typical result of the image processing is shown in figure 2(a), where the extracted interface is superimposed onto a top view of the instantaneous interfacial pattern.

Figure 2. Results from an experiment in the D56/H56 latex cell (see table 1) with $Q=500~\text{ml}~\text{min}^{-1}$ . (a) Interfacial pattern in top view with results from the image processing overlaid. The solid curve shows the extracted position of the gas–liquid interface. The dashed and dot-dashed curves connect the finger tips and bases, respectively. The dotted curve shows the circle of radius $\bar{R}$ . (b) Enlargement of the square region in (a), showing points with locally maximum (crosses) and minimum (circles) radii. The line segments illustrate the criterion for determining if a minimum is considered to be the base of a finger, see text for more details. (c) Data extracted from the experiment, showing half the finger length, $\hat{R}$ , with error bars indicating the standard deviation $\pm \unicode[STIX]{x1D6FF}\hat{R}$ in the pattern, as a function of the mean radius $\bar{R}$ .

Once the interface of an instantaneous pattern had been identified, the number $N(t)$ of fingers was determined. This required a criterion to distinguish the local minima of the interfacial radius associated with the base of a fully formed finger from those arising from a slight indentation at the tip of a wider finger that was only just beginning to split, and hence could be ignored (see e.g. figures 1–3). We used the criterion that a minimum $\boldsymbol{x}_{j}$ is counted if it lay closer to the centre than the line segment connecting the points $\boldsymbol{x}_{j+\unicode[STIX]{x1D706}M}$ and $\boldsymbol{x}_{j-\unicode[STIX]{x1D706}M}$ , which are positioned $\unicode[STIX]{x1D706}M$ points along the circumference to either side of $\boldsymbol{x}_{j}$ , as shown in figure 2(b). We used $\unicode[STIX]{x1D706}=0.015$ in the experiments with latex sheets and $\unicode[STIX]{x1D706}=0.0315$ for the experiments with polypropylene sheets. Hence, the point $\boldsymbol{x}_{i}$ in figure 2(b) was counted as the base of a finger, whereas the point $\boldsymbol{x}_{j}$ was not. We confirmed that this method worked well over the whole range of experiments performed with the same elastic sheet. Alternative criteria based on Fourier decomposition of the patterns or local curvatures of the interface were found to be less robust.

The local minima of the interfacial radius that were deemed to indicate the presence of a finger were then connected with a spline to yield an inner envelope of the pattern (the dot-dashed curve in figure 2 a). A corresponding outer envelope was constructed by fitting a spline through every local maximum; here we included multiple maxima on the same finger.

From the outer and inner envelopes of the pattern, the corresponding radial coordinates, $R_{outer}^{i}$ and $R_{inner}^{i}$ (for $1\leqslant i\leqslant M$ ) were found at the $M$ sampling points, and the average radius of the interface was calculated as

(2.1) $$\begin{eqnarray}\bar{R}(t)=\frac{1}{M}\mathop{\sum }_{i=1}^{M}\frac{R_{outer}^{i}(t)+R_{inner}^{i}(t)}{2}.\end{eqnarray}$$

Throughout this paper, we will plot the evolution of the instability as a function of $\bar{R}(t)$ as a proxy for time $t$ . This is because the interface undergoes rapid initial growth followed by a much slower evolution, so that a nonlinear scale would be required if $t$ was used as the independent variable. Moreover, data could only be collected while the interface was inside the field of view of the camera, so each data series ends at approximately the same value of $\bar{R}$ but at different values of $t$ .

The amplitude of the perturbations to the interface is quantified using the average and standard deviation,

(2.2a,b ) $$\begin{eqnarray}\hat{R}(t)=\frac{1}{M}\mathop{\sum }_{i=1}^{M}\frac{R_{outer}^{i}(t)-R_{inner}^{i}(t)}{2},\quad \unicode[STIX]{x1D6FF}\hat{R}(t)=\sqrt{\frac{1}{M}\mathop{\sum }_{i=1}^{M}\left(\frac{R_{outer}^{i}(t)-R_{inner}^{i}(t)}{2}\right)^{2}-\hat{R}(t)^{2}}.\end{eqnarray}$$

We will refer to $\hat{R}$ as the ‘perturbation amplitude’ or the ‘finger length’ (although the average distance from the tip of a finger to its base is $2\hat{R}$ ). The standard deviation $\unicode[STIX]{x1D6FF}\hat{R}$ is often a significant fraction of $\hat{R}$ (see e.g. figures 1–3 and 11).

Finally, the time-series data for $\bar{R}$ , $\hat{R}$ and $\unicode[STIX]{x1D6FF}\hat{R}$ were all smoothed using the same Savitzky–Golay filter used to smooth the raw interface. A typical result is shown in figure 2(c). The error bars show one standard deviation $\unicode[STIX]{x1D6FF}\hat{R}$ above and below the average value $\hat{R}$ . Similarly, experimental results compiled in later figures are shown as mean points with error bars indicating one standard deviation to either side.

As shown in figure 2(c), the perturbation amplitude $\hat{R}$ initially increases rapidly before reaching a maximum and subsequently decreases upon further expansion of the interface. This behaviour is observed in all experiments. We define $\hat{R}_{peak}$ to be the maximal value of the perturbation amplitude $\hat{R}$ reached in each experiment (i.e. its value before the instability starts decaying) and $\bar{R}_{peak}$ and $N_{peak}$ to be the values of the average radius $\bar{R}$ and the number $N$ of fingers at this point.

The variability in finger length observed within a single pattern and shown with error bars in figure 2(c) is comparable to the variability between different experimental runs performed under the same control parameters. This can be seen in figure 3(a), which compares results for $\hat{R}$ in three identical experiments, both in the D56/H56 latex cell and in the polypropylene cell. Both $\hat{R}_{peak}$ and $\bar{R}=\bar{R}_{peak}$ vary significantly between runs as shown in figure 3(c). This is in contrast with the time evolution of the average radius $\bar{R}(t)$ between different experimental runs which shows far less variability (see Pihler-Puzović et al. (Reference Pihler-Puzović, Juel, Peng, Lister and Heil2015)).

Figure 3. Comparison of three experimental runs (see table 1) in the D56/H56 latex cell (left) and the polypropylene cell (right) at $Q=500~\text{ml}~\text{min}^{-1}$ . (a) The finger length $\hat{R}$ as a function of the average radius $\bar{R}$ . The circles indicate the maximum finger length $\hat{R}_{peak}$ during the pattern evolution. (b) Evolution of the growth rate $\unicode[STIX]{x1D70E}$ obtained from the smoothed data. The insets are enlargements of the region near $\unicode[STIX]{x1D70E}=0$ . (c) The corresponding top views of the interface for the points 1–6 and the number of fingers as determined by the criterion described in § 2.2.

The fixed field of view of the camera used to record the interface in experiments with latex sheets posed a problem for determining $\hat{R}_{peak}$ because any maximum in $\hat{R}$ appearing after $\bar{R}\approx 70~\text{mm}$ could not be detected. In some cases, such as run 1 in figure 3(a), $\hat{R}$ reached a clear maximum and started decaying well before $\bar{R}=70~\text{mm}$ . In other cases, such as run 3, the values of $\hat{R}$ appear to plateau around $\hat{R}_{peak}$ , but a well defined global maximum was not seen implying that decay occurred for values of $\bar{R}>70$  mm. In these cases, the plateau value provided a good estimate of $\hat{R}_{peak}$ , but the value of $\bar{R}$ at which the maximum occurred could not be detected. In the experiments with polypropylene sheets, $\hat{R}$ always reached a maximum value within the field of view.

We characterize the instantaneous growth rate of the instability as

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(t)=\frac{1}{\hat{R}}\frac{\text{d}\hat{R}}{\text{d}t}.\end{eqnarray}$$

The plots of $\unicode[STIX]{x1D70E}$ in figure 3(b) were obtained by differentiation of the data in figure 3(a) and indicate the typical variability between different experimental runs performed under the same control parameters. The results for $\unicode[STIX]{x1D70E}$ form the basis of our comparison with the linear stability analysis.

3.1 The model

Following Pihler-Puzović et al. (Reference Pihler-Puzović, Périllat, Russell, Juel and Heil2013), Pihler-Puzović, Juel & Heil (Reference Pihler-Puzović, Juel and Heil2014), Peng et al. (Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015) and Pihler-Puzović et al. (Reference Pihler-Puzović, Juel, Peng, Lister and Heil2015), we employ horizontal polar coordinates ( $r,\unicode[STIX]{x1D703}$ ) centred at the injection point and model the deflection of the elastic sheet using the Föppl–von-Kármán equations

(3.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{4}\unicode[STIX]{x1D719}+\frac{Ed}{2}\{h,h\}=0,\quad p=\frac{Ed^{3}}{12(1-\unicode[STIX]{x1D708}^{2})}\unicode[STIX]{x1D6FB}^{4}h-\{\unicode[STIX]{x1D719},h\},\end{eqnarray}$$

where $h(r,\unicode[STIX]{x1D703},t)$ is the cell depth, $\unicode[STIX]{x1D719}(r,\unicode[STIX]{x1D703},t)$ is an Airy stress function (related to the components of the stress tensor $\unicode[STIX]{x1D70E}_{ij}$ via $\unicode[STIX]{x1D70E}_{rr}=(1/r)(\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}r)+(1/r^{2})(\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{2})$ , $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D703}\unicode[STIX]{x1D703}}=(\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}r^{2})$ and $\unicode[STIX]{x1D70E}_{r\unicode[STIX]{x1D703}}=-(\unicode[STIX]{x2202}/\unicode[STIX]{x2202}r)((1/r)(\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}))$ ), $p$ is the net upward pressure acting on the sheet and

(3.2) $$\begin{eqnarray}\{f,g\}=\frac{1}{r^{2}}\left(\frac{\unicode[STIX]{x2202}^{2}f}{\unicode[STIX]{x2202}r^{2}}\frac{\unicode[STIX]{x2202}^{2}g}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{2}}+\frac{\unicode[STIX]{x2202}^{2}g}{\unicode[STIX]{x2202}r^{2}}\frac{\unicode[STIX]{x2202}^{2}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{2}}\right)+\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}r}\frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}r}\right)-2\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}g}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)\end{eqnarray}$$

is the Monge–Ampère bracket.

We describe the motion of the viscous fluid occupying the narrow gap underneath the deformed elastic sheet using the lubrication approximation, which yields the equations for the vertically averaged velocity $\boldsymbol{u}(r,\unicode[STIX]{x1D703},t)$ and the fluid pressure as

(3.3a,b ) $$\begin{eqnarray}\boldsymbol{u}=-\frac{h^{2}}{12\unicode[STIX]{x1D707}}\unicode[STIX]{x1D735}p,\quad \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(h\boldsymbol{u})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\left(\frac{h^{3}}{12\unicode[STIX]{x1D707}}\unicode[STIX]{x1D735}p\right)\quad \text{in }r>R(\unicode[STIX]{x1D703},t),\end{eqnarray}$$

where $R(\unicode[STIX]{x1D703},t)$ is the bubble radius.

In the region $r<R(\unicode[STIX]{x1D703},t)$ , the gas bubble has a uniform pressure $p_{g}(t)$ . We note that since the liquid wets the walls of the cell, thin films of liquid are left behind the advancing bubble tip. We assume that the films have no dynamical effect in this region, so that the pressure acting on the elastic sheet is given by

(3.4) $$\begin{eqnarray}p=p_{g}(t)\quad \text{in }r<R(\unicode[STIX]{x1D703},t).\end{eqnarray}$$

At $r=R(\unicode[STIX]{x1D703},t)$ , we impose the kinematic and dynamic conditions

(3.5a,b ) $$\begin{eqnarray}(1-f_{1})U_{n}=\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n},\quad [p]_{-}^{+}=-\unicode[STIX]{x1D6FE}\left(\frac{\unicode[STIX]{x03C0}}{4}\left[\frac{1}{R}-\frac{1}{R^{2}}\frac{\unicode[STIX]{x2202}^{2}R}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{2}}\right]+\frac{2}{h}f_{2}\right),\end{eqnarray}$$

where $U_{n}$ is the speed of the gas–liquid interface in the direction of the unit normal, $\boldsymbol{n}$ , to that interface. The functions $f_{1}$ and $f_{2}$ model the effect of the films that are deposited on the upper and lower walls. Comparisons against full free-surface Navier–Stokes simulations by Peng et al. (Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015), Pihler-Puzović et al. (Reference Pihler-Puzović, Juel, Peng, Lister and Heil2015) showed that, for the capillary numbers encountered in our experiments, the presence of these films has a significant effect on the axisymmetric base flow. The choices

(3.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle f_{1}(Ca)=\frac{Ca^{2/3}}{0.76+2.16Ca^{2/3}},\\ \displaystyle f_{2}(Ca)=1+\frac{Ca^{2/3}}{0.26+1.48Ca^{2/3}}+1.59Ca,\quad \text{for }Ca=\frac{\unicode[STIX]{x1D707}U_{n}}{\unicode[STIX]{x1D6FE}},\end{array}\right\}\end{eqnarray}$$

which were obtained by a fit to Reinelt & Saffman’s (Reference Reinelt and Saffman1985) computational data for the width of the finger penetrating a viscous fluid between parallel plates and the pressure drop across the tip of that finger as a function of the capillary number, produced results that were in excellent agreement with solutions of the full Navier–Stokes equations. We neglect the small corrections to the expressions (3.6) due to the walls of the cell being non-parallel, and refer to papers by Halpern & Jensen (Reference Halpern and Jensen2002) and Jensen et al. (Reference Jensen, Horsburgh, Halpern and Gaver2002) which study the asymptotic structure of the flow in similar geometries with non-parallel walls. The effect of the films deposited on the upper and lower walls can be neglected by taking $f_{1}=0$ and $f_{2}=1$ . Alternative models which also account for the effect of the films in rigid cells can be found in Martyushev & Birzina (Reference Martyushev and Birzina2008), Anjos & Miranda (Reference Anjos and Miranda2013) and Jackson et al. (Reference Jackson, Stevens, Giddings and Power2015).

In the experiments, the sheet rests freely on the viscous layer and there is no significant deformation far ahead of the gas–liquid interface. We mimic this behaviour in the outer boundary conditions of the elastic sheet, but use simplified conditions for the fluid, by imposing

(3.7a,b ) $$\begin{eqnarray}h=h_{0},\quad \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}r}=\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}r}+\frac{1}{r^{2}}\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}^{2}}=-\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\right)=\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}r}=0\quad \text{at }r=R_{out},\end{eqnarray}$$

with $R_{out}=200~\text{mm}$ . We checked that changing $R_{out}$ from $150~\text{mm}$ (edge of the fluid layer in experiments) to $300~\text{mm}$ (corresponding to the diagonal of the elastic sheet) changed the data reported in § 4 by less than $2\,\%$ .

Initially, the cell is unperturbed and, as in the experiments, contains a small bubble of radius $R_{init}=5~\text{mm}$ , so we imposed the initial conditions:

(3.8a,b ) $$\begin{eqnarray}h(r,\unicode[STIX]{x1D703},t=0)=h_{0},\quad R(r,\unicode[STIX]{x1D703},t=0)=R_{init}.\end{eqnarray}$$

Finally, mass conservation requires that

(3.9) $$\begin{eqnarray}\int _{0}^{R_{out}}\int _{0}^{2\unicode[STIX]{x03C0}}(h-h_{0})r\,\text{d}\unicode[STIX]{x1D703}\,\text{d}r=Qt.\end{eqnarray}$$

3.2 Linear stability analysis

We investigate the early states of the fingering instability using linear stability analysis. The base state is axisymmetric, so we employ a Fourier decomposition of the linear perturbations in the azimuthal direction and consider each mode separately. Each quantity $f$ is thus decomposed into a sum

(3.10) $$\begin{eqnarray}f(r,\unicode[STIX]{x1D703},t)=\bar{f}(r,t)+\unicode[STIX]{x1D716}\,\text{Re}[\hat{f}(r,t)\text{e}^{\text{i}n\unicode[STIX]{x1D703}}]+O(\unicode[STIX]{x1D716}^{2})\end{eqnarray}$$

of an axisymmetric but time-dependent base solution $\bar{f}$ and a perturbation with small amplitude $\unicode[STIX]{x1D716}\ll 1$ and azimuthal wavenumber $n$ . The resulting equations are given in appendix B. At $O(1)$ we recover the governing equations for the axisymmetric base state (B 1), which were analysed by Peng et al. (Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015), and at $O(\unicode[STIX]{x1D716})$ , we obtain the governing equations for the perturbations (B 2), which we shall analyse here. Henceforth, we use primes and dots to denote differentiation with respect to $r$ and $t$ , respectively.

The evolution of the perturbation was followed from the initial conditions

(3.11a,b ) $$\begin{eqnarray}\hat{R}=1,\quad {\hat{h}}=0,\end{eqnarray}$$

which corresponds to an initial perturbation to the interface position only.

We solve the governing equations, (B 1) for the base flow and (B 2) for each linear perturbation with wavenumber $1\leqslant n\leqslant 50$ , simultaneously using a fully implicit backward-Euler finite-difference scheme with adaptive spatial grid and adaptive time stepping (for details see appendix B). We define $\unicode[STIX]{x1D70E}$ for each azimuthal wavenumber $n$ as in (2.3). We also compare the instantaneous growth rates to the growth rates obtained using the ‘frozen-time’ approximation (the eigenvalue analysis of (B 2), conducted at each time step using the numerical solutions to (B 1)).

In figure 4(a) we show the time evolution of a representative base-state height profile $\bar{h}$ and interfacial position. The injected gas deflects the sheet upwards and causes the gas–liquid interface to spread outwards. Ahead of the interface, the displaced liquid accumulates in a wedge whose size increases with time. In figure 4(a) below the axis $\bar{h}=0$ the symbols show the interface observed in experiments, with each point corresponding to one of the 12 repeated experimental runs and the horizontal bars indicating the range $\bar{R}\pm \hat{R}$ observed in them. The model can be seen to capture the mean position of the interface in the experiments satisfactorily.

Figure 4. Results of the linear stability analysis for the D56/H56 latex cell with $Q=400~\text{ml}~\text{min}^{-1}$ . (a) Instantaneous base-state profiles $\bar{h}$ (solid curves), and the corresponding approximate wedge profiles (dashed curves) at $t=0.36$ , $2.38$ and $7.38~\text{s}$ , respectively. The solid vertical lines indicate the position of the interface, and points with the horizontal error bars below indicate the position $\bar{R}\pm \hat{R}$ observed in each of the 12 experiments. (b,c) The perturbations to the height and pressure profiles for a wavenumber $n=25$ , normalized by the maximum absolute value. In (b), the dotted curves show the results from the ‘frozen-time’ approximation. In (c), $\hat{p}$ is only plotted in the range $r>\bar{R}$ , and the dashed curves show the profiles from the rigid-wedge approximation. (d) The evolution of amplitudes $\hat{R}$ for a range of different wavenumbers, comparing the results from the full analysis (solid curves) to those obtained when using the rigid-wedge approximation (dashed curves). The maximum for each wavenumber is marked by a vertical bar.

To analyse the characteristic features of the perturbations, we focus on the mode that is predicted to grow to the largest amplitude $\hat{R}$ in figure 4, namely $n=25$ . The height perturbations ${\hat{h}}$ are shown in figure 4(b) using solid lines, for the same three times as in figure 4(a). They have multiple undulations in the vicinity of the wedge region, but are otherwise very small. A more detailed inspection of the data shows that the perturbations to the fluid pressure and to the position of the air–liquid interface are much larger than the perturbations to the deflection of the sheet. Specifically, $({\hat{h}}/\bar{h})/(\hat{R}/\bar{R})$ never rises above a value of $O(10^{-6})$ in the time interval shown in this figure. In contrast, the corresponding pressure perturbations $\hat{p}$ , illustrated in figure 4(c) using solid lines, are much more significant, with $(\hat{p}/\bar{p})/(\hat{R}/\bar{R})$ being in the range $O(10^{2}{-}10^{3})$ over the same time interval. We will exploit this observation in § 3.3 below.

In figure 4(b), we also show the results obtained from the ‘frozen-time’ approximation using dotted curves. We find that within the experimental field of view the maximum instantaneous growth rate $\unicode[STIX]{x1D70E}$ obtained using this approximation is no more than 1 % larger than that obtained from the direct calculations for all wavenumbers between 5 and 50.

In figure 4(d), we show the logarithm of the perturbation amplitude $\hat{R}$ as a function of the mean radius $\bar{R}$ for wavenumbers between 10 and 40 in increments of 5. A large number of modes can be seen to have comparable (positive) growth rates and the linear analysis predicts all of these modes to grow to a maximum amplitude before decaying. Which one of these modes will actually become dominant in a specific experiment depends on the relative initial amplitudes of the perturbations and on nonlinear effects that start to play a role when the amplitude of the perturbation becomes sufficiently large, an issue we will discuss in more detail below. It is, however, interesting to observe that in figure 4(d) a single mode ( $n=25$ ) maintains the largest instantaneous growth rate throughout the system’s evolution and therefore also grows to the largest amplitude. This behaviour was found to be typical for all the cases considered in this study (see table 1).

3.3 The rigid-wedge approximation and the physical mechanism of fingering suppression

The results in figure 4 show that the most rapidly growing modes have relatively large wavenumber, $n\geqslant 10$ , and hence comparatively short azimuthal (lateral) length scales. As discussed in § 3.2, we also find that the pressure perturbations $\hat{p}$ only generate a small perturbation ${\hat{h}}$ to the sheet (i.e. gap height), suggesting that ${\hat{h}}$ can be neglected. Moreover, since the perturbations were found to be mostly confined to the wedge region (see figure 4 b,c), we approximate the perturbation flow ahead of the interface as a flow inside a wedge-shaped domain with straight converging walls. These are precisely the ‘rigid-lid’ and ‘wedge’ approximations, first introduced by Peng & Lister (Reference Peng and Lister2018) for the idealized rectangular elastic-walled Hele-Shaw cell. Using these approximations, we will now derive an analytical expression for the instantaneous growth rate $\unicode[STIX]{x1D70E}$ of the different modes in terms of the (numerically calculated) base-state quantities at the interface. This will then allow us to identify the physical mechanisms by which the presence of the elastic sheet weakens (or even suppresses) the fingering instability.

In order to simplify the notation, we denote the base-state height, slope and (compressive) longitudinal velocity gradient at the interface (subscript ‘ $i$ ’) by

(3.12a-c ) $$\begin{eqnarray}h_{i}=\bar{h},\quad \unicode[STIX]{x1D6FC}_{i}=-\bar{h}^{\prime },\quad c_{i}=-\bar{u}^{\prime }\quad \text{at }r=\bar{R}^{+}.\end{eqnarray}$$

By assuming ${\hat{h}}=0$ (the ‘rigid-lid’ approximation), substituting $\dot{\hat{R}}=\unicode[STIX]{x1D70E}\hat{R}$ and using (B 1d ), the interfacial conditions (B 2e ), (B 2f ) become

(3.13a ) $$\begin{eqnarray}\displaystyle & \displaystyle (1-f_{1}-\overline{Ca}f_{1}^{\prime })\unicode[STIX]{x1D70E}=\frac{\hat{u} }{\hat{R}}-c_{i}, & \displaystyle\end{eqnarray}$$

(3.13b ) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\hat{p}}{\hat{R}}=\frac{12\unicode[STIX]{x1D707}\dot{\bar{R}}}{h_{i}^{2}}(1-f_{1})-\unicode[STIX]{x1D6FE}\left[\frac{\unicode[STIX]{x03C0}}{4}\frac{n^{2}-1}{\bar{R}^{2}}+\frac{2f_{2}}{h_{i}^{2}}\unicode[STIX]{x1D6FC}_{i}\right]-\frac{2f_{2}^{\prime }}{h_{i}}\unicode[STIX]{x1D707}\unicode[STIX]{x1D70E}, & \displaystyle\end{eqnarray}$$

$\hat{u}$

$\hat{p}$

$r=\bar{R}^{+}$

The dynamic condition (3.13b ) describes how the pressure perturbation is generated. The first term on the right-hand side describes how the base viscous pressure drop associated with the mean advance of the interface results in a relative increase in pressure at the tip of the fingers. This is the driving force of the instability, as the increased pressure drives a perturbation flow that causes growth of the fingers. The term includes a stabilizing correction due to the wetting films left behind the interface, as these films reduce the amount of fluid that needs to be displaced to achieve a given interfacial speed $\dot{\bar{R}}$ and hence reduce the base viscous pressure drop. There are also two stabilizing surface-tension terms in (3.13b ). The first term is due to the perturbations to the in-plane curvature; this is the main restoring effect in the classical case of viscous fingering in a rigid-walled Hele-Shaw cell. The second term is due to surface tension acting on the vertical curvature of the interface as observed in a radial cross-section of the cell. Due to the converging geometry of the cell, the tip of the finger sits in a narrower gap and hence has a larger vertical curvature. The result is a lowering of the pressure at the tip, which has a restoring effect. This is the ‘taper’ mechanism described by Al-Housseiny et al. (Reference Al-Housseiny, Tsai and Stone2012).

In the kinematic condition (3.13a ), the first term on the right-hand side describes how perturbations to the interface grow or decay due to advection of the interface by the perturbation velocity, which is driven by the pressure perturbation from (3.13b ). The second term describes how spatial variations in the base flow affect the advection of the finger tip. The simulations show that typically $c_{i}=-\bar{u}^{\prime }>0$ , so this is a stabilizing effect and we call it ‘kinematic compression’. (In Juel et al. (Reference Juel, Pihler-Puzović and Heil2018), the term ‘geometry’ is used instead.) The correction factors in the parentheses on the left-hand side account for the effects of the thin wetting films left behind the tip.

In order to solve for $\unicode[STIX]{x1D70E}$ in (3.13), we need to calculate the ratio $\hat{u} /\hat{p}=-h_{i}^{2}\hat{p}^{\prime }/(12\unicode[STIX]{x1D707}\hat{p})$ by solving the perturbation lubrication equation (B 2c ) (while neglecting ${\hat{h}}$ ). We introduce a local variable $x=\bar{R}+h_{i}/\unicode[STIX]{x1D6FC}_{i}-r$ , and approximate the base-state height profile as a wedge $\bar{h}=\unicode[STIX]{x1D6FC}_{i}x$ with slope $\unicode[STIX]{x1D6FC}_{i}$ and height $h_{i}$ at the interface $r=\bar{R}$ . Assuming that the wedge is short compared with the radius, $x\ll \bar{R}$ , allows us to simplify the lubrication equation further by neglecting the terms arising from the radial geometry to obtain

(3.14) $$\begin{eqnarray}0=\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}\left(x^{3}\frac{\unicode[STIX]{x2202}\hat{p}}{\unicode[STIX]{x2202}x}\right)-\frac{n^{2}}{\bar{R}^{2}}x^{3}\hat{p}\quad \text{in }0<x<\frac{h_{i}}{\unicode[STIX]{x1D6FC}_{i}},\quad \hat{p}^{\prime }=0\text{ at }x=0,\end{eqnarray}$$

where the boundary condition is analogous to (B 2i ). The solution is $\hat{p}\propto \text{I}_{1}(nx/\bar{R})/x$ , where $\text{I}_{1}$ is the order- $1$ modified Bessel function of the first kind. Following Peng & Lister (Reference Peng and Lister2018), we rescale the ratio $\hat{u} /\hat{p}$ (evaluated at $r=\bar{R}^{+}$ or equivalently $x=h_{i}/\unicode[STIX]{x1D6FC}_{i}$ ) to define the ‘admittance’

(3.15) $$\begin{eqnarray}Y=\frac{\bar{R}}{n}\frac{12\unicode[STIX]{x1D707}}{h_{i}^{2}}\frac{\hat{u} }{\hat{p}}=\frac{\bar{R}}{n}\frac{1}{\hat{p}}\frac{\unicode[STIX]{x2202}\hat{p}}{\unicode[STIX]{x2202}x}=\frac{I_{1}^{\prime }(s)}{I_{1}(s)}-\frac{1}{s},\quad \text{where }s=\frac{nh_{i}}{\unicode[STIX]{x1D6FC}_{i}\bar{R}}.\end{eqnarray}$$

The scaling is chosen so that $Y=1$ for an infinite cell with rigid and parallel walls, while the restricted flow in a converging wedge gives a smaller value of $\hat{u} /\hat{p}$ and hence $Y<1$ .

Using the admittance, we can eliminate the unknowns $\hat{p}/\hat{R}$ and $\hat{u} /\hat{R}$ from (3.13) to obtain

(3.16) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(t)={\displaystyle \frac{\dot{\hat{R}}}{\hat{R}}}={\displaystyle \frac{n\dot{\bar{R}}}{\bar{R}}}{\displaystyle \frac{\displaystyle Y\left[1-f_{1}-{\displaystyle \frac{1}{12\overline{Ca}}}\left({\displaystyle \frac{\unicode[STIX]{x03C0}}{4}}{\displaystyle \frac{h_{i}^{2}}{\bar{R}^{2}}}(n^{2}-1)+2f_{2}\unicode[STIX]{x1D6FC}_{i}\right)\right]-{\displaystyle \frac{\bar{R}c_{i}}{n\dot{\bar{R}}}}}{\displaystyle 1-f_{1}-\overline{Ca}f_{1}^{\prime }+{\displaystyle \frac{2n\,Yf_{2}^{\prime }h_{i}}{12\bar{R}}}}}.\end{eqnarray}$$

The evolution of the perturbation amplitude is then given by

(3.17) $$\begin{eqnarray}\hat{R}(t)=\exp \left[\int _{0}^{t}\unicode[STIX]{x1D70E}(t^{\prime })\,\,\text{d}t^{\prime }\right].\end{eqnarray}$$

The denominator in (3.16) is always positive, so the stability of the perturbation is determined by the sign of the numerator. The perturbation decays, i.e.  $\unicode[STIX]{x1D70E}<0$ , provided that

(3.18) $$\begin{eqnarray}f_{1}(\overline{Ca})+\frac{2f_{2}(\overline{Ca})}{12\overline{Ca}}\unicode[STIX]{x1D6FC}_{i}+\frac{\unicode[STIX]{x03C0}/4}{12\overline{Ca}}\frac{h_{i}^{2}}{\bar{R}^{2}}(n^{2}-1)+\frac{\bar{R}c_{i}}{n\dot{\bar{R}}Y}>1.\end{eqnarray}$$

These four terms represent the stabilizing effects of films deposited behind the interface, taper, horizontal surface tension and kinematic compression, respectively, and they give a quantitative measure of the stabilization due to each physical mechanism.

Using the base-state equations (B 1b ) and (B 1d ), we can further decompose the longitudinal compression into

(3.19) $$\begin{eqnarray}c_{i}=-\bar{u}^{\prime }(\bar{R}^{+})=\frac{\dot{\bar{R}}f_{1}\unicode[STIX]{x1D6FC}_{i}}{h_{i}}+\frac{(1-f_{1})\dot{\bar{R}}}{\bar{R}}+\frac{\text{d}h_{i}/\text{d}t}{h_{i}}.\end{eqnarray}$$

These three terms represent three mechanisms that stretch the base flow near the interface in the azimuthal and vertical directions, and cause the base flow to compress longitudinally through mass conservation. In order to understand these three terms, we consider the local frame of reference moving with the interface at velocity $\dot{\bar{R}}$ so that the interface is stationary.

The first term in (3.19) is the same as in the unidirectional steady peeling case (Peng & Lister Reference Peng and Lister2018), and describes how the flow towards the interface, with velocity $\bar{u}-\dot{\bar{R}}=-f_{1}\dot{\bar{R}}$ , experiences a vertically widening cell at the relative rate $\unicode[STIX]{x1D6FC}_{i}/h_{i}$ per unit length in the longitudinal direction. The second term is due to the radial geometry of the cell, and is the relative expansion rate $\dot{\bar{R}}/\bar{R}$ of the circumference multiplied by the correction factor $(1-f_{1})$ to account for the deposited films. Finally, the third term is the relative rate of change of the interfacial gap height (in the frame of reference moving with the interface). Typically, the height is (slowly) increasing with time, so this represents a stretching in the vertical direction. All three terms in (3.19) are positive and therefore contribute to the suppression of the instability.

For $n\gg 1$ , the expression (3.16) for $\unicode[STIX]{x1D70E}$ reduces to the expression found by Peng & Lister (Reference Peng and Lister2018) if we identify the lateral wavenumber as $k=n/\bar{R}$ . The only difference is in the compression term involving the velocity gradient $c_{i}=-\bar{u}^{\prime }$ . The kinematic compression was also found to be stabilizing in the cell studied by Peng & Lister (Reference Peng and Lister2018), but in their two-dimensional time-independent system the peeling process did not involve the additional mechanisms provided by the second and third terms on right-hand side of (3.19). Hence, we conclude that the present radial time-evolving system is more stable than the one studied by Peng & Lister (Reference Peng and Lister2018).

Our numerical results (e.g. figure 4) show that, in the regime under consideration, the perturbation pressure profiles obtained using the full calculation and the rigid-wedge approximation agree well overall, particularly near the interface (see figure 4 c). Consequently there is excellent agreement between the models regarding the evolution of the perturbation amplitudes $\hat{R}$ and growth rate $\unicode[STIX]{x1D70E}$ (figure 4 d). Thus, for simplicity, we focus henceforth on results from the rigid-wedge approximation only.