3.1. Time evolution of a thin layer

Figure 2(a) shows snapshots from a numerical solution of the thin-film equations (2.23)–(2.29) with $B=0.05$ and $A=0.2$. For any viscoplastic simulation, the initial perturbation must be sufficiently large in order to trigger instability because a sufficiently large pressure gradient must be created to overcome the yield stress and make the fluid yield. Here, $A$ is large enough that the fluid in the centre of the domain yields, but $Y=0$ near the boundaries so the fluid there is initially rigid (figure 2a, $t=0$). There is an initial period during which there is minimal deformation in $H$, but $Y$ deforms significantly and, by $t=20$, $Y>0$ for all $z\in (0,L)$ indicating that the whole layer is yielded. This initial yielding period causes a delay in the growth of the instability, as can be seen when $\max _zH(z,t)$ is compared with the same quantity from a Newtonian ($B=0$) simulation in figure 2(b). After the initial yielding period, there is a period of significant deformation of the free surface. This coincides with a peak in $\max _zY(z,t)$ and a peak in the wall shear stress, $|\tau _w|$, at around $t=100$. This indicates that there is significant shear during this period, with more of the layer exhibiting shear-dominated flow and the region of plug-like flow becoming smaller.

At late times, the layer relaxes slowly towards a final marginally yielded static shape in which $Y\rightarrow 0$ across the whole layer (figure 2(a), $t=600$). Figure 2(b) shows that $\max _zY$ decays towards zero at a rate proportional to $t^{-1}$. From (2.30), we note that as $Y\rightarrow 0$, $|\tau _w|\rightarrow B$ across the entire layer, which indicates that even when the layer reaches its final static shape, viscous and capillary effects apply a uniform stress on the tube wall equal to the yield stress. Figure 2(b) also indicates that $\max _z|\tau _w|$ decays towards $B$ at a rate proportional to $t^{-1}$, in contrast to the Newtonian result that the peak wall shear stress decays towards zero at a rate proportional to $t^{-1/4}$ (Jones & Wilson Reference Jones and Wilson1978; Hammond Reference Hammond1983). The late time shape of the layer consists of a large collar of fluid around $z=L$ and a small collar around $z=0$. Figure 2(b) shows that the final value of $\max _zH$ is lower for the $B=0.05$ solution compared with the Newtonian solution. This is because, unlike in the Newtonian evolution, not all of the fluid drains into the large collar at late times. Instead, some is trapped in the small collar so the peak height of the layer is decreased. Thus, the yield stress inhibits the growth of the instability. In §§ 3.2 and 3.3 we quantify how increasing $B$ affects the final marginally yielded static shape of the layer.

In the early time period of gradual yielding, capillary-wave-like structures form ahead of the travelling yield surfaces (where $Y$ makes contact with zero). They can be identified by observing the structure of $Y$ and the pressure gradient $p_z$ (figure 2c). There is a jump discontinuity and a change of sign in $p_z$ between each of these structures, indicating that the direction of flow reverses. At each point that $p_z$ passes through zero, we also have $Y=0$. We expect there to exist additional waves with smaller wavelengths ahead of those observed in figure 2(c), but the numerical method can only resolve the largest few since it is limited by the size of the finite difference grid spacing. The observed structures resemble closely the capillary waves identified by Balmforth et al. (Reference Balmforth, Ghadge and Myers2007b) and Jalaal & Balmforth (Reference Jalaal and Balmforth2016), and studied in detail by Jalaal et al. (Reference Jalaal, Stoeber and Balmforth2021) in the context of spreading viscoplastic droplets. They are a feature common to surface-tension-driven viscoplastic flows, occurring when a yield surface advances into a region of unyielded fluid. Unlike in previously studied flows, here the capillary waves, in general, only exist transiently, during the early time period of gradual yielding until the whole layer yields. However, we find that for some values of $A$ and $B$ (mostly very large $A$), one or more of the jump discontinuities in $p_z$ which develop can persist for the whole evolution. In § 3.4 we discuss how this phenomenon can affect the final static shape of the layer, and present criteria for it to occur. Until then, we focus on the case that the flow is unidirectional ($Y>0$ for all $z\in (0,L)$) after the early time capillary waves have passed.

3.2. Late time asymptotics for the thin-film evolution equation

To analyse the late time dynamics of the layer, we look for a solution in which $Y=\mathcal {Y}\rightarrow 0$ as $t\rightarrow \infty$. Numerical simulations suggest that $Y$ decays like $t^{-1}$ (figure 2b), so we make the expansions

(3.1a,b)\begin{equation} H = H_0(z;B) + \frac{H_1(z;B)}{Bt}+\cdots, \quad Y = \mathcal{Y} = \frac{Y_1(z;B)}{Bt}+\cdots, \mbox{as}\ t\rightarrow\infty.\end{equation}

