2.1. Governing equations and boundary conditions

The flow of the fluid is governed by the incompressible, axisymmetric Navier–Stokes equations

(2.1a) \begin{align} \bar {\rho }(\bar {u}_{\bar {t}}+\bar {u}\bar {u}_{\bar {r}}+\bar {w}\bar {u}_{\bar {z}})&=-\bar {p}_{\bar {r}}+\bar {\mu }\left [\frac {1}{\bar {r}}\partial _{\bar {r}}(\bar {r}\bar {u}_{\bar {r}})+\bar {u}_{\bar {z}\bar {z}}-\frac {\bar {u}}{\bar {r}^2}\right ]\!, \end{align}

(2.1b) \begin{align} \bar {\rho }(\bar {w}_{\bar {t}}+\bar {u}\bar {w}_{\bar {r}}+\bar {w}\bar {w}_{\bar {z}})&=-\bar {p}_{\bar {z}}+\bar {\mu }\left [\frac {1}{\bar {r}}\partial _{\bar {r}}(\bar {r}\bar {w}_{\bar {r}})+\bar {w}_{\bar {zz}}\right ]-\bar {\rho }\bar {g}, \end{align}

(2.1c) \begin{align} \frac {1}{\bar {r}}\partial _{\bar {r}}(\bar {r}\bar {u})+\bar {w}_{\bar {z}}&=0, \end{align}

where $(\bar {u},\bar {w})$ are the velocity components in the $(\bar {r},\bar {z})$ directions, respectively, and $\bar {t}$ is time. Pressure is denoted $\bar {p}$ ; $\bar {\mu }$ , $\bar {\rho }$ and $\bar {g}$ denote viscosity, density and acceleration due to gravity, respectively. Letter subscripts denote partial derivatives, and overbars denote dimensional quantities.

At the tube wall $\bar {r}=\bar {a}(\bar {z},\bar {t})$ , there is a kinematic boundary condition $\bar {u}=\bar {a}_{\bar {t}}+\bar {w}\bar {a}_{\bar {z}}$ . Only radial disturbances to the tube will be considered here; thus in the case of no-slip boundary conditions, one would prescribe $\bar {w}=0$ and $\bar {u}=\bar {a}_{\bar {t}}$ at the wall. Here, we explore the impact of slip on the film flow using a Navier slip boundary condition with slip length $\bar {\varLambda }$ ; due to only radial wall deflection being considered, the resulting boundary conditions at the tube wall are

(2.2a) \begin{align} \bar {w}&=-\bar {\varLambda }\, \bar {w}_{\bar {r}}, \end{align}

(2.2b) \begin{align} \bar {u}&=\bar {a}_{\bar {t}}-\bar {\varLambda }\, \bar {w}_{\bar {r}}\bar {a}_{\bar {z}}. \end{align}

At the free surface, $\bar {r}=\bar {R}(\bar {z},\bar {t})$ , there are three boundary conditions: a kinematic boundary condition

(2.3) \begin{equation} \bar {u}=\bar {R}_{\bar {t}}+\bar {w}\bar {R}_{\bar {z}}; \end{equation}

continuity of tangential stress

(2.4) \begin{equation} \bar {\mu }\Big [(\bar {w}_{\bar {r}}+\bar {u}_{\bar {z}})\Big(1-\bar {R}^2_{\bar {z}}\Big)+2(\bar {u}_{\bar {r}}-\bar {w}_{\bar {z}})\bar {R}_{\bar {z}}\Big ]=\bar {\tau }^{(g)}, \end{equation}

with $\bar {\tau }^{(g)}=\bar {\mu }^{(g)} [(\bar {w}^{(g)}_{\bar {r}}+\bar {u}^{(g)}_{\bar {z}})(1-\bar {R}^2_{\bar {z}})+2(\bar {u}^{(g)}_{\bar {r}}-\bar {w}^{(g)}_{\bar {z}})\bar {R}_{\bar {z}} ]$ denoting the tangential stress exerted by the gas flow on the film’s free surface; and jump in normal stress (according to the Young–Laplace law)

(2.5) \begin{align}& 2\bar {\mu }\Big(\overline {u}_{\overline {r}}+\overline {w}_{\overline {z}}\overline {R}^2_{\overline {z}}\Big)+\bar {\mu }\overline {R}_{\overline {z}}\left(\overline {w}_{\overline {r}}+\overline {u}_{\overline {z}}\right)-2\bar {\mu }^{(g)}\Big(\overline {u}_{\overline {r}}^{(g)}+\overline {w}_{\overline {z}}^{(g)}\overline {R}^2_{\overline {z}}\Big)-\bar {\mu }^{(g)}\overline {R}_{\overline {z}}\Big(\overline {w}^{(g)}_{\overline {r}}+\overline {u}^{(g)}_{\overline {z}}\Big)\nonumber\\ &\quad= \Big(\overline {p}-\overline {p}^{(g)}\Big)\Big(1+\overline {R}_{\overline {z}}^2\Big)+\bar {\sigma }\Big(1+\overline {R}^2_{\overline {z}}\Big)\left (\frac {1}{\overline {R}\Big(1+\overline {R}^2_{\overline {z}}\Big)^{1/2}}-\frac {\overline {R}_{\overline {zz}}}{\Big(1+\overline {R}^2_{\overline {z}}\Big)^{3/2}}\right )\!, \end{align}

with $\bar {\sigma }$ the surface tension and superscripts of $(g)$ denoting variables in the core gas flow.

The tube is assumed to consist of a thin, impermeable wall with density $\bar {\rho }_w$ , damping coefficient $\bar {\alpha }$ , circumferential tension $\bar {T}_{\theta }$ , longitudinal tension $\bar {T}_l$ , thickness $\bar {h}_0$ and dimensionless Poisson ratio $\gamma$ . Requiring equilibrium of normal forces at the wall yields the equation

(2.6) \begin{equation} \bar {\rho }\bar {h}_0\bar {\alpha } \bar {a}_{\bar {t}}\boldsymbol{\hat {r}}_{{w}}\boldsymbol{\cdot }\boldsymbol{\hat {n}}_{{w}}+\frac {\bar {T}_{\theta }}{\bar {a}}\Big (1+\bar {a}_{\bar {z}}^2\Big )^{-\frac {1}{2}}-\bar {T}_{\bar {\ell }}\bar {a}_{\bar {zz}}\Big (1+\bar {a}_{\bar {z}}^2\Big )^{-\frac {3}{2}}=-\boldsymbol{\hat {n}}_{{w}}\boldsymbol{\cdot }\boldsymbol{\bar {T}}\boldsymbol{\cdot }\boldsymbol{\hat {n}}_{{w}}, \end{equation}

with $\boldsymbol{\bar {T}}=-\bar {p}\boldsymbol{I}+\bar {\mu } (\boldsymbol{\nabla }\boldsymbol{\bar {u}}+\boldsymbol{\nabla }\boldsymbol{\bar {u}}^T )$ and where $\boldsymbol{\hat {r}}_w$ and $\boldsymbol{\hat {n}}_w$ are the unit radial and normal vectors at the tube wall. Equation (2.6) may be expressed as the following evolution equation for $\bar {a}$ :

