Two-interface and thin-filament approximation in Hele-Shaw channel flow

Abstract

When a viscous fluid partially fills a Hele-Shaw channel, and is pushed by a pressure difference, the fluid interface is unstable due to the Saffman–Taylor instability. We consider the evolution of a fluid region of finite extent, bounded between two interfaces, in the limit that the interfaces are close, that is, when the fluid region is a thin liquid filament separating two gases of different pressure. In this limit, we derive a second-order ‘thin-filament’ model that describes the normal velocity of the filament centreline, and evolution of the filament thickness, as functions of the thickness, centreline curvature and their derivatives. We show that the second-order terms in this model, that include the effect of transverse flow along the filament, are necessary to regularise the instability. Numerical simulation of the thin-filament model is shown to be in accordance with level-set computations of the complete two-interface model. Solutions ultimately evolve to form a bubble of rapidly increasing radius and decreasing thickness.

1. Introduction

For a standard model of Hele-Shaw flow in a rectilinear channel, consisting of a semi-infinite inviscid fluid region and a semi-infinite viscous fluid region separated by a single interface, the interface exhibits the Saffman–Taylor instability when the inviscid fluid displaces the viscous, and is stable when the viscous fluid displaces the inviscid (Saffman & Taylor 1958). In the absence of surface tension, exact solution methods exist that exhibit either finite-time singularity or long-time finger formation (Howison 1986a,b); however, the problem is ill-posed. A regularisation such as surface tension is needed to stabilise sufficiently large-wavenumber perturbations. In the presence of surface tension, solutions generically tend to a linearly stable travelling wave solution known as the Saffman–Taylor finger, although for very small surface tension, finite-amplitude noise in experiments and numerical simulations can lead to tip-splitting (Casademunt 2004).

The traditional Saffman–Taylor instability is normally studied on the assumption that the viscous fluid region extends infinitely far along the channel, so that there is only a single interface to consider. However, in reality the fluid region will only have finite extent, so that there are in fact two interfaces: one in which the viscous fluid is being displaced by the driving inviscid fluid and one in which the viscous fluid is the one advancing (see figure 1). In this case, the force driving the fluid region is the pressure difference between the two inviscid fluids on either end of the fluid region. In the absence of surface tension, some classes of exact solutions to the two-interface problem have been found through use of special functions, in both channel geometries (Feigenbaum, Procaccia & Davidovich 2001; Crowdy & Tanveer 2004) as well as annular geometries (Crowdy 2002; Dallaston & McCue 2012). The exact solutions in these studies, however, do not exhibit the interfaces becoming closely separated. Richardson (1982, 1996) similarly finds exact solutions for two-interface channel flow, some of which exhibit topological changes when two cusps form on different parts of the interface at the same point in space (none of these solutions correspond to the two semi-infinite regions of different pressure meeting, however). Farmer & Howison (2006) consider an approximate model where the two interfaces in an unbounded Hele-Shaw cell are very close, resulting in a thin filament of viscous liquid, and construct exact solutions to this approximate model (see § 4.1). All of these zero-surface-tension models are ill-posed.

Figure 1. A finite fluid region in a Hele-Shaw channel. The fluid is driven in the positive $y$ direction by a pressure difference.

The two-interface Hele-Shaw model in the limit that the interfaces are close is of great importance as, even if the fluid region is not initially thin, the effect of the Saffman–Taylor instability will result in a thin fluid region or filament developing (see, for example, the experimental results in Ward & White (2011), Morrow, De Cock & McCue (2023) and Cardoso & Woods (1995), although the latter study focuses on a much less viscous thin annular region in a more viscous ambient fluid). This formation of a thin filament precedes the fluid ‘bursting’, at which point the two inviscid regions meet and the pressure rapidly equalises. The breakup of a thin viscous filament in a Hele-Shaw cell (with surface tension but in the absence of a driving pressure difference) has also been examined by the use of the lubrication approximation (Constantin et al. 1993; Dupont et al. 1993; Goldstein, Pesci & Shelley 1993, 1995; Almgren 1996; Almgren, Bertozzi & Brenner 1996; Goldstein, Pesci & Shelley 1998). The complicated, but self-similar breakup behaviour of the filament in particular (where the filament thickness goes to zero at a finite time and point in space) is detailed in Almgren et al. (1996).

In this article we consider two-interface Hele-Shaw flow in a channel, including the effects of both driving pressure difference and surface tension, with particular focus on the case in which the fluid interfaces become close together. In § 2 we describe the full two-interface model and its stability. In § 3 we derive a simplified model (the thin-filament approximation) that applies when the two interfaces are close together, by applying the lubrication approximation (up to second order) in a coordinate system following the filament centreline. The application of lubrication theory in an evolving, curvlinear coordinate system is similar to that applied in Van De Fliert, Howell & Ockendon (1995) and Howell (2003), although taking such an approximation to higher order is not standard. Our approximation also represents a generalisation of the model by Farmer & Howison (2006) by including both the effects of surface tension as well as higher order terms, which include the effect of flow along the filament. In § 4 we demonstrate the importance of including these higher-order terms by showing (both numerically and by constructing similarity solutions) that the leading-order problem generically blows up in finite time, even in the presence of surface tension. We also find quasi-travelling wave solutions, which may be thought of as the analogue of Saffman–Taylor fingers, which generalise the ‘grim reaper’ exact solution found by Farmer & Howison (2006). In § 5 we compute numerical solutions of the thin-filament model including the higher-order regularising terms. Solutions of the original two-interface model, found using a level-set method, are seen to approach the thin-filament model results as the resolution of the level-set method is increased. We show the general behaviour of our thin-filament model is not to tend toward a quasi-travelling wave solution, but instead to develop a rapidly expanding ‘bubble’ of circular shape and decreasing thickness.

2. Formulation of the two-interface Hele-Shaw flow problem

2.1. Hele-Shaw flow equations

We consider flow in a Hele-Shaw channel of non-dimensional width $2{\rm \pi}$ in the $x$ direction. The fluid region $\varOmega$ is bounded above and below by interfaces $\partial \varOmega _{U}$ and $\partial \varOmega _{L}$, respectively. A non-dimensional pressure difference $P$ acts to push the fluid region in the positive $y$ direction, while surface tension acts on both interfaces (see figure 1). The standard governing equations in non-dimensional form are (e.g. Morrow et al. 2021)

(2.1a)$$\begin{gather} \nabla^2\phi = 0,\quad \boldsymbol x \in \varOmega , \end{gather}$$

(2.1b)$$\begin{gather}v_{n} = \frac{\partial \phi}{\partial n},\quad \boldsymbol x \in \partial\varOmega_{L}, \partial\varOmega_{U}, \end{gather}$$

(2.1c)$$\begin{gather}\phi ={-}P - \sigma\kappa_{L},\quad \boldsymbol x \in \partial\varOmega_{L}, \end{gather}$$

(2.1d)$$\begin{gather}\phi = \sigma\kappa_{U},\quad \boldsymbol x \in \partial\varOmega_{U}, \end{gather}$$

with $\phi$ the velocity potential, $v_{n}$ the normal velocity of each interface, $\kappa _{U}$ and $\kappa _{L}$ the curvatures of each interface, $P$ is the non-dimensional pressure difference and $\sigma$ the non-dimensional surface tension. In (2.1) the normal vector $\hat {\boldsymbol n}$ of each interface is defined so that a flat interface has normal in the positive $y$ direction, and the signs of curvatures $\kappa$ are chosen such that positive $\kappa$ implies a concave interface in the positive $y$ direction (for example, both interfaces in figure 1 have negative curvature at $x=0$ and positive curvature at $x=\pm {\rm \pi}$). This choice of sign convention for both interfaces will simplify the thin-filament derivation in the next section.

For definiteness, we consider no-flux boundary conditions on the channel walls:

(2.2a,b)\begin{equation} \frac{\partial \phi}{\partial x} = 0,\quad y ={\pm} {\rm \pi}, \end{equation}

although many of our results, in particular the thin-filament model derived in § 3, do not rely on these conditions, and could equally apply to an infinitely long fluid region in an unbounded cell, or an annular fluid region represented by a closed curve (in that case, the constant pressure difference $P$ would physically require injection or extraction of air into the interior of the cell).

The system (2.1) is derived from the dimensional system by introducing the length scale $[x]$, the pressure scale $[p]$ and defining the dimensionless velocity potential $\phi$ in terms of the pressure $p$:

(2.3a–c)\begin{equation} [x] = \frac{\ell}{2{\rm \pi}},\quad [p] = \frac{12\mu [x]^{2}}{b^2[t]},\quad \phi = \frac{p_U-p}{[p]}, \end{equation}

so that

(2.4a,b)\begin{equation} P = \frac{p_U-p_L}{[p]},\quad \sigma = \frac{\hat\gamma}{[p][x]}, \end{equation}

where $p_U$ and $p_L$ are the upper and lower inviscid fluid pressures, $b$ is the plate separation, $\ell$ is the channel width, $\mu$ is the viscosity, $\hat \gamma$ is the surface tension and $[t]$ is an arbitrary time scale. While one of the parameters $P$ or $\sigma$ (if non-zero) could be scaled out by choosing an appropriate time scale, retaining both parameters aids in understanding the effects of the terms that arise in our later analysis.

2.2. Linear stability

The two-interface Hele-Shaw configuration has an exact base state where both interfaces are horizontal, and the fluid region moves upward at constant velocity $v_{n} = P/h_0$, where $h_0$ is the distance between the two interfaces. To examine the stability of this configuration, we write the system (2.1) in Cartesian ($x,y$) coordinates, and define the upper and lower interfaces as $y = f_{U}(x,t)$ and $y = f_{L}(x,t)$, respectively. In addition to Laplace's equation for $\phi$, the boundary conditions (2.1b)–(2.1d) are

