Surface-tension-driven evolution of a viscoplastic liquid coating the interior of a cylindrical tube

Abstract

One mechanism for airway closure in the lung is the surface-tension-driven instability of the mucus layer which lines the airway wall. We study the instability of an axisymmetric layer of viscoplastic Bingham liquid coating the interior of a rigid tube, which is a simple model for an airway that takes into account the yield stress of mucus. An evolution equation for the thickness of the liquid layer is derived using long-wave theory, from which we also derive a simpler thin-film evolution equation. In the thin-film case we show that two branches of marginally yielded static solutions of the evolution equation can be used to both predict the size of the initial perturbation required to trigger instability and quantify how increasing the capillary Bingham number (a parameter measuring yield stress relative to surface tension) reduces the final deformation of the layer. Using numerical solutions of the long-wave evolution equation, we quantify how the critical layer thickness required to form a liquid plug in the tube increases as the capillary Bingham number is increased. We discuss the significance of these findings for modelling airway closure in obstructive conditions such as cystic fibrosis, where the mucus layer is often thicker and has a higher yield stress.

1. Introduction

The surface-tension-driven instability of a liquid layer lining the interior of a cylindrical tube is of physiological importance since, when it occurs in a lung airway, it can cause obstruction or airway closure by redistributing the liquid lining the airway, potentially leading to formation of a liquid plug. The liquid that lines the lungs' airways consists primarily of mucus, a non-Newtonian fluid exhibiting various rheological properties such as shear thinning and viscoelasticity (Hill et al. 2022). Importantly, mucus also has a yield stress, which is significantly increased in diseases such as cystic fibrosis (CF) and chronic obstructive pulmonary disease (COPD) compared with its typical value in healthy lungs (Patarin et al. 2020). Increased prevalence of airway obstruction by mucus plugging is also a key symptom of CF and COPD (Mall 2016). Motivated by this application, we study the effect of viscoplastic liquid rheology on the evolution of a layer coating the interior of a cylindrical tube, which is a simple model for an airway. Additionally, there are numerous potential applications of this class of flow in engineering and industry, as highlighted by Craster & Matar (2009) for thin-film and coating flows, and Balmforth, Frigaard & Ovarlez (2014) for thin-film and free-surface viscoplastic flows.

The surface-tension-driven flow of a viscous film coating a rigid circular cylinder has been well studied in the case that the liquid layer is Newtonian. Goren (1962) identified that a flat layer can be linearly unstable. Everett & Haynes (1972) found and analysed capillary-static configurations of a volume of liquid inside a tube, which were either annular collars of fluid or liquid plugs. A nonlinear evolution equation using thin-film theory was first derived and solved by Hammond (1983), who found that an initially flat layer evolves into a configuration with large quasi-static annular collars separated by thin films which slowly drain into the collars. Lister et al. (2006) studied the long-time dynamics of the thin-film system, finding that at very long times collars can translate along the tube, potentially consuming other collars in the process, provided other physical effects do not intervene first, while Xu & Jensen (2017) showed how collars can be pinned by wall roughness. Hammond's theory was extended by Gauglitz & Radke (1988) to predict plug formation by retaining certain higher-order terms in the thin-film theory, notably the exact free-surface curvature. This approach provides a composite approximation to the evolution of layers with thickness comparable to the tube radius, which accurately determines capillary-static effects whilst approximating the dynamics well where the layer is thin. They identified a critical average layer thickness, approximately $12\,\%$ of the tube radius, required for a liquid plug to form during the evolution. A similar composite approximation was compared with full two-dimensional numerical simulations by Johnson et al. (1991); whilst their quasi-one-dimensional theory could predict when a plug would form, it could not capture the genuinely two-dimensional dynamics which occur around coalescence. Otis et al. (1993) derived a similar reduced-order model by making a long-wave assumption when simplifying the governing equations.

At this point, we clarify the distinction between thin-film and long-wave approaches to deriving reduced-order evolution equations: in thin-film theory, it is assumed that the thickness of the layer is much smaller than the radius of the tube, while in long-wave theory, it is assumed that the tube radius is much smaller than the characteristic axial length scale of the flow but the layer is not necessarily thin compared with the radius. Making the thin-film assumption results in an evolution equation with the same mobility function as would appear in the planar case, and the curvature of the cylindrical geometry is felt only through the linearised free-surface curvature. In long-wave theory, additional terms appear in the mobility which arise due to the curvature of the geometry, and the full expression for the free-surface curvature is generally retained. Thin-film models benefit from their relative simplicity, but long-wave models capture the effects of the curved geometry more accurately (Camassa & Ogrosky 2015), allowing the dynamics leading to plug formation to be described.