(2.7) \begin{align} &\bar {\rho }_w\bar {h}_0\bar {\alpha } \bar {a}_{\bar {t}}\Big(1+\bar {a}_{\bar {z}}^2\Big)^{-1/2}+\frac {\bar {T}_0}{\bar {a}}\Big(1+\bar {a}_{\bar {z}}^2\Big)^{-1/2}+\frac {\bar {E}\bar {h}_0}{1-\gamma ^2}\frac {\bar {a}-\bar {a}_0}{\bar {a}^2}\Big(1+\bar {a}_{\bar {z}}^2\Big)^{-1/2}\nonumber \\ &\quad -\bar {T}_l\bar {a}_{\bar {z}\bar {z}}\Big(1+\bar {a}_{\bar {z}}^2\Big)^{-3/2}={\bar {p}-\bar {p}_{\textit{ext}}-\frac {2\bar {\mu }}{1+\bar {a}_{\bar {z}}^2}\Big [\bar {u}_{\bar {r}}+\bar {w}_{\bar {z}}\bar {a}_{\bar {z}}^2- (\bar {u}_{\bar {z}}+\bar {w}_{\bar {r}})\bar {a}_{\bar {z}}\Big ]}, \end{align}

where we have used the linear version ( $n=1$ ) of the strain law of Elad, Foux & Kivity (Reference Elad, Foux and Kivity1988; Elad, Kamm & Shapiro Reference Elad, Kamm and Shapiro1988) to set the circumferential tension $\bar {T}_{\theta }=\bar {T}_0+\bar {E}\bar {h}_0 (\bar {a}-\bar {a}_0)/[(1-\gamma )^2\bar {a}]$ with $\bar{E}$ the modulus of elasticity (similar to the approach of Halpern & Grotberg (Reference Halpern and Grotberg1992)). The first term in (2.7) represents wall damping. The circumferential tension term with coefficient $\bar {E}\bar {h}_0/(1-\gamma ^2)$ (which will be the primary one retained in what follows) provides a restoring force, acting to return the tube wall to its mean radius $\bar {a}_0$ . The longitudinal tension term (a diffusion-type term) provides a restoring force in the axial direction. The strength of each of these terms depends on the structural properties of the airways themselves as well as the intrapleural pressure and the surrounding soft parenchymal tissue – which itself exhibits viscoelastic properties – to which the airways are tethered. The magnitude of each of these terms can also be expected to be a function of airway generation and patient health; while precise values are not known for every airway, it is known that healthy (diseased) airways tend to be more compliant (stiff) and structurally stable (unstable) (Maghsoudi-Ganjeh, Sattari & Eskandari Reference Maghsoudi-Ganjeh, Sattari and Eskandari2021).

There is a steady solution with constant $\bar {a}=\bar {a}_0$ and $\bar {R}=\bar {R}_0$ given by

(2.8a) \begin{align}& \bar {w}_0=\frac {1}{4}\left (\frac {\bar {p}_{\bar {z}}^{(g)}+\bar {\rho }\bar {g}}{\bar {\mu }}\right )\left (\bar {r}^2-\bar {a}^2-2\bar {R}_0^2\ln \frac {\bar {r}}{\bar {a}}-\frac {2\bar {\varLambda }}{\bar {a}}\big (\bar {a}^2-\bar {R}_0^2\big )\right )+\bar {R}_0\bar {\tau }^{(g)}\left (\ln \frac {\bar {r}}{\bar {a}}-\frac {\bar {\varLambda }}{\bar {a}}\right )\!, \end{align}

(2.8b) \begin{align}&\qquad\qquad\qquad \bar {u}_0=0,\quad \quad \bar {p}_0=\bar {p}^{(g)}-\frac {\bar {\sigma }}{\bar {R}_0},\quad \quad \bar {p}_{\textit{ext}}=\bar {p}_0-\frac {\bar {T}_0}{\bar {a}_0}. \end{align}

2.2. Model derivation

A long-wave asymptotic model for the evolution of the free surface and tube wall will be derived next. Equations (2.1)–(2.7) may be made dimensionless using the following reference scales:

(2.9) \begin{equation} r=\frac {\bar {r}}{\bar {R}_0}, \quad z=\frac {\bar {z}}{\bar {\kappa }}, \quad u=\frac {\bar {u}}{\bar {U}_0}, \quad w=\frac {\bar {w}}{\bar {W}_0}, \quad t=\frac {\bar {t}\bar {W}_0}{\bar {\kappa }}, \quad p=\frac {\epsilon \bar {p}{\bar {R}_0}}{\bar {\mu }\bar {W}_0}, \quad \tau =\frac {\bar {\tau }\bar {R}_0}{\bar {\mu }\bar {W}_0}, \end{equation}

where $\bar {\kappa }$ is a typical wavelength of the perturbed free surface, $\epsilon =\bar {R}_0/\bar {\kappa }$ is an aspect ratio and where $\bar {W}_0=\bar {\sigma }/\bar {\mu }$ and $\bar {U}_0=\epsilon \bar {W}_0$ are reference velocity scales.

A long-wave approximation, $\epsilon \ll 1$ , will be employed to develop the model equation. Given that the focus here is on highly viscous films, a further assumption that the Reynolds number ${\textit{Re}}=\bar {\rho }\bar {W}_0\bar {R}_0/\bar {\mu }$ is $O(\epsilon )$ will be made. After substituting (2.9) into (2.1) and truncating at $O(\epsilon )$ , the dimensionless governing equations are given by

(2.10) \begin{align} p_r=0,\quad \quad p_z=\frac {1}{r}\partial _{r}(rw_r)-{\textit{Bo}},\quad \quad \frac {1}{r}\partial _r(ru)+w_z=0, \end{align}

where ${\textit{Bo}}=\bar {\rho }\bar {g}\bar {R}_0^2/\bar {\sigma }$ . At the tube wall, $r=a(z,t)$ , the truncated dimensionless boundary conditions are

(2.11) \begin{align} w&=-\varLambda w_r,\quad \quad u=a_t-\varLambda w_ra_z, \end{align}

where $\varLambda =\bar {\varLambda }/\bar {R}_0$ . At the free surface, $r=R(z,t)$ , the boundary conditions are

(2.12) \begin{align} u=R_t+\textit{wR}_z,\quad \quad w_r&=\tau ^{(g)},\quad \quad p-p^{(g)}=-\frac {\epsilon }{R}+\epsilon ^3R_{zz}. \end{align}

Note that one term of $O(\epsilon ^3)$ has been included in (2.12). Despite not being strictly valid to retain, this term is commonly kept in film flow models of this type as it has been shown to provide the correct cutoff wavenumber in linear stability analysis and produces growth rates in good agreement with the full governing equations.

The model equation for the dimensionless free surface $R(z,t)$ may be found by integrating the continuity equation from (2.1c ) across the fluid layer from $R(z,t)$ to $a(z,t)$ and applying the necessary boundary conditions from (2.11) and (2.12) to produce

(2.13) \begin{align} \frac {1}{2}\frac {\partial }{\partial t}\big(a^2-R^2\big)+\frac {\partial }{\partial z}\int _{R(z,t)}^{a(z,t)}rw\,{\rm d}r=0. \end{align}

The velocity profile $w$ may be found by solving the dimensionless $w$ -momentum equation subject to the remaining boundary conditions in (2.11) and (2.12), resulting in

(2.14) \begin{align} w&=\frac {1}{4}\left (p_z+{\textit{Bo}}\right )\left [r^2-a^2-2R^2\ln \frac {r}{a}-\frac {2\varLambda }{a}\big(a^2-R^2\big)\right ]+R\tau ^{(g)}\left (\ln \frac {r}{a}-\frac {\varLambda }{a}\right ); \end{align}

