2.1 Governing equations

Our physical model is similar to the one investigated by Peng et al. (Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015), who discuss several of the assumptions made below in more detail. The horizontal length scales of the system are assumed to be much longer than the vertical length scales, so that vertically integrated quantities can be employed which are functions of the horizontal position $\boldsymbol{x}=(x,y)$ and time $t$ only. We use $\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{y})$ to denote the horizontal gradient operator, a prime to denote the longitudinal derivative $\unicode[STIX]{x2202}_{x}$ , in particular, and a dot to denote a time derivative.

The elastic sheet is modelled either as a membrane under isotropic tension or as an elastic beam, which we refer to as the ‘peeling-by-pulling’ and ‘peeling-by-bending’ cases, respectively, following Lister et al. (Reference Lister, Peng and Neufeld2013). (These are two limiting cases of the Föppl–von Kármán equations which describe a linearly elastic sheet with small strain and small slope, as discussed in more detail by e.g. Lister et al. (Reference Lister, Peng and Neufeld2013) and Peng et al. (Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015).) Thus, the net upward pressure $p(x,y,t)$ on the sheet is related to its deflection $h-h_{0}$ by either

(2.1a,b ) $$\begin{eqnarray}p=-T\unicode[STIX]{x1D6FB}^{2}h\quad \text{(pulling)}\quad \text{or}\quad p=B\unicode[STIX]{x1D6FB}^{4}h\quad \text{(bending)},\end{eqnarray}$$

where $T$ is the isotropic tension in the sheet or $B$ is its bending stiffness (which is related to the Young’s modulus $E$ , thickness $d$ and Poisson’s ratio $\unicode[STIX]{x1D708}$ by $B=Ed^{3}/12(1-\unicode[STIX]{x1D708}^{2})$ ). (We investigate the non-isotropic case where the spanwise component of tension is reduced towards zero in appendix C.)

We neglect the effects of gravity on both the sheet and the fluid, and measure all pressures relative to atmospheric pressure. In the liquid region ahead of the meniscus, lubrication theory applies and hence the pressure is vertically uniform to leading order and equal to $p$ . Gradients in this pressure drive a flow with a parabolic Poiseuille profile, and hence the depth-averaged velocity and the time evolution of the height profile are given by the lubrication equations (Reynolds Reference Reynolds1886)

(2.2a,b ) $$\begin{eqnarray}\boldsymbol{u}=(u,v)=-\frac{h^{2}}{12\unicode[STIX]{x1D707}}\unicode[STIX]{x1D735}p,\quad {\dot{h}}=-\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 }x>R(y,t).\end{eqnarray}$$

(It can be shown that horizontal motion of the sheet due to the viscous shear stress is negligible compared to the motion of the fluid.)

In the gas region behind the meniscus, the pressure in the gas (relative to atmospheric) is assumed to be zero. Physically, this could be due to the cell being open to the atmosphere at the left end. Alternatively, if the region under consideration forms part of a larger-scale system such as an elastic-walled Hele-Shaw cell with gas injection (Lister et al. Reference Lister, Peng and Neufeld2013), then the separation of length scales implies that the pressure in the gas is negligible at leading order. The top and bottom surfaces are coated by thin films of liquid, but these films have negligible vertical variation in pressure and negligible pressure jump across the very slightly curved gas–liquid interface. Hence, to a good approximation,

(2.3) $$\begin{eqnarray}p=0\quad \text{in }x<R(y,t).\end{eqnarray}$$

The liquid and gas regions are connected by a region near the meniscus in which the horizontal length scales of the flow are comparable with the vertical length scales. We assume that the meniscus is nearly perpendicular to the longitudinal $x$ -direction, which is the main direction of flow, and that variations in the $y$ -direction occur on a length scale much larger than the vertical scale. Hence, the flow in this region is locally similar to that near the meniscus of a bubble advancing in a Hele-Shaw cell with parallel and rigid walls.

Following Peng et al. (Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015), we approximate the meniscus region as a sharp interface in the lubrication model by parametrizing the flow in this region using two functions $f_{1}$ and $f_{2}$ , which respectively describe the amount of fluid deposited into the thin films behind the meniscus (as seen in figure 1) and the apparent pressure drop across the meniscus region due to capillary and viscous effects. The dependence of these functions on the local instantaneous capillary number $\widetilde{Ca}(y,t)=\unicode[STIX]{x1D707}{\dot{R}}/\unicode[STIX]{x1D6FE}$ was taken by Peng et al. (Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015) to be

(2.4a,b ) $$\begin{eqnarray}f_{1}(\widetilde{Ca})=\frac{\widetilde{Ca}^{2/3}}{0.76+2.16\widetilde{Ca}^{2/3}},\quad f_{2}(\widetilde{Ca})=1+\frac{\widetilde{Ca}^{2/3}}{0.26+1.48\widetilde{Ca}^{2/3}}+1.59\widetilde{Ca},\end{eqnarray}$$

as an approximation to the asymptotic results by Bretherton (Reference Bretherton1961) and the numerical results by Reinelt & Saffman (Reference Reinelt and Saffman1985). The matching conditions at the meniscus then become