Various physical effects that can modify the evolution of the coating layer have previously been incorporated into models. Halpern, Fujioka & Grotberg (2010) used a long-wave evolution equation to model the effect of viscoelasticity, showing that the critical layer thickness for plug formation is not changed but the time to form a plug can be shortened by increasing the Weissenberg number (a parameter proportional to the relaxation time of the fluid). Romanò et al. (2019) used a volume-of-fluid method to model the pre- and post-coalescence phases of Newtonian plug formation, and recently extended this work to include the effect of viscoelasticity (Romanò et al. 2021). They found that post-coalescence bi-frontal plug growth can induce significant stresses on the tube wall, and that viscoelasticity can induce additional wall stress due to the occurrence of an elastic instability. Erken et al. (2022) studied the instability of a two-layer coating on the interior of a tube as a model for a lung airway which takes into account the periciliary liquid layer that lies beneath the mucus layer and is generally less viscous than mucus. They found that plug formation occurs more quickly in a two-layer model due to the lubricating effect of the base layer, but that the combined critical thickness of the layers required for plug formation can be significantly larger than the single-layer result of Gauglitz & Radke (1988). Halpern & Grotberg (1992) modelled the evolution of a liquid layer coating an elastic tube and subsequently extended the model to include the effect of insoluble surfactants (Halpern & Grotberg 1993). They found that the presence of surfactant can significantly increase the critical layer thickness required to form a plug and can delay plug formation when it does occur, while decreasing the wall stiffness has the opposite effect, decreasing both the critical layer thickness for plug formation and the closure time. Heil, Hazel & Smith (2008) showed that the volume of liquid required for a plug to form is significantly decreased if there is non-axisymmetric collapse of the elastic tube wall. Halpern & Grotberg (2003) developed a thin-film model that included the effect of an oscillating air flow in the centre of the tube, showing that at certain frequencies of oscillation, air flow can suppress deformation of the liquid layer. Camassa, Ogrosky & Olander (2014) developed a long-wave model for gravity-driven flow, and identified families of travelling-wave solutions which they found can be used to predict the critical thickness for plug formation as a function of the Bond number. Camassa, Ogrosky & Olander (2017) also found travelling-wave solutions for the case of flow driven by air flow in the centre of the tube and recently Ogrosky (2021) extended the long-wave model to include the combined effects of gravity, air flow and surfactant.

Turning to viscoplastic flows with applications in airway modelling, Craster & Matar (2000) modelled surfactant-driven flow on a single-layer or two-layer film of viscoplastic or Newtonian fluids, showing that, at least in the single-layer case, yield stress decreases spreading rates and can cause the layer to become frozen in a non-trivial static shape. Modelling of propagation of viscoplastic liquid plugs in tubes and channels has shown that increasing the yield stress increases the stress applied to the wall and increases the thickness of the layer of liquid left behind as a plug propagates (Zamankhan et al. 2012; Zamankhan, Takayama & Grotberg 2018). Rupture of viscoplastic liquid plugs has also been modelled both experimentally (Hu et al. 2015) and numerically (Hu, Romanò & Grotberg 2020), showing that increased yield stress can inhibit plug rupture because a larger pressure drop is required across the plug to make it yield. Recently, Bahrani et al. (2022) proposed a model of elastoviscoplastic plugs, which they validated against experimental results, showing that increased yield stress slows the propagation of a plug but can speed up its rupture since the trailing film thickness is increased. The distribution of mucus throughout a whole lung has also been studied using a viscoplastic model for mucus, showing that the yield stress and the strength of air flow (in this case modelling air flow induced by chest physiotherapy) influence the mucus layer thickness in each airway generation (Mauroy et al. 2011, 2015).

The exposition of viscoplastic thin-film theory by Balmforth & Craster (1999) has provided the basis for various studies of canonical viscoplastic free-surface flows. In their theory, there are regions of plug-like flow near the free surface, which have the same structure as the ‘pseudo-plugs’ first identified in a bounded annular flow by Walton & Bittleston (1991). Viscoplastic thin-film theory was used by Balmforth et al. (2000) to study axisymmetrically spreading gravity currents and the work was recently extended to model droplets spreading under surface tension as well as gravity (Jalaal, Stoeber & Balmforth 2021), showing that after spreading, the fluid is frozen in a non-trivial static shape in which the hydrostatic or capillary pressure is balanced by resistance from the yield stress. Jalaal et al. (2021) also compared results for the final shape and size of the droplets computed using thin-film theory to results from computational fluid dynamics (CFD) simulations showing good agreement except when the capillary Bingham number was very large. Gravity-driven flow down inclined planes has been well studied using viscoplastic thin-film theory, as reviewed by Balmforth et al. (2007a). Balmforth, Ghadge & Myers (2007b) investigated the surface-tension-driven fingering instability of a film travelling down an inclined plane, finding that increasing the Bingham number (which measures yield stress relative to viscous stress) slows growth of the linear instability and, when it is above a critical value, instability is fully suppressed. Jalaal & Balmforth (2016) also used thin-film theory to model the steady propagation of a bubble through a tube filled with viscoplastic fluid, and compared their results to CFD simulations, showing that thin-film theory predicts the behaviour accurately when the liquid film is thin but less well when the Bingham number is increased and the film is thicker. Viscoplastic flows are often solved using regularisation of the constitutive equation; Frigaard & Nouar (2005) review popular regularisation approaches. Jalaal (2016) introduced a regularisation specifically designed for thin-film flows which we describe in § 2.4 and use when solving our evolution equations numerically.

With this study, we aim to quantify the effect of viscoplastic liquid rheology on the surface-tension-driven Rayleigh–Plateau instability of a layer coating the interior of a cylindrical tube. This flow has not previously been studied in the case that the liquid layer is viscoplastic. We ask how the yield stress affects the dynamics during the evolution of the layer and the critical layer thickness required to form a plug. To answer these questions, we derive an evolution equation using long-wave theory, with the detailed flow structure inspired by the viscoplastic thin-film theory of Balmforth & Craster (1999), but we include additional terms arising from the cylindrical geometry which are neglected in the thin-film approximation. We then show how this model reduces, in the appropriate limit, to a thin-film evolution equation, analogous to the Newtonian version derived by Hammond (1983). Other complicating effects are neglected in the model so that the effect of the viscoplastic rheology can be examined in isolation: the tube is rigid, the flow is axisymmetric, surface tension is constant and the air in the centre of the tube is passive and inviscid. We use the Bingham model for the liquid layer since it is the simplest viscoplastic rheology, without the potentially complicating effects of, for example, shear thinning, elasticity or thixotropy. We compute marginally yielded static solutions to the thin-film evolution equation, and show how they can be used to predict both the size of perturbation to a flat layer required to trigger instability, and the final shape of the layer when there is instability. The thin-film theory cannot, however, predict formation of liquid plugs. By solving the long-wave evolution equation numerically, we examine the critical layer thickness required for plug formation and the time taken for plugs to form, quantifying how both can be increased by increasing the capillary Bingham number.