this is equivalent to (2.15) in Schwitzerlett et al. (Reference Schwitzerlett, Ogrosky and Topaloglu2023) up to a choice of velocity scales except that here the tube radius $a(z,t)$ is no longer a constant.

There are several ways in which one could incorporate the effects of pressure-driven airflow on the free surface. For small-amplitude disturbances, it is typical to specify a constant shear stress at the free surface; here, in order to study the evolution of disturbances beyond the linear regime, a variable shear stress will be considered. We elect to use what has been referred to in prior studies as a ‘locally Poiseuille’ approximation in which the airflow – presumed to be much faster than the highly viscous liquid flow – essentially experiences the free surface as a rigid wall. This allows for a decoupled treatment of the airflow in which the free-surface stresses $p_z^{(g)}$ and $\tau ^{(g)}$ may be estimated and incorporated into (2.14). Briefly, the air is assumed to flow at a constant volume flux $\bar {Q}^{(g)}$ through the core region (though a constant pressure gradient could also be incorporated). Modelling the (assumed laminar) airflow using Poiseuille flow with the long-wave assumption that the free surface varies slowly in $z$ produces the dimensionless local estimates of the stresses

(2.15) \begin{align} p_z^{(g)}=-\frac {4{\textit{Ca}}^{(g)}}{R^4},\quad \quad \tau ^{(g)}=\frac {2{\textit{Ca}}^{(g)}}{R^3}, \end{align}

where ${\textit{Ca}}^{(g)}$ is the capillary number, defined as ${\textit{Ca}}^{(g)}=2\bar {\mu }^{(g)}\bar {Q}^{(g)}\!/\!\pi \bar {R}_0^2\bar {\sigma }$ . Additional details of these calculations can be found in, e.g. Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012) and Schwitzerlett et al. (Reference Schwitzerlett, Ogrosky and Topaloglu2023). Note that fully coupled dynamics could also be used as in, e.g. Dietze & Ruyer-Quil (Reference Dietze and Ruyer-Quil2015).

Substituting (2.14) into (2.13) results in the final governing equation for the film’s free surface

(2.16) \begin{align} \frac {1}{2}\big(a^2-R^2\big)_t=\left [{\textit{Ca}}^{(g)}f_1(R;a)+\frac {1}{16}\big ({\textit{Bo}}\,R^2+R_z+R^2R_{zzz}\big )f_2(R,a)\right ]_z, \end{align}

where

(2.17a) \begin{align} f_1(R;a)&=\frac {1}{4R^4}\left [-a^4+2a^2R^2-R^4-4\varLambda a\big(a^2-R^2\big)\right ]\!, \end{align}

(2.17b) \begin{align} {f_2(R,a)}&=\frac {a^4}{R^2}+3R^2-4a^2+4R^2\ln \left (\frac {a}{R}\right )+\frac {4\varLambda }{aR^2}\big(a^2-R^2\big)^2, \end{align}

and where we have returned to the original aspect ratio by rescaling $z$ and $t$ by $\epsilon$ .

For the tube wall, substitution of (2.9) into (2.7) and truncating at $O(\epsilon )$ (and again returning to the original aspect ratio) gives

(2.18) \begin{align} \psi a_t=T_la_{zz}-\frac {1}{\varGamma }\frac {a-a_0}{a^2}+T_0\frac {a-a_0}{a_0a}+1-\frac {1}{R}+R_{zz}, \end{align}

where

(2.19) \begin{align} \psi =\frac {\bar {\rho }_w\bar {h}_0\bar {\alpha }\bar {R}_0}{\bar {\mu }},\quad T_0=\frac {\bar {T}_0}{\bar {\sigma }},\quad T_l=\frac {\bar {T}_l}{\bar {\sigma }},\quad \varGamma =\frac {(1-\gamma ^2)\bar {\sigma }}{\bar {E}\bar {h}_0}. \end{align}

The parameter $\psi$ is a damping parameter proportional to the product of the wall-to-fluid mass ratio and damping ratio, and $\varGamma$ is the ratio of surface tension forces to elastic forces. Both the damping and longitudinal tension terms are of $O(\epsilon ^2)$ (as well as the curvature term $R_{zz}$ ) but are retained here as in Halpern & Grotberg (Reference Halpern and Grotberg1992) and Zhou et al. (Reference Zhou, Peng, Zhang and Zhuge2016) in order to assess the impacts of flexibility on the film flow dynamics. Also note that the pressure $p_{\textit{ext}}$ is assumed to satisfy $(p-p_{\textit{ext}})/\epsilon =1-({1}/{R})+R_{zz}+( {T_0}/{a_0})$ .

Equations (2.16) and (2.19) comprise the model. In the case of no base flow ( ${\textit{Ca}}^{(g)}={\textit{Bo}}=0$ ) and no slip ( $\varLambda =0$ ), the model could be termed a long-wave version of the thin-film model studied by Halpern & Grotberg (Reference Halpern and Grotberg1992); we will explore the impact of base flow and slip on both thin and moderately thick films. In the case that ${\textit{Ca}}^{(g)}=\varLambda =0$ , the problem set-up is a Newtonian version of the problem studied by Zhou et al. (Reference Zhou, Peng, Zhang and Zhuge2016) (viscoelastic film flow down a flexible tube); their model started from the same equation for tube motion (2.7) but was derived using an integral boundary layer method. In the case of a rigid tube, so that $a=a_0$ (which could be accomplished by, e.g. taking $\psi \rightarrow \infty$ or $\varGamma \rightarrow 0$ ), the single-partial differential equation model of Schwitzerlett et al. (Reference Schwitzerlett, Ogrosky and Topaloglu2023) is recovered.

2.3. Parameter values

Here, we briefly discuss the range of parameter values which will be considered. There are a total of fifteen parameters to be specified: eight describing the tube: $\bar {T}_l$ , $\bar {T}_0$ , $\bar {E}$ , $\bar {\alpha }$ , $\bar {\rho }_w$ , $\bar {h}_0$ , $\bar {a}_0$ , $\bar {\gamma }$ ; three parameters describing the fluid: $\bar {\mu }$ , $\bar {\sigma }$ , $\bar {\rho }$ ; two parameters describing the airflow: $\bar {\mu }^{(g)}$ , $\bar {Q}^{(g)}$ ; and two corresponding to the film thickness and boundary condition: $\bar {R}_0$ and $\bar {\varLambda }$ , respectively. The final model contains eight dimensionless parameters – ${\textit{Ca}}^{(g)}$ , ${\textit{Bo}}$ , $a_0$ , $\psi$ , $T_0$ , $T_l$ , $\varGamma$ , $\varLambda$ – that govern the film flow; this reduction in number of parameters is not only the result of non-dimensionalisation but also simplifying assumptions (e.g. small $Re$ , ‘locally Poiseuille’ airflow) and the way in which the original parameters were grouped in the tube wall equation (2.7). Our parameter choices will largely fall in a range roughly corresponding to those considered by Halpern & Grotberg (Reference Halpern and Grotberg1992) and Zhou et al. (Reference Zhou, Peng, Zhang and Zhuge2016). Parameter values that motivated Halpern & Grotberg (Reference Halpern and Grotberg1992) for the tube wall were, e.g. $\bar {E}=60\,000$ dyn cm⁻ $^2$ , $\bar {\rho }_w=1$ g cm^– $^3$ , $\bar {T}_l=25$ dyn cm^–1, $\gamma =0.5$ , $\bar {T}_0=0$ dyn/cm, $\bar {\alpha }=10$  s $^{-1}$ , $\bar {h}_0=0.0025$ cm; for the film, we will take $\bar {\rho }=1$ g cm^– $^3$ and surface tension $\bar {\sigma }=20$ dyn cm^–1. Similar to Halpern & Grotberg (Reference Halpern and Grotberg1992), a range of dimensionless parameter values – which is motivated by the above choices but also recognises the wide variety of tube and fluid properties possible – will be considered: $T_l\in [0,1]$ , $T_0=0$ and $\varGamma \in [0,1]$ .