(2.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle u={\dot{R}}(1-f_{1}(\unicode[STIX]{x1D707}{\dot{R}}/\unicode[STIX]{x1D6FE}))\quad \text{at }x=R^{+}, & \displaystyle\end{eqnarray}$$

(2.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle [p]_{-}^{+}=\frac{\unicode[STIX]{x03C0}}{4}\unicode[STIX]{x1D6FE}\unicode[STIX]{x2202}_{y}^{2}R-\frac{2\unicode[STIX]{x1D6FE}f_{2}(\unicode[STIX]{x1D707}{\dot{R}}/\unicode[STIX]{x1D6FE})}{h}\quad \text{at }x=R, & \displaystyle\end{eqnarray}$$

$[p]_{-}^{+}$

$u$

${\dot{R}}$

$f_{1}$

$f_{2}$

Although the pressure, which is the second or fourth derivative of $h$ in (2.1), has a discontinuity at the interface, lower derivatives of $h$ do not. Hence, we have the conditions

(2.6) $$\begin{eqnarray}\left\{\begin{array}{@{}ll@{}}h\text{ and its first derivative are continuous}\quad & \text{(pulling)},\\ h\text{ and its first, second and third derivatives are continuous}\quad & \text{(bending)},\end{array}\right.\quad \text{at }x=R.\end{eqnarray}$$

Far ahead of the peeling region, the sheet and liquid are undisturbed and hence

(2.7) $$\begin{eqnarray}h\rightarrow h_{0}\quad \text{as }x\rightarrow \infty .\end{eqnarray}$$

At the left-hand end, we assume that the elastic sheet is being peeled off either at a fixed slope $\unicode[STIX]{x1D6FC}$ (peeling-by-pulling) or with a fixed curvature $\unicode[STIX]{x1D705}$ (peeling-by-bending),

(2.8a,b ) $$\begin{eqnarray}h^{\prime }\rightarrow -\unicode[STIX]{x1D6FC}\quad \text{(pulling)},\quad h^{\prime \prime }\rightarrow \unicode[STIX]{x1D705}\quad \text{(bending)}\quad \text{as }x\rightarrow -\infty ,\end{eqnarray}$$

where primes denote derivatives with respect to $x$ . These conditions were shown by Lister et al. (Reference Lister, Peng and Neufeld2013) to be the appropriate matching conditions when the peeling region forms the inner boundary-layer region in a matched asymptotic expansion.

From now on, we will, for simplicity, work in a frame of reference moving at constant speed $U$ to the right. The evolution equation (2.2b ) becomes

(2.9) $$\begin{eqnarray}{\dot{h}}-Uh^{\prime }=-\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 }x>R(y,t)-Ut,\end{eqnarray}$$

while the remaining equations (2.1) and (2.3)–(2.8) remain the same but with the interfacial position $R$ replaced by $R-Ut$ . (However, we keep $u$ to denote the original velocity in the stationary frame, rather than transforming it into the moving frame.)

2.2 Base flow

We first seek a $y$ -independent base solution describing steady peeling at the given speed $U$ . Owing to translational invariance, we can choose the position of the interface to be at $x=0$ in the moving coordinate system (i.e.  $x=Ut$ in the original frame of reference). Using overbars to denote the steady solution, the elastic laws (2.1) become

(2.10a,b ) $$\begin{eqnarray}\bar{p}=-T\bar{h}^{\prime \prime }\quad \text{(pulling)}\quad \text{or}\quad \bar{p}=B\bar{h}^{\prime \prime \prime \prime }\quad \text{(bending)}.\end{eqnarray}$$

Combining (2.3), (2.8) and (2.10) yields analytical solutions for the gas region,

(2.11a,b ) $$\begin{eqnarray}\bar{h}=-\unicode[STIX]{x1D6FC}x+A_{1}\quad \text{(pulling)}\quad \text{or}\quad \bar{h}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D705}x^{2}+A_{2}x+A_{3}\quad \text{(bending)}\quad \text{in }x<0,\end{eqnarray}$$

where $A_{1,2,3}$ are unknown constants. Hence, only the liquid region will need to be solved numerically.

Integration of the evolution equation (2.9) with ${\dot{h}}=0$ , and use of the far-field condition

(2.12) $$\begin{eqnarray}\bar{h}\rightarrow h_{0}\quad \text{as }x\rightarrow \infty\end{eqnarray}$$

from (2.7), yields

(2.13) $$\begin{eqnarray}\bar{u}=-\frac{\bar{h}^{2}}{12\unicode[STIX]{x1D707}}\bar{p}^{\prime }=U\frac{\bar{h}-h_{0}}{\bar{h}}\quad \text{in }x>0.\end{eqnarray}$$

This equation describes how the leftward flux $\bar{h}(U-\bar{u})$ in the travelling-wave frame is everywhere equal to its value $Uh_{0}$ in the far field.

The kinematic and dynamic interfacial conditions (2.5) simplify to

(2.14a,b ) $$\begin{eqnarray}\bar{h}=\frac{h_{0}}{f_{1}},\quad \bar{p}=-\frac{2\unicode[STIX]{x1D6FE}f_{2}}{\bar{h}}\quad \text{at }x=0^{+}.\end{eqnarray}$$

In this equation and henceforth, the functions $f_{1}$ and $f_{2}$ and their derivatives are evaluated at the base value $Ca=\unicode[STIX]{x1D707}U/\unicode[STIX]{x1D6FE}$ unless otherwise stated. The first condition (2.14a ) is readily interpreted as requiring that the flux $Uf_{1}\bar{h}$ of fluid being deposited into thin films behind the interface must be equal to the amount $Uh_{0}$ of fluid arriving from the far field. If no such films are formed (i.e.  $f_{1}=0$ ), then a steady base state would not exist.

Combining the continuity conditions (2.6) with the analytical solutions (2.11) for the gas region yields

(2.15a,b ) $$\begin{eqnarray}\bar{h}^{\prime }=-\unicode[STIX]{x1D6FC}\quad \text{(pulling)},\quad \bar{h}^{\prime \prime }=\unicode[STIX]{x1D705},~\bar{h}^{\prime \prime \prime }=0\quad \text{(bending)}\quad \text{at }x=0^{+}.\end{eqnarray}$$

Typically, one would consider $\unicode[STIX]{x1D6FC}$ or $\unicode[STIX]{x1D705}$ as given and seek to determine the unknown peeling speed $U$ , but for convenience we will instead assume that $U$ is known (either because it is imposed directly or because the base state has been solved for already), so that two equations in (2.15) simply yield $\unicode[STIX]{x1D6FC}$ or $\unicode[STIX]{x1D705}$ in terms of $U$ and the other parameters of the problem, and only the equation $\bar{h}^{\prime \prime \prime }=0$ for peeling-by-bending is a condition to impose.

There are two independent non-dimensional parameters in this problem, which we choose to be the capillary number and the aspect ratio,

(2.16a,b ) $$\begin{eqnarray}Ca=\frac{\unicode[STIX]{x1D707}U}{\unicode[STIX]{x1D6FE}},\quad {\mathcal{A}}=\left\{\begin{array}{@{}ll@{}}(T/12\unicode[STIX]{x1D707}U)^{1/3}\quad & \text{(pulling)},\\ (B/12\unicode[STIX]{x1D707}Uh_{0}^{2})^{1/5}\quad & \text{(bending)}.\end{array}\right.\end{eqnarray}$$

The choice of ${\mathcal{A}}$ is obtained from balancing all terms in the lubrication equation (2.13), assuming that $h\sim h_{0}$ and $x\sim {\mathcal{A}}h_{0}$ . A full non-dimensionalization and the resulting governing equations are given in appendix A, but we retain the use of dimensional quantities in the main text in order to highlight the physical mechanisms more clearly.

Strictly speaking, our governing equations apply to the asymptotic limit ${\mathcal{A}}\gg 1$ , but we will consider all values of ${\mathcal{A}}\geqslant 1$ , with the understanding that for moderate values of ${\mathcal{A}}$ the results are only qualitatively correct. (Equivalently, we later consider slopes $\unicode[STIX]{x1D6FC}_{i}\leqslant 1$ when in reality $\unicode[STIX]{x1D6FC}_{i}\ll 1$ is necessary for asymptotic accuracy, and in the linear stability analysis we strictly need $Ca\ll 1$ for the main wavelengths to be long compared with the vertical length scale, but we consider all values $Ca\leqslant 1$ .)

The base-state governing equations (2.10)–(2.15) are straightforward to solve numerically. In fact, both cases have been solved before (McEwan & Taylor Reference McEwan and Taylor1966; Balmforth et al. Reference Balmforth, Craster and Hewitt2014), albeit with different choices of $f_{1}$ and $f_{2}$ , which do not change the qualitative behaviour of the system. A few cases are shown in figure 2. The base-state solutions have height $h_{0}/f_{1}$ at the interface, as imposed by (2.14a ) due to conservation of flux. Ahead of the interface, the height decreases towards $h_{0}$ , resulting in a wedge-like shape.

Figure 2.

Base-state solutions satisfying (2.10)–(2.15) for various values of the capillary number $Ca$ and the aspect ratio ${\mathcal{A}}$ defined in (2.16).

Near the tip of the wedge, there are peeling-wave oscillations with length scale $O({\mathcal{A}}h_{0})$ due to the interplay of the lubrication pressure and the elastic effects (Lister et al. Reference Lister, Peng and Neufeld2013), but these are not a significant effect for the stability calculations here. For peeling-by-bending, if ${\mathcal{A}}$ or $Ca$ is too small then no steady-state base peeling solution exists (Peng et al. Reference Peng, Pihler-Puzović, Juel, Heil and Lister2015), as can be seen, for example, in figure 5. For each given value of $Ca$ , as ${\mathcal{A}}\rightarrow \infty$ the solution profile tends to a limiting shape which corresponds to replacing $1/{\mathcal{A}}$ by $0$ in the pressure-jump condition (A 4), since that is the only occurrence of ${\mathcal{A}}$ in the non-dimensional base-state equations (A 2)–(A 4).

As we proceed with the linear stability analysis, two quantities from the base-state solution will be of particular importance. These are the height $h_{i}$ and slope $\unicode[STIX]{x1D6FC}_{i}$ of the sheet at the interface:

(2.17a ) $$\begin{eqnarray}h_{i}=\bar{h}(0)=\frac{h_{0}}{f_{1}},\quad \unicode[STIX]{x1D6FC}_{i}=-\bar{h}^{\prime }(0)>0.\end{eqnarray}$$

We will also make use of the fact that the longitudinal velocity gradient $\bar{u}^{\prime }$ can be simplified using (2.13) to give

(2.17b ) $$\begin{eqnarray}\bar{u}^{\prime }=\frac{Uh_{0}\bar{h}^{\prime }}{\bar{h}^{2}}=-\frac{Uh_{0}\unicode[STIX]{x1D6FC}_{i}}{{h_{i}}^{2}}<0\quad \text{at }x=0,\end{eqnarray}$$

which means that the longitudinal flow (in the $x$ -direction) is compressive near the interface, due to the flow being stretched vertically as it approaches the meniscus in a widening gap.

2.3 Linearized governing equations for small perturbations

We now derive the equations governing small perturbations to the base states calculated above. Since the base state is independent of $y$ and $t$ , we can consider separately perturbations with different wavenumbers $k$ in the $y$ -direction and different temporal growth rates $\unicode[STIX]{x1D70E}$ (which we seek to determine as functions of $k$ ). Hence, we decompose each quantity as a sum of the base-state solution and a perturbation with amplitude $\unicode[STIX]{x1D716}\ll 1$ , real wavenumber $k$ and complex growth rate $\unicode[STIX]{x1D70E}$ :

(2.18a ) $$\begin{eqnarray}\displaystyle & \displaystyle h=\bar{h}(x)+\unicode[STIX]{x1D716}\,\text{Re}[{\hat{h}}(x)\text{e}^{\text{i}ky+\unicode[STIX]{x1D70E}t}]+O(\unicode[STIX]{x1D716}^{2}),\quad p=\bar{p}(x)+\unicode[STIX]{x1D716}\,\text{Re}[\hat{p}(x)\text{e}^{\text{i}ky+\unicode[STIX]{x1D70E}t}]+O(\unicode[STIX]{x1D716}^{2}),\qquad & \displaystyle\end{eqnarray}$$

(2.18b ) $$\begin{eqnarray}\displaystyle & \displaystyle u=\bar{u}(x)+\unicode[STIX]{x1D716}\,\text{Re}[\hat{u} (x)\text{e}^{\text{i}ky+\unicode[STIX]{x1D70E}t}]+O(\unicode[STIX]{x1D716}^{2}),\quad R=Ut+\unicode[STIX]{x1D716}\,\text{Re}[\hat{R}\text{e}^{\text{i}ky+\unicode[STIX]{x1D70E}t}]+O(\unicode[STIX]{x1D716}^{2}),\qquad & \displaystyle\end{eqnarray}$$

$\text{Re}$

$O(\unicode[STIX]{x1D716}^{2})$

The elastic equations (2.1) become

(2.19a,b ) $$\begin{eqnarray}\hat{p}=T(k^{2}-\unicode[STIX]{x2202}_{x}^{2}){\hat{h}}\quad \text{(tension)},\quad \hat{p}=B(k^{2}-\unicode[STIX]{x2202}_{x}^{2})^{2}{\hat{h}}\quad \text{(bending)}.\end{eqnarray}$$

In the gas region, (2.3) becomes simply $\hat{p}=0$ , which gives the analytical solution

(2.20a,b ) $$\begin{eqnarray}{\hat{h}}=A_{4}\text{e}^{kx}\quad \text{(pulling)},\quad {\hat{h}}=(A_{5}+A_{6}x)\text{e}^{kx}\quad \text{(bending)}\quad \text{ in }x<0,\end{eqnarray}$$

where $A_{4,5,6}$ are as yet unknown constants. (We have imposed the condition that ${\hat{h}}\rightarrow 0$ as $x\rightarrow -\infty$ in order to exclude exponentially growing solutions and thus satisfy (2.8).) The conditions (2.6) for continuity of $h$ and its derivatives at the interface $x=0$ take the form

(2.21a ) $$\begin{eqnarray}\displaystyle & \displaystyle [{\hat{h}}+\bar{h}^{\prime }\hat{R}]_{-}^{+}=[{\hat{h}}^{\prime }+\bar{h}^{\prime \prime }\hat{R}]_{-}^{+}=0\quad \text{(pulling)}, & \displaystyle\end{eqnarray}$$

(2.21b ) $$\begin{eqnarray}\displaystyle & \displaystyle [{\hat{h}}+\bar{h}^{\prime }\hat{R}]_{-}^{+}=[{\hat{h}}^{\prime }+\bar{h}^{\prime \prime }\hat{R}]_{-}^{+}=[{\hat{h}}^{\prime \prime }+\bar{h}^{\prime \prime \prime }\hat{R}]_{-}^{+}=[{\hat{h}}^{\prime \prime \prime }+\bar{h}^{\prime \prime \prime \prime }\hat{R}]_{-}^{+}=0\quad \text{(bending)},\qquad & \displaystyle\end{eqnarray}$$

(2.22a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\hat{h}}^{\prime }-k{\hat{h}}+\bar{h}^{\prime \prime }\hat{R}=0\quad \text{at }x=0^{+}\text{ (pulling)}, & \displaystyle\end{eqnarray}$$

(2.22b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\hat{h}}^{\prime \prime }-2k{\hat{h}}^{\prime }+k^{2}{\hat{h}}={\hat{h}}^{\prime \prime \prime }-2k{\hat{h}}^{\prime \prime }+k^{2}{\hat{h}}^{\prime }+\bar{h}^{\prime \prime \prime \prime }\hat{R}=0\quad \text{at }x=0^{+}\text{ (bending)},\quad & \displaystyle\end{eqnarray}$$

(2.23) $$\begin{eqnarray}{\hat{h}}\rightarrow 0\quad \text{as }x\rightarrow \infty .\end{eqnarray}$$

In the liquid region itself, the lubrication equations in the travelling-wave frame, (2.2a ) and (2.9), become

(2.24a,b ) $$\begin{eqnarray}\hat{u} =-\frac{\bar{h}^{2}}{12\unicode[STIX]{x1D707}}\hat{p}^{\prime }-\frac{2\bar{h}\bar{p}^{\prime }}{12\unicode[STIX]{x1D707}}{\hat{h}},\quad \unicode[STIX]{x1D70E}{\hat{h}}=\left(\frac{\bar{h}^{3}}{12\unicode[STIX]{x1D707}}\hat{p}^{\prime }+\frac{3\bar{h}^{2}\bar{p}^{\prime }}{12\unicode[STIX]{x1D707}}{\hat{h}}+U{\hat{h}}\right)^{\prime }-k^{2}\frac{\bar{h}^{3}}{12\unicode[STIX]{x1D707}}\hat{p}\quad \text{in }x>0.\end{eqnarray}$$

At the interface, the linearized kinematic condition from (2.5a ) is

(2.25a ) $$\begin{eqnarray}\hat{u} +\bar{u}^{\prime }\hat{R}=U(1-f_{1}-f_{1}^{\prime }Ca)\unicode[STIX]{x1D70E}\hat{R}\quad \text{at }x=0^{+}.\end{eqnarray}$$

On the right-hand side, the term $f_{1}$ describes the deposited liquid films, while the term $f_{1}^{\prime }Ca$ captures their variation due to the interfacial velocity being perturbed. However, the sum of these two terms is always less than $1$ , so the parenthesis remains positive in all cases.

On the left-hand side of (2.25a ), the first term represents the perturbation flow, which advects the interface and is the main cause of finger growth or decay. The second term describes how longitudinal variations in the base-flow velocity cause the tips of invading air fingers to be advected by a different velocity compared with the base of the fingers. Since $\bar{u}^{\prime }(0)<0$ (2.17b ), we call this effect ‘kinematic compression’, and it helps stabilize the system by advecting the curved interface back towards its mean position.

The linearized dynamic condition is

(2.25b ) $$\begin{eqnarray}\hat{p}+\bar{p}^{\prime }\hat{R}=-\frac{\unicode[STIX]{x03C0}}{4}\unicode[STIX]{x1D6FE}k^{2}\hat{R}+\frac{2\unicode[STIX]{x1D6FE}f_{2}}{\bar{h}^{2}}({\hat{h}}+\bar{h}^{\prime }\hat{R})-\frac{2\unicode[STIX]{x1D6FE}f_{2}^{\prime }Ca}{\bar{h}}\frac{\unicode[STIX]{x1D70E}\hat{R}}{U}\quad \text{at }x=0^{+},\end{eqnarray}$$

and we interpret it as quantifying the perturbation pressure $\hat{p}$ at the interface (the first term on the left-hand side) in response to various physical effects. The second term on the left-hand side describes how the lack of viscous pressure drop in the air fingers (compared with in the liquid between the fingers) tends to lead to the pressure in the liquid being larger ahead of the fingers than between the fingers. This is the fundamental driving mechanism of the viscous-fingering instability. In the absence of other effects, it would drive a spanwise flow away from the finger tips and cause them to grow.

On the right-hand side of (2.25b ), the first term describes the stabilizing effect of horizontal surface tension. The middle pair of terms encapsulate the effects of vertical surface tension (as modified by the viscous effects captured by $f_{2}$ ) which cause the liquid pressure ahead of the interface to decrease due to the increase in vertical curvature when the cell height at the interface decreases. The first of this pair of terms describes deflection of the sheet and the second describes the ‘taper mechanism’ (Al-Housseiny et al. Reference Al-Housseiny, Tsai and Stone2012). Finally, the last term on the right-hand side of (2.25b ) describes the dependence of the dynamic condition (2.5b ) on the meniscus speed via  $f_{2}$ .

For each value of the wavenumber $k$ , equations (2.19) and (2.22)–(2.25) form an eigenvalue problem for the infinite-dimensional eigenvector $(\hat{R},{\hat{h}}(x))$ and the associated eigenvalue $\unicode[STIX]{x1D70E}$ . We are interested in the dominant eigenmode, i.e. the one whose growth rate has the largest real part and hence grows fastest or decays slowest. Henceforth, we only discuss the dominant solution, and thus consider $\unicode[STIX]{x1D70E}$ as a function of the wavenumber $k$ (for a given base state). If the real part of $\unicode[STIX]{x1D70E}$ is positive, then the dominant eigenmode is growing and the system is linearly unstable. If the real part is negative, then the system is linearly stable to all perturbations with this wavenumber. If the system is linearly stable to perturbations with all wavenumbers $k$ , then it is linearly stable.

3.1 Review of classical case – rigid parallel cell

We briefly review the classical case of viscous fingering, where the upper boundary is a rigid horizontal plate. In this case, we have a fixed uniform height $\bar{h}=h_{0}$ with ${\hat{h}}=0$ , instead of the elastic equations (2.1). It follows from flux conservation (2.2b ) that $\bar{u}=-h_{0}^{2}\bar{p}^{\prime }/12\unicode[STIX]{x1D707}$ is uniform, and the interfacial conditions (2.5) then yield the base solution

(3.1a,b ) $$\begin{eqnarray}\bar{u}=U(1-f_{1}),\quad \bar{p}=-\frac{2\unicode[STIX]{x1D6FE}f_{2}}{h_{i}}-\frac{12\unicode[STIX]{x1D707}U}{{h_{i}}^{2}}(1-f_{1})x\quad \text{in }x>0,\end{eqnarray}$$

where we have used $h_{i}=h_{0}$ for consistency with later equations.

The perturbation pressure is found by solving the lubrication equation (2.24b ) with ${\hat{h}}=0$ , which now expresses flux conservation, together with the dynamic meniscus condition (2.25b ) and the far-field decay condition (2.23). This yields

(3.2) $$\begin{eqnarray}\hat{p}=\hat{R}\left[\frac{12\unicode[STIX]{x1D707}U}{{h_{i}}^{2}}(1-f_{1})-\frac{\unicode[STIX]{x03C0}}{4}\unicode[STIX]{x1D6FE}k^{2}-\frac{2\unicode[STIX]{x1D6FE}f_{2}^{\prime }Ca}{h_{i}}\frac{\unicode[STIX]{x1D70E}}{U}\right]\text{e}^{-kx}\quad \text{in }x>0.\end{eqnarray}$$

Substitution into the kinematic condition (2.25a ) yields the growth rate

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D70E}=kU\frac{1-f_{1}-\displaystyle {\displaystyle \frac{\unicode[STIX]{x03C0}/4}{12Ca}}k^{2}{h_{i}}^{2}}{1-f_{1}-f_{1}^{\prime }Ca+\displaystyle {\displaystyle \frac{2f_{2}^{\prime }kh_{i}}{12}}}.\end{eqnarray}$$

We can trace the various parts of (3.3) back to the physical mechanisms discussed in § 2.3. The last two terms in the denominator capture the dependence of the meniscus conditions (2.5) on the speed of the meniscus via $f_{1}$ and $f_{2}$ , but they do not affect the stability of the system (i.e. the sign of $\unicode[STIX]{x1D70E}$ ) since the denominator is always positive. The stability is thus determined by the sign of the numerator. The first part $(1-f_{1})$ of the numerator describes the viscous effect that drives the instability, while the second part proportional to $k^{2}$ describes the stabilizing effect of horizontal surface tension. If the wavenumber $k$ exceeds a critical wavenumber $k_{c}=[48Ca(1-f_{1})/\unicode[STIX]{x03C0}]^{1/2}/h_{i}$ , then $\unicode[STIX]{x1D70E}<0$ and the perturbation decays; conversely, if $0<k<k_{c}$ then $\unicode[STIX]{x1D70E}>0$ .

We conclude that in a rigid parallel cell there is always a band $0<k<k_{c}$ of unstable wavenumbers and the system is unstable regardless of the value of $Ca$ . However, for $Ca\ll 1$ , corresponding to strong surface tension, we have $k_{c}\sim \sqrt{Ca}/h_{i}$ so the unstable wavelengths are very long and the maximum growth rate $\unicode[STIX]{x1D70E}\sim \sqrt{Ca}\,U/h_{i}$ is small.

Based on this analysis, we define viscous-fingering scales for the wavenumber, growth-rate and pressure perturbations,

(3.4a ) $$\begin{eqnarray}k_{vf}=\frac{\sqrt{Ca}}{h_{i}},\quad \unicode[STIX]{x1D70E}_{vf}=\frac{\sqrt{Ca}\,U}{h_{i}},\quad \hat{p}_{vf}=\frac{12\unicode[STIX]{x1D707}U}{{h_{i}}^{2}}\hat{R}.\end{eqnarray}$$

We also introduce a scale for the height perturbations,

(3.4b ) $$\begin{eqnarray}{\hat{h}}_{vf}=\frac{\hat{p}_{vf}}{Tk^{2}}\quad \text{(pulling)}\quad \text{or}\quad {\hat{h}}_{vf}=\frac{\hat{p}_{vf}}{Bk^{4}}\quad \text{(bending)},\end{eqnarray}$$

based on a balance in (2.19) with $\hat{p}\sim \hat{p}_{vf}$ . We will use these scales to non-dimensionalize the numerical results and to confirm when the instability is of viscous-fingering type.

3.2 Numerical results

Returning to the elastic case, we use a numerical method to investigate the behaviour of (2.19)–(2.25) for various values of the capillary number $Ca=\unicode[STIX]{x1D707}U/\unicode[STIX]{x1D6FE}$ and the peeling aspect ratio ${\mathcal{A}}=(T/12\unicode[STIX]{x1D707}U)^{1/3}$ or $(B/12\unicode[STIX]{x1D707}Uh_{0}^{2})^{1/5}$ . The equations are discretized using a finite-difference scheme and the dominant linear mode is found using a combination of the Arnoldi algorithm (which is used to find the maximal eigenvalue of a linear operator) and the Newton–Raphson method (which is used to refine the result and for numerical continuation to other parameter values) as described in appendix A.

Figure 3.

Structure of the pressure and height perturbations of the dominant mode for $Ca=1$ , ${\mathcal{A}}=10$ , and three values of the spanwise wavenumber $k$ , rescaled using the viscous-fingering scales (3.4). In the top row, the translational mode (3.5) is overlaid as dotted curves. Bottom row: base-state profile showing the extent of the liquid wedge.

We begin by focusing on the particular values $Ca=1$ and ${\mathcal{A}}=10$ . Figure 3 shows the structure of the dominant mode for three different values of the spanwise wavenumber $k$ , and the base-state height profile for comparison. For both peeling-by-pulling and peeling-by-bending, the mode with largest wavenumber decays quickly away from the interface, on the length scale $1/k$ , and becomes negligible before it reaches the tip of the liquid wedge. The modes with smaller wavenumber instead vary on the length scale $L_{w}$ of the wedge, and reach past the wedge tip into the peeling-wave region. For all $k$ , the pressure perturbation scales like $\hat{p}_{vf}$ , which is to be expected from a balance on the left-hand side of (2.25b ). The height perturbations scale like ${\hat{h}}_{vf}$ for large $k$ , when all of the terms balance in (2.19), but for small $k$ , the $k^{2}$ terms become negligible compared with the $\unicode[STIX]{x2202}_{x}^{2}\sim 1/L_{w}^{2}$ terms and hence ${\hat{h}}\ll {\hat{h}}_{vf}$ .

Figure 4.

Growth rate $\unicode[STIX]{x1D70E}$ of the dominant mode for each wavenumber $k$ , for ${\mathcal{A}}=10$ and $Ca=10^{0},10^{-0.5},10^{-1},10^{-1.5},10^{-2}$ , rescaled using the viscous-fingering scales (3.4). For the dashed parts of the curves at small $k$ , the dominant growth rate is complex and only the real part is plotted. The kinks in the curves for small $k$ correspond to different modes being dominant.

Figure 5.

Contours of the rescaled maximal growth rate $\unicode[STIX]{x1D70E}_{m}/\unicode[STIX]{x1D70E}_{vf}=\unicode[STIX]{x1D70E}_{m}/(U\sqrt{Ca}/h_{i})$ as a function of the parameters $Ca$ and ${\mathcal{A}}$ . In the shaded region, the system is stable ( $\unicode[STIX]{x1D70E}\leqslant 0$ for all $k$ ). In the unshaded region, the contour spacing is $0.1$ , with labels indicating the values of some contours. In (b), for sufficiently small $Ca$ and ${\mathcal{A}}$ , no steady base-state solution exists. The results from the rigid-lid and wedge approximations from § 4 are shown for comparison.

The growth rates $\unicode[STIX]{x1D70E}$ of the dominant mode for various values of $k$ are shown in figure 4. We have rescaled $k$ and $\unicode[STIX]{x1D70E}$ using the viscous-fingering scales $k_{vf}=\sqrt{Ca}/h_{i}$ and $\unicode[STIX]{x1D70E}_{vf}=U\sqrt{Ca}/h_{i}$ from (3.4a ). For not too small $Ca$ , the system has growing modes for a finite range of wavenumbers $k$ . For $k$ larger than this range, the growth rate is negative due to the stabilizing effects of surface tension, like in the classical rigid case. For $k$ smaller than this range, there are kinks in the curve, which indicates that different modes are dominant for different ranges of $k$ , with each kink being a transition from one mode being dominant to another.

Several of these small- $k$ modes have complex $\unicode[STIX]{x1D70E}$ and hence oscillatory behaviour, in which case we plot only the real part of $\unicode[STIX]{x1D70E}$ in figure 4. Since all of the oscillating modes are decaying, we will not consider them further, and the modes discussed in detail in this paper all have purely real growth rates.

For $k=0$ the dominant perturbation is a translational mode with $\unicode[STIX]{x1D70E}=0$ , which can be found analytically by observing that $\bar{h}(x-\unicode[STIX]{x1D716}\hat{R})$ is a steady solution of the original equations (in the travelling-wave frame) for any $\unicode[STIX]{x1D716}\hat{R}$ , so there exists a linear perturbation

(3.5a-c ) $$\begin{eqnarray}{\hat{h}}(x)=\lim _{\unicode[STIX]{x1D716}\rightarrow 0}\frac{\bar{h}(x-\unicode[STIX]{x1D716}\hat{R})-\bar{h}(x)}{\unicode[STIX]{x1D716}}=-\bar{h}^{\prime }(x)\hat{R},\quad \hat{p}=-\bar{p}^{\prime }\hat{R},\quad \hat{u} =-\bar{u}^{\prime }\hat{R},\end{eqnarray}$$

that satisfies the perturbation equations exactly. Indeed, the numerical results for $k=0.001k_{vf}$ in figure 3 agree excellently with the profiles (3.5) shown with dotted curves.

For small $k>0$ , a long-wavelength expansion can be performed about the translational mode (3.5); see appendix B for details. In particular, the leading-order corrections to the translation come from the $O(k)$ terms in the interfacial matching conditions (2.22) for the elastic sheet, and result in a negative $O(k)$ correction to $\unicode[STIX]{x1D70E}$ . The physical interpretation is that at long wavelengths each longitudinal strip of the elastic sheet is independent of the others and has the steady-peeling structure at leading order. Where the interface is slightly ahead of its mean position, the peeling solution has advanced further and hence the elastic sheet is slightly higher. Elasticity acts on the resultant spanwise corrugation of the sheet in the gas region to impose a reduction in the peeling slope or curvature at the advanced parts of the interface, which slows down their further advance. (The classical viscous and capillary effects of § 3.1 only come in at  $O(k^{2})$ .)

Hence, unlike in the classical case (§ 3.1), $\unicode[STIX]{x1D70E}$ is negative for sufficiently small $k$ , and hence there is always a lower bound to the range of unstable wavenumbers $k$ (as well as the upper bound provided by stabilizing capillary effects). The unstable range depends on the system parameters, and as $Ca$ decreases, corresponding to surface tension becoming stronger, eventually a critical value (which depends on the value of ${\mathcal{A}}$ ) is reached below which all perturbations decay and hence the system is stable to linear perturbations.

Given a choice of parameters $(Ca,{\mathcal{A}})$ , we can calculate the dominant growth rate $\unicode[STIX]{x1D70E}$ for a range of wavenumbers $k$ to find the maximum growth rate $\unicode[STIX]{x1D70E}_{m}$ , attained at the maximizing wavenumber $k_{m}$ . Figure 5 shows how $\unicode[STIX]{x1D70E}_{m}$ depends on $Ca$ and ${\mathcal{A}}$ . In the shaded region (small $Ca$ or ${\mathcal{A}}$ ), the system is stable, and the maximal growth rate (for any parameter value) is $\unicode[STIX]{x1D70E}_{m}=0$ , attained by the translational mode (3.5). In the unshaded region (large $Ca$ or ${\mathcal{A}}$ ), the system is unstable and figure 5 shows contours of $\unicode[STIX]{x1D70E}_{m}$ rescaled by the viscous-fingering scaling (3.4a ).

Figure 5 shows that the maximal growth rate $\unicode[STIX]{x1D70E}_{m}=O(\unicode[STIX]{x1D70E}_{vf})$ in the unstable region (and we have also verified that the associated wavenumber $k_{m}=O(k_{vf})$ ), which indicates that the instability is of viscous-fingering type. Reducing $Ca$ (e.g. by increasing surface tension) stabilizes the system as expected, like in the classical rigid case. We can give the following heuristic argument for why increasing ${\mathcal{A}}$ has a destabilizing effect. As ${\mathcal{A}}$ increases, the base peeling solution becomes longer and thinner, reducing the effect of the taper mechanism discussed in § 2.3 and becoming geometrically more similar to the rigid cell. Hence, the mechanisms that help stabilize the elastic system are likely to be weaker, making the system more unstable. We will identify these mechanisms and give a more detailed explanation in § 4 below.

4.1 Governing equations for the rigid-lid approximation

We define the rescaled wavenumber

(4.1) $$\begin{eqnarray}K=kh_{i}/\unicode[STIX]{x1D6FC}_{i},\end{eqnarray}$$

and show in § D.1 that the height perturbations ${\hat{h}}$ can be neglected both in the lubrication equation (2.24) and in the dynamical interfacial condition (2.25b ) if $K\gg 1$ . We will make these simplifications, but allow more moderate values $K\gtrsim 1$ in general.

The lubrication equation (2.24) simplifies to

(4.2a,b ) $$\begin{eqnarray}\hat{u} =-\frac{\bar{h}^{2}}{12\unicode[STIX]{x1D707}}\hat{p}^{\prime },\quad 0=(\bar{h^{3}}\hat{p}^{\prime })^{\prime }-k^{2}\bar{h}^{3}\hat{p}\quad \text{in }x>0.\end{eqnarray}$$

Equation (4.2b ) expresses flux conservation in terms of the perturbation pressure $\hat{p}$ but, unlike the classical case (§ 3.1) where $\bar{h}$ is a constant, $\hat{p}$ must be found numerically. However, since (4.2) is linear, homogeneous and second order, and there is a homogeneous boundary condition $\hat{p}\rightarrow 0$ as $x\rightarrow \infty$ , the solution is determined up to a multiplicative constant. We can calculate a numerical solution $\hat{p}=\hat{p}_{rig}(x)$ for the arbitrary normalization $\hat{p}_{rig}(0)=1$ , and the actual solution is then given simply by

(4.3) $$\begin{eqnarray}\hat{p}(x)=\hat{p}(0)\hat{p}_{rig}(x).\end{eqnarray}$$

In particular, we can express the perturbation velocity (4.2a ) at the interface as

(4.4) $$\begin{eqnarray}\hat{u} (0)=\frac{{h_{i}}^{2}}{12\unicode[STIX]{x1D707}}kY(Ca,{\mathcal{A}},k)\,\hat{p}(0),\quad \text{where }Y(Ca,{\mathcal{A}},k)=-\frac{{\hat{p}_{rig}}^{\prime }(0)}{k\hat{p}_{rig}(0)}.\end{eqnarray}$$

The quantity $Y$ , which we call the ‘admittance’, depends only on the base-state profile $\bar{h}(x)$ and the wavenumber $k$ , and captures how the geometry of the cell affects the perturbation flow. For a rigid parallel cell ( $\bar{h}=h_{0}$ ), the perturbation pressure profile would be exponential like (3.2), $\hat{p}_{rig}(x)=\exp (-kx)$ , and we would obtain $Y=1$ . For the non-parallel peeling geometry, the converging cell restricts the flow and hence a given pressure perturbation results in a weaker perturbation flow, and we find that $Y<1$ .

For perturbations that decay rapidly away from the meniscus ( $K\gg 1$ ), we can use a Taylor expansion of $\bar{h}$ about the meniscus position:

(4.5) $$\begin{eqnarray}\bar{h}(x)\approx \bar{h}(0)+\bar{h}^{\prime }(0)x=\unicode[STIX]{x1D6FC}_{i}s,\quad \text{where }s=h_{i}/\unicode[STIX]{x1D6FC}_{i}-x.\end{eqnarray}$$

We call this the ‘wedge approximation’, as it is equivalent to approximating the base-state profile as being a triangular wedge. Inserting this approximation into (4.2) and solving for $\hat{p}$ with the condition that $\hat{p}$ does not diverge as $s\rightarrow 0$ yields

(4.6) $$\begin{eqnarray}\hat{p}\propto \frac{\text{I}_{1}(ks)}{ks},\end{eqnarray}$$

where $\text{I}_{1}$ is the modified Bessel function of the first kind with order $1$ (cf. Grenfell-Shaw & Woods Reference Grenfell-Shaw and Woods2017). This result for $\hat{p}$ gives the admittance as

(4.7a,b ) $$\begin{eqnarray}Y_{w}(K)=\frac{\text{I}_{1}^{\prime }(K)}{\text{I}_{1}(K)}-\frac{1}{K},\quad Y_{w}(K)\sim 1-\frac{3}{2K}+\frac{3}{8K^{2}}\rightarrow 1\quad \text{as }K\rightarrow \infty ,\end{eqnarray}$$

where $K=kh_{i}/\unicode[STIX]{x1D6FC}_{i}$ as before.

Finally, neglecting ${\hat{h}}$ in the interfacial conditions (2.25), and combining with (4.4), yields the main result for the growth rate:

(4.8) $$\begin{eqnarray}\unicode[STIX]{x1D70E}=kU\frac{\displaystyle Y\left(1-f_{1}-{\displaystyle \frac{\unicode[STIX]{x03C0}/4}{12Ca}}k^{2}{h_{i}}^{2}-{\displaystyle \frac{2f_{2}}{12Ca}}\unicode[STIX]{x1D6FC}_{i}\right)-{\displaystyle \frac{-\bar{u}^{\prime }(0)}{kU}}}{\displaystyle 1-f_{1}-f_{1}^{\prime }Ca+{\displaystyle \frac{2f_{2}^{\prime }Ykh_{i}}{12}}}.\end{eqnarray}$$

(We have refrained from eliminating $h_{i}$ and $-\bar{u}^{\prime }(0)$ using (2.17) in order to highlight the physical interpretation of the various terms.)

As in the classical rigid case (§ 3.1), the denominator of (4.8) contains terms that capture the dependence of the meniscus conditions (2.5) on the speed of the meniscus via $f_{1}$ and $f_{2}$ , but it remains positive and hence does not affect the stability of the system. In the numerator, we recognize the first three terms as describing the viscous driving force and horizontal surface tension. The fourth term captures the stabilizing effect of the taper mechanism as discussed in § 2.3. These first four terms describe the pressure perturbation generated by a displacement of the interface, and are multiplied by the admittance $Y$ which quantifies how the flow response to the pressure perturbation is reduced by the converging geometry of the cell. Finally, the fifth term captures the effect of kinematic compression, which is stabilizing since $\bar{u}^{\prime }<0$ .

4.2 Comparison with full linear calculation

Figure 6 shows perturbation pressure profiles from the full linear calculation (two of which also shown in figure 3) with the rigid-lid and wedge approximations overlaid. We find that for the largest value of $k$ the approximations are excellent throughout the domain. For smaller values of $k$ , the approximations are excellent near the meniscus but do not fare as well towards the tip of the wedge. However, the approximations can still be expected to give a reasonable estimate of the growth rate $\unicode[STIX]{x1D70E}$ , since it is determined mainly by the behaviour near the meniscus. For even smaller values of $k$ (not shown), the approximations fail, as is to be expected when the wavelength of the perturbation is greater than the wedge length.

Figure 6.

Profiles of the pressure perturbations of the dominant mode for $Ca=1$ , ${\mathcal{A}}=10$ and three values of $k$ , rescaled using the viscous-fingering scales (3.4a ), comparing results from the full linear calculation and the rigid-lid and wedge approximations (cf. figure 3). The exponential profiles for the classical rigid case are shown rescaled so that they agree with the other results at the interface $x=0$ . Bottom row: base-state height profile and the wedge-approximation profile (4.5).

For comparison, the exponential profiles $\hat{p}\propto \text{e}^{-kx}$ (with an arbitrary overall factor) for a rigid parallel cell are also shown in figure 6. The exponential profile agrees with the others only if $k$ is very large, corresponding to the perturbation decaying on such a short length scale that the converging geometry of the cell has no effect. For smaller values of $k$ , there is a distinct difference between the profiles, and the rigid-lid and wedge approximations constitute a genuine improvement over treating the cell as approximately parallel.

Figure 7.

Growth rate $\unicode[STIX]{x1D70E}$ of the dominant mode for each wavenumber $k$ , for ${\mathcal{A}}=10$ and $Ca=10^{0},10^{-0.5},10^{-1},10^{-1.5},10^{-2}$ , comparing results from the full linear calculation and the rigid-lid and wedge approximations (cf. figure 4).

Figure 7 shows the growth rate $\unicode[STIX]{x1D70E}$ as a function of wavenumber $k$ for ${\mathcal{A}}=10$ and various values of $Ca$ , similarly to figure 4. We find that the rigid-lid and wedge approximations are in excellent agreement with the full linear results over the full range of wavenumbers relevant for the viscous-fingering instability. Only for much smaller $k$ , where $\unicode[STIX]{x1D70E}<0$ , do they deviate from the full linear results, and predict that $\unicode[STIX]{x1D70E}$ tends to a negative value as $k\rightarrow 0$ when in reality $\unicode[STIX]{x1D70E}\rightarrow 0$ . The contour plot of the rescaled maximal growth rate $\unicode[STIX]{x1D70E}_{m}/\unicode[STIX]{x1D70E}_{vf}$ shown in figure 5 includes the results from the rigid-lid and wedge approximations, and we observe that the excellent agreement extends to the entire parameter space. We conclude that the rigid-lid and wedge approximations have completely captured the physics of the viscous-fingering instability in the elastic-walled Hele-Shaw cell.

4.3 Comparison between rigid-lid and wedge approximations

The rigid-lid and wedge approximations differ only in the calculation of the admittance $Y$ : the rigid-lid approximation uses the actual base-state profile $\bar{h}(x)$ , whereas the wedge approximation goes further and uses a simple triangular wedge profile (4.5) with the same height and slope at the interface.

Figure 8.

Admittance $Y$ calculated from the rigid-lid approximation for ${\mathcal{A}}=10$ and $Ca=10^{0},10^{-0.5},10^{-1}$ , and from the wedge approximation (4.7). In (b), results for a parabolic wedge (E 1) are also shown.

Figure 8 shows a comparison between the admittance calculated from the rigid-lid approximation (for ${\mathcal{A}}=10$ and various values of $Ca$ ) and from the wedge approximation $Y=Y_{w}(K)$ (which is a function of $K$ alone). We find excellent agreement between the two for large $K=kh_{i}/\unicode[STIX]{x1D6FC}_{i}$ (i.e.  $1/k\ll L_{w}$ ), as expected, since the wedge approximation is derived from a Taylor expansion of the rigid-lid profile. As $K$ decreases towards 1, the approximations start to differ since $1/k\sim L_{w}$ . For all $K$ , as $Ca$ decreases (and the wedge becomes larger), the rigid-lid results approach those for the wedge, which suggests that the base-state profile is becoming more triangular. For peeling-by-pulling, the agreement is excellent. For peeling-by-bending (figure 8 b), slightly better results are obtained by expanding to second order and approximating the wedge as parabolic; see appendix E for a detailed discussion.

Under the wedge approximation, the expression (4.8) for the growth rate $\unicode[STIX]{x1D70E}$ does not depend directly on ${\mathcal{A}}$ or on the elasticity model for the lid. Instead, these affect the base-state profile and hence the slope $\unicode[STIX]{x1D6FC}_{i}$ , which appears in (4.8) explicitly as well as via the admittance for the wedge approximation (4.7) and via the compression term (2.17b ). Apart from an overall dimensional factor $U/h_{i}$ (or equivalently $U/h_{0}$ , since $h_{i}=h_{0}/f_{1}(Ca)$ ), the growth rate $\unicode[STIX]{x1D70E}$ thus only depends on the capillary number $Ca$ , the slope $\unicode[STIX]{x1D6FC}_{i}$ at the interface and the rescaled wavenumber $kh_{i}$ . Maximizing $\unicode[STIX]{x1D70E}$ over all $kh_{i}$ yields the maximal growth rate $\unicode[STIX]{x1D70E}_{m}$ , which apart from the overall factor $U/h_{i}$ now only depends on $Ca$ and $\unicode[STIX]{x1D6FC}_{i}$ .

Figure 9.

Contours of the rescaled maximal growth rate $\unicode[STIX]{x1D70E}_{m}/\unicode[STIX]{x1D70E}_{vf}=\unicode[STIX]{x1D70E}_{m}/(U\sqrt{Ca}/h_{i})$ (cf. figure 5) as a function of capillary number $Ca$ and slope $\unicode[STIX]{x1D6FC}_{i}$ for peeling-by-pulling and peeling-by-bending under the full linear calculation, and for the universal wedge approximation. The stability boundary $\unicode[STIX]{x1D6FC}_{i}=6Ca$ (4.10) for $Ca\ll 1$ is also shown.

Figure 9 shows contours of the rescaled growth rate as a function of $Ca$ and $\unicode[STIX]{x1D6FC}_{i}$ for the full linear calculations and the wedge approximation. We notice that the results for both peeling-by-pulling and peeling-by-bending collapse excellently onto the universal result obtained from the wedge approximation. This confirms that the instabilities in this system behave like viscous-fingering instabilities in a rigid triangular wedge, and that the aspect ratio ${\mathcal{A}}$ and the elasticity of the lid only affect the instability via the slope of the lid at the meniscus.

4.4 Analysis of contributions to stability

Having verified the rigid-lid and wedge approximations against the full linear calculations, we can now analyse the expression (4.8) for the growth rate to determine the relative importance of the various physical effects. The condition $\unicode[STIX]{x1D70E}<0$ for stability of the mode with a given wavenumber $k$ is equivalent to the condition

(4.9) $$\begin{eqnarray}f_{1}+\frac{2f_{2}}{12Ca}\unicode[STIX]{x1D6FC}_{i}+\frac{\unicode[STIX]{x03C0}/4}{12Ca}k^{2}{h_{i}}^{2}+\frac{-\bar{u}^{\prime }(0)}{kUY}>1.\end{eqnarray}$$

The four terms on the left-hand side are positive, so we interpret them as the ‘contributions to stability’ of four distinct physical mechanisms.

We recognize the first and third terms in (4.9) from the classical rigid case. They represent the effects of wetting films reducing the driving viscous pressure drop and of the horizontal surface tension acting to straighten out the interface. The second term represents the taper mechanism providing a larger resisting pressure drop across the meniscus on fingers penetrating deeper into the converging cell. The fourth term represents kinematic compression in the base flow providing a restoring advection, and this effect is strengthened by the reduction in admittance $Y$ due to the converging cell.

Figure 10.

The contribution to stability of the four terms in (4.9), stacked from bottom to top, for $Ca=1$ and ${\mathcal{A}}=10$ , as a function of wavenumber $k$ , under the rigid-lid and wedge approximations.

We plot the four contributions to stability in (4.9) using stacked-area graphs to highlight both their sum and the relative importance of each term. For example, figure 10 shows the contributions as a function of the wavenumber $k$ for the case $Ca=1$ and ${\mathcal{A}}=10$ . The wavenumbers for which the sum is less than $1$ are unstable (cf. figure 7). Since the base state is fixed, the terms representing films and taper have a constant contribution, while the horizontal surface tension stabilizes larger wavenumbers and kinematic compression stabilizes smaller wavenumbers.

Figure 11.

Properties at marginal stability (i.e. the boundary between stable and unstable in parameter space) as functions of $Ca$ , for peeling-by-pulling and peeling-by-bending using the rigid-lid approximation, and using the wedge approximation. (a) The wavenumber $k$ of the marginal mode, rescaled using the viscous-fingering scale $k_{vf}=\sqrt{Ca}/h_{i}$ (3.4a ). The power law $k/k_{vf}=O(Ca^{7/18})$ is predicted by (D 4b ). (b) The contribution to stability of the four terms in (4.9) for the marginal mode, stacked from bottom to top.

The curve of marginal stability, i.e. the stability boundary, is important since its location in figures 5 and 9 determines whether fingering occurs in any given physical system. On the curve, the sum of the four contributions to stability in (4.9) is larger than $1$ for all $k$ except the wavenumber of the crucial ‘marginal mode’, for which the sum equals $1$ and the growth rate is exactly $\unicode[STIX]{x1D70E}=0$ .

For each value of $Ca$ , we calculate (under both the rigid-lid or wedge approximation, which are in excellent agreement) the marginally stable base state, the associated marginal mode, and the contributions to stability of the various physical mechanisms for this mode. The wavenumber $k$ of the marginal mode, shown in figure 11(a), is $O(k_{vf})$ for this range of $Ca$ , but $k/k_{vf}$ decreases (very) slowly as $Ca\rightarrow 0$ . The contributions to stability (4.9), shown in figure 11(b), sum exactly to $1$ as expected, and for $Ca\gtrsim 0.1$ all four terms are of comparable magnitude. We conclude that all four mechanisms play an important role in stabilizing the system for moderate $Ca$ , and neglecting any one of them would yield a significantly different prediction for the stability boundary.

In the limit $Ca\rightarrow 0$ , the contribution to stability of the taper term dominates in figure 11(b). Hence, the condition (4.9) for the system to be stable to perturbations of all wavenumbers reduces to

(4.10) $$\begin{eqnarray}\unicode[STIX]{x1D6FC}_{i}>\frac{12Ca}{2f_{2}}\approx 6Ca\end{eqnarray}$$

(cf. Al-Housseiny et al. Reference Al-Housseiny, Tsai and Stone2012), which agrees well with the full linear results in figure 11(b). A more detailed analysis of the limit $Ca\rightarrow 0$ for the unstable region, given in § D.2, shows that the wavenumber $k_{m}$ of the mode with maximal growth rate $\unicode[STIX]{x1D70E}_{m}$ satisfies $k_{m}\sim k_{vf}$ away from the stability boundary but decreases to $k_{m}\sim Ca^{7/18}k_{vf}$ (D 4b ) on the stability boundary, which is in excellent agreement with the results in figure 11(a). The analysis also shows that $K=k_{m}h_{i}/\unicode[STIX]{x1D6FC}_{i}\gg 1$ , which justifies the use of the rigid-lid approximation, and that the asymptotic corrections (D 6) to the leading-order result (4.10) are of relative order $O(Ca^{2/3})$ (due to the effects of films and the correction to $f_{2}$ ) and $O(Ca^{7/9})$ (due to the effects of horizontal surface tension and kinematic compression).