The rest of the paper will be organised as follows. In § 2 we formulate the models, presenting the long-wave evolution equation in § 2.2, and the thin-film equation in § 2.3. A brief discussion of the methods for solving these equations is given in § 2.4. Results for thin layers are presented in § 3. We discuss a representative numerical solution of the thin-film evolution equation in § 3.1, and we examine the behaviour of the layer at long times in § 3.2. We compute and analyse static solutions of the thin-film evolution equation in § 3.3, and investigate the dependence of the evolution on the initial conditions and the capillary Bingham number in § 3.4. Results for layers with finite thickness are presented in § 4. We discuss an example numerical solution of the long-wave equation in § 4.1 and examine the dependence of the evolution on the capillary Bingham number, layer thickness and initial conditions in § 4.2, including discussion of the critical layer thickness for plug formation and the time taken to form a plug. A summary of the results, and a discussion of their significance for modelling lung airways is given in § 5.

2. Model formulation

2.1. The Stokes system

We consider a rigid circular cylinder of radius $a$ coated on the inside by a layer of Bingham fluid. The rest of the tube is filled with a gas which is assumed inviscid with spatially uniform pressure. The geometry is illustrated in figure 1. We consider only the flow in the liquid layer. The flow is assumed to be axisymmetric, and is described by cylindrical coordinates $(r^*,z^*)$. The air–liquid interface is located at $r^* = R^*(z^*,t^*)=a-H^*(z^*,t^*)$. The fluid velocity in the film is $(u^*(r^*,z^*,t^*),w^*(r^*,z^*,t^*))$, where $u^*$ and $w^*$ are measured in the positive $r^*$ and $z^*$ directions, respectively. The non-zero components of the shear-rate tensor, $\boldsymbol {\dot \gamma }^* = \boldsymbol {\nabla } \boldsymbol {u}^*+{\boldsymbol {\nabla } \boldsymbol {u}^*}^\mathrm {T}$, are therefore

(2.1a–d)\begin{equation} \dot\gamma^*_{rr} = 2 \partial^*_ru^*, \quad \dot\gamma^*_{rz} = \partial^*_rw^* + \partial^*_zu^*, \quad \dot\gamma^*_{\theta\theta} = 2\frac{u^*}{r^*}, \quad \dot\gamma^*_{zz} = 2\partial^*_zw^*.\end{equation}

The liquid is assumed to be incompressible and have no inertia, so the flow is governed by the Stokes equations,

(2.2a)$$\begin{gather} 0 = \partial^*_zw^* + \frac{1}{r^*}\partial_r^*(r^*u^*), \end{gather}$$

(2.2b)$$\begin{gather}0 ={-}\partial^*_rp^* + \frac{1}{r^*}\partial_r^*(r^*\tau^*_{rr}) + \partial^*_z\tau^*_{rz} - \frac{\tau^*_{\theta\theta}}{r^*}, \end{gather}$$

(2.2c)$$\begin{gather}0 ={-}\partial^*_zp^* + \frac{1}{r^*}\partial_r^*(r^*\tau^*_{rz}) + \partial^*_z\tau^*_{zz}, \end{gather}$$

where $\boldsymbol {\tau }^*(r^*,z^*,t^*)$ is the stress tensor and $p^*(r^*,z^*,t^*)$ is pressure measured relative to the gas pressure. The Bingham fluid constitutive relation is

(2.3)\begin{equation} \left. \begin{array}{ll@{}} \displaystyle \tau_{ij}^* = \left(\eta + \dfrac{\tau_Y}{\dot\gamma^*}\right)\dot\gamma_{ij}^* & \mbox{if}\ \tau^* > \tau_Y,\\ \dot\gamma_{ij}^* = 0 & \mbox{if}\ \tau^*\leq\tau_Y, \end{array}\right\}\end{equation}

where $\eta$ is a viscosity, $\tau _Y$ is the yield stress, and $\dot \gamma ^*$ and $\tau ^*$ are the second invariants of shear rate and stress, respectively. The second invariant of a tensor $\boldsymbol{\mathsf{T}}$ is defined as ${\mathsf{T}} = \sqrt {{{\mathsf{T}}_{ij}{\mathsf{T}}_{ij}}/{2}}$.

Figure 1. Sketches of the geometry and flow structure in the long-wave model. In non-dimensionalised variables the free surface is at $r=R$ and the layer thickness is $\epsilon H$. The surface $r=\varPsi$ separates a region of shear-dominated flow near the cylinder wall and a region of plug-like flow near the free surface. We define $\varPsi \equiv \min (1,\psi )$ with $\psi$ defined in (2.13) and $\varPsi =1$ corresponding to regions of unyielded fluid. The thin-film model has qualitatively the same flow structure.

The boundary conditions are as follows. There is no slip and no penetration at the cylinder wall,

(2.4)\begin{equation} u^* = w^* = 0 \ \mbox{on}\ r^* = a.\end{equation}

The kinematic boundary condition at the free surface is

(2.5)\begin{equation} \partial^*_tR^* + w^* \partial^*_zR^* = u^* \ \mbox{on}\ r^*=R^*(z^*,t^*). \end{equation}

The gas phase applies no shear stress to the liquid film so