As the focus here is on situations where gravity may play a non-negligible role in the film’s evolution, we will consider a range of tube radii and film thicknesses: $\bar {a}_0\in [0.025,0.5]$ cm, $\bar {R}_0\in [0.0225,0.4]$ cm; here the lower bounds correspond to values considered by Halpern & Grotberg (Reference Halpern and Grotberg1992) appropriate for the bronchioles in which gravity may be neglected, and the upper end of the range reasonable for upper airways. As the viscosity of mucus is a property that varies widely (especially in the presence of shear stress created by airflow considered here), a wide range of viscosities may be appropriate, e.g. $\bar {\mu }\in [0.01,100]$ P. Given these uncertainties, a wide range of dimensionless parameters will be used: $a_0\in [1.05,1.5]$ , ${\textit{Bo}}\in [0,10]$ , $\psi \in [0,10^5]$ , although values outside these ranges may be considered briefly in order to make a point about the model solutions.

For the airflow, $\bar {\mu }^{(g)}=1.81\times 10^{-4}$ ; in previous experiments designed to model upper airways, volume fluxes of $\bar {Q}^{(g)}\in [10\,1000]$ cm $^3$ s^–1 were used Kim et al. (Reference Kim, Rodriguez, Eldridge and Sackner1986), Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012), Kim, Iglesias & Sackner (Reference Kim, Iglesias and Sackner1987). For lower airways this value may be expected to be much smaller. In an effort to capture the range of values possible, we will consider ${\textit{Ca}}^{(g)}\in [0,1]$ . Finally, for slip length, (Brochard-Wyart et al. Reference Brochard-Wyart, de Gennes, Hervert and Redon1994) note that a slip length of $1{-}10\,\unicode{x03BC}$ m may be appropriate for silicone oils in experiments like those of Camassa et al. (Reference Camassa, Forest, Lee, Ogrosky and Olander2012), leading to values of $\varLambda$ as high as 0.01; other modelling studies with slip have investigated higher dimensionless slip lengths (e.g. Liu & Ding (Reference Liu and Ding2017), who study up to $\varLambda \approx 0.1$ using the scalings here). We note that apparent wall slip depends not only on fluid and polymer properties but on properties of the tube wall, and that while slip boundary conditions have been found to aid in addressing discrepancies between models and experiments in a wide variety of problems, there is some debate over the validity of their use in some settings (see, e.g. Lauga et al. Reference Lauga, Brenner, Stone, Tropea, Yarin and Foss2005; Bäumchen & Jacobs Reference Bäumchen and Jacobs2009; Kavokine, Netz & Bocquet Reference Kavokine, Netz and Bocquet2021 and references therein). Motivated by these previous studies, we will consider $\varLambda \in [0,0.1]$ .

A set of reference parameter values is given in table 1 for § 3, and a second set is given in table 2 for § 4.3.

Table 1. Reference values used in § 3.

Table 2. Standard case values used in § 4.3.

3.1. No base flow

In the case of no base flow ( ${\textit{Ca}}^{(g)}={\textit{Bo}}=0$ ), the growth rate without damping is

(3.4) \begin{align} i\omega =\Big[g_5\big(k^4-k^2\big)\big(1+k^2T_l\chi a_0^2\big)\Big]\Big /\Big[\chi a_0^3\big(1-k^2\big)-1-k^2T_l\chi a_0^2\Big]. \end{align}

Note that in the rigid tube case ( $\chi =0$ ) these growth rates are $i\omega =g_5(k^2-k^4)$ , leading to the wavenumber of maximum growth rate $k_{\textit{max}}=1/\sqrt {2}$ . In the other extreme, as noted by Halpern & Grotberg (Reference Halpern and Grotberg1992), if $\chi$ is sufficiently large, here $\chi \gt 1/a_0^3$ , then there is an infinite growth rate for finite $k$ given by

(3.5) \begin{align} k_{\infty }=\left (1-\frac {1}{\chi a_0^3}\right )\bigg /\left (1+\frac {T_l}{a_0}\right )\!. \end{align}

Figure 2. (a) Value of $\varGamma =1/a_0^3$ (solid line) and smallest $\varGamma$ for which $k_{\textit{max}}=0$ (dashed line) as a function of $a_0$ . The $\times$ symbols correspond to growth rate plots in (b). (b)–(c) Growth rates for a variety of $\varGamma$ values; (b) $\varLambda =0$ ; (c) $\varLambda =0.02$ . Unless otherwise mentioned, parameter values are the reference values in table 1; thick red line denotes reference value growth rates. Thick black line denotes the rigid tube case $(\varGamma =0$ ). Dotted black line denotes maximum growth rate with $\varLambda =0$ as $\varGamma$ varies.

In the case of damping, with $\psi \gt 0$ , there are two solutions to (3.1); one of these is unstable for $0\lt k\lt 1$ , with $k_c=1$ the cutoff wavenumber. Note that at $k=0$ , these solutions are

(3.6) \begin{align} i\omega =0,\,\,\frac {1}{\psi }\left (\frac {1}{\chi a_0^2}-a_0\right )\!, \end{align}

so that for $\chi \gt 1/a_0^3$ the $k=0$ mode is also unstable. This can lead to the ‘compliant collapse’ scenario observed by Halpern & Grotberg (Reference Halpern and Grotberg1992) in simulations of a thin-film model, in which the Plateau–Rayleigh instability becomes subdominant to this wall collapse in driving the dynamics. In their thin-film model, $\varGamma \gt 1$ was required for the $k=0$ mode to be unstable, although they noted that, in nonlinear simulations, smaller values of $\varGamma$ were sufficient to trigger this compliant collapse.

Here, film thickness is also seen to lower the value of $\varGamma$ required for the $k=0$ mode to be linearly unstable. Figure 2(a) shows the region of $a_0$ – $\varGamma$ space (i) for which the growth rate at $k=0$ is positive, and (ii) for which the growth of $\varGamma$ is largest at $k=0$ , i.e. $k_{\textit{max}}=0$ . Figure 2(b) shows the growth rates for $a_0=1.1$ and several values of $\varGamma$ corresponding to the $\times$ symbols in (a). Increasing $\varGamma$ increases growth rates for all $k$ and decreases the wavenumber of maximum growth rate $k_{\textit{max}}$ which eventually reaches 0 for sufficiently large $\varGamma$ . Note that the curve with $\varGamma =0$ represents the rigid tube growth rate curve, while the red curve represents the parameter values in table 1.

Figure 2(c) shows growth rates with identical parameters to (b) but with $\varLambda =0.02$ . While slip leaves the $k=0$ growth rates unchanged, it has the impact of promoting the Plateau–Rayleigh instability, leading to larger values of $k_{\textit{max}}$ for all $\varGamma$ shown here.