(2.5a)$$\begin{gather} \phi|_{y=f_{L}} ={-}P - \sigma \frac{(f_{L})_{xx}}{(1+(f_{L})_x^2)^{3/2}},\quad \phi|_{y=f_{U}} = \sigma \frac{(f_{U})_{xx}}{(1+(f_{U})_x^2)^{3/2}}, \end{gather}$$

(2.5b)$$\begin{gather}(f_{L})_t = \phi_y|_{y=f_{L}} - \phi_x|_{y=f_{L}}(f_{L})_x,\quad (f_{U})_t = \phi_y|_{y=f_{U}} - \phi_x|_{y=f_{U}}(f_{U})_x \end{gather}$$

(for notational compactness we use variable subscripts to refer to partial derivatives throughout this article; text or numerical subscripts, however, do not refer to partial derivatives). The base state is (up to arbitrary translation in $y$) represented by $f_{L} = (P/h_0)t - h_0$, $f_{U} = (P/h_0)t$ and $\phi = (P/h_0) \hat y$, where $\hat y = (y-Pt/h_0)$ is the coordinate in the travelling frame.

To examine linear stability then, we impose perturbations on $f_{L}$, $f_{U}$:

(2.6)\begin{equation} \left.\begin{gathered} f_{L}(x,t) = \frac{P}{h_0}t - h_0 + f_{{L}1}(x,t), \\ f_{U}(x,t) = \frac{P}{h_0}t + f_{{U}1}(x,t), \\ \phi(x,y,t) = \frac{P}{h_0}\hat y + \phi_1(x,\hat y,t). \end{gathered}\right\} \end{equation}

On substitution into the boundary conditions:

(2.7)\begin{equation} \left.\begin{gathered} f_{{L}1t} = \phi_{1y},\quad \phi_1 ={-}\sigma f_{{L}1xx} - \frac{P}{h_0}f_{{L}1},\quad \hat y ={-}h_0,\\ f_{{U}1t} = \phi_{1y},\quad \phi_1 = \sigma f_{{U}1xx} - \frac{P}{h_0}f_{{U}1},\quad \hat y = 0. \end{gathered}\right\} \end{equation}

A perturbation with wavenumber $k$ in the $x$ direction takes the form

(2.8a–c)\begin{align} f_{{U}1} = A\cos(kx)\,\mathrm e^{\lambda t},\quad f_{{L}1} = B\cos(kx)\,\mathrm e^{\lambda t},\quad \phi_1 = (c_1\,\mathrm e^{k\hat y} + c_2\,\mathrm e^{{-}k\hat y})\cos(kx)\,\mathrm e^{\lambda t}, \end{align}

where the pressure boundary conditions give $c_1$ and $c_2$ in terms of the interfacial amplitudes $A$ and $B$, and the kinematic conditions result in an eigenvalue problem for $\lambda$:

(2.9)\begin{equation} \lambda \begin{bmatrix} A \\ B \end{bmatrix} =\frac{k}{\sinh(kh_0)} \begin{bmatrix} -\cosh(kh_0)(\sigma k^2 + P/h_0) & (-\sigma k^2 + P/h_0) \\ -(\sigma k^2 + P/h_0) & \cosh(kh_0)(-\sigma k^2 + P/h_0) \end{bmatrix} \begin{bmatrix} A \\ B \end{bmatrix}. \end{equation}

The eigenvalues of this system are thus

(2.10)\begin{equation} \lambda ={-}\sigma k^3 \coth(kh_0) \pm k\sqrt{\left(\frac{P}{h_0}\right)^2 + \frac{\sigma^2 k^4}{\sinh^2(kh_0)}}. \end{equation}

The growth rate curves ($\lambda = \lambda (k)$) are shown in figure 2 (solid line) for parameter values $\sigma = 0.1, h_0 = 0.2$. One eigenvalue is negative for all wavenumbers $k$, while the other is positive for a finite band of wavenumbers that ranges from zero up to a critical wavenumber. The presence of surface tension regularises the system by stabilising the wavenumbers for large $k$. We emphasise that the system is unstable no matter the direction of the pressure gradient (regardless of the sign of $P$), as each direction involves one of the interfaces moving in the unstable direction according to the Saffman–Taylor instability.

Figure 2. Growth rate $\lambda$ of perturbations of given wavenumber $k$ for the two-interface Hele-Shaw model (2.10), the filament model (3.23) (second order in the lubrication parameter) and the leading-order version of this model, which lacks terms corresponding to transverse flow (4.2). Each model has the parameter values $P=1$, $h_0 = 0.2$ and $\sigma = 0.1$. All systems have two eigenvalues for each wavenumber. For the Hele-Shaw model and the filament model, the most unstable eigenvalue is stabilised at a cut-off wavenumber $k$, and agree closely. The leading-order filament model (4.2) is unstable for all wavenumbers, with the positive eigenvalue tending to a constant as $k\to \infty$.

Our linear stability analysis here for a finite viscous fluid evolving in a Hele-Shaw channel is analogous to that presented for a radial geometry with two interfaces (Morrow et al. 2023). Generalisations and alterations to include different viscous fluids on either side of the interfaces have been conducted for both linear and weakly nonlinear frameworks (Gin & Daripa 2015; Anjos & Li 2020).

We close this section by noting some limiting behaviour of (2.10). For large $h_0$, in order to keep the speed of the base state $O(1)$, we need to keep the driving pressure difference $P=O(h_0)$. In that case, $\lambda \sim -\sigma k^3\pm kP/h_0$ as $h_0\rightarrow \infty$. This limit agrees with the well-studied single-interface problem (an infinite body of viscous fluid), with the plus (minus) sign associated with the unstable (stable) direction of flow. Of a more particular interest here is the other limit of $h_0\ll 1$. Again, supposing that $P=O(h_0)$ in order to keep the interface speed $O(1)$, we have

(2.11)$$\begin{gather} \lambda \sim \left[-\sigma\left(\frac{k}{h_0}\right)^2\pm\left(\frac{k}{h_0}\right) \sqrt{\left(\frac{P}{h_0}\right)^2 + \sigma^2\left(\frac{k}{h_0}\right)^2}\right]h_0 \quad\mbox{as}\ h_0\rightarrow 0,\quad k=O(h_0), \end{gather}$$

(2.12)$$\begin{gather}\lambda \sim \left\{\left[\frac{1}{2\sigma}\left(\frac{P}{h_0}\right)^2- \frac{\sigma k^4}{2}\right]h_0,\ -2\sigma\left(\frac{k}{h_0}\right)^2h_0\right\} \quad\mbox{as}\ h_0\rightarrow 0,\quad k=O(1), \end{gather}$$

(2.13)$$\begin{gather}\lambda \sim{-}\frac{\sigma(\cosh(kh_0)\mp 1)}{\sinh(kh_0)} \frac{1}{h_0^3}\quad\mbox{as}\ h_0\rightarrow 0,\quad k=O(1/h_0), \end{gather}$$

where for (2.12) $k=O(1)$ means strictly of order one, while for (2.13) $k=O(1/h_0)$ means strictly of order $1/h_0$.

3. Thin-filament approximation

3.1. Model in intrinsic coordinate system

In this section we derive an approximation of the two-interface flow by considering the thickness of the fluid region (that is, the distance between the two interfaces) to be small. This approximation is valid when the separation of the two interfaces is smaller than the radius of curvature of either interface, but still larger than the plate separation (if the interface separation is of the same order as the plate separation, the depth averaging that leads to two-dimensional Hele-Shaw flow (2.1) is no longer valid).

We start by converting the governing equations and boundary conditions into an intrinsic coordinate system $(x,y) \mapsto (s,n)$ in which the filament is represented by a centreline curve $\bar {\boldsymbol x}(s,t) = (\bar x, \bar y)$, with $s$ the arclength coordinate and $n$ the coordinate normal to this centreline; thus

(3.1)\begin{equation} \boldsymbol x = \bar{\boldsymbol x} + n\boldsymbol n \end{equation}

(e.g. Van De Fliert et al. 1995; Howell 2003). The coordinate system is represented schematically in figure 3. Here $\boldsymbol n = (-\bar y_s, \bar x_s)$ is the centreline normal (again we use the variable subscript $s$ to represent the $s$ partial derivative, for brevity). In this coordinate system we specify the upper and lower interfaces to occur at $n = \pm h(s,t)/2$, respectively, so that the normal thickness is $h$. The upper interface is thus given by the curve

(3.2)\begin{equation} \boldsymbol x_{{U}} = \bar{\boldsymbol x} + \frac{h}{2} \boldsymbol n, \end{equation}

and a (non-unit) normal to this interface (as opposed to the centreline) is then

(3.3)\begin{equation} \boldsymbol n_{{U}} = \left(1 - \frac{h\kappa}{2}\right) \boldsymbol n - \frac{h_s}{2}\boldsymbol s, \end{equation}

where $\boldsymbol s = (\bar x_s, \bar y_s)$ is the centreline unit tangent vector and $\kappa = \bar x_s \bar y_{ss} - \bar y_s \bar x_{ss}$ is the centreline curvature.

Figure 3. A schematic of the coordinate system (3.1) used to derive the second-order lubrication model (3.19), which describes the velocity $v_{n}$ of the centreline and the evolution of the film thickness $h$, in a direction normal to the centreline. In the derivation of the model, it is important to distinguish between the centreline normal $\boldsymbol n$ and normals to the interface (e.g. $\boldsymbol n_{{U}}$ on the upper interface), and velocity components in these respective directions.