(2.6)\begin{equation} -p^*n_i + \tau^*_{ij}n_j = \sigma\kappa^* n_i \ \mbox{on}\ r^*=R^*(z^*,t^*), \end{equation}

where $n_i$ are the components of the normal to the free surface, $\sigma$ is the constant value of surface tension and

(2.7)\begin{equation} \kappa^* = \frac{1}{\sqrt{1+(\partial^*_zR^*)^2}} \left[\frac{1}{R^*}-\frac{\partial^*_{zz}R^*}{1+(\partial^*_zR^*)^2}\right] \end{equation}

is the free-surface curvature. At the side boundaries, we impose symmetry boundary conditions,

(2.8)\begin{equation} \partial^*_zR^*=\tau^*_{rz}=w^*=0 \ \mbox{at}\ z=\{0,L^*\}. \end{equation}

Rather than solve the full Stokes problem defined above, we will derive reduced-order models using long-wave and thin-film theories, presented in §§ 2.2 and 2.3, respectively.

Surface tension at the air–liquid interface introduces an associated energy, proportional to the surface area of the air–liquid interface,

(2.9)\begin{equation} E^* \equiv \sigma\int_0^{L^*}2{\rm \pi} R^*\sqrt{1+(\partial^*_zR^*)^2}\,\mathrm{d}z^*.\end{equation}

In Appendix A.1 we show that the Stokes equations and boundary conditions (2.2)–(2.8) imply that

(2.10)\begin{equation} \partial_t^*E^* ={-}\int_V(\eta(\dot\gamma^*)^2+\tau_Y\dot\gamma^*)\,\mathrm{d}V\leq 0, \end{equation}

where $V$ is the volume of the layer, so the interfacial energy is always decreasing. In the Newtonian problem the final shape that the layer reaches after its evolution can be found by solving for shapes which locally minimise interfacial energy (Everett & Haynes 1972). However, in the viscoplastic problem this is not necessarily the case because we expect that the yield stress may freeze some or all of the layer before it has reached a minimal energy state. Analysing the final static shapes of viscoplastic layers will form a large part of our discussion, particularly in the thin-film case (cf. § 3.3).

2.2. The long-wave model

We non-dimensionalise the governing equations and boundary conditions (2.2)–(2.8) by defining

(2.11)\begin{equation} \left. \begin{array}{c@{}} \displaystyle ( r, z) = \left(\dfrac{r^*}{a},\dfrac{z^*}{a}\right),\quad (u,\hat{w}) = \dfrac{\eta}{\sigma}(u^*,w^*), \quad R = \dfrac{R^*}{a}, \quad\boldsymbol{ \tau} = \dfrac{a}{\sigma}\boldsymbol{\tau}^*, \\ \displaystyle \boldsymbol{{\dot\gamma}}=\dfrac{\eta a}{\sigma}\boldsymbol{\dot\gamma}^*, \quad \hat{t}= \dfrac{\sigma}{a\eta}t^*,\quad \hat{p}= \dfrac{a}{\sigma}p^*, \quad \hat\kappa = a\kappa^*, \quad L =\dfrac{L^*}{a}, \end{array} \right\}\end{equation}

where hats are used to distinguish $\hat {w}$, $\hat {t}$, $\hat {p}$ and $\hat {\kappa }$ from the scaled, thin-film quantities $w$, $t$, $p$ and $\kappa$ which we will define in § 2.3. After non-dimensionalising, we consider (2.2)–(2.8) in a long-wave limit by introducing a characteristic axial length scale for the flow, $a/\delta$, where $\delta \ll 1$. There is no assumption at this point that the liquid layer is thin. We then rescale the system using the small aspect ratio, $\delta$, defining

(2.12a–e)\begin{equation} \bar{z}\equiv \delta z, \quad \bar{u}\equiv\frac{u}{\delta}, \quad \bar{t}\equiv \delta\hat{t}, \quad \bar{p}\equiv\delta \hat{p}, \quad \bar{L}\equiv \delta L, \end{equation}

with other variables remaining unstretched. The resulting scaled, dimensionless equations are given in Appendix B.

To derive the long-wave evolution equation, we propose a similar flow structure in the long-wave model to that of the thin-film theory of Balmforth & Craster (1999). Where the fluid is yielded, the flow is separated into a shear-dominated region, $\psi \leq r\leq 1$, adjacent to the no-slip boundary, where the shear stress is large compared with the normal stresses, and a region, $R\leq r < {\psi }$, adjacent to the free surface, where we say the flow is ‘plug-like’ as the axial velocity, $\hat {w}$, is independent of $r$ (figure 1). We make separate expansions for the velocities and stresses in the shear-dominated and plug-like regions. In the plug-like region the shear stress is below the yield stress, but the fluid is still yielded because the normal stresses are large enough that the total stress exceeds the yield stress. We solve at leading order in $\delta$ in the shear-dominated and plug-like regions, match the solutions from the two regions together, and finally arrive at the evolution equation which we state below. The full derivation is given in Appendix B.

Following the approach of, e.g. Camassa et al. (2012), we write the evolution equation in terms of the unscaled variables (2.11) rather than the scaled variables (2.12a–e), so $\delta$ does not appear in the equations. However, the limit in which the theory is formally valid remains $\delta \ll 1$.

From now on, we will use subscripts to denote derivatives. We determine the surface between the shear-dominated and plug-like regions to be

(2.13)\begin{equation} {\psi}(z,\hat{t}) = \frac{ \hat{B}}{|\hat{p}_z|} \left(1+\sqrt{1+\left(\frac{|\hat{p}_z| {R}}{ \hat{B}}\right)^2}\right),\quad\mbox{where}\ \hat{B}\equiv\frac{\tau_Ya}{\sigma}\end{equation}