The $t^{-1}$ rate of decay of $H$ towards a steady state is also consistent with the numerical results in figure 2. For this analysis, we assume that the capillary pressure is monotonic, or equivalently, the pressure gradient is one-signed, $-H_{0,z}-H_{0,zzz}<0$ for all $z\in [0,L]$. This is equivalent to assuming unidirectional flow at late times. Substituting the expansions (3.1a,b) into the definition of $\mathcal {Y}$ (2.25) and the evolution equation (2.27) gives

(3.2a)$$\begin{gather} H_0(H_{0,z}+H_{0,zzz}) = B, \end{gather}$$

(3.2b)$$\begin{gather}B(H_1-Y_1) + H_0^2(H_{1,z}+H_{1,zzz}) = 0, \end{gather}$$

(3.2c)$$\begin{gather}H_1 = \tfrac{1}{2}(Y_1^2)_z. \end{gather}$$

The boundary conditions (2.28) imply that

(3.3)\begin{equation} H_{0,z}=H_{1,z}=Y_1=0\ \mbox{at}\ z=\{0,L\},\end{equation}

and mass conservation implies that

(3.4a,b)\begin{equation} \int_0^LH_0\,\mathrm{d}z=\frac{L}{\rm \pi}, \quad \int_0^LH_1\,\mathrm{d}z = 0.\end{equation}

From (3.2a), note that $H=H_0(z;B)$ is a static solution of the evolution equation (2.27) in which $\mathcal {Y}=0$ uniformly. This is not a capillary-static solution in which the pressure is everywhere uniform; instead it is a state in which the layer is uniformly marginally yielded. From (2.30) we note that, in the static solution, $\tau _w=B$ uniformly, indicating that there is a stress being applied in the positive $z$-direction, but it is resisted by the yield stress, preventing flow. The functions $H_1(z;B)$ and $Y_1(z;B)$ quantify the rate at which the layer approaches the static solution at late times.

The equations (3.2a) and (3.2c) were solved subject to (3.3) and (3.4a,b) using a boundary value problem solver in Matlab. A solution for $B=0.05$ is shown in figure 3 with comparison to the final snapshot of the numerical simulation from § 3.1. Figure 3(a) shows the agreement is very good between $H_0$ and the late time shape of the layer from the numerical solution. Figures 3(b) and 3(c) show that $H_1$ and $Y_1$ approximate well the rates of decay of $H$ towards $H_0$, and $Y$ towards zero, respectively. This indicates that the expansion (3.1a,b) accurately describes the late time dynamics of the evolution and confirms the $O(t^{-1})$ rate of decay in $Y$ determined empirically in figure 2(b).

Figure 3.

Late time asymptotic solutions for $B=0.05$, compared with the final snapshot of the numerical simulation with $A=0.2$ and $B=0.05$ at $t=10^4$. (a) The static solution, $H_0$, compared with the layer height, $H(z,t=10^4)$, from the simulation. (b) Function $H_1$ compared with $Bt[H(z,t=10^4)-H_0]$ from the simulation, which represents the rate of decay of $H$ towards $H_0$. (c) Function $Y_1$ compared with $BtY(z,t=10^4)$ from the simulation, which represents the rate of decay of $Y$ towards zero.

3.3. Static solutions

The marginally yielded static solutions, $H_0(z;B)$, can predict the final state of the layer. To investigate the dependence of $H_0$ on $B$, we solve (3.2a) with (3.3) and (3.4a), varying $B$. We find that there exists a value $B_*\approx 0.163$ such that, for all $0\leq B< B_*$, exactly two solutions for $H_0$ exist and, for $B>B_*$, no solutions exist. There is a bifurcation at $B=B_*$. Figure 4(a) shows $\max _zH_0(z;B)$ for all the solutions, which always coincides with $H_0(z=L;B)$. Figure 4(b) shows several example solutions, some lying on the upper branch in figure 4(a) and some lying on the lower branch.

Figure 4.

Solutions $H_0(z;B)$ of (3.2a), (3.3) and (3.4a), which are static solutions of the evolution equation (2.27). (a) Bifurcation diagram showing $\max _zH_0$ for all values of $B$ such that solutions exist. Dotted black lines are the asymptotic approximations (3.6), (3.8) and (3.9). The arrows indicate time evolution as determined by the near-bifurcation asymptotic analysis (Appendix C). The five dots on (a) correspond to the example solutions shown in (b), at $B=0.05$, $B=0.12$ and $B=B_{*}$.