The velocity normal to the interface (not the centreline) $v_{{U}}$ is equal to the normal derivative of the potential via the kinematic condition (2.1b):

(3.4)\begin{align} v_{{U}} &= \boldsymbol{\nabla} \phi \boldsymbol{\cdot} \frac{\boldsymbol n_{{U}}}{|\boldsymbol n_{{U}}|} \nonumber\\ &= \frac{1}{|\boldsymbol n_{{U}}|}\left[\frac{\partial \phi}{\partial n}\boldsymbol n + \frac{1}{1-n\kappa} \frac{\partial \phi}{\partial s}\boldsymbol s\right] \boldsymbol{\cdot} \left[\left(1 - \frac{h\kappa}{2}\right) \boldsymbol n - \frac{h_s}{2}\boldsymbol s\right] \nonumber\\ &= \frac{1}{|\boldsymbol n_{{U}}|}\left[\left(1-\frac{h\kappa}{2}\right)\frac{\partial \phi}{\partial n} - \frac{h_s}{2(1-h\kappa/2)} \frac{\partial \phi}{\partial s}\right]. \end{align}

The component of the upper interface velocity in the centreline-normal direction, $v_{{nU}}$, is then

(3.5a)\begin{equation} v_{{nU}} = \frac{|\boldsymbol n_{{U}}| v_{n}}{\boldsymbol n_{{U}} \boldsymbol{\cdot} \boldsymbol n} = \frac{\partial \phi}{\partial n} - \frac{h_s}{2(1-h\kappa/2)^2} \frac{\partial \phi}{\partial s},\quad n = \frac{h}{2}. \end{equation}

Similarly, on the lower interface, the velocity in the centreline-normal direction is

(3.5b)\begin{equation} v_{{nL}} = \frac{\partial \phi}{\partial n} + \frac{h_s}{2(1+h\kappa/2)^2} \frac{\partial \phi}{\partial s},\quad n ={-}\frac{h}{2}. \end{equation}

From the interface velocities (3.5) we compute the centreline velocity and thickness evolution. Let $\eta$ be a Lagrangian coordinate (such that $\eta$ is constant at a point on the curve $\boldsymbol x$ moving in the direction normal to the centreline); then

(3.6)\begin{equation} \frac{{\rm D}\boldsymbol x_{{U}}}{{\rm D} t} = \frac{{\rm D}\bar{\boldsymbol x}}{{\rm D} t} + \frac{h}{2} \frac{{\rm D}\boldsymbol n}{{\rm D} t} + \frac{1}{2}\frac{{\rm D} h}{{\rm D} t}\boldsymbol n. \end{equation}

Here the notation $\textrm {D}/\textrm {D} t$ represents the Lagrangian derivative (holding $\eta$ constant). Taking the centreline-normal component of each term in this equation, and noting that $\textrm {D}\boldsymbol n/\textrm {D} t$ is orthogonal to $\boldsymbol n$, we find

(3.7a,b)\begin{equation} v_{{nU}} = v_{n} + \frac{1}{2}\frac{{\rm D} h}{{\rm D} t},\quad v_{{nL}} = v_{n} - \frac{1}{2}\frac{{\rm D} h}{{\rm D} t}, \end{equation}

with the result for the lower interface derived in very similar fashion to that for the upper. Arranging for $v_{n}$ and $Dh/Dt$, and substituting (3.5), we then find expressions for the normal velocity and thickness evolution in the centreline-normal direction:

(3.8a,b)\begin{align} v_{n} =\left\langle \frac{\partial \phi}{\partial n} \right\rangle - \frac{h_s}{4} \left[\frac{1}{(1-n\kappa)^2}\frac{\partial \phi}{\partial s}\right]_{{-}h/2}^{h/2},\quad \frac{Dh}{Dt} = \left[\frac{\partial \phi}{\partial n}\right]_{h/2}^{h/2} - h_s\left\langle \frac{1}{(1-n\kappa)^2}\frac{\partial \phi}{\partial s}\right\rangle ,\end{align}

where the notation $\langle\,f \rangle = (f|_{n=-h/2} + f|_{n=h/2})/2$ represents a value averaged between the two interfaces.

As is standard in free-boundary lubrication problems, it is useful to be able to express the thickness evolution in a conservative form. To do so, we first note that Laplace's equation (2.1a), under the coordinate transformation (3.1), becomes

(3.9)\begin{equation} \frac{1}{(1-n\kappa)^2}\frac{\partial^2\phi}{\partial s^2} + \frac{\partial^2\phi}{\partial n^2}+ \frac{n\kappa_s}{(1-n\kappa)^3} \frac{\partial\phi}{\partial s} - \frac{\kappa}{1-n\kappa}\frac{\partial \phi}{\partial n} = 0 \end{equation}

which is equivalent to

(3.10)\begin{equation} \frac{\partial }{\partial n}\left[(1-n\kappa)\frac{\partial \phi}{\partial n}\right] + \frac{\partial }{\partial s} \left[\frac{1}{1-n\kappa}\frac{\partial \phi}{\partial s}\right] = 0. \end{equation}

Integrating over $n$ gives

(3.11)\begin{equation} \left[\frac{\partial \phi}{\partial n}\right]_{{-}h/2}^{h/2} - h\kappa\left\langle \frac{\partial \phi}{\partial n} \right\rangle ={-}\int_{{-}h/2}^{h/2} \frac{\partial }{\partial s}\left[ \frac{1}{1-n\kappa}\frac{\partial \phi}{\partial s}\right]{\rm d} n \end{equation}

so that, using (3.8a,b),

(3.12)\begin{align} \frac{{\rm D} h}{{\rm D} t} - \kappa h v_{n} &={-}\int_{{-}h/2}^{h/2} \frac{\partial }{\partial s} \left[\frac{\phi_s}{1-n\kappa}\right]{\rm d} n - h_s\left\langle\frac{\phi_s}{(1-n\kappa)^2}\right\rangle + \frac{hh_s\kappa}{4}\left[\frac{\phi_s}{(1-n\kappa)^2}\right]_{{-}h/2}^{h/2} \nonumber\\ &={-}\int_{{-}h/2}^{h/2} \frac{\partial }{\partial s}\left[\frac{\phi_s}{1-n\kappa}\right]{\rm d} n - \frac{h_s}{2}\left.\frac{\phi_s}{1-n\kappa}\right|_{n=h/2} - \frac{h_s}{2} \left.\frac{\phi_s}{1-n\kappa}\right|_{n={-}h/2}. \end{align}

Thus,

(3.13)\begin{equation} \frac{{\rm D} h}{{\rm D} t} = \kappa h v_{n} -\frac{\partial }{\partial s}\int_{{-}h/2}^{h/2} \frac{1}{1-n\kappa}\frac{\partial \phi}{\partial s}{\rm d} n. \end{equation}

The first term on the right-hand side of (3.13) represents the effect of dilation or compression of the film as its length changes, while the second term represents the contribution of flux along the filament.

The expressions (3.8a,b) and (3.13) are exact, but we have not yet determined the potential $\phi$. We now formally take the lubrication approximation by substituting $n = \epsilon N$, $h = \epsilon H$ and $t = \epsilon T$, taking $\epsilon \ll 1$. In this expansion we assume the curvature $\kappa = O(1)$ (or smaller), and the surface tension $\sigma = O(1)$. Unlike most standard lubrication models where the leading-order behaviour in $\epsilon$ is all that is required, in our case we will see that terms at order $\epsilon ^2$ are necessary to fully regularise the instability, and so we expand to this order.

Writing $\phi = \phi _0 + \epsilon \phi _1 + \epsilon ^2\phi _2 + \cdots$ and substituting into Laplace's equation (3.9), we have

(3.14a–c)\begin{equation} \frac{\partial ^2\phi_0}{\partial N^2} = 0,\quad \frac{\partial ^2\phi_1}{\partial N^2} = \kappa\frac{\partial \phi_0}{\partial N},\quad \frac{\partial ^2\phi_2}{\partial N^2} = \kappa\frac{\partial \phi_1}{\partial N} - \kappa^2N\frac{\partial \phi_0}{\partial N} - \frac{\partial ^2\phi_0}{\partial s^2}, \end{equation}

with pressure conditions from (2.1c) and (2.1d):

(3.15)\begin{equation} \left.\begin{gathered} \phi_0 = \sigma\kappa, \quad \phi_1 = \sigma K_1, \quad \phi_2 = \sigma K_2,\quad N= \frac{H}{2}, \\ \phi_0 ={-} P - \sigma\kappa, \quad \phi_1 = \sigma K_1, \quad \phi_2 ={-}\sigma K_2,\quad N={-}\frac{H}{2}. \end{gathered}\right\} \end{equation}

Here $\kappa$ is again the centreline curvature, and $K_1$ and $K_2$ are order $\epsilon$ and $\epsilon ^2$ corrections to the curvature on the upper interface, respectively:

(3.16a,b)\begin{equation} K_1 = \frac{H\kappa^2}{2} + \frac{H_{ss}}{2},\quad K_2 = \frac{H^2\kappa^3}{4} + \frac{HH_s\kappa_s}{4} + \frac{HH_{ss}\kappa}{2} + \frac{H_s^2\kappa}{8} \end{equation}

(the corrections on the lower interface are the same, with appropriate changes in sign). Solving for each component of the potential $\phi$ term by term gives

