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. Reference Morrow, Moroney, Dallaston and McCue2021)

(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. Reference Morrow, De Cock and McCue2023). 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 Reference Gin and Daripa2015; Anjos & Li Reference Anjos and Li2020).

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.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. Reference Van De Fliert, Howell and Ockendon1995; Howell Reference Howell2003). 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. Reference Constantin, Dupont, Goldstein, Kadanoff, Shelley and Zhou1993; Dupont et al. Reference Dupont, Goldstein, Kadanoff and Zhou1993; Goldstein et al. Reference Goldstein, Pesci and Shelley1993, Reference Goldstein, Pesci and Shelley1995; Almgren Reference Almgren1996; Almgren et al. Reference Almgren, Bertozzi and Brenner1996; Goldstein et al. Reference Goldstein, Pesci and Shelley1998). 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.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 (Reference Farmer and Howison2006). 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 (Reference Farmer and Howison2006) 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 (Reference Farmer and Howison2006) 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 (Reference Farmer and Howison2006)). 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 (Reference Farmer and Howison2006), 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 (Reference Farmer and Howison2006) (also relevant for the zero-surface-tension Saffman–Taylor finger (Saffman & Taylor Reference Saffman and Taylor1958) and for curve shortening flow (Angenent Reference Angenent1991)), 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.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. (Reference Morrow, Moroney, Dallaston and McCue2021, Reference Morrow, De Cock and McCue2023), 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. (Reference Morrow, Moroney, Dallaston and McCue2021).

We solve (2.1a) for the velocity potential $\phi$ via a finite-difference stencil. Following Chen et al. (Reference Chen, Merriman, Osher and Smereka1997), 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. (Reference Moroney, Lusmore, McCue and McElwain2017) 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).