is a capillary Bingham number, which measures yield stress relative to capillary stress. We use the hat notation to distinguish $\hat {B}$ from the thin-film version, $B$, which we will introduce in § 2.3. The capillary pressure is proportional to the free-surface curvature, and is given by

(2.14)\begin{equation} \hat{p} ={-}\hat{\kappa} ={-}\frac{1}{ R \sqrt{1+R_z^2}}\left(1-\frac{RR_{zz}}{ 1+R_z^2}\right).\end{equation}

We have retained the full expression for capillary pressure (2.14) including all terms which are higher order in $\delta$. Although this is not strictly consistent with the asymptotic analysis, it allows the evolution equation to describe capillary-static effects accurately and is a widely used device (e.g. Gauglitz & Radke 1988). The axial velocity is defined separately in the shear-dominated and plug-like regions,

(2.15)\begin{align} \hat{w} = \left\{ \begin{array}{@{}ll} \frac{1}{2}\hat{p}_z[\frac{1}{2}( {r}^2-1)- {\varPsi}^2\log (r)] + \hat{B}\, {\rm sgn}(\hat{p}_z)[ {\varPsi}\log{(r)}+1- r], & R\leq r < \varPsi, \\ \frac{1}{2}\hat{p}_z[\frac{1}{2}( \varPsi^2-1)- {\varPsi}^2\log (\varPsi)] + \hat{B}\, {\rm sgn}(\hat{p}_z)[ {\varPsi}\log{(\varPsi)}+1- \varPsi], & \varPsi \leq r \leq 1, \end{array} \right. \end{align}

where $\varPsi (z,\hat {t})\equiv \min (1,\psi )$. The function $\varPsi$ is defined so that (2.15) applies to the whole layer, including regions of unyielded fluid. Where $\varPsi =1$, the fluid is unyielded, so there is no motion and $\hat {w}=0$. The axial flux, $\hat {Q}$, is calculated by radially integrating $\hat {w}$.

The long-wave evolution equation is

(2.16)\begin{equation} R_{\hat{t}}+\frac{1}{ {R}}\hat{Q}_z=0, \quad \mbox{where}\ \hat{Q}= \frac{\hat{p}_z}{16}f_1(R,\varPsi)+\frac{\hat{B}}{12}\, {\rm sgn}(\hat{p}_z) f_2(R,\varPsi),\end{equation}

with non-negative functions $f_1$ and $f_2$ (see figure 9 below) given by

(2.17a)$$\begin{gather} f_1(R,\varPsi)\equiv (1- {\varPsi}^2)^2-2 {R}^2(1- {\varPsi}^2 +2 {\varPsi}^2\log{ {\varPsi}}), \end{gather}$$

(2.17b)$$\begin{gather}f_2(R,\varPsi)\equiv 2-3 {\varPsi}+ {\varPsi}^3-6 {R}^2( {\varPsi} \log{ {\varPsi}}+1- \varPsi). \end{gather}$$

The boundary conditions at the sides of the domain are

(2.18)\begin{equation} R_z = \hat{Q} = 0 \ \mbox{at}\ z = \{0,L\}.\end{equation}

The initial conditions which we impose when solving (2.16) are

(2.19)\begin{equation} R(z,t=0) = \sqrt{(1-\epsilon)^2-\epsilon^2A^2/2}-\epsilon A\cos \left(\frac{{\rm \pi} z}{L}\right),\end{equation}

which corresponds to a flat layer perturbed by a single Fourier mode with wavelength $2L$. The constant $A$ is the perturbation amplitude and the constant $\epsilon$ is the ratio of average layer thickness to tube radius when $A=0$. The constant term in (2.19) is chosen so that the total volume of the layer is independent of $A$ for a given $\epsilon$.

We derive an expression for the shear stress in $\varPsi \leq r\leq 1$ (B11), which when evaluated at $r=1$ gives the stress exerted on the tube wall,

(2.20)\begin{equation} \hat{\tau}_w \equiv \tfrac{1}{2}\hat{p}_z(1-\varPsi^2) + \hat{B}\, {\rm sgn}(\hat{p}_z)\varPsi.\end{equation}

Note that (2.20) only holds in regions where the fluid is yielded (where $\varPsi <1$). In unyielded regions, the stress is not defined by the constitutive relation (2.3) but we do know that the wall stress must be bounded by the yield stress, so $|\hat \tau _w| \leq \hat {B}$ where $\varPsi =1$.

From (2.9), the dimensionless interfacial energy is

(2.21)\begin{equation} E \equiv 2{\rm \pi}\int_0^LR\sqrt{1+R_z^2}\,\mathrm{d}z.\end{equation}

In Appendix A.2 we deduce directly from (2.14)–(2.18) that $E_{\hat {t}}\leq 0$. Hence, the result (2.10) is preserved in the long-wave theory.

2.3. The thin-film model

We now derive the analogous evolution equation for a thin film, no longer requiring $\delta \ll 1$, but assuming that $|1-R|\ll 1$. The thin-film evolution equation can be derived using the approach of Balmforth & Craster (1999), but here we derive it by taking a thin-film limit of the long-wave system (2.13)–(2.18). The thin-film approximation acts to flatten the geometry, so that terms in the long-wave mobility (2.17) which arise due to the curvature of the geometry are negligible. The effect of the cylindrical geometry is then only felt through the free-surface curvature, which is linearised.

We consider a layer with characteristic thickness $\epsilon a$, and now let $\epsilon \ll 1$. We rescale time, defining $t\equiv \epsilon ^3\hat {t}$, then define the dimensionless film thickness, $H(z,t)$, which satisfies

(2.22)\begin{equation} R(z,\hat{t}) = 1-\epsilon H(z,t).\end{equation}

It is convenient to define a radial coordinate, $y$, measured from the no-slip boundary, which satisfies $r=1-\epsilon y$. Then the free surface is located at $y=H$. Substituting (2.22) into (2.14) gives $\hat {p}=-\hat {\kappa }=-1-\epsilon (H+H_{zz})+O(\epsilon ^2)$. We define the thin-film curvature and capillary pressure as $\kappa \equiv ({\hat {\kappa }-1})/{\epsilon }$ and $p \equiv ({1+\hat {p}})/{\epsilon }$, then linearise in $\epsilon$, so that the pressure gradient driving the flow is

(2.23)\begin{equation} p_z={-}\kappa_z ={-}H_z-H_{zzz}.\end{equation}

We define the thin-film capillary Bingham number as

(2.24)\begin{equation} B \equiv \frac{\hat{B}}{\epsilon^2} = \frac{\tau_Ya}{\sigma\epsilon^2}. \end{equation}

Expanding (2.13), we find that $\psi = 1 - \epsilon \mathcal {Y}+O(\epsilon ^2)$, where

(2.25)\begin{equation} \mathcal{Y}\equiv H-\frac{B}{|p_z|}.\end{equation}

As before, we augment this definition so that it holds in yielded and unyielded regions: we define $Y\equiv \max (0,\mathcal {Y})$, which obeys $\varPsi =1-\epsilon Y+O(\epsilon ^2)$, with $Y=0$ corresponding to regions of unyielded fluid. Where $Y>0$, the fluid is yielded, and the flow structure is qualitatively the same as in the long-wave model. The value of $Y$ then indicates the boundary between the shear-dominated and plug-like regions of flow. The axial velocity (2.15) becomes $\hat {w} = \epsilon ^3w + O(\epsilon ^4)$, where the thin-film axial velocity is

(2.26)\begin{equation} w = \left\{ \begin{array}{@{}ll} \frac{1}{2}p_z y(y-2Y), & 0\leq y < Y, \\ -\frac{1}{2}p_zY^2, & Y \leq y \leq H. \end{array} \right. \end{equation}

Finally, substituting (2.22)–(2.25) into (2.16) and (2.17), and linearising in $\epsilon$, gives

(2.27)\begin{equation} H_t + \tfrac{1}{6}[p_zY^2(Y-3H)]_z =0, \quad \mbox{where}\ Y = \max(0,\mathcal{Y}) .\end{equation}

Equation (2.27), with definitions (2.23) and (2.25), is the thin-film evolution equation.

In the thin-film limit the boundary conditions (2.18) become

(2.28)\begin{equation} H_z = Q = 0 \ \mbox{at}\ z = \{0,L\},\end{equation}

where we define the thin-film flux as $Q\equiv p_zY^2(Y-3H)/6$. Note that in the Newtonian problem (e.g. Hammond 1983), enforcing zero flux at $z=\{0,L\}$ is equivalent to enforcing zero third derivative, $R_{zzz}=0$ or $H_{zzz}=0$. Here, the zero flux conditions (2.18) and (2.28) are preferable because the third derivatives, $R_{zzz}$ and $H_{zzz}$, generally become discontinuous at $z=\{0,L\}$ during the evolution. This also occurs at any interior points where the direction of flow changes (cf. figure 2c). This is an inconsistency in the theory which could be resolved by finding a solution in the inner region (likely of axial length $O(\epsilon )$ in the thin-film system or $O(\delta )$ in the long-wave system) around each of these points and matching these to the global outer solution which we compute. Following the approach of, e.g. Balmforth et al. (2000), we do not solve in these inner regions and assume that the solution which we compute captures the global dynamics of the layer sufficiently accurately.

Figure 2. Time evolution of a thin film with $B=0.05$ and $A=0.2$. (a) Snapshots of layer height $H(z,t)$ and the internal surface $Y(z,t)$ (supplementary movie 1 is available at https://doi.org/10.1017/jfm.2022.479 for the full evolution of $H$ and $Y$ up to $t=1000$). Axial velocity, $w$, as defined in (2.26), is also shown. The plug-like region lies between $y=Y$ and $y=H$, showing significant transient deformation. (b) Time evolution of $\max _zH$ compared with the same quantity for a Newtonian ($B=0$) simulation, and time evolution of $\max _zY$ and $\max _z|\tau _w|-B$ for the $B=0.05$ simulation. (c) Example of the capillary-wave-like structures that are observed ahead of the travelling yield surface at early times, similar to those discussed in Jalaal et al. (2021). Here $Y$ and $p_z$ are shown near to the travelling yield surfaces at $t=1$. The sign changes in $p_z$ indicate reversals in the direction of flow between these structures.

After linearising in $\epsilon$ and combining with (2.22), the initial condition (2.19) becomes

(2.29)\begin{equation} H(z,0) = 1 - A\cos{\left(\frac{{\rm \pi} z}{L}\right)},\end{equation}

which is the initial condition we will use when solving (2.27). We define the thin-film wall shear stress, $\tau _w\equiv \hat {\tau }_w/\epsilon ^2$, then linearise in $\epsilon$ to get

(2.30)\begin{equation} \tau_w = p_zY + B\, {\rm sgn}(p_z) = p_zH,\end{equation}

where we used (2.25) in the second equality. Note that (2.30) only holds in regions where the fluid is yielded ($Y>0$), but we have the bound $|\tau _w|\leq B$ in unyielded regions ($Y=0$). The interfacial energy (2.21), when expanded in powers of $\epsilon$, becomes

(2.31)\begin{equation} E = 2{\rm \pi} L-\epsilon V_0+{\rm \pi}\epsilon^2\int_0^L (H_z^2-H^2)\,\mathrm{d}z+O(\epsilon^3),\end{equation}

where $V_0$ is the (constant) total volume of the layer. We show in Appendix A.3 that (2.23), (2.27) and (2.28) imply that $E_t\leq 0$ for $\epsilon \ll 1$. Hence, the thin-film approximation also preserves the result (2.10).

2.4. Solution methods

When solving both the thin-film and long-wave equations, we choose the domain length to be $L=\sqrt {2}{\rm \pi}$. This length corresponds to the half-wavelength of the most unstable mode in the Newtonian linear stability analysis (Hammond 1983). The instability in the viscoplastic problem is inherently nonlinear, but this choice for $L$ allows direct comparison to previous literature on the Newtonian and viscoelastic versions of the problem (e.g. Gauglitz & Radke 1988; Halpern et al. 2010). We found that small changes in $L$ do not qualitatively affect our results so $L=\sqrt {2}{\rm \pi}$ can be considered a representative domain length. The form of perturbation in the initial conditions, (2.19) or (2.29), then corresponds to the single unstable Fourier mode that exists in the domain. In the long-wave theory $\delta$ is defined as the ratio of tube radius $a$ to a typical axial length scale. If that axial length scale is taken to be the wavelength of the initial disturbance, then $\delta =1/(2L)=1/(2\sqrt {2}{\rm \pi} )$. Shorter wavelength structures also develop in the thin-film and long-wave simulations (cf. capillary waves discussed in § 3.1), testing the validity of the long-wave theory. Pending validation by computations of the full problem (2.1a–d)–(2.8), we anticipate that our results provide a good approximation to the true behaviour, with additional accuracy gained from retaining the exact expression for $\hat {\kappa }$ in (2.14). We do not observe any significantly different behaviour when applying periodic boundary conditions at the sides of the domain compared with the boundary conditions (2.18) or (2.28), validating our use of the latter in all the results presented.

When solving the systems (2.13)–(2.18) or (2.23)–(2.28) numerically, we use the method proposed by Jalaal (2016). The evolution equations are regularised by redefining $Y\equiv \max (Y_{min},\mathcal {Y})$ and $\varPsi \equiv \min (\varPsi _{max},\psi )$, where $Y_{min},1-\varPsi _{max}\ll 1$. We choose $Y_{min}=1-\varPsi _{max}=10^{-6}$ after confirming this is small enough that the results are not sensitive to the precise value of $Y_{min}$. For example, for the simulation presented in figure 6 below, the absolute errors in $\max _zH(z,t=120)$ and $t_p$ (the time to form a plug) are bounded above by $800Y_{min}^2$ and $40000Y_{min}^2$, respectively, for all $10^{-6}\leq Y_{min}\leq 10^{-3}$, which we find to be typical of convergence rates in simulations. Where $Y=Y_{min}$ or $\varPsi =\varPsi _{max}$, there is a very weak regularisation-induced flow, but since $Y_{min}$ is chosen small enough for this flow to be negligible, we treat these regions as unyielded, treating $Y=Y_{min}$ or $\varPsi =\varPsi _{max}$ as equivalent to $Y=0$ or $\varPsi =1$. The regularised equations are solved using the method of lines: the spatial derivatives are approximated using second-order centred finite differences and the resulting system of ordinary differential equations (ODEs) is solved through time using a stiff solver in Matlab. We have confirmed that the number of spatial grid points used is large enough that the precise value does not affect our results.

3. Results: thin-film theory

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 1978; Hammond 1983). 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. (2007b) and Jalaal & Balmforth (2016), and studied in detail by Jalaal et al. (2021) 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 1983), 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. Results: long-wave theory

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 (1988) 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. (2010), 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. 1991), 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 (1988). 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.