(3.17a–c)\begin{equation} \left.\begin{gathered} \phi_0 = V_0N - \frac{P}{2},\quad \phi_1 = \frac{\kappa V_0}{2} \left(N^2 - \frac{H^2}{4}\right) + \sigma K_1, \\ \phi_2 = \frac{2\kappa^2V_0 - V_{0ss}}{6}\left(N^3 - \frac{H^2}{4}N\right) + \frac{2\sigma N}{H}K_2, \end{gathered}\right\} \end{equation}

where

(3.18a)\begin{equation} V_0 = \frac{P + 2\sigma \kappa}{H} \end{equation}

is the leading-order centreline velocity. Substitution of these expressions into the expansions of (3.8a,b) and (3.13) results in expressions for the rescaled normal velocity $V_{n} = \epsilon ^{-1}v_{n}$ and $\textrm {D} H/rD T = \textrm {D} h/\textrm {D} t$:

(3.18b)\begin{equation} V_{n} = V_0 + \epsilon^2\left[\frac{2H^2\kappa^2V_0 - H^2V_{0ss}}{12} - \frac{HH_sV_{0s}}{4} + \frac{2\sigma K_2}{H} \right] \end{equation}

and

(3.18c)\begin{equation} \frac{{\rm D} H}{{\rm D} T}= \kappa H V_{n} - \epsilon^2\frac{\partial }{\partial s} \left[-\frac{H^3\kappa_sV_0}{12} - \frac{H^2 H_s\kappa V_0}{4} + H\sigma K_{1s} \right].\end{equation}

In (3.18b) and (3.18c), the order of terms in the lubrication parameter $\epsilon$ is explicit. We can of course formally remove the scaling by returning to $h, t$ variables, resulting in

(3.19a)$$\begin{gather} v_{n} = v_0 + \left(\frac{2h^2\kappa^2v_0 - h^2v_{0ss}}{12} - \frac{hh_s v_{0s}}{4} + \frac{2\sigma \kappa_2}{h} \right) , \end{gather}$$

(3.19b)$$\begin{gather}\frac{{\rm D} h}{{\rm D} t} = \kappa h v_{n} - \frac{\partial }{\partial s}\left[-\frac{h^3\kappa_s v_0}{12} - \frac{h^2h_s\kappa v_0}{4} + \sigma h \kappa_{1s}\right] , \end{gather}$$

with $v_0 = (P + 2\sigma \kappa )/h$ the leading-order velocity and $\kappa _1 = \epsilon ^{-1} K_1$, $\kappa _2 = \epsilon ^{-2}K_2$ the rescaled versions of (3.16a,b):

(3.19c)\begin{equation} \kappa_1 = \frac{h\kappa^2}{2} + \frac{h_{ss}}{2},\quad \kappa_2 = \frac{h^2\kappa^3}{4} + \frac{hh_s\kappa_s}{4} + \frac{hh_{ss}\kappa}{2} + \frac{h_s^2\kappa}{8}. \end{equation}

The higher-order terms have been included in the above derivation as they are needed to fully regularise the instability. This will become apparent when we examine the behaviour of the leading-order system in § 4. In particular, the higher spatial derivative of thickness $h$ resulting from the interface curvature correction term $\kappa _{1}$ plays a crucial role. In the case where driving pressure $P$ and centreline curvature $\kappa$ both vanish, (3.19) reduces to the well-known thin-film equation

(3.20)\begin{equation} \frac{\partial h}{\partial t} ={-}\frac{\sigma}{2}\frac{\partial }{\partial s}\left[h\frac{\partial ^3h}{\partial s^3}\right], \end{equation}

where $h = h(s,t)$ and arclength $s$ becomes time-independent. This equation has been studied extensively in the context of droplet breakup in Hele-Shaw flow (Constantin et al. 1993; Dupont et al. 1993; Goldstein et al. 1993, 1995; Almgren 1996; Almgren et al. 1996; Goldstein et al. 1998). In our case, this fourth-order spatial derivative term stabilises high-wavenumber perturbations, as is shown in the following stability analysis of the filament model.

3.2. Stability

As with the full two-interface model (2.1), the thin-filament approximation (3.19) has an exact solution comprising a straight filament of uniform thickness $h_0$ moving upward with speed $P/h_0$. To test the stability of this straight filament, we write the centreline $\bar {\boldsymbol x}(s,t)$ in Cartesian coordinates as $y = f(x,t)$, write thickness $h$ as a function of $x$ and $t$ and perturb the straight filament:

(3.21a,b)\begin{equation} h(x,t) = h_0 + \tilde h\exp({\mathrm i k x + \lambda t}),\quad f(x,t) = \frac{P}{h_0}t + \tilde f\exp({\mathrm i k x + \lambda t}). \end{equation}

In the linear approximation, $\textrm {D}/\textrm {D} t = \partial /\partial t + O(f_x^2)$, $s = x + O(f_x^2)$ and $\kappa = f_{xx} + O(f_x^2)$. On substituting these expansions into (3.19) we obtain the eigenvalue problem

(3.22)\begin{equation} \lambda \begin{bmatrix} \tilde f \\ \tilde h \end{bmatrix} =\begin{bmatrix} -2\sigma k^2/h_0 & -P/h_0^2 + Pk^2/12 \\ -Pk^2 - Ph_0^2 k^4 /12 & -\sigma h_0 k^4/2 \end{bmatrix} \begin{bmatrix} \tilde f \\ \tilde h \end{bmatrix}. \end{equation}

The eigenvalues are thus

(3.23)\begin{equation} \lambda ={-}\left(\frac{\sigma k^2}{h_0} + \frac{\sigma h_0 k^4}{4}\right) \pm \sqrt{\left(\frac{\sigma k^2}{h_0} + \frac{\sigma h_0 k^4}{4}\right)^2 + \frac{P^2k^2}{h_0^2} - \sigma^2 k^6 - \frac{P^2h_0^2 k^6}{144}}. \end{equation}

We note that the linear stability of the filament model (3.19) does not agree exactly with the full model, even for zero surface tension. When $\sigma =0$, the error in the filament model is due to the final term $P^2h_0^2 k^6/144$, which arises from the product of the two higher-order correction terms in the off-diagonal terms in the determinant. This term is (implicitly) of order $\epsilon ^4$, which is of higher order than the model is accurate to. If further terms were taken in the lubrication expansion, they would contribute further order $\epsilon ^4$ terms which would presumably cancel with the term present here.

For non-zero surface tension $\sigma$, (3.23) has the same limiting behaviours (2.11) and (2.12) as the full model, but differs slightly from (2.13). In the latter case, both (2.10) and (3.23) give $\lambda =O(1/h_0^3)$ as $h_0\rightarrow 0$ for $k=O(1/h_0)$, but with a different prefactor. Thus, our thin-filament approximation recovers the same leading-order linear stability behaviour as the full model for small and moderate wavenumbers, which is all we can expect from such a lubrication model. For very large wavenumbers (i.e. very small wavelengths of perturbation), with $k\gg 1/h_0$, while the scalings between the eigenvalues of the full and filament model may differ, these modes of perturbation decay very quickly and therefore any differences in the models are of no practical consequence.

In figure 2 we compare the eigenvalues of the thin-filament approximation (3.23) against those of the full problem (2.10). For even moderately small thickness ($h_0 = 0.2$) and small surface tension ($\sigma = 0.1$), the agreement is excellent in the region in which eigenvalues are positive; the difference for negative values occurs for much larger $k$ than depicted. These linear stability results do indicate, however, that the approximation will not be valid if surface tension is zero (or much smaller than the filament thickness), which is unsurprising given the implicit assumption that $\sigma = O(1)$ in the lubrication parameter $\epsilon$.

4. Properties of the leading-order model

In this section we consider the properties of the leading-order version of (3.19), that is, neglecting terms of order $\epsilon ^2$:

(4.1a,b)\begin{equation} v_{n} = v_0 = \frac{P + 2\sigma \kappa}{h},\quad \frac{{\rm D} h}{{\rm D} t} = \kappa h v_0, \end{equation}

where $\kappa$ is the centreline curvature. In the leading-order model (4.1a,b), there is no fluid flow tangent to the centreline; the change in filament thickness is due purely to the stretching or compression of the filament.

Potential issues with the leading-order model can first be observed in the linear stability analysis of a flat filament. The stability analysis of the leading-order model (4.1a,b) is the same as for the higher-order model with the higher-order terms absent; the eigenvalues of this leading-order system are thus

(4.2)\begin{equation} \lambda ={-}\frac{\sigma k^2}{h_0} \pm \sqrt{\left(\frac{\sigma k^2}{h_0}\right)^2 + \frac{P^2k^2}{h_0^2}}. \end{equation}

For small $k$, (4.2) coincides with the leading-order behaviour (2.11), while for $k=O(1)$, (4.2) is no longer a reasonable approximation for the full model. Indeed, unlike in (3.23), where both eigenvalues become negative for sufficiently large wavenumber $k$, in (4.2), one eigenvalue tends to a positive, constant value as $k\to \infty$:

(4.3a,b)\begin{equation} \lambda \sim \frac{P^2}{2\sigma h_0},\quad k \to \infty \end{equation}

(see figure 2). Thus the system is (significantly) unstable to perturbations at arbitrarily small spatial scales, even in the presence of surface tension. While this analysis does not suggest the leading-order problem is technically ill-posed if $\sigma > 0$ (the eigenvalues do not become arbitrarily large), this large-$k$ behaviour strongly suggests that singularities in curvature will generically form (see e.g. Dallaston & McCue 2014). We confirm the self-similar formation of curvature singularities in § 4.2.