How does the magnitude of wall damping parameter $\psi$ affect the growth rates? Figure 3(a) shows the growth rate curves for various values of $\psi$ , with $\varGamma =0.5$ . Increasing $\psi$ decreases the growth rates and increases $k_{\textit{max}}$ . In the limit $\psi \rightarrow \infty$ , the growth rates approach those of the rigid tube. In the limit $\psi \rightarrow 0$ , the growth rates approach some curve with reduced $k_{\textit{max}}$ and a maximum growth rate roughly 40 % larger than that of the rigid tube case for these parameter values. Note that, if $\chi \gt 1/a_0^3$ , $k_{\textit{max}}$ may be again reduced all the way to zero (not shown here). The impact of $T_l$ is shown in figure 3(b). Increasing $T_l$ has a similar impact on growth rates as increasing $\psi$ does, with the rigid tube approached as $T_l\rightarrow \infty$ .

In rigid tubes, slip has been shown previously to increase growth rates in film flow models by allowing liquid to flow more readily into a growing wave, promoting the Plateau–Rayleigh instability (e.g. Halpern & Wei Reference Halpern and Wei2017; Liu & Ding Reference Liu and Ding2017), as seen earlier in figure 2(c). This is also the case for fixed $\varGamma =0.5$ in figure 3(c), in which increasing the slip length $\varLambda$ results in increased growth rates without modifying $k_{\textit{max}}$ . In the case where $\varGamma \gt 1/a_0^3$ , slip does modify $k_{\textit{max}}$ (not shown here), with the value of $k_{\textit{max}}$ closer to its rigid tube value of $1/\sqrt {2}$ the larger the values of slip lengths.

Figure 3. Growth rates for various parameter values. Unless otherwise mentioned, parameter values are the reference values in table 1; thick red lines denote reference value growth rates. Thick black lines denote the rigid tube case $(\varGamma =0$ ). Growth rates are shown for various (a) $\psi$ , (b) $T_l$ and (c) $\varLambda$ .

3.2. Base flow

Next, the impact of base flow, i.e. with either ${\textit{Ca}}^{(g)}\gt 0$ or ${\textit{Bo}}\gt 0$ , on the linear stability is examined. Note that, in the absence of damping, base flow does not change the growth rates, as positive values for ${\textit{Ca}}^{(g)}$ and ${\textit{Bo}}$ only modify the imaginary part of the solution $i\omega =-\tilde {c}/\tilde {b}$ , creating non-zero phase speed.

Figure 4(a) shows the impact of varying ${\textit{Bo}}$ on growth rates in the presence of damping ( $\psi \gt 0$ ) for $\chi \lt 1/a_0^3$ . Increasing ${\textit{Bo}}$ decreases the maximum growth rates, while increasing growth rates for the smallest values of $k$ . The wavenumber of maximum growth rate $k_{\textit{max}}$ increases with increasing ${\textit{Bo}}$ .

For $\chi \gt 1/a_0^3$ , varying ${\textit{Bo}}$ has a more complicated impact on growth rates. As shown in figure 4(b), for sufficiently large ${\textit{Bo}}$ , there is a second unstable mode (e.g. ${\textit{Bo}}=0.075$ ); if ${\textit{Bo}}$ is large enough, these two branches of the growth rate curve meet and pinch off (e.g. ${\textit{Bo}}=0.1$ ). Also note that the value of $k_{\textit{max}}$ changes abruptly from a value near 0.5 to $k_{\textit{max}}=0$ as ${\textit{Bo}}$ increases from 0 to 0.075.

Figure 4. (a)–(b) Growth rates are shown for various ${\textit{Bo}}$ : (a) $\varGamma =0.5$ ; (b) $\varGamma =1$ . (c)–(d) Phase speed for various values of ${\textit{Bo}}$ and ${\textit{Ca}}^{(g)}$ ; thick red lines denote speeds with parameter values in table 1. (e) Profile of free surface and wall for $k$ and ${\textit{Bo}}$ values corresponding to dots in (a), (c). Note that the magnitude of disturbances is arbitrary, although relative magnitude is meaningful, determined by the eigenvalues and eigenvectors. Unless otherwise mentioned, parameter values are those of table 1.

The presence of base flow also introduces a non-zero phase speed, shown in figure 4(c,d), for various values of ${\textit{Bo}}$ and ${\textit{Ca}}^{(g)}$ , respectively. In figure 4(c), the phase speed shows dependence on $k$ only for very small $k$ . As expected, increasing ${\textit{Bo}}$ results in waves falling down the tube faster. Similarly, as ${\textit{Ca}}^{(g)}$ increases in figure 4(d), the phase speed increases as free-surface waves move up the tube. Again only weak dependence on $k$ is seen.

The profile of the wall and free surface are shown in figure 4(e) for solutions corresponding to the dots in (a), (c). In all cases the disturbance to the free surface is larger than the disturbance to the wall. Increasing ${\textit{Bo}}$ slows the growth in the free-surface disturbance and, more noticeably, in the wall disturbance. A phase shift between free surface and wall is introduced by increasing ${\textit{Bo}}$ , with the free-surface disturbance leading the wall disturbance.

We close this section by mentioning that increasing the film thickness increases growth rates.

4.1. Numerical methods

Equations (2.16) and (2.19) are solved on a domain of length $L$ using periodic boundary conditions. A pseudospectral method was used, in which spatial derivatives were calculated in Fourier space and nonlinearities were calculated in physical space. An explicit second-order predictor–corrector method for time integration was used. As the quantity $a^2-R^2$ is a conserved quantity, conservation of film volume was monitored throughout each simulation; when the change in film volume between successive time steps exceeded some chosen threshold (set to be 10 $^{-8}$ or smaller in all simulations), the update to the free surface and tube wall were recalculated with the time step decreased by a factor of two. Initial conditions consisted of a constant solution plus a series of small-amplitude perturbations (with random amplitude and phase shift) with varying wavenumber

(4.1) \begin{align} R=1+\sum _{k=1}^N\tilde {R}_k\cos (2\pi kz/L+c_{1,k}),\quad a=a_0+\sum _{k=1}^N\tilde {a}_k\cos (2\pi kz/L+c_{2,k}). \end{align}

Travelling-wave solutions were also found by seeking steady solutions in a moving reference frame with speed $c$ , i.e. seeking functions $Q(Z)=R(z-ct)$ and $b(Z)=a(z-ct)$ , where $Z=z-ct$ . Substituting this ansatz into the model equations results in a fourth-order ordinary differential equation (ODE) and a second-order ODE. The equation for $Q$ may be integrated once, yielding a third-order equation (with a constant of integration $K$ ); the resulting equations are