5. Discussion

To summarise, we have quantified how viscoplastic rheology can either inhibit growth of, or fully suppress, the surface-tension-driven instability of a layer of liquid coating the interior of a cylindrical tube. We found that for both thin layers and layers with finite thickness, the final shape after evolution depends sensitively on the capillary Bingham number, $B$, as well as on the initial amplitude of perturbation, $A$. Using thin-film theory, we showed that when $A$ is below a critical value, which depends on $B$, there is no unstable growth because the fluid does not yield. When there is unstable growth, the final shape of the layer either coincides with the marginally yielded static solution, $H_0(z;B)$, from the upper branch of figure 4(a), or the final shape has non-monotonic pressure. Figure 5(a) shows that the quantity $1-H_0(z=0;B)$ from the static solutions accurately predicts both the minimum $A$ required to trigger instability, and the large set of $A$ and $B$ for which the final shape of the layer is $H_0(z;B)$. By solving the long-wave evolution, we quantified how the critical layer thickness, $\epsilon _c$, required to form a liquid plug is increased by increasing $B$. Figure 7 shows that $\epsilon _c$ can be increased significantly beyond the Newtonian value of $\epsilon _c\approx 0.12$ found by Gauglitz & Radke (1988), primarily because, when $B$ is sufficiently large, there is no yielding so no unstable growth. For $0.12\leq \epsilon \leq 0.13$, it is possible for there to be no unstable growth, unstable growth leading to plug formation, or unstable growth with no plug formation, depending on the values of $A$ and $B$ (figure 8). For $\epsilon >0.13$, if $A$ is large enough to trigger instability, a plug will form.