4.1. The unregularised (zero surface tension) leading-order model

In the case $\sigma = 0$, the leading-order model (4.1a,b) reduces to that considered by Farmer & Howison (2006). In this case the problem is indeed ill-posed; the eigenvalues (4.2) are $\lambda = \pm (P/h_0)k$, as are the eigenvalues of the full two-interface problem (2.10), so one eigenvalue is arbitrarily large as $k\to \infty$. Solutions that exhibit this ill-posedness may be constructed by showing that the centreline is given by level curves of a harmonic function. We summarise the approach of Farmer & Howison (2006) and further examine this approach here.

Assuming the filament centreline $(x(\eta,t),y(\eta,t))$ and the thickness $h(\eta,t)$ are parametrised by a Lagrangian coordinate $\eta$, then

(4.4a,b)\begin{equation} \frac{\partial x}{\partial t}={-}\frac{P y_\eta}{h\sqrt{x_\eta^2 + y_\eta^2}},\quad \frac{\partial y}{\partial t}=\frac{P x_\eta}{h\sqrt{x_\eta^2 + y_\eta^2}}. \end{equation}

For $\sigma =0$, the evolution of $h$ (4.1a,b) is equivalent to conservation of mass between a point $\eta$ and reference point $\eta _0$. That is, we may define an area function $A(\eta )$, such that

(4.5)\begin{equation} A(\eta)=\int_{\eta_0}^\eta h(\bar{\eta},0)\,\mathrm{d}\bar{\eta} =\int_{\eta_0}^\eta h(\bar{\eta},t)\sqrt{x_\eta^2 + y_\eta^2}\,\mathrm{d}\bar{\eta} \end{equation}

is constant in time on a point on the centreline moving with its normal velocity. Choosing to scale time such that $P=1$, we arrive at the following:

(4.6a,b)\begin{equation} \frac{\partial x}{\partial t}={-}\frac{\partial y}{\partial A},\quad \frac{\partial y}{\partial t}=\frac{\partial x}{\partial A}. \end{equation}

These are Cauchy–Riemann equations relating $(x,y)$ to $(A,t)$. This line of argument was used by Farmer & Howison (2006) to demonstrate that $w = A + \mathrm i t$ must be an analytic function of the complex spatial variable $z=x+\mathrm {i}y$; thus, for a given time $t$, the centreline is the level curve of the harmonic function $t(x,y)$.

Given the definition of $A$ (4.5), the thickness $h$ may be calculated by

(4.7)\begin{equation} h = \frac{A_xx_\eta + A_y y_\eta}{\sqrt{x_\eta^2 + y_\eta^2}} = \frac{A_x t_y - A_y t_x}{\sqrt{t_x^2 + t_y^2}} = \sqrt{t_x^2 +t_y^2} = |w'(z)|. \end{equation}

This thickness will go to zero at a critical point $z_c$ where $w'(z_c) = 0$. As this is also a point where the conformal map between $w$ and $z$ ceases to be smooth, we expect to see a singularity in the curvature in the centreline there. The preimage of the straight line $w = A + \mathrm i t_c$ that passes through the critical point is the centreline in the $z$ plane; assuming $w''(z_c)\neq 0$ the centreline must therefore have a corner with an angle of ${\rm \pi} /2$ at $z_c$ (and not a cusp, as suggested by Farmer & Howison (2006)). If the initial condition is such that $w''(z_c)$ also vanishes but the third derivative is non-zero, the corner angle is ${\rm \pi} /3$, and so on.

As an example, consider an initial condition with the centreline on $y=0$ and initial thickness given by $h(x,0) = \delta [1-a\cos x]$, with $\delta > 0$ and $0 < a < 1$. This initial condition corresponds to an initially horizontal filament that is thinner near $x=0$. We thus have $A = \delta [x - a\sin x]$ at $t=0$ (determined by choosing our reference point $\eta _0$ to lie on $x=0$), and, analytically continuing into the complex plane,

(4.8)\begin{equation} A + \mathrm i t = \delta(z - a\sin z). \end{equation}

Taking the imaginary part, we find that the centreline location is given implicitly by

(4.9)\begin{equation} t = \delta(y - a\cos x\sinh y). \end{equation}

The critical point occurs for $z_c = \mathrm i\cosh ^{-1}(1/a)$ and time $t_c = \delta (\cosh ^{-1}(1/a) - \sqrt {1-a^2})$. The centreline profiles of this solution, along with the upper and lower interfaces (found by adding and subtracting half the thickness $h$ (4.7), respectively, in the normal direction) are plotted in figure 4(a), showing the formation of the ${\rm \pi} /2$ angle as $t\to t_c^-$.

Figure 4. Exact solutions to the thin-filament equation in the absence of surface tension. (a) The solution (4.9) that evolves from an initial condition that results in corner formation with the generic angle of ${\rm \pi} /2$ and (b) the solution (4.11) that evolves from an initial condition that results in corner formation with a non-generic angle of ${\rm \pi} /3$. The parameter values are $a=0.1$, $\delta =0.2$. In both of these examples, the centrelines are plotted in solid blue while the black dashed lines are the upper and lower interfaces (found by adding or subtracting half the thickness $h$ in the normal direction). (c) The profiles and (d) the thickness for the example (4.12) with $\delta =0.2$ and $a=0.5$, which forms a zero-angle cusp with infinite thickness.

As an example of an initial condition that leads to a non-${\rm \pi} /2$ angle, consider

(4.10)\begin{equation} h(x,0) = \delta[1 - a\cos x]^2 \end{equation}

which on integrating results in

(4.11)\begin{equation} A + \mathrm i t = \delta \left[\left(1 + \frac{a^2}{2}\right)z + \frac{a^2}{4}\sin 2z - 2a\sin z\right]. \end{equation}

Again, the centreline is determined by taking the contours of the imaginary part of this function. In this case, since $w'(z_c) = w''(z_c) = 0$, the corner angle at the critical point where $h\to 0$ is ${\rm \pi} /3$. Profiles of this solution are shown in figure 4(b).

The previous two exact solutions show only a subset of the possible singular behaviours of the ill-posed zero-surface-tension model; as centrelines are given by level curves of a harmonic function $t = \Im (w(z))$, with thickness given by (4.7), a variety of other solution behaviours are possible. One such possibility is a zero-angle cusp that results from a square root singularity, an example of which is the complex function

(4.12)\begin{equation} A + \mathrm i t = w(z) = \delta (z - a\mathrm i\sqrt{1-\mathrm e^{-\mathrm i z}}). \end{equation}

This function has a square root singularity at $z=0$, as well as the appropriate periodicity. The level curve that passes through this singularity occurs at $t=0$, so as $t\to 0^{-}$, a zero-angle, backward-facing cusp forms on the interface. The thickness $h$, which is given by $|w'(z)|$, becomes infinite at the cusp, so that the interface at the cusp becomes stationary. The profiles and thickness corresponding to (4.12) are shown in figures 4(c) and 4(d), respectively, for the parameter values of $\delta = 0.2$ and $a=0.5$.

4.2. Singularity formation in the leading-order model with surface tension

We now further consider the leading-order model (4.1a,b) in the presence of surface tension. As eigenvalues (4.2) are positive for arbitrarily large wavenumber, we expect the possibility of singularities in curvature, and this does indeed generically occur, as we both demonstrate through the construction of similarity solutions and corroborate with numerical simulations.

To establish the existence of curvature singularities we perform a self-similar analysis. This analysis is simplest to carry out in Cartesian coordinates. Let the thickness $h = h(x,t)$ be a function of $x$ and $t$, and define the centreline to be given by the function $y = f(x,t)$; we then have

(4.13a,b)\begin{equation} v_0 = \frac{1}{\sqrt{1+f_x^2}} \frac{\partial f}{\partial t},\quad \frac{{\rm D} h}{{\rm D} t} = \frac{\partial h}{\partial t} - \frac{v_0\, f_x}{\sqrt{1+f_x^2}} \frac{\partial h}{\partial x}, \end{equation}

where $\partial /\partial t$ is the time derivative with $x$ held constant. Using (4.1a,b) then results in

(4.14a,b)\begin{equation} \frac{\partial f}{\partial t} = \frac{\sqrt{1+f_x^2}}{h}\left(P + 2\sigma \frac{f_{xx}}{(1+f_x^2)^{3/2}}\right),\quad \frac{\partial h}{\partial t} = \left(\frac{f_x h_x}{1+f_x^2} + \frac{f_{xx}h}{(1+f_x^2)^2}\right) \frac{\partial f}{\partial t}. \end{equation}

Further simplification is achieved by defining the filament thickness in the $y$ direction ${\bar h}(x,t) = h\sqrt {1+f_x^2}$ (as opposed to the thickness $h$ in the centreline-normal direction), which results in the system of equations

(4.15a,b)\begin{align} \frac{\partial f}{\partial t} = \frac{1+f_x^2}{{\bar h}}\left(P + 2\sigma \frac{f_{xx}}{(1+f_x^2)^{3/2}}\right),\quad \frac{\partial {\bar h}}{\partial t} = \frac{\partial }{\partial x}\left[f_x\left(P + 2\sigma\frac{f_{xx}}{(1+f_x^2)^{3/2}} \right) \right]. \end{align}

By introducing ${\bar h}$, (4.15a,b) has been written in conservative form, which highlights the fact that mass is indeed conserved in this system.

Assume a singularity occurs at a time $t=t_c$ at $x=x_c$, and let