(4.2a) \begin{align} K&=\frac {c}{2}\big (b^2-Q^2\big )+{\textit{Ca}}^{(g)}f_1(Q,b)+\frac {1}{16}\big ({\textit{Bo}}\,Q^2+Q'+Q^2Q^{\prime\prime\prime}\big )f_2(Q,b), \\[-12pt] \nonumber \end{align}

(4.2b) \begin{align} -c\psi a'&=T_la^{\prime\prime}-\frac {1}{\varGamma }\frac {a-a_0}{a^2}+T_0\frac {a-a_0}{a_0a}+\left (1-\frac {1}{Q}+Q^{\prime\prime}\right )\!. \\[9pt] \nonumber \end{align}

Equation (4.2) were solved by first identifying an equilibrium solution $Q=Q_0$ and $b=b_0$ . By varying the wave speed, a Hopf bifurcation in the solution family may be identified. The numerical continuation package XPP/AUTO was used for this task; XPP is a numerical package for analysing dynamical systems (Ermentrout Reference Ermentrout2002) that provides a convenient way to access AUTO (Doedel et al. Reference Doedel, Champneys, Dercole, Fairgrieve, Kuznetsov, Oldeman, Paffenroth, Sandstede, Wang and Zhang2008). Moving onto a branch of periodic solutions from the Hopf bifurcation, families of travelling-wave solutions may be found while varying any model parameter, including period size. The mean value of $b^2-Q^2$ was held fixed throughout these continuations. Details of this numerical approach can be found in Camassa et al. (Reference Camassa, Marzuola, Ogrosky and Vaughn2016), including the addition (and subsequent removal) of a small viscosity-like term to the model. This aids in finding Hopf bifurcations by avoiding zero-Hopf bifurcations which can be more difficult to identify numerically; this viscosity term was set to zero after moving onto a branch of periodic solutions.

4.2. No base flow

For the case of no base flow, i.e. $Ca^{(g)}=Bo=0$ , simulations were first conducted for a period of $L\approx 2.86\pi$ . This choice for the domain matches that used by Halpern & Grotberg (Reference Halpern and Grotberg1992) as linear stability analysis shows it to be the most unstable wave. As has been discussed in previous studies, e.g. Hammond (Reference Hammond1983), Halpern & Grotberg (Reference Halpern and Grotberg1992), when the film is thinner than some critical thickness, the growth in amplitude of a wave saturates with the free surface approaching some steady-state profile. When the film is thicker than this critical thickness, enough fluid drains into the wave from the surrounding substrate to allow the wave to undergo accelerated growth with $\min _z R\rightarrow 0$ in finite time.

Figure 5 shows the values of $\min _zR(z,t)$ and $\min _za(z,t)$ for three simulations; one for a rigid tube ( $\varGamma =0$ , $\varLambda =0$ ), one for a flexible tube ( $\varGamma =0.5$ , $\varLambda =0$ ) and one for a flexible tube with slip ( $\varGamma =0.5$ , $\varLambda =0.02$ ). The beginning stages of this accelerated growth are evident for the flexible tube solutions, while in the rigid tube case the film thickness parameter $a_0=1.1$ is (barely) sufficiently small to prevent the formation of a plug, with $\min _zR(z,t)$ approaching a constant value.

Figure 5. Evolution of $\min _z R(z,t)$ and $\min _z a(z,t)$ for three simulations with $a_0=1.1$ , $\psi =16,000$ , ${\textit{Bo}}={\textit{Ca}}^{(g)}=T_l=T_0=0$ .

Snapshots of the free surface are shown for these three simulations in figure 6. In the second case, flexibility enhances the Plateau–Rayleigh instability, resulting in faster growth and the eventual formation of a plug. The inclusion of slip further promotes plug formation, with a plug forming quicker than in the no-slip case. The tube wall deformation is slightly less in the presence of slip.

Figure 6. (a)–(c) Snapshots of the free-surface profile and tube wall corresponding to the $\varGamma =0$ , $\varLambda =0$ case from figure 5. (d)–( f) Snapshots corresponding to the $\varGamma =0.5$ , $\varLambda =0$ case. (g)–(i) Snapshots corresponding to the $\varGamma =0.5$ , $\varLambda =0.02$ case.

Figure 7. (a) Closure time $t_c$ as a function of $\varGamma$ for various $\psi$ and $\varLambda$ ; $a_0=1.1$ , ${\textit{Ca}}^{(g)}={\textit{Bo}}=T_l=T_0=0$ . (b) Closure time as a function of $a_0$ for various $\varGamma$ , $\psi$ and $\varLambda$ ; ${\textit{Ca}}^{(g)}={\textit{Bo}}=T_l=T_0=0$ . (c) Ratio of closure time with slip to closure time without slip for two values of $\psi$ as a function of $\varGamma$ ; marker symbols correspond to those in (a). (d) Ratio of closure time with slip and/or flexibility to closure time of rigid tube with no slip; marker symbols correspond to those in (b).

Figure 7(a) shows the closure time (here approximated as the time needed for $\min R$ to reach 0.55) for a variety of $\varGamma$ , $\psi$ and $\varLambda$ values. The chosen values of $\psi$ and $\varGamma$ produce the same phenomenon observed by Halpern & Grotberg (Reference Halpern and Grotberg1992) in which the closure time decreases with decreasing $\varGamma$ ; for large damping ( $\psi =$ 16 000), this decrease is particularly pronounced for values of $\varGamma$ between 0.5 and 1, which corresponds to $\varGamma$ crossing the $1/a_0^3$ threshold. The impact of slip is seen as well; for small damping and large $\varGamma$ , the impact of slip is minimal as flexibility provides the main enhancement to the Plateau–Rayleigh instability, leading to plug formation. For large damping or for small $\varGamma$ – i.e. for tubes that are ‘weakly’ flexible – the contribution of slip is more significant. For each value of $\psi$ , figure 7(c) shows the ratio of closure times found with slip to no-slip. For ‘weakly’ flexible tubes, slip can cut the closure time by a factor of 4 or more. We note that a threshold value is a common way in thin-film modelling studies to estimate the closure time (see, e.g. Cassidy et al. (Reference Cassidy, Halpern, Ressler and Grotberg1999)) as the model cannot be run all the way to $\min R=0$ due to the inverse powers and logarithms of $R$ in (2.16). While the chosen threshold value of $\min {R}=0.55$ results in a slight underestimate of the closure time, extended simulations run for some cases show that every time $\min {R}$ decreases below 0.55, a plug is formed very soon afterwards.

How does the film thickness parameter affect the closure time? Figure 7(b) shows the closure time for a variety of $a$ values and several combinations of $\varLambda$ and $\varGamma$ . It appears that for thin films, closure time is more sensitive to the change in $\varGamma$ (flexibility), while for thicker films, the closure time is more sensitive to the change in $\varLambda$ (slip). This is confirmed in figure 7(d), which shows the ratio of each of the three curves with markers in panel (b) (namely rigid tube with slip, flexible tube with slip and flexible tube without slip) to the closure time found for a rigid tube with no slip. In all cases, the greatest drop in closure time occurs when both slip and flexibility are present, but the relative impact of these two traits depends on the film thickness.

4.3. Base flow

The case of base flow due to gravity is examined next. In the case of a rigid tube, it has been shown that base flow due to gravity, although not contributing to linear instability, provides a nonlinearly stabilising force that inhibits plug formation and increases the critical thickness required for plug formation (Dietze et al. Reference Dietze, Lavalle and Ruyer-Quil2020; Ogrosky, Reference Ogrosky2021a ; Camassa, Ogrosky & Olander Reference Camassa, Ogrosky and Olander2024). How do tube flexibility and slip impact the evolution of falling films?

Figure 8. (a) Time snapshots of $h=a_0-R$ in a solution to model equations (2.16) and (2.19) with domain length $L=24\pi$ and with standard base flow parameter values in table 2; gravity acts right to left. Snapshots shown every $\Delta t\approx 7.2$ in a frame of reference moving with the linearly most unstable wavenumber. (b) Snapshot in a tube corresponding to final snapshot in (a), zoomed in around largest wave; plug forms prior to next snapshot.

Simulations for a variety of parameter value combinations were run; the ‘standard’ case parameter values used throughout this section are shown in table 2. Figure 8(a) shows a waterfall plot of the free surface with these values; as it turns out, these values correspond to a critical thickness very near $a_0=1.12$ . The film settles into a series of waves which fall down the tube; wave coalescence occurs between several of the waves, with the largest of the waves eventually becoming large enough that during a wave merger, $R\rightarrow 0$ in finite time, indicating the formation of a plug.

How does each parameter – such as $\varGamma$ , $\psi$ , $\varLambda$ and $T_l$ – impact whether plugs form, and the time until a plug is formed? Model solutions were found for the standard case in table 2, as well as cases in which a single parameter value is changed from those standard values. Figure 9 shows the evolution of $\min _zR(z,t)$ and $\min _za(z,t)$ for each case with a domain length of $4\pi$ , which is of sufficient width to produce a single wave in each case. Identical initial conditions were used in each simulation.

The only solution in which a plug did not form was the rigid tube ( $\varGamma =0$ ) solution; $a=a_0$ is a constant in this case, and $\min _zR(z,t)$ approaches a fixed value. Of the flexible tube simulations, the standard case took the greatest amount of time to produce a plug; increasing $T_l$ to one had little impact on the solution. Increasing $\varGamma$ to one – past the value $1/a_0^3$ – resulted in a plug forming faster due to the $k=0$ mode of the wall decreasing. Increasing slip reduces the plug formation time considerably, while resulting in the smallest wall deviations. Decreasing damping results in the greatest deviation of $\min _za(z,t)$ from $a_0$ , and plugs form fastest in this case.

Figure 9. Evolution of $\min _z R(z,t)$ and $\min _z a(z,t)$ for simulations with base flow and domain $L=4\pi$ . Standard case (thick cyan line) refers to values in table 2. All other cases have identical parameter values except those indicated by text labels.

In order to assess the impact of varying each parameter on the critical thickness required for plug formation, simulations could be repeatedly run for each parameter value with a variety of film thicknesses known to be near the critical thickness; the critical thickness would then be in a range between the largest thickness for which plugs did not form and the smallest thickness for which plugs formed. Based on figure 9, $a_0=1.15$ appears to be below the critical thickness for the rigid tube case, but above the critical thickness for all other cases considered here.

Figure 10. Travelling-wave solution families with period $L=4\pi$ . Standard case (thick cyan line) refers to values in table 2. All other cases have identical parameter values except those indicated by text labels. Colours correspond to parameter values in figure 9.

A second way to assess the impact each parameter has on this critical thickness is to examine turning points in families of travelling-wave solutions, as in, e.g. Camassa et al. (Reference Camassa, Ogrosky and Olander2014), Ding et al. (Reference Ding, Liu, Liu and Yang2019), Dietze et al. (Reference Dietze, Lavalle and Ruyer-Quil2020) and Dietze (Reference Dietze2024). Figure 10 shows families of travelling-wave solutions for each of the cases shown in figure 9 (except $\varGamma =1$ ) as a function of $a_0$ . Each point on a solid curve represents a travelling-wave solution with $\min _ZR(Z)$ ; the corresponding value of $\min _Za(Z)$ lies on the corresponding dashed curve. Each family contains a turning point at some value of $a_0$ , in which a lower branch and upper branch merge. Upper branch solutions in the rigid tube case are unstable with a single large positive real eigenvalue (Camassa et al. Reference Camassa, Marzuola, Ogrosky and Vaughn2016); waves with these amplitudes in all cases are never seen in solutions to the evolution equations, and are not explored further here. As has been shown in previous studies the value of $a_{0,c}$ corresponding to each turning point provides a proxy for the critical thickness required for plugs to form. Note that $a_0=1.15$ lies to the left of this turning-point thickness for the rigid tube, but to the right of all other cases, consistent with figure 9; it is only slightly to the right of the standard and $T_l=1$ cases, helping to explain the long time for plug formation to occur in those cases.

There is no deflection in the wall for the rigid tube case, as shown by the straight line for $\min _Za(Z)$ in the $\varGamma =0$ case. The largest deflection again occurs for the small-damping case $\psi =16$ . The amplitude of the turning-point wave is much less in all flexible tube cases than in the rigid tube case.

Figure 11. Turning-point thickness dependence on (a) $\psi$ , (b) $\varLambda$ and (c) $\varGamma$ . All other parameter values correspond to the standard case in table 2. The $\times$ symbols denote parameter values in figure 9. (c) Shading denotes region where $\varGamma \gt 1/a_0^3$ .

One benefit of this approach to identifying the critical thickness is that it is possible to explore the dependence of $a_{0,c}$ on parameter values very quickly. Figure 11 shows this dependence of $a_{0,c}$ on $\psi$ , $\varLambda$ and $\varGamma$ . This proxy thickness $a_{0,c}$ decreases with increasing slip ( $\varLambda$ ) and increasing flexibility ( $\varGamma$ ). As was noted in Schwitzerlett et al. (Reference Schwitzerlett, Ogrosky and Topaloglu2023), there appears to be a limiting value of $a_{0,c}\gt 1$ as $\varLambda$ gets arbitrarily large. The dependence of $a_{0,c}$ on $\psi$ is non-monotonic, with some minimal value for $a_{0,c}$ obtained at finite $\psi$ (roughly $\psi \approx 44$ in figure 11 a).

One caveat is that in the case of large $\varGamma$ , in which $\varGamma \gt 1/a_0^3$ , the travelling-wave solutions seen in figure 11(c) are not observed in solutions to the evolution equations. Instead, growth of the $k=0$ mode overtakes all other features, as described by Halpern & Grotberg (Reference Halpern and Grotberg1992) for the ‘compliant collapse’ case. These turning points thus cannot necessarily be taken as a proxy for the critical thickness.

Figure 12. Travelling-wave solutions corresponding to $\times$ symbols in figure 11.

The waves corresponding to each turning point in figure 10 (also denoted by $\times$ symbols in figure 11) are shown in figure 12. Note the shape of the wall profile in each case: the rigid tube case in (a) has no displacement, while the standard case (b) and slip case ( f) show very little deviation from $a_0$ other than a small $k=0$ perturbation. In the large $\varGamma$ case, there is marked deviation in the $k=0$ mode of $a$ , along with asymmetric displacement of $a$ in the wave support region. In the small-damping case (e), the wall exhibits an asymmetric ripple in the wall at the leading edge of the wave; for moderate damping (g) a combination of the asymmetry seen in (c) and (e) is visible. The free surface in (a), (b), (c) and ( f) displays a capillary ripple on the leading edge of the wave, as is typical in falling film flows. In (e) and (g), however, there is no such ripple visible in the free surface, although one is visible in the wall profile. The free surface (wall) profiles are plotted on top of one another in figures 12(d) and 12(h); the reduced amplitude (relative to the rigid tube case) is visible in all cases, and the $k=0$ mode deviations are apparent in (h).

Figure 13. Turning-point thickness dependence on ${\textit{Bo}}$ for various values of $\psi$ and $\varGamma$ .

How does the value of ${\textit{Bo}}$ impact $a_{0,c}$ in a flexible tube? In a rigid tube, the presence of base flow due to gravity plays a nonlinearly stabilising role, suppressing plug formation and increasing $a_{0,c}$ .

The impact of ${\textit{Bo}}$ on $a_{0,c}$ in a flexible tube is explored in figure 13. As with the rigid tube case, increasing ${\textit{Bo}}$ results in larger $a_{0,c}$ , inhibiting plug formation. In the case of strong damping ( $\psi =16\,000$ ), the impact of increasing $\varGamma$ has a moderate impact for small ${\textit{Bo}}$ ( ${\textit{Bo}}\lessapprox 0.02$ for $\varGamma =0.5$ ), perhaps the least impact for a ‘transition’ region $0.02\lessapprox {\textit{Bo}}\lessapprox 0.2$ (for $\varGamma =0.5$ ) and the most significant impact for higher ${\textit{Bo}}$ ( ${\textit{Bo}}\gtrapprox 0.2$ ). For weak damping, there still appears to be a similar type of transition region, $0.6\lessapprox {\textit{Bo}}\lessapprox 2$ . For both small or large damping, the values of $a_{0,c}$ are nearly identical in the case of either very small or very large ${\textit{Bo}}$ ; to say the same thing another way, the damping has the most noticeable effect on $a_{0,c}$ for intermediate values of ${\textit{Bo}}$ .

It was shown in Ogrosky (Reference Ogrosky2021a ) and Camassa et al. (Reference Camassa, Ogrosky and Olander2024) that the critical thickness can be made arbitrarily thick by increasing the value of ${\textit{Bo}}$ far enough. This is not the case for flexible tubes; there appears to be some $a_{0,c,\infty }$ that $a_{0,c}$ approaches as ${\textit{Bo}}$ gets arbitrarily large, with this thickness (horizontal asymptote) decreasing with $\varGamma$ .

To summarise, figures 11 and 13 demonstrate the efficiency of this alternate way of identifying the critical thickness required for plug formation to occur. In the traditional approach, one would repeatedly run simulations for a single set of parameter values to sufficiently narrow down the region in which $a_{0,c}$ lies. This would then need to be repeated, varying each parameter one at a time in small increments, to provide figures like figures 11 and 13. The alternate approach provides a shortcut; through numerical continuation of turning points, each of the curves in these plots can be found in a matter of seconds or minutes.

One caveat, however, is that all of the simulations above were conducted with a domain of $L=4\pi$ . How accurate is the critical thickness identified by turning points in predicting plug formation in larger domains, in which plugs may not only form through uninhibited growth of a single wave, but also through the merger of two waves?

Figure 14 shows the values of $\min _zR(z,t)$ and $\min _za(z,t)$ for the same cases as in figure 9, but with domain length $L=24\pi$ and thickness $a_0=1.12$ . Each simulation again began with identical initial conditions. In the rigid tube case ( $\varGamma =0$ ), plugs were not seen to form in the simulation, with wave growth saturating and wave trains progressing down the tube.

Figure 14. Evolution of $\min _z R(z,t)$ and $\min _z a(z,t)$ for simulations with base flow. Standard case (cyan) parameters are those of table 2 except that $a_0=1.12$ ; domain length $L=24\pi$ . Other cases have all parameters identical to the standard case, with the exception of $a_0=1.12$ and the parameter value indicated by the text.

On the other extreme, the case of decreasing damping ( $\psi =1$ ) provides the only case in which a plug formed through growth of an individual wave, rather than through wave mergers. For the remaining cases, wave mergers occur in the simulation, with most mergers marked by a sharp drop in $\min _zR(z,t)$ , followed by a slow partial rebounding of $\min _zR(z,t)$ to some new quasi-steady value (seen most clearly in the first few mergers in the standard case). But for all cases except the rigid tube and small-damping case, a plug eventually formed through the coalescence of two waves, with the standard case taking the longest amount of time for a plug to form. This possibility of plug formation through wave mergers in longer tubes results in a decrease in the film thickness required for plugs to form.

For both the $\psi =16$ and $\varGamma =1$ cases, the tube wall pinched in noticeably (i.e. $\min _za(z,t)$ decreased noticeably prior to plug formation), demonstrating how increased flexibility and decreased damping promote the Plateau–Rayleigh instability through tube constriction. Increasing slip ( $\varLambda =0.02$ ) and longitudinal tension ( $T_l=1$ ) also promoted faster plug formation than the standard case, with slip noticeably promoting growth in the waves at early times, consistent with the increased growth rates that slip produces; in both these cases, wall deformations remained minimal.

It is interesting to note that the standard and $T_l=1$ cases are virtually identical until $t\approx375$ , at which point a wave merger occurs in the $T_l=1$ case, but does not for the standard case. Note that the $T_l$ case leads to significantly shorter closure time than the standard case, in contrast to the $L=4\pi$ simulations. Also note that the $\varGamma =1$ case again corresponds to a case of ‘compliant collapse’ in the tube wall, seen by the significantly reduced value of $\min _za(z,t)$ .

Figure 15. Evolution of (a) average volume flux $\tilde {Q}(z,t)=({1}/{L})\int _0^LQ(z,t)\,{\rm d}z$ and (b) $\min _z R(z,t)$ for simulations with base flow due to core flow; domain length $L=24\pi$ , ${\textit{Ca}}^{(g)}=0.00625$ , $a_0=1.1$ , ${\textit{Bo}}=T_l=T_0=0$ .

Before concluding, we briefly explore the case of core-driven film flow in which ${\textit{Ca}}^{(g)}\gt 0$ . Much of the dynamics is similar to that of ${\textit{Bo}}\gt 0$ : long-wave growth saturates in a series of waves which interact and propagate in the direction of $Q^{(g)}$ . One notable difference: in the ‘locally Poiseuille’ model used here plug formation does not occur in simulations due to the assumed constant volume flux $Q^{(g)}$ and assumed large viscosity gradient in the model derivation.

How do tube flexibility and slip impact the transport of the film along the tube? The volume flux

(4.3) \begin{align} Q(z,t)&=2\pi \Big [{\textit{Ca}}^{(g)}f_1(R;a)+\frac {1}{16}\big ({\textit{Bo}}\,R^2+R_z+R^2R_{zzz}\big )f_2(R,a)\Big ]\!, \end{align}

may be averaged over the domain

(4.4) \begin{align} \tilde {Q}(t)=\frac {1}{L}\int _0^LQ(z,t)\,{\rm d}z. \end{align}

Figure 15(a) shows $\tilde {Q}(t)$ for six simulations with various combinations of $\varGamma$ , $\psi$ and $\varLambda$ ; figure 15(b) shows the value of $\min _zR(z,t)$ . In each case, much of the film transport comes from the leading-order velocity profile corresponding to perfect core–annular flow with a constant free surface; the value of this leading-order flux can be seen in the initial period of the simulation ( $t\lessapprox 100$ ). Slip plays a large role in promoting film transport as can be expected due to its impact on this leading-order velocity profile.

In all cases, the growth of waves also contributes to the flux, enhancing film transport. While tube flexibility does not impact the leading-order velocity, it does enhance film transport through larger wave growth and higher-amplitude waves. Table 3 shows the percent increase in film transport over the duration of each simulation in figure 15, calculated by averaging $\tilde {Q}$ from $t=500$ to the end of the simulation, dividing this by the mean of $\tilde {Q}$ from $t=0$ to $t=50$ , and subtracting 1.

Table 3. Per cent increase in film flux $\tilde {Q}$ from beginning to end of simulation, calculated by dividing the mean of $\tilde {Q}$ from $t=500$ to $t=623$ by the mean of $\tilde {Q}$ from $t=0$ to $t=50$ in figure 15 and subtracting 1.