The upper-branch solutions are significantly deformed layers with a large collar around $z=L$ and a small collar around $z=0$. These shapes are approached by the evolving layer at late times (figure 3). The upper branch in figure 4(a) quantifies the decrease in the size of the large collar formed by the layer as $B$ is increased. This decrease can be significant: $\max _zH_0$ for the upper branch ranges from the Newtonian value, ${\max _zH_0(z;0)=2\sqrt {2}\approx 2.83}$, down to $\max _zH_0(z;B_*)\approx 1.98$. This indicates that increased yield stress can significantly inhibit deformation of the film.

The lower-branch solutions are near flat for small $B$, becoming more deformed as $B$ is increased. In numerical simulations an unstable layer evolves away from a near-flat configuration towards a strongly deformed (upper-branch) static shape. This suggests that the upper-branch solutions are stable and the lower-branch solutions are unstable. We confirm the stability of the two branches using asymptotic analysis near to the bifurcation, $B\approx B_*$, presented in full in Appendix C. We find that $H(z,t)\sim H_0(z;B_*)+\mu \mathcal {A}(T)\phi _1(z)$ as $\mu \equiv \sqrt {B_*-B}\rightarrow 0$, where $\phi _1(z)$ is a solution to a linear ODE, $T=\mu ^3t$ is a slow time scale, and $\mathcal {A}(T)$ is an amplitude function which solves an ODE of the form

(3.5)\begin{equation} \mathcal{A}_T = C_0(\mathcal{A}^2-\mathcal{A}_0^2)^2,\end{equation}

where $C_0$ and $\mathcal {A}_0$ are positive constants. Equation (3.5) has two fixed points, $\mathcal {A}=\pm \mathcal {A}_0$, and we compute $\mathcal {A}_0\approx 2.20$. Solutions evolve away from the negative fixed point towards the positive one (figure 10b below), indicating that the positive fixed point is stable and the negative one unstable. These fixed points correspond to two static solutions for $H$ which are the upper- and lower-branch solutions, respectively. With this stability result, we identify the bifurcation at $B=B_*$ as a saddle-node bifurcation. We also approximate the location of the branches in figure 4(a) by

(3.6)\begin{equation} \max_zH_0\sim H_0(L,B_*)\pm\mu\mathcal{A}_0\phi(L)\approx 1.98\pm2.20\mu \ \mbox{as} \ \mu\equiv\sqrt{B_*-B}\rightarrow0.\end{equation}

Since the lower-branch solutions satisfy $\mathcal {Y}=0$, they are marginal states between rigid layers ($\mathcal {Y}\leq 0$) and fully yielded layers ($\mathcal {Y}>0$). We expect that if a layer is initially more deformed than the lower-branch solution it will be yielded and unstable, but if it is less deformed initially it is likely to be rigid and thus stabilised. We provide evidence from numerical simulations to confirm this in § 3.4, where we show that the lower-branch static solutions correspond almost exactly to the minimum amplitude of initial perturbation required to trigger unstable growth.

Asymptotic analysis for small $B$ shows that the lower-branch solutions have the regular expansion,

(3.7)\begin{equation} H_0(z;B) = 1 + B\left[x-\frac{L}{2}-\sin{x}+ \frac{(1-\cos{L})\cos{x}}{\sin{L}}\right]+\cdots \ \mbox{as} \ B\rightarrow0.\end{equation}

Taking the maximum value of (3.7) gives

(3.8)\begin{equation} \max_{z}H_0= 1 + B\left[\frac{L}{2}+\cot{L}-\mathrm{cosec}\,{L}\right] +\cdots\ \mbox{as} \ B\rightarrow0,\end{equation}

which approximates the lower branch in figure 4(a). To approximate the upper-branch solutions for small $B$, we construct a solution by matched asymptotic expansions. The analysis is presented in full in Appendix D and illustrated in figure 11 below. In addition to the large collar around $z=L$ and the small collar around $z=0$, we identify a third, inner region located around $z=L-{\rm \pi}$ where $H_0\sim O(B^2)$. In contrast to the Newtonian problem, where the inner region is described by an ODE of the form $H^3H_{zzz}=1$ (Hammond Reference Hammond1983), here the relevant ODE is of the form $HH_{zzz}=1$ (see (D 5)). The difference arises because in the Newtonian problem there is constant flux across the inner region during the late time draining regime, while here there must be constant stress across the inner region since $H_0$ is marginally yielded. After expanding and solving for $H_0(z,B)$ in each of the three regions and matching the solutions, a composite approximation to $H_0(z,B)$ is found. This also provides an approximation to the maximum value,

(3.9)\begin{equation} \max_zH_0 = \frac{2L}{\rm \pi}+4a_1B^{1/2}+\frac{4-{\rm \pi}^2}{L}B +\cdots\ \mbox{as} \ B\rightarrow0,\end{equation}