One application of our results is to modelling mucus flow and airway closure in lungs. We used a thin-film capillary Bingham number, $B\equiv a\tau _Y/\sigma \epsilon ^2$, to measure the relative strength of the yield stress in the flow. To estimate $B$ for a $12$th generation healthy airway, we propose the following typical values: airway radius $a=0.4$ mm (Hsia, Hyde & Weibel 2016), surface tension $\sigma =30$ mNm (Chen et al. 2019) and mucus yield stress $\tau _Y=0.27$ Pa (Patarin et al. 2020). For a mucus layer with thickness $\epsilon =0.125$, this gives $B\approx 0.2$. Figure 8(b) shows that, for $B= 0.2$, we expect plug formation to be possible, but only for $A\gtrsim 0.6$, i.e. only if the mucus layer is significantly deformed initially. Patarin et al. (2020) measured the yield stress in CF mucus to be $6.34$ Pa, which would correspond to $B\approx 5$ when $\epsilon =0.125$. Figure 8(b) shows that $B=5$ is well inside the yield-stabilised region for all $A$, suggesting that airway closure would not occur via this mechanism for these parameter values. However, other key symptoms of CF are increased volume of mucus in airways and surfactant deficiency (Tiddens et al. 2010), which would correspond to increases in $\epsilon$ and $\sigma$ and, hence, a potentially significant decrease in $B$, making plug formation more likely to occur. Thus, the net effect of CF symptoms on the likelihood of airway closure by this mechanism is not obvious. Experiments and numerical modelling of plug rupture (Hu et al. 2015, 2020) suggest that, once an airway does close, increased yield stress makes airway reopening more difficult, which would contribute to the increased prevalence of plugged airways in CF. Our results also suggest that airway closure could be triggered if $B$ is suddenly decreased, which could be caused by applying certain therapies which are commonly used in CF, such as mucolytics which decrease the mucus yield stress (Patarin et al. 2020) or expectorants which increase the volume of liquid (Donaldson et al. 2006). However, detailed modelling of the effect of such therapies would be required to confirm this conjecture.