(4.16a–c)\begin{equation} f \sim f_0(t) + (t_c-t)^\alpha F(\xi),\quad {\bar h} \sim (t_c-t)^\beta H(\xi),\quad \xi = \frac{x-x_c}{(t_c-t)^\gamma}, \end{equation}

where the similarity exponents $\alpha,\beta,\gamma$ are to be determined. Furthermore, we assume the profile and thickness are symmetric in $x$. Assuming that $\alpha > 1$, the dominant term in the velocity $f_t$ is $\dot f_0 = \dot f_0(t_c)$. Thus, on balancing terms we find $\beta = -1$ and $\alpha = 2\gamma -1$, with $\gamma$ being undetermined (a second-kind self-similarity). Given $\alpha > \gamma$ (which we will check for consistency after the fact), the dominant terms in (4.15a,b) become

(4.17a,b)\begin{equation} \dot f_0 = 2\sigma \frac{F''}{H},\quad H + \gamma\xi H' = 2\sigma[F'F'']'. \end{equation}

These equations can be further scaled to remove $\dot f_0$ and $\sigma$. Let $F = (\,\dot f_0)^{-1}\hat F$ and $H = 2\sigma (\,\dot f_0)^{-2} H$, then

(4.18a,b)\begin{equation} \hat H = \hat F'',\quad \hat H + \gamma\eta \hat H' = [\hat F'\hat F'']'. \end{equation}

Let $u = \hat F'$ and eliminate $\hat F$, so that $u$ satisfies the second-order equation

(4.19)\begin{equation} u'' = \frac{u'(u'-1)}{\gamma\xi - u}. \end{equation}

For symmetry we require $u$ to be odd in $\xi$; thus, when $\xi =0$, $u=0$, which is therefore a singular point of the ordinary differential equation. By expanding near this point we find that the similarity exponent $\gamma$ is uniquely specified by the first odd power in the expansion greater than unity:

(4.20)\begin{equation} u \sim \xi - C\xi^n,\quad \xi \to 0,\enspace n = 3, 5, \ldots, \end{equation}

where $C$ is an arbitrary scaling factor. Expanding near the singular point results in $\gamma = n/(n-1) > 1$, which is consistent with the assumptions made on the similarity exponents above.

Given $\gamma = n/(n-1)$, the equation for $u$ can be solved implicitly by letting $u$ be the independent variable, which allows us to construct parametric solutions for $\hat F$ and $\hat G$:

(4.21a–c)\begin{equation} \xi = u + Cu^n,\quad \hat H = \frac{1}{1+nCu^{n-1}},\quad \hat F = \frac{u^2}{2} + \frac{n}{n+1} C u^{n+1}. \end{equation}

While $n$ can be any odd number $\geq 3$, dependent on the initial condition, the most generic case will be $n=3$, in which case $\gamma = 3/2$ and $\alpha =2$. The constant $C$ is arbitrary, due to scale invariance in (4.15a,b) in the limit that curvature dominates over the driving pressure $P$. In a given case, the scale $C$, velocity $\dot f_0$ and exponent $n$ will all depend on the initial condition of the time-dependent problem.

We provide numerical evidence for the curvature singularity formation by numerically solving the system (4.15a,b). This computation is performed in MATLAB using finite-difference discretisation along with MATLAB's ode15s algorithm for time stepping. Parameters $P = 1$, $\sigma = 0.5$ and an initial condition of $f(x,0) = -\cos (x)$ and $h(x,0) = 1$ is chosen in order to start with high curvature near $x=0$, where the singularity will occur.

In figure 5 we plot the results of this numerical computation. The singularity time $t_c$ is estimated by fitting a straight line through $h_{max}(t)^{-1}=(\max _x h(x,t))^{-1}$, which occurs at $x=0$, and should (according to the analysis above) go to zero linearly. The centreline velocity $\dot f_0(t)$ at the maximum thickness is observed to tend to a non-zero constant, from which $\dot f_0$ is estimated. Scaling the profiles near the singularity time into similarity variables $\xi, \hat F, \hat H$, we observe collapse. The exact similarity profiles (4.21a–c), with a suitable fitted value of the arbitrary constant $C$, match well with the numerical profiles.

Figure 5. Curvature singularity formation in the leading-order Hele-Shaw filament model (4.15a,b) for initial condition $f(x,0) = -\cos (x)$, and ${\bar h}(x,0) = 1$. (a) Centreline profiles $f(x,t)$ approaching a curvature singularity and (b) filament thickness becoming unbounded at singularity time of $t_c \approx 0.77$. These profiles (solid lines) collapse onto similarity profiles (c) $\hat F(\xi )$ and (d) $\hat H(\xi )$, which asymptotically match the exact similarity solution (4.21a–c) (dotted lines) as $t\to t_c$. The scaling factor $C\approx 0.3$ is fitted to the profiles. (e) The singularity time (large circle) is found by fitting a straight line approximation to the reciprocal of the maximum thickness near the singularity time, while (f) the speed of the centreline $\dot f_0 = f_t(0, t_c)$ (large circle) at the singularity is similarly estimated. The five smaller points in (e,f) are the times at which the scaled profiles are plotted in (c,d).

The curvature singularity exhibited by this model is weaker than a corner singularity, in that as $\alpha > \gamma$, the first derivative goes to zero in the neighbourhood of the singularity, even while the curvature becomes unbounded. It is also interesting to note that the singularity is not dependent on (and does not require) the driving pressure $P$; its cause is the presence of surface tension pulling regions of high positive curvature inward, concentrating the thickness at a single point. It is thus only present as the leading-order model (4.1a,b) has a surface-tension-dependent velocity, but no regularising term that penalises higher derivatives of thickness, as appears in the second-order model in (3.19b).

4.3. Quasi-travelling wave solutions

The system (4.15a,b), while insufficiently regularised to simulate the full dynamics of a thin filament, does exhibit a quasi-travelling wave solution of the form (in Cartesian variables)

(4.22a,b)\begin{equation} f ={-}B\log(t_0-t) + F(x),\quad {\bar h}(x,t) = (t_0-t)H(x), \end{equation}

where $t_0$ is a finite time. A solution of this form is similar to, but not exactly, a travelling wave, as the centreline has a fixed shape but moves to infinity (with speed unbounded) as $t\to t_0$, while the thickness linearly decreases to zero. The parameter $B$ is the analogue of a travelling wave speed. Solutions to (4.15a,b) of this form generalise the ‘grim reaper’ solutions to the zero-surface-tension problem found in Farmer & Howison (2006), and are also candidates for asymptotically valid solutions to the higher-order system (3.19), since as the thickness goes to zero, we would expect the higher-order terms to vanish at a faster rate than the leading-order terms. We will see in § 5 that these quasi-travelling waves do not appear to be attractors, but will compute them here for completeness.

On substitution of the ansatz (4.22a,b) into the leading-order model in Cartesian form (4.15a,b), and scaling the variables according to

(4.23a–c)\begin{equation} H = \frac{P}{B}\hat H,\quad F = B\hat F,\quad x = B\hat x , \end{equation}

we obtain