where $a_1$ is a constant which depends on $L$. For $L=\sqrt {2}{\rm \pi}$, we compute $a_1\approx -0.105$. Equation (3.9) approximates the upper branch in figure 4(a).

At the saddle-node bifurcation, $B=B_*$, the two static solutions annihilate each other, so it is only possible for the layer to select the upper-branch solution if $B< B_*$. Note that in this section we have only computed the static solutions that have monotonic pressure, but other static solutions may exist. We show in the next section that even when $B< B_*$, the layer may select a different final static shape, depending on the initial conditions of the layer, and that it is possible to have some yielding and unstable growth for some $B>B_*$ if the initial perturbation is sufficiently large.

3.4. Dependence on capillary Bingham number and initial conditions

The final shape of the layer can depend on its initial conditions as well as $B$. To investigate this dependence, we solve the thin-film equations (2.23)–(2.29) numerically for a range of values of $B$ and a range of initial perturbation amplitudes $A$. We run simulations on a regularly spaced grid of points in the range $(0\leq B\leq 0.5,0\leq A\leq 0.99)$. All simulations are run to a fixed, long time, which we choose to be $t=1000$. In figure 5(a) the final maximum height, $\max _z H(z,t=1000)$, is plotted for each simulation.

Figure 5.

(a) Data from numerical solutions of the thin-film evolution equation at various $B$ and $A$. Each dot corresponds to a simulation with the colour indicating $\max H(z,t=1000)$. The data are linearly interpolated to produce the black contour lines, which are evenly spaced. When $A$ is small and $B$ is large, the layer is yield stabilised: no unstable growth occurs. The critical amplitude for any yielding to occur, $A_m(B)$ (dashed red), defined in (3.10), is a strict lower bound on the boundary of the yield-stabilsed region. When there is growth, the final shape either has monotonic or non-monotonic pressure, depending on $A$ and $B$. The quantity $1-H_0(z=0,B)$ (magenta) from the static solutions (figure 4) predicts the boundary of the monotonic pressure region. The maximum value of $B$ along the magenta curve is $B_*$. Two large black dots indicate the location of the example solutions shown in (b) and (c), with $H$, $70Y$ (scaled for clarity) and $p_z$ plotted at $t=1000$. Plot (b) shows a solution in the monotonic pressure region ($p_z\leq 0$) and (c) shows a solution in the non-monotonic pressure region ($p_z$ changes sign once in the domain).

Figure 5(a) shows that there are three qualitatively different possible outcomes for an evolving layer, depending on the values of $B$ and $A$. First, there is a region for small $A$ and large $B$ with $\max _zH(z,t=1000)=1+A$, so the final maximum height is equal to the initial maximum height. In these cases, the initial perturbation does not generate a large enough pressure gradient to make the layer yield, so the yield stress entirely suppresses unstable growth. This region also extends to all of $B>0.5$ for all $0\leq A<1$. We call this the yield-stabilised region. We can seek a bound for the yield-stabilised region using the minimum amplitude, $A=A_m(B)$, such that the whole layer is initially unyielded if and only if $A< A_m$. We identify that $A_m$ is the value of $A$ such that the initial condition (2.29) makes $\mathcal {Y}$ non-negative at exactly one point in the domain, and, thus, we find $A_m$ is given implicitly by

(3.10)\begin{equation} \frac{\sqrt{1+8A_m^2}-1}{(4A_m^2+\sqrt{1+8A_m^2}-1)^{3/2}} = \frac{k|1-k^2|}{2^{5/2}B},\end{equation}

where $k=L/{\rm \pi} =\sqrt {2}/2$. The curve $A=A_m(B)$ provides a strict lower bound on the boundary of the yield-stabilised region (figure 5a). A small number of yield-stabilised simulations have $A>A_m(B)$: in these simulations the fluid initially yields a small amount in the centre of the domain but rigidifies before the sides of the domain yield, so there is no growth in $\max _zH$.

In the second possible outcome, there is unstable growth and the layer evolves towards the upper-branch static solution with monotonic pressure which we computed in § 3.3. This occurs for a large set of $B$ and $A$ as indicated in figure 5(a), with $\max _zH(z,t=1000)$ being independent of $A$ in this region since all simulations approach the same final shape for a given $B$. Figure 5(b) shows the final shape for a simulation with $A=0.4$, $B=0.1$, which is within the monotonic pressure region. The pressure gradient $p_z$ is non-positive meaning the flow is unidirectional in the positive $z$-direction at $t=1000$ when the simulation is stopped. Here $Y$ is positive but close to zero for all $z\in (0,L)$ at $t=1000$, so the layer is almost rigid and $H$ is very close to the upper-branch static solution. Within the monotonic pressure region in figure 5(a), the value of $\max H(z,t=1000)$ decreases as $B$ is increased, consistent with the decrease shown in the upper branch in figure 4(a).