We have also shown that yield stress can delay plug formation when it does occur (figures 6c, 7 and 8c). If we set $\eta =10^{-2}\mathrm {Pa}\,\mathrm {s}$ then the dimensional plugging time is $t^*_p=(a\eta /\sigma \epsilon ^3)t_p\mathrm {s}\approx 0.07t_p\mathrm {s}$ for $\epsilon =0.125$. For the Newtonian simulation in figure 6(c), $t_p\approx 70$, corresponding to $t^*_p\approx 5\ \mathrm {s}$, which is approximately the length of a breathing cycle. The $B=0.05$ simulation in figure 6(c) takes about twice as long to form a plug, so $t^*_p$ is likely to be longer than a breathing cycle meaning airway closure is less likely to occur via this mechanism. If we relate $\eta$ to the measured viscosity of mucus (Lai et al. 2009), $\eta =10^{-2}\mathrm {Pa}\,\mathrm {s}$ is a feasible value but it could also be significantly larger, meaning plug formation for these simulations could be on the scale of minutes or hours instead of seconds. The layer thickness also strongly influences $t_p$ (figure 7), and also $t^*_p$ depends inversely on $\epsilon ^3$, so a modest increase in layer thickness can significantly decrease the time taken for plug formation to occur.

Our results also suggest that the shear stress exerted on the tube wall during the pre-coalescence phase of plug formation can be decreased by yield stress (figure 6b). This has physiological significance because a large shear stress exerted on an airway wall may cause epithelial cell damage (Huh et al. 2007). However, we expect that the wall shear stress is likely to be much larger in the post-coalescence phase, as is the case when the liquid is Newtonian (Romanò et al. 2019), so we cannot make conclusions about the effect of yield stress on wall shear stress during the entire closure process. Additionally, when the layer is too thin to form a plug, the wall shear stress at late times approaches the yield stress (figure 2b), so in these cases it increases as yield stress is increased.

We have focused on investigating the effect of viscoplastic rheology, so other physical effects which are relevant to airway modelling have been neglected. Various extensions to our work could be made to investigate how viscoplastic effects interact with, for example, shear stress induced by air flow, elastic tube walls or surfactant, all of which have been studied in the case that the liquid is Newtonian (Halpern & Grotberg 1992, 1993, 2003). Additionally, in order to isolate the effects of the viscoplastic rheology, we have not incorporated shear thinning or viscoelastic rheologies, which are also known to be exhibited by mucus (Hill et al. 2022). It remains an interesting open question how the addition of other rheological properties would affect the dynamics of a viscoplastic layer as studied here.

There are some limitations to the thin-film and long-wave theories that we have used to derive reduced-order models. Thin-film theory cannot predict the formation of liquid plugs, and the quasi-one-dimensional long-wave theory cannot capture the fully two-dimensional dynamics which develop as a liquid plug is forming (requiring simulations to be stopped just before coalescence). Viscoplastic thin-film theory is known to break down at points where the direction of flow changes and the pressure gradient has a jump discontinuity (Balmforth et al. 2000), and we observed this same behaviour in our long-wave simulations. Additionally, the long-wave theory is strictly valid for $\delta \equiv a/L\ll 1$ but we solved the evolution equation in a finite domain with a small but finite value of $\delta$. We have not solved the full axisymmetric Stokes problem here, which could be used to validate the long-wave model.

Our model predicts that viscoplastic rheology can significantly alter the evolution of a layer coating a cylindrical tube. When the layer is thin, key aspects of the dynamics and the final shape of the layer can be understood by studying marginally yielded static solutions. When the layer has finite thickness, the critical thickness required to form a liquid plug can depend strongly on the capillary Bingham number. These results have implications for modelling real-world problems where the coating liquid has a yield stress, such as models of airway closure, particularly in the context of diseases which alter mucus rheology.