(4.24)\begin{equation} -\left[\hat F'\left(1 + \hat\sigma \frac{\hat F''}{(1+\hat F'^2)^{3/2}}\right)\right]' = (1+\hat F'^2)\left(1 + \hat\sigma \frac{\hat F''}{(1+\hat F'^2)^{3/2}}\right),\quad \hat\sigma = \frac{2\sigma}{BP}. \end{equation}

This equation may be solved numerically directly, but to further simplify we cast it into the following curvature–angle formulation. Let $\theta = -\tan ^{-1} \hat F'$ be the angle between the tangent to the centreline and the $x$ axis (counting positive for negative $F'$), and scaled curvature $K = B^{-1}\kappa = \hat F''/(1+\hat F'^2)^{3/2}$ (see figure 6a). We then obtain the first-order equation

(4.25)\begin{equation} \frac{\mathrm dK}{\mathrm d\theta} ={-}\frac{1 + \hat\sigma K}{\hat\sigma K \sin\theta} \left(1+\frac{K}{\cos\theta}\right). \end{equation}

For a semi-infinite curve, the appropriate interval is $-{\rm \pi} /2 < \theta < {\rm \pi}/2$, where the nose at $\theta =0$ is a singular point, at which we require $K(0) = -1$. In the limit of the tail of the travelling wave $\theta \to {\rm \pi}/2$, we must have $K \sim \theta -{\rm \pi} /2 \to 0$.

Figure 6. (a) Depiction of the curvature–angle coordinate system used to solve the quasi-travelling wave solutions (4.25). (b) Quasi-travelling wave solutions to (4.25), for $\epsilon =0$ and $\epsilon \to \infty$ (dashed curves), and $\hat \sigma \in \{0.1, 0.5, 1, 10\}$ (solid curves). The arrows indicate direction of increasing $\hat \sigma$.

We solve (4.25) numerically for different values of $\hat \sigma$, and reconstruct the $\hat x, \hat F$ coordinates, using

(4.26a,b)\begin{equation} \frac{\mathrm d \hat x}{\mathrm d\theta} ={-}\frac{\cos\theta}{K},\quad \frac{\mathrm d\hat F}{\mathrm d\theta} = \frac{\sin\theta}{K}. \end{equation}

The $\hat \sigma =0$ and $\hat \sigma \to \infty$ limits are amenable to exact solutions. If $\hat \sigma =0$ then $K = K_0 = -\cos (\theta )$, which is the ‘grim reaper’ solution found by Farmer & Howison (2006) (also relevant for the zero-surface-tension Saffman–Taylor finger (Saffman & Taylor 1958) and for curve shortening flow (Angenent 1991)), equivalent to a centreline profile $\hat F = \ln (\cos \hat x)$. On the other hand, as $\hat \sigma \to \infty$ the equation for $K$ becomes linear, and has exact solution

(4.27)\begin{equation} K = K_\infty ={-}\cot\theta\log\left(\frac{\cos(\theta/2) + \sin(\theta/2)}{\cos(\theta/2) - \sin(\theta/2)}\right). \end{equation}

Plots of these quasi-travelling wave solutions are shown in figure 6(b). The effect of $\hat \sigma$, and thus of surface tension, is weak, as all solutions are bounded between the $\hat \sigma =0$ and $\hat \sigma \to \infty$ limits. The scaled curvature at the nose is required to be $-1$ in all cases, and thus returning to the unscaled system, we expect the curvature at the nose to be $-1/B$, that is, inversely proportional to the wave speed parameter.

The corresponding scaled thickness in the $y$ direction is $\bar H = \sec ^2\theta (1+\hat \sigma \kappa )$. The time-dependent thickness in the centreline-normal direction is thus

(4.28)\begin{equation} h = \frac{P(t_0-t)}{B\cos\theta} (1+\hat\sigma \kappa). \end{equation}

At a fixed time $t_0$ then, the thickness is unbounded in the tails as $\theta \to \pm {\rm \pi}/2$. If, however, we consider a fixed value of $y$ (in the non-travelling coordinate system) as the singular time $t_0$ is approached, then from (4.22a,b) and (4.26a,b), the corresponding value of $\theta$ behaves as

(4.29)\begin{equation} {\rm \pi}/2 - \theta \sim (t_0-t)\,\mathrm e^{y/B},\quad t\to t_0. \end{equation}

Thus $h = O(\mathrm e^{-y/B})$ is bounded for a fixed value of $y$ as $t\to t_0$. The exponential decay in $y$ of the thickness also means that a quasi-travelling wave solution does not require infinite mass to form, even as its length is unbounded in the $y$ direction as $t\to t_0$.

In § 5 we simulate the second-order thin-filament model, including with initial conditions close to a quasi-travelling wave solution. These simulations demonstrate that these quasi-travelling wave solutions do not appear to be stable.

5. Numerical simulation and comparison

In this section we describe numerical simulations of the second-order thin-filament model (3.19), and compare against simulations of the full two-interface problem (2.1).

Our solution method for the thin-filament model (3.19) is a front tracking method, whereby the centreline is represented by $N$ points $\boldsymbol x_j = (x_j, y_j),\ j=1,\ldots, N$. The thickness also has a value $h_j$ at each point. At a given time, the normal vector, curvature and arclength derivatives of $h$ and related quantities in (3.19) are calculated using a central finite-difference scheme. The points are then moved in time in the normal direction with velocity given by (3.19a) (consistent with (3.19b) being the Lagrangian evolution equation for $h$) using MATLAB's ode15s. Since moving points with normal velocity results in highly unevenly spaced points on the centreline, we remesh onto an evenly spaced grid when the ratio between minimum and maximum node spacing drops below a threshold value. We typically use 2000–5000 points (increasing this resolution does not have an appreciable effect in the simulations reported here) and remesh when the minimum-to-maximum node spacing is less than $0.8$.

5.1. Validation against solution of the two-interface problem

To validate solutions of the thin-filament model (3.19) it is valuable to compare it with solutions of the full two-interface system (2.1), which is a more challenging numerical problem. We solve (2.1) using the numerical framework proposed by Morrow et al. (2021, 2023), which we briefly summarise here. The framework is based on the level-set method, in which we represent each interface, $f_{L}$ and $f_{U}$, as the zero level set of the associated level-set functions $\psi _{{L}}$ and $\psi _{{U}}$. Each level-set function is chosen such that the viscous fluid will occupy the region where both $\psi _{{L}}$ and $\psi _{{U}} > 0$; otherwise the region is filled with inviscid fluid. Both level-set functions are updated according to

(5.1a,b)\begin{equation} \frac{\partial \psi_{{L}}}{\partial t} + f_{L} | \boldsymbol{\nabla} \psi_{{L}}| = 0 \quad \textrm{and} \quad \frac{\partial \psi_{{U}}}{\partial t} + f_{U} | \boldsymbol{\nabla} \psi_{{U}}| = 0, \end{equation}

where

(5.2a,b)\begin{equation} f_{L} = \boldsymbol{\nabla} \phi \boldsymbol{\cdot} \boldsymbol{n}_{L} \quad \textrm{and} \quad f_{U} = \boldsymbol{\nabla} \phi \boldsymbol{\cdot} \boldsymbol{n}_{U}, \end{equation}

and $\boldsymbol {n}_{L} = \boldsymbol {\nabla } \psi _{{L}} / |\boldsymbol {\nabla } \psi _{{L}}|$ and $\boldsymbol {n}_{U} = \boldsymbol {\nabla } \psi _{{U}} / |\boldsymbol {\nabla } \psi _{{U}}|$ are the unit (outward) normals. We discretise the spatial derivatives in (5.1a,b) using a second-order essentially non-oscillatory scheme and integrate in time using second-order total-variation-diminishing Runge–Kutta with time step $\Delta t = 10^{-5}$. We perform simulations on the computational domain $-{\rm \pi} \le x \le {\rm \pi}$ and $-0.5 \le y \le 4$, which is discretised into equally spaced nodes with mesh size $\Delta x = 2 {\rm \pi}/ 750$. Simulations are concluded when the minimum distance between the two interfaces is less than $3 \Delta x$. Further, we periodically perform reinitialisation to maintain both $\psi _{{L}}$ and $\psi _{{U}}$ as signed distance functions. The details of this reinitialisation procedure are described in Morrow et al. (2021).

We solve (2.1a) for the velocity potential $\phi$ via a finite-difference stencil. Following Chen et al. (1997), we modify the stencil at nodes adjacent to either interface, corresponding to where $\psi _{{L}}$ or $\psi _{{U}}$ changes sign, by imposing a ghost node on the interface to incorporate the appropriate dynamic boundary conditions (2.1c) and (2.1d). Here, $\kappa _{L} = \boldsymbol {\nabla } \boldsymbol{\cdot } \boldsymbol {n}_{L}$ and $\kappa _{U} = \boldsymbol {\nabla } \boldsymbol{\cdot } \boldsymbol {n}_{U}$. By solving for $\phi$, we can compute $f_{L}$ and $f_{U}$ from (5.2a,b). These choices of $f_{L}$ and $f_{U}$ satisfy the kinematic boundary conditions (2.1b) and give a continuous expression for $f_{L}$ and $f_{U}$ in the region occupied by the viscous fluid $\boldsymbol {x} \in \varOmega$. However, to solve (5.1a,b), we require expressions for $f_{L}$ and $f_{U}$ over the entire computational domain. To extend our expressions for $f_{L}$ and $f_{U}$ into the gas regions, we follow Moroney et al. (2017) by solving the biharmonic equation

(5.3a,b)\begin{equation} \nabla^4 f_{L} = 0 \quad \textrm{and} \quad \nabla^4 f_{U} = 0 \quad \textrm{in}\ \boldsymbol{x} \in \mathbb{R}^2 \backslash \varOmega. \end{equation}

By doing so, we obtain smooth continuous normal velocities over the entire computational domain, allowing us to solve (5.1a,b) for $\psi _{{L}}$ and $\psi _{{U}}$.

In figure 7 we compare the results of the filament model and full two-interface model, with initial conditions

(5.4a,b)\begin{equation} f_{U}(x, 0) = 1/24 - 0.00375 {\rm sech}^2 x,\quad f_{L}(x, 0) ={-}1/24 + 0.00375 {\rm sech}^2 x. \end{equation}

These initial conditions correspond to an initial centreline location of $y=0$ and an initial thickness $h(x,0) = 1/12 - 0.0075\textrm {sech}^2 x$, which is an almost flat filament with a small thinner region near the centre of the channel, $x=0$. As demonstrated in figure 7(a), the agreement between the two methods is initially very good, and only start to disagree quantitatively when the filament thins near the central region, leading to a large increase in velocity. This is to be expected, as the level-set method will become inaccurate when the filament thickness becomes too thin to capture accurately using a level-set function on a discretised mesh. To demonstrate this point, in figure 7(b) we plot the minimum thickness $h_{min}(t)$ against time, for level-set computations of increasing resolution. We see that the level-set result does approach that of the filament model as the resolution is increased (of course, we do not expect the results to be identical due to the approximation made in deriving the thin-filament model (3.19), but this error is much smaller than the achievable discretisation error in the level-set simulation). Ultimately in the level-set simulations, the thickness saturates at a value dependent on the resolution. This limitation means that the level-set simulations ultimately underestimate the speed at which the filament advances when it becomes very thin, as can be observed in the profiles shown in figure 7(a).

Figure 7. (a) A comparison between the interfaces predicted from numerical solution of the filament model and those from numerical solution of the full problem using the level-set method described in § 5.1. The initial condition is given by (5.4a,b) with pressure and surface tension parameters $P=1$ and $\sigma = 0.1$, respectively. For clarity, only the interfaces of the filament model, found from the actual variables of centreline and thickness, are plotted. (b) A comparison of the minimum thickness $h_{min}$ (which occurs at the nose $x=0$) between the filament model and level-set simulations, as the resolution of the level-set simulation is increased. For smaller $\Delta x$, the level-set method tends to the thin-filament result, while the filament thickness in the level-set method saturates at a value that depends on the numerical resolution. The points marked on the filament curve correspond to the times at which the profiles are plotted in (a).

5.2. Numerical results for late times

Here we use the front-tracking scheme to solve the thin-filament model for later times, in the regime where the thickness becomes too small for the level-set method to accurately resolve. We use the initial condition

(5.5a,b)\begin{equation} y = 0,\quad h = 0.2[1-0.1\cos(x)], \end{equation}

which is the same as the initial condition of the exact solution (4.9) depicted in figure 4(b). In order to observe the effect of different surface tensions $\sigma$, we run simulations for $\sigma = 0.1$ (figure 8), as well as $\sigma = 0.05$ and $\sigma = 0.2$ (figure 9).

Figure 8. Evolution of a filament with initial position $y = 0$ and thickness $h = 0.2(1-0.1\cos (x))$, with $P=1$ and surface tension $\sigma = 0.1$. (a) The centreline profiles over time (solid lines) and the thickness (dashed lines), and (b) the thickness against $x$, showing the initially thinner part of the filament bulge outward into a ‘bubble’, while the bulk of the fluid is driven out to the edges of the filament. Ultimately the profile intersects with the channel boundary at $x=\pm {\rm \pi}$. The (c) maximum normal velocity $v_{n\,max}(t)$ and (d) minimum thickness $h_{min}(t)$ appear to become unbounded and go to zero, respectively. (e) The curvature at the nose ($x=0$) over time, initially increases in magnitude, then rapidly heads toward zero as the bubble expands. (f) The curvature scaled by the curvature at the nose, showing the curvature tending to a constant over the bubble. The circles marked in (c–e) correspond to the times at which profiles are plotted in (a,b,f), and the arrows in (a,b,f) indicate the direction of increasing time.

Figure 9. Evolution of a filament with initial position $y = 0$ and thickness $h = 0.2(1-0.1\cos (x))$, for (a,b) $\sigma = 0.05$ and (c,d) $\sigma = 0.2$. The formation of a circular bubble occurs in a similar fashion to what happens for $\sigma = 0.1$, with the major effect being the width of the neck region that precedes the circular bubble.

In the $\sigma =0.1$ simulation (figure 8), the filament bulges outward in the centre near $x=0$, where the thin filament is initially thinner (and so the filament moves faster). This bulge becomes a ‘bubble’ that rapidly expands in radius while the thickness rapidly decreases. The majority of the fluid is pushed into the outer regions of the channel. At a finite time, the bubble intersects the channel walls at $x = \pm {\rm \pi}$. Unlike a fluid-filled Hele-Shaw cell, there is nothing to prevent the filament reaching the channel walls; this does not correspond to singular behaviour in the mathematical model, but physically represents an area of gas of lower pressure being trapped by the filament.

To further understand this behaviour, we note that (3.19) has as a solution a perfectly circular bubble of radius $R(t)$, which evolves according to

(5.6a,b)\begin{equation} R(t) = \frac{2\sigma}{P} + \left(R(0) - \frac{2\sigma}{P}\right) \exp({(P/c)t}),\quad h = \frac{c}{R}. \end{equation}

Here $c$ is a constant that depends on the initial thickness. In such a solution, the radius (if initially greater than $2\sigma /P$) grows exponentially, but does not exhibit a finite-time singularity. It may be that this is the behaviour that is governing late stages of evolution of the filament depicted in this section. In figure 8(e) we plot the curvature $\kappa _{nose}$ at the nose, or front, of the bubble (where $x=0$), while in figure 8(f) we plot $\kappa /\kappa _{nose}$. The curvature initially grows in magnitude (in our convention the curvature is negative in the bubble) when the bulge initially grows, but then rapidly decreases in magnitude at the time when the bubble expands outwards. When this happens, the curvature tends to a uniform state in the bubble region ($\kappa /\kappa _{nose} \to 1$), implying the convergence to a circular shape. While we have only plotted the interface up to the time at which it intersects the channel wall, mathematically the simulation continues for longer, and the bubble becomes more circular in shape. Ultimately, the thickness becomes so small and the velocity so large that the numerical method no longer converges.

This behaviour is also seen for smaller and larger surface tension values $\sigma = 0.05$ and $0.2$, as depicted in figure 9. The notable effect of different surface tension is that the bulging region that forms the ‘neck’ of the bubble is smaller or larger in width as the surface tension is smaller or larger.

5.3. Initial condition near a quasi-travelling wave

In addition to the previous near-flat initial conditions, we demonstrate the instability of the quasi-travelling wave solutions considered in § 4.3 by choosing an initial condition close to one such solution. For simplicity we use the ‘grim reaper’ exact solution for zero surface tension $K_0 = -\cos \theta$, and a scale factor $B=1$, for which

(5.7a,b)\begin{equation} y = \log(\cos x),\quad h = \frac{b}{\cos x} \end{equation}

(the exact solution is a good approximation for small $\sigma$, and although we do not show the results here, we have observed the same behaviour from numerically constructed initial conditions by solving (4.25) for $\hat \sigma > 0$). In terms of the quasi-travelling wave solution, the factor $b$ in the thickness $h$ is arbitrary as it corresponds to translation in time. As the travelling wave is semi-infinite, we must of course approximate it by a finite curve that respects the required boundary conditions at $x=\pm {\rm \pi}$. We thus choose the closely related ‘hairclip’ curve, defined by

(5.8a,b)\begin{equation} y = \sinh^{{-}1}(\alpha\cos x),\quad h = \frac{b}{\cos\theta} = b \sqrt{1 + \left(\frac{\alpha\sin x}{\cosh y}\right)^2},\quad \alpha \gg 1. \end{equation}

The larger the parameter $\alpha$ is made, the more elongated the initial finger becomes. We note that the thickness will become large on the sides of the initial finger, representing the blow-up in thickness of the quasi-travelling wave in the limit $\theta \to {\rm \pi}/2$ discussed in § 4.3.

We depict the result of such an initial condition in figure 10, for the values $\sigma =0.05$, $\alpha =20$, $b=0.1$. The regions with large thickness on the sides of the finger correspond to the exponential growth in thickness required for the travelling wave solution as $y\to -\infty$ (see § 4.3); the filament in this region does not evolve to a great extent over the simulation. At the front, while the finger does initially propagate forwards, the nose bulges outward and becomes more circular in shape, leading to the same late-time behaviour as seen for the more general initial conditions.

Figure 10. The evolution of a filament that starts near the $\sigma =0$ limit of the quasi-travelling wave described in § 4.3 (the precise initial condition is (5.8a,b)). (a) The quasi-travelling wave appears unstable, with the front evolving toward the expanding circle seen for general initial conditions. (b) The curvature at the nose of the filament ($x=0$) starts at $-1$ but decreases over time as the filament expands. The circles marked in (b) correspond to the profiles plotted in (a).

6. Discussion

In this paper we have developed a simplified but highly accurate second-order lubrication flow model (3.19) that describes two-interface Hele-Shaw flow very well in regions where the thickness of the fluid region becomes small. Due to the instability of one of the interfaces, this limit is one that is generally reached, even if initially the thickness is not very small. By examining the generic singular behaviour of the leading-order model, even in the presence of surface tension (§ 4), we have shown why a second-order model is necessary to represent the dynamics of the full problem. In particular, the fourth-order spatial derivative term arising from the difference in curvature between the upper and lower interfaces is necessary to regularise the system. Although we have included all terms formally of order $\epsilon ^2$ in our model and simulations, we have also observed that the leading-order model along with the addition of only this regularising term (proportional to $[hh_{sss}]_s$) will behave in a qualitatively similar manner (with small quantitative differences). We have not shown these results here for the sake of brevity.

Here we note some clear differences between the instability of a thin filament, and the classical Saffman–Taylor instability in a semi-infinite fluid region that results inthe Saffman–Taylor finger with (in the small-surface-tension limit) width half that of the channel. One difference is that the thin-filament model does not feel the effects of the channel wall away from the thick neck regions, and thus the rapidly expanding bubble may intersect the channel walls at finite time with no breakdown in the mathematical model. Physically speaking, this phenomenon would correspond to trapping a part of the lower-pressure gas inside the fluid region. In addition, as the walls have no strong effect, the orientation of the filament motion to be mainly in the positive $y$ direction is a somewhat artificial consequence of the initial condition. For this reason it may be more natural to consider the fluid in an unbounded Hele-Shaw cell, with the length scale set by the initial thickness.

The specifics of the late stages of evolution of the filament (depicted in § 5), wherein the thickness becomes very small and the velocity correspondingly large, are not resolved. In order to understand the dynamics at late times of this system, an analysis of the stability of the quasi-travelling wave solutions, and the stability of the expanding circular bubble (5.6a,b) to non-radially symmetric perturbations, would be very valuable. We observe that in our model, the filament does not appear to exhibit finite-time ‘bursting’ behaviour, that is, the thickness does not go to zero at an isolated finite point in space and time. In the case of self-similar breakup in the manner of the unforced lubrication equation described in Almgren et al. (1996), the thickness goes to zero at a point where the curvature becomes infinite, while in our solutions depicted in § 5, the thickness becomes small in the expanding circular region where the curvature is also decreasing, while the curvature is largest in magnitude in the neck regions, where the thickness does not decrease. Bursting in physical systems is likely to require fully three-dimensional effects to explain (that is, when the filament thickness becomes of the same order as the separation between plates in the Hele-Shaw apparatus). Once these two length scales are comparable, the filament model and indeed the two-dimensional Hele-Shaw model are no longer valid. The dynamics of the filament model studied in this work indicates that such a regime will generically be approached very rapidly in the form of the expanding bubble depicted in § 5.