The final possible outcome involves unstable growth of the layer leading to the final static shape having non-monotonic pressure. Figure 5(a) shows that this occurs mainly when $A$ is very large. In these simulations, $\max _zH(z,t=1000)$ depends on both $B$ and $A$ and the final shape is not predicted by the upper-branch static solutions in figure 4(a). An example of a late time shape with non-monotonic pressure is shown in figure 5(c) from a simulation with $A=0.95$, $B=0.1$. Comparing this with figure 5(b), there is less fluid in the collar near $z=0$ and so the collar near $z=L$ is larger, even though $B$ is the same. In figure 5(c) there is exactly one point in the domain where the sign of $p_z$ changes, which corresponds to a point where the direction of flow changes. For most values of $B$ and $A$ in the non-monotonic pressure region, the final shape selected by the layer has exactly one sign change in $p_z$ but, for $A$ very close to $1$, we found that there can be more. The sign changes in $p_z$ develop during the early time yielding period, caused by the presence of capillary waves near the yield surfaces (figure 2c). Most of the sign changes exist only during this early time period, but we see here that one or more can persist and affect the late time dynamics, suggesting that in these cases the early time capillary waves can influence the entire evolution. In general, the location of the sign change(s) in $p_z$ affects the final shape of $H$, and the location of the sign change(s) depends on the initial shape of the layer.

There is a sharp boundary in the numerical data between the yield-stabilised and monotonic pressure regions, across which $\max _zH(z,t=1000)$ jumps significantly. There is also a clear boundary between the monotonic and non-monotonic pressure regions, indicated by where the contours of $\max _zH(z,t=1000)$ begin to curve. Figure 5(a) shows that we can predict the locations of both boundaries using the quantity $1-H_0(z=0;B)$, where $H_0(z;B)$ are the static solutions with monotonic pressure computed in § 3.3.

First, we give an explanation for why $1-H_0(z=0;B)$ coincides with the boundary between the yield-stabilised and monotonic pressure regions. Consider a simulation with $A$ just above the boundary, so the fluid yields just enough to trigger instability (e.g. figure 2). After the initial period of gradual yielding, $\mathcal {Y}$ is very close to zero but positive everywhere. At this point the layer's shape is very close to the lower-branch static solution with the same $B$, which has $\mathcal {Y}=0$ everywhere. Since the fluid at $z=0$ is unyielded for most of the initial period of the evolution in which the layer gradually yields, the height at $z=0$ does not change significantly in this period. Hence, the height at $z=0$ must have initially been very close to $H_0(z=0;B)$, the height of the lower-branch static solution at that point. The initial height at $z=0$ is $H(0,0)=1-A$, so the boundary between the yield-stabilised and monotonic pressure regions is predicted by $A = 1-H_0(z=0;B)$, and the lower-branch static solutions effectively correspond to the minimum amplitude of perturbation required to trigger instability.

Secondly, we give an explanation for why $1-H_0(z=0;B)$ coincides with the boundary between the monotonic and non-monotonic pressure regions in figure 5(a). To do this, we determine a condition for a final solution with monotonic or non-monotonic pressure to be selected during the evolution. If the initial height of the layer at $z=0$ is smaller than the height of the corresponding (i.e. for the same $B$) upper-branch static solution at $z=0$, then there must be fluid flow towards $z=0$ in the negative $z$-direction during the evolution if this solution is to be selected. This would mean that the flow at late times must not be unidirectional, and, hence, the pressure of the final shape must be non-monotonic. (The transient early time capillary waves in these simulations create flow reversal but it is very weak so can be neglected in this argument.) So, if $H(z=0,t=0)< H_0(z=0;B)$, for a given value of $B$, then the layer will select a final shape with non-monotonic pressure. Noting again that $H(0,0)=1-A$, we see that $A>1-H_0(z=0;B)$ is an equivalent condition for the layer to select a final shape with non-monotonic pressure.

The results in figure 5(a) are specific to the sinusoidal form of initial perturbation used in (2.29). However, we also ran simulations with initial conditions of the form $H(z,0)=1+A\tanh {(2z-L)}$ and the results for $\max _z H(z,t=1000)$ were qualitatively, and largely quantitatively, the same. The quantity $1-H(z=0;B)$ was still found to accurately bound the monotonic pressure region.

In this section we have illustrated the complex dependence of the evolution of a thin layer on $B$ and $A$, and shown how the static solutions with monotonic curvature can provide insight into the dynamics of the layer. However, thin-film theory cannot capture the full range of possible dynamics for the system because the volume of fluid in a thin layer is too small to form a liquid plug in the tube. We now address this by using long-wave theory to model layers with finite thickness.

4.1. Time evolution of a layer with finite thickness

Figure 6(a) shows snapshots from a numerical solution of the long-wave equations (2.13)–(2.19) with film thickness $\epsilon =0.125$, capillary Bingham number $B=0.05$ and initial perturbation amplitude $A=0.2$. For ease of comparison with the thin-film results, we describe solutions in terms of thin-film parameters and variables: instead of $\hat {B}$ we use $B=\hat {B}/\epsilon ^2$, instead of $\hat {t}$ we use $t=\epsilon ^3\hat {t}$, instead of $R(z,\hat {t})$ we use $H(z,t)$, instead of $\varPsi (z,\hat {t})$ we use $\hat {Y}(z,t)\equiv (1-\varPsi )/\epsilon$ and instead of $\hat {\tau }_w$ we use $\tilde {\tau }_w\equiv \hat {\tau }_w/\epsilon ^2$.

Figure 6.

(a) Snapshots of a numerical simulation with $\epsilon =0.125$, $B=0.05$ and $A=0.2$ (supplementary movie 2 is available at https://doi.org/10.1017/jfm.2022.479 for the full evolution of $H$ and $\hat {Y}$). Axial velocity, $\hat {w}$, as defined in (2.15), is also shown. (b) Time evolution of $\max _zH$, $\max _z\hat {Y}$ and $\max _z|\tilde {\tau }_w|$ for the same simulation, with plugging occurring at $t_p\approx 140$. The simulation is stopped when $\max _zH=0.7/\epsilon$ as it is clear that a plug will form. The evolution of $\max _zH$ and $\max _z|\tilde {\tau }_w|$ for a Newtonian ($B=0$) simulation with the same $\epsilon$ and $A$ shows plug formation occurring earlier, $t_p\approx 70$. (c) Snapshots of the $B=0.05$ simulation at evenly spaced time points between $t=50$ (darkest lines) and $t=130$ (lightest lines). Inset shows $\hat {Y}$ becoming equal to zero across progressively more of the domain, indicating that the small collar of fluid around $z=0$ is rigidifying. There is a small jump in $\hat {p}_z$ where $\hat {Y}$ makes contact with zero, but still $\hat{p}_z\leq 0$ everywhere.

The early time behaviour is qualitatively the same as in the thin-film simulations. There is a delay to the growth as the fluid gradually yields (figure 6b) and capillary waves develop, which are qualitatively the same as in the thin-film case. There is then a peak in $\max _z\hat {Y}$ at around $t\approx 30$, coinciding with significant deformation of the layer. Figure 6(b) also shows that there is a small associated peak in the wall shear stress, $\max _z|{\tilde \tau }_w|$, which occurs slightly later, around $t\approx 40$. The layer then evolves towards a shape with a large collar near $z=L$ and a smaller collar near $z=0$ (figure 6a, $t=70$).

Gauglitz & Radke (Reference Gauglitz and Radke1988) identified the critical thickness required to form a plug in their Newtonian simulations as $\epsilon =0.12$. Since $\epsilon =0.125$ is larger than this critical value, we expect plug formation may be possible in this simulation. Indeed, figure 6(b) shows that at around $t=140$, $\max _zH$ begins to rapidly increase towards the centre of the tube, which is located at $1/\epsilon =8$. The long-wave theory cannot model coalescence so, following the approach of, e.g. Halpern et al. (Reference Halpern, Fujioka and Grotberg2010), we stop simulations when $\max _zH=0.7/\epsilon$, but when we do run the simulation further, $\max _zH$ rapidly approaches $1/\epsilon$, so it is clear that a plug will form. We denote the time taken to form a plug as $t_p$, and use the time at which we stop the simulation as a proxy for $t_p$. Figure 6(b) shows that the plugging time, $t_p\approx 140$, is significantly longer than the plugging time for a Newtonian ($B=0$) simulation, $t_p\approx 70$. This is partly due to the delay caused by the initial yielding period, and partly because the rheology slows down the subsequent period of growth. Throughout the evolution, $\max _z|\tilde {\tau }_w|$ is larger for the $B=0$ simulation than for $B=0.05$, suggesting that the yield stress decreases the wall shear stress during this pre-coalescence phase of plug formation. Note that $\max _z\hat {Y}$ and $\max _z|\tilde {\tau }_w|$ increase rapidly around $t\approx t_p$, suggesting that the fluid in the large lobe is strongly yielded and the wall shear stress increases as plug formation occurs. However, we cannot expect the theory to remain accurate during this period since the assumption that radial velocities are weak no longer holds. As in the Newtonian problem (Johnson et al. Reference Johnson, Kamm, Ho, Shapiro and Pedley1991), fully two-dimensional theory is required to capture the coalescence phase of the evolution.

A new phenomenon which we have observed only in long-wave simulations is that the small collar of fluid which forms near $z=0$ can rigidify during the evolution. At $t=30$ (figure 6a), the layer is fully yielded with $\hat {Y}>0$ for all $z\in (0,L)$. Figure 6(c) shows that a yield surface (where $\hat {Y}$ makes contact with zero) then travels from $z=0$ at around $t=50$ to around $z=1.8$ at $t=130$, indicating that almost the entire small collar has rigidified by this point. The speed at which the yield surface travels through the domain decreases slightly through the evolution. We observe a small jump in $\hat {p}_z$ at the point where $\hat {Y}$ becomes zero, but $\hat {p}_z$ remains non-positive everywhere including in the rigid region. We generally observe this rigidification of the small collar in our simulations whenever there is enough fluid to form a plug, i.e. $\epsilon \gtrsim 0.12$. For thinner layers, $\epsilon \lesssim 0.11$, we generally do not observe this phenomenon; instead, $\hat {Y}$ decays to zero from above everywhere in the domain, qualitatively the same behaviour as we observed in thin-film simulations (e.g. figure 2).

4.2. Dependence on capillary Bingham number, layer thickness and initial conditions

We investigate the dependence of the long-wave evolution on the capillary Bingham number, layer thickness and initial conditions by running large numbers of numerical simulations varying $B$, $\epsilon$ and $A$. This allows us to examine the dependence of the critical thickness required to form a plug, $\epsilon _c$, on $B$ and $A$. In figures 7 and 8 we plot the final maximum heights of the layers from the numerical simulations. We run all simulations to $t=1000$, or if a plug begins to form before this time, the simulation is stopped and the stopping time is identified as $t_p$.

Figure 7.

Data from numerical solutions of the long-wave evolution equation for various values of $B$ and $\epsilon$, with $A=0.25$. Coloured dots correspond to simulations which did not plug, with the colour indicating $\max _zH(z,t=1000)$. Grey crosses correspond to simulations which are stopped early due to a plug forming. The critical thickness $\epsilon _c$ required for plug formation can be identified as the boundary of the plugging region. The plugging time, $t_{p}$, is indicated by the grey contours. The maximum $B$ for any yielding to occur, $B_m(\epsilon,0.25)$ (dashed red), provides a strict upper bound on the yield-stabilised region.

Figure 8.

Data from numerical solutions of the long-wave evolution equation for various values of $B$ and $A$, with (a) $\epsilon =0.12$, (b) $\epsilon = 0.125$ and (c) $\epsilon =0.13$. Each dot corresponds to a solution with the colour indicating $\max _zH(z,t=1000)$. The data are interpolated linearly to produce the black contour lines, which are evenly spaced. Grey points correspond to simulations which are stopped early due to a plug forming. In (c) grey contours in the plugging region indicate the plugging time, $t_p$. The critical amplitude for any yielding to occur, $A_m(B,\epsilon )$ (dashed red), provides a strict lower bound on the boundary of the yield-stabilised region. Unlike in the thin-film case (figure 5a), the quantity $1-H_0(z=0;B,\epsilon )$ (solid magenta), where $H_0$ are the static solutions computed in Appendix E, does not provide useful information on the final shape of these layers.

In figure 7 we vary $\epsilon$ and $B$ between simulations while the initial perturbation amplitude is fixed at $A=0.25$. There are three distinct regions in the data, corresponding to three qualitatively different outcomes for the layer. Firstly, when $B$ is sufficiently large, there is no growth so we say the layer is yield stabilised. As in the thin-film case in § 3.4, we can find the minimum amplitude for yielding, $A=A_m(B,\epsilon )$, such that the layer is initially fully rigid if and only if $A\leq A_m$. This defines a corresponding capillary Bingham number, $B=B_m(\epsilon,A)$, such that for a given $A$ and $\epsilon$, the layer is fully rigid if and only if $B\geq B_m$. We compute $B_m(\epsilon,A)$ numerically by finding $B=B_m$ such that the initial condition (2.19) makes $\psi =1$ at exactly one point in the domain, where $\psi$ is defined in (2.13). Figure 7 shows that $B=B_m(\epsilon,A=0.25)$ provides a strict upper bound on the boundary of the yield-stabilised region. The second possible outcome for the layer is that there is unstable growth but a liquid plug does not form, and instead the final shape at $t=1000$ is a two-collar configuration, like in the thin-film simulations. Figure 7 shows that this occurs for roughly $0\leq B\lesssim 0.1$ and $\epsilon \lesssim 0.12$, and the figure also indicates how the final peak height of the large collar, $\max _zH(z,t=1000)$, depends on both $\epsilon$ and $B$.

In the third possible outcome, the layer forms a liquid plug, which occurs when the layer is sufficiently thick, and $B$ is not too large for it to be yield stabilised. The boundary of the plugging region in figure 7 corresponds to the critical thickness required for plug formation to occur, $\epsilon _c$, which can be seen to depend strongly on $B$. For very small $B$, the critical value is around $\epsilon _c\approx 0.12$, which is consistent with the results of Gauglitz & Radke (Reference Gauglitz and Radke1988). As $B$ is increased, $\epsilon _c$ first increases slowly, to somewhere in the range $0.125<\epsilon _c<0.1275$ when $B=0.11$, with two-collar final shapes being observed at $\epsilon =0.125$ when $B=0.11$. This suggests that for a few values of $\epsilon$ around $\epsilon \approx 0.125$, the layer can exhibit plugging, two-collar or yield-stabilised behaviour, depending on the value of $B$. For $B\geq 0.12$, the plugging region is bounded by the yield-stabilised region and $\epsilon _c$ increases rapidly as $B$ is increased. Figure 7 also shows how the plugging time, $t_p$, decreases as $\epsilon$ is increased and increases as $B$ is increased. There is a rapid increase in $t_p$ near to the boundary of the plugging region, suggesting that we have located the boundary accurately by running simulations to $t=1000$. Simulations in the plugging region which are near the boundary spend a long time in a near-static two-collar shape (e.g. figure 6a, $t=70$) before eventually transitioning to form a plug.

In figure 8 we investigate the dependence of the evolution on $B$ and $A$, for $\epsilon =0.12, 0.125, 0.13$. We can again identify a yield-stabilised region, a two-collar region and a plugging region for each value of $\epsilon$. The boundary of the yield-stabilised region does not change significantly between $\epsilon =0.12$ and $\epsilon =0.13$, which is consistent with the results in figure 7. This boundary corresponds to the minimum amplitude of perturbation required to trigger instability, which can be seen to strongly depend on $B$. Again, we put a strict lower bound on the boundary of this region using $A_m(B,\epsilon )$, the minimum amplitude for any yielding to occur. The size of the two-collar region decreases quickly as $\epsilon$ is increased, and it has almost entirely disappeared when $\epsilon =0.13$. This suggests that for $\epsilon \geq 0.13$, as long as $A$ is large enough to trigger growth, a plug is guaranteed to form. When the two-collar region does exist (figure 8a,b), the location of its boundary with the plugging region depends on $A$ and $B$. When $A$ is small, this boundary is largely independent of $A$, but when $A$ is large, the plugging region extends to higher $B$. We propose that this is because highly deformed initial conditions place most of the fluid near $z=L$, so, compared with simulations with small $A$, less fluid is trapped in the small collar near $z=0$, making the large collar larger and more unstable to plug formation. The results in figure 8 show that the boundary of the plugging region depends on $A$ as well as on $B$ and $\epsilon$. Hence, the critical thickness for plug formation, $\epsilon _c$, must also depend on $A$ as well as $B$. The plugging time also depends on both $B$ and $A$, as indicated in figure 8(c), with $t_p$ increasing when $B$ is increased and decreasing when $A$ is increased. As in figure 7, $t_p$ increases rapidly near to the boundary of the plugging region.

In the thin-film problem, the quantity $1-H_0(z=0;B)$ from the static solutions could be used to predict the final shape of the layer in a large number of cases (figure 5a). We have also computed static solutions, $H_0(z;B,\epsilon )$, for the long-wave problem by solving $\psi =1$ and assuming monotonic pressure (see Appendix E). The quantity $1-H_0(z=0;B,\epsilon )$ is plotted for $\epsilon =0.12,0.125,0.13$ in figure 8. It predicts the threshold amplitude required to trigger instability for small $B$, but underestimates it for larger $B$, and the prediction becomes less accurate as $\epsilon$ is increased. The upper branch of the curve does not appear to be correlated with the numerical results in any way. In numerical simulations with $\epsilon \geq 0.12$, we generally observe that rigidification of the small collar near $z=0$ (figure 6b) occurs whether a plug forms or not. When this rigidification happens, the final static shape of the layer does not satisfy $\psi =1$ everywhere, so the layer selects different static solutions than those we have computed. The static solutions we have computed may predict the final shapes of the layer when $\epsilon$ is smaller, but they cannot do so when there is rigidification of the small collar or when there is plug formation.

The results in figures 7 and 8 are specific to the sinusoidal form of initial perturbation used in (2.19). However, we also ran simulations with initial conditions of the form $H(z,0)=C_0+A\tanh {(2z-L)}$, where $C_0$ is a constant chosen so that the total volume of fluid is independent of $A$, and found the same qualitative behaviour and only minor quantitative differences. This suggests that the observed behaviour is not strongly dependent on the exact form of initial perturbation.