A long-wave model for a falling upper convected Maxwell film inside a tube

Abstract

A long-wave asymptotic model is developed for a viscoelastic falling film along the inside of a tube; viscoelasticity is incorporated using an upper convected Maxwell model. The dynamics of the resulting model in the inertialess limit is determined by three parameters: Bond number Bo, Weissenberg number We and a film thickness parameter $a$. The free surface is unstable to long waves due to the Plateau–Rayleigh instability; linear stability analysis of the model equation quantifies the degree to which viscoelasticity increases both the rate and wavenumber of maximum growth of instability. Elasticity also affects the classification of instabilities as absolute or convective, with elasticity promoting absolute instability. Numerical solutions of the nonlinear evolution equation demonstrate that elasticity promotes plug formation by reducing the critical film thickness required for plugs to form. Turning points in travelling wave solution families may be used as a proxy for this critical thickness in the model. By continuation of these turning points, it is demonstrated that in contrast to Newtonian films in the inertialess limit, in which plug formation may be suppressed for a film of any thickness so long as the base flow is strong enough relative to surface tension, elasticity introduces a maximum critical thickness past which plug formation occurs regardless of the base flow strength. Attention is also paid to the trade-off of the competing effects introduced by increasing We (which increases growth rate and promotes plug formation) and increasing Bo (which decreases growth rate and inhibits plug formation) simultaneously.

1. Introduction

The flow of highly viscous films along the inside or outside of tubes arises in numerous applications in engineering and the sciences (Oron, Davis & Bankoff 1997; Grotberg & Jensen 2004; Craster & Matar 2009), and many modelling, experimental and computational studies have been conducted to better understand these flows. In contrast to films falling along an inclined plane, in which surface tension acts to stabilize any free-surface instabilities, these flows are unstable to long-wave disturbances as surface tension plays both a stabilizing and destabilizing role due to the curvature of the free surface (see, e.g. Goren 1962; Yih 1967; Hickox 1971). This Plateau–Rayleigh instability is the same mechanism at work in the breakup of liquid jets; depending on the set-up, other instability mechanisms can be at play as well, e.g. the Kapitza instability (Kapitza 1948) which arises due to inertial effects. As the focus of the current study is highly viscous films, the Plateau–Rayleigh instability is primarily in mind here.

Experiments show that there are multiple dynamical outcomes for falling film flows inside a tube. If the film is thinner than some critical thickness, the instability growth will saturate resulting in a free surface consisting of travelling wave trains. However, if the film is thicker than this critical thickness, wave growth can accelerate with the wave crest approaching the centre of the tube in finite time, resulting in the formation of a liquid plug. This critical thickness depends on the strength of the base flow relative to the effects of surface tension (measured by, e.g. the Bond number Bo). The presence of a base flow, here driven by gravity, provides a nonlinearly stabilizing mechanism to halt wave growth; hence the critical thickness is larger with stronger base flow.

Lubrication theory provides a modelling framework to study these flows; a brief review of the modelling studies most closely related to the current set-up is now given. Hammond (1983) derived a thin-film model for the free-surface evolution of a film coating a tube in the absence of any base flow due to gravity or core flow. Frenkel (1992) derived a thin-film model valid for a falling film along the inside or outside of a tube, and numerically studied the free surface evolution and travelling wave solutions (Kerchman & Frenkel 1994). Kalliadasis & Chang (1994) used self-similar solutions to the model of Frenkel (1992) to find an expression for the blow-up of solitary wave solutions which provides a condition for plug or droplet formation. Craster & Matar (2006), Kliakhandler, Davis & Bankoff (2001), Ji et al. (2019), Camassa & Lee (2006), Camassa, Ogrosky & Olander (2014) and others all developed what will be termed here as ‘long-wave’ models for the same set-up, with the first three focused on flow outside of a tube, and the last two focused on flow inside a tube. The designation ‘long-wave’ is used here to distinguish these models – in which the film is not necessarily assumed to be thin relative to the tube radius – from those in which the thinness of the film is assumed small and directly exploited in the development of the model. A weighted-residual integral boundary-layer (WRIBL) model was developed by Dietze & Ruyer-Quil (2015) for the same set-up, also accounting for the presence of core flow.

The long-wave and WRIBL models have been shown to provide good quantitative agreement with experimental outcomes, with the WRIBL models capable of describing flows with small to moderate Reynolds number particularly well. The linear growth and propagation of free-surface disturbances seen in experiments, as well as the relationship between volume flux and film thickness, has been well captured by both types of models (Smolka, North & Guerra 2008; Camassa et al. 2014). When linear stability analysis is conducted from a spatiotemporal viewpoint, instabilities may be classified as absolute or convective (borrowing from the nomenclature of the study of plasmas). Duprat et al. (2007) used the model of Craster & Matar (2006) to show that this classification quantitatively matched experiments on the exterior of a tube in which instability growth was or was not visible near the inlet, respectively; this was subsequently extended for the case of flows along the interior of a tube by Camassa et al. (2014). Jensen (2000) studied the conditions under which symmetric unduloids might form plugs. Further studies using long-wave and WRIBL models have highlighted the role gravity plays in breaking wave symmetry and modifying the critical thickness required for plugs to form; numerical simulations of these models predict this critical thickness well (Camassa et al. 2014; Dietze, Lavalle & Ruyer-Quil 2020).

This critical thickness is also well captured by families of travelling wave solutions found for long-wave and WRIBL models. Such families, when found for constant values of Bo and varying film thickness, e.g. have a turning point at some thickness, where two branches of travelling wave solution families parametrized by thickness merge; for larger thicknesses, no such solutions are found. This turning-point thickness has been shown to be a good proxy for the critical thickness required for plug formation to occur in both types of models over a variety of parameters (Camassa et al. 2014, 2016; Ding et al. 2019; Dietze et al. 2020; Schwitzerlett, Ogrosky & Topaloglu 2023).

In some applications – e.g. the flow of airway surface liquid lining human airways (Chen et al. 2019) – the film is non-Newtonian. A variety of constitutive equations have been proposed to incorporate the effects of elasticity, with two of the simplest and most common being the Oldroyd-B and upper convected Maxwell (UCM) models. The linear stability of a falling viscoelastic film has been studied by Zhou et al. (2014), who conducted linear stability analysis of a UCM film in the presence of shear and surfactant. They showed that elasticity is destabilizing for a falling film in the absence of shear, but that elasticity can have a stabilizing or destabilizing effect in the presence of shear. Zhou et al. (2016) developed an integral boundary layer model for a viscoelastic film falling down a flexible tube to study the impact of elasticity on the Kapitza and Plateau–Rayleigh instabilities; travelling wave solutions and transient numerical solutions demonstrated the impact elasticity has on wave speed, profile and plug formation. Kang & Chen (1995) developed an asymptotic model for the flow of an Oldroyd-B film down an inclined plane and considered both elastic instability and instability arising from inertia; as with other previous studies, they found that elasticity increases linear instability growth rates, and the weakly nonlinear growth of instabilities was studied. Khayat & Kim (2006) modelled the flow of an Oldroyd-B film flowing over axisymmetric substrates, studying the effects of both inertia and elasticity, and showed the role of viscoelasticity in increasing the accumulation of the film at the inlet. Halpern, Fujioka & Grotberg (2010) developed a lubrication-theory-type model for a viscoelastic film inside a tube in the absence of any base flow, showing that growth rates, plug formation and wall stresses were promoted by elasticity. The impact of elasticity on wall stresses, particularly before and after plug formation or rupture, was further studied numerically by Romano et al. (2021), who used an axisymmetric finite volume method to investigate the wall stresses for Oldroyd-B and FENE-CR films; they showed that elasticity results in a second peak in stress gradients arising after plug formation, with this second peak large enough to cause damage to epithelial cells in human airways. Patne (2021) studied the linear stability of a viscoelastic, shear-thinning film over a flexible wall with surfactant at the free surface in the presence of shear from airflow; liquid elastic modes, excited by the first normal stress difference across the free surface, and solid elastic modes, triggered by the shear, can undergo resonance leading to an amplified growth rate. Patne (2024) studied the impacts of air temperature on a UCM film sheared by turbulent airflow inside a tube, demonstrating the stabilizing effects of warm air on free-surface instabilities. The role of viscoelasticity on mucociliary clearance in lungs has been studied by Choudhury et al. (2023), who used the slip model of Bottier et al. (2017) and Bottier (2017) to model the beating of cilia on clearance of an Oldroyd-B film and also explored the impact of shear-thinning on mucociliary clearance. The impact of plug formation on mucus clearance and flow resistance in the lungs was studied by Fujioka et al. (2016), who developed a semiempirical formula for flow resistance caused by plugs, and Bahrani et al. (2022), who considered elastoviscoplastic plugs. Erken et al. (2023) incorporated elastoviscoplastic properties of mucus for healthy patients and patients with a variety of airway diseases in an airway closure model. The effect of surfactant on a viscoplastic layer was explored in detail by Shemilt et al. (2022, 2023), who showed that both surfactant and yield stress suppress plug formation.

Experimentally, Kim et al. (1986) studied air-driven transport of a viscoelastic film up a tube, showing that viscoelasticity greatly reduced the pressure drop across the tube; viscoelasticity was also found to increase the speed of film transport in vertical tubes with (asymmetric) oscillatory airflow (Kim, Iglesias & Sackner 1987) and with constant upward airflow (Kim et al. 1986). Boulogne, Pauchard & Giorgiutti-Dauphiné (2012) studied flow on a fibre, showing that a shear-thinning film was thinner than its Newtonian counterpart for fixed volume flux, and that growth rate increased. They also point out that the polymers responsible for non-Newtonian rheology also created an effective surface tension lower than the original. Boulogne et al. (2013) discuss how polymers can even suppress the Plateau–Rayleigh instability for viscoelastic films on a fibre. For films inside a tube, experiments by Olander (2020) demonstrated that the addition of a small concentration of polymers may lower the surface tension of a fluid film while producing negligible elastic properties.

In the current study, we develop an asymptotic model for the free-surface evolution of a falling viscoelastic film using a UCM model for the constitutive equations. This model will be used to determine the impact of elasticity on the linear instability of a film, with an analytical expression obtained for the classification of instabilities as absolute or convective. The nonlinear evolution of the free surface will be explored numerically, and travelling wave solution families are explored via numerical continuation of the limit points that serve as proxy for the critical film thickness in model simulations. This continuation succinctly shows the reduction in critical thickness that elasticity induces in the model.

The rest of the paper is organized as follows. The model is developed in § 2, and parameter values corresponding to previous experiments are discussed. Linear stability analysis from both a temporal and spatiotemporal viewpoint is conducted in § 3. Solutions to the model, including those for travelling waves, are found numerically and discussed in § 4. Conclusions are given in § 5.

2. Long-wave model

In this section, a long-wave asymptotic model is derived for a highly viscous viscoelastic falling film coating the interior of a vertical rigid tube.

2.1. Governing equations and UCM model

The governing equations are the incompressible, axisymmetric Navier–Stokes equations in cylindrical coordinates,

(2.1a)$$\begin{gather} \bar{\rho}(\partial_{\bar{t}}\bar{u}+\bar{u}\partial_{\bar{r}}\bar{u}+ \bar{w}\partial_{\bar{z}}\bar{u})={-}\partial_{\bar{r}}\bar{p}+\frac{1}{\bar{r}} \partial_{\bar{r}}(\bar{r}\bar{\sigma}_{{r}{r}})+\partial_{\bar{z}} \bar{\sigma}_{rz}{-\frac{\bar{\sigma}_{{\theta}{\theta}}}{\bar{r}}}, \end{gather}$$

(2.1b)$$\begin{gather}\bar{\rho}(\partial_{\bar{t}}\bar{w}+\bar{u}\partial_{\bar{r}}\bar{w}+\bar{w}\partial_{\bar{z}}\bar{w})={-}\partial_{\bar{z}}\bar{p}+\frac{1}{\bar{r}}\partial_{\bar{r}}(\bar{r}\bar{\sigma}_{rz})+ \partial_{\bar{z}}\bar{\sigma}_{zz}+\bar{\rho} \bar{g}, \end{gather}$$

(2.1c)$$\begin{gather}\frac{1}{\bar{r}}\partial_{\bar{r}}(\bar{r}\bar{u})+\partial_{\bar{z}}\bar{w}=0, \end{gather}$$

where $\bar {\rho }$ is the density of the fluid, $\bar {g}$ is acceleration due to gravity, $\bar {p}$ is pressure, $\bar {t}$ is time, and $\bar {\sigma }$ is the deviatoric stress tensor with components $\bar {\sigma }_{rr}$, $\bar {\sigma }_{rz}$ and $\bar {\sigma }_{zz}$. Overbars denote dimensional quantities and partial derivatives with respect to, e.g. $\bar {r}$, are denoted $\partial _{\bar {r}}$. The coordinates are $(\bar {r},\bar {\theta },\bar {z})$ with associated velocity components $(\bar {u},\bar {v},\bar {w})$; see figure 1. Here, the axial coordinate $\bar {z}$ increases in the downward direction along the tube. The flow is assumed to be axisymmetric, so the $\bar {\theta }$ and $\bar {v}$ components are taken to be zero.

Figure 1. Sketch of the flow configuration and variable definitions.

A model for the deviatoric stress components is needed to close the system. Many models have been proposed and seen widespread use; here, we choose to use the UCM model. Under the UCM model, each of the stress components $\bar {\sigma }_{ij}$ are assumed to obey the following constitutive equations:

(2.2a) \begin{gather} \bar{\sigma}_{rr}+\bar{\lambda}[\partial_{\bar{t}}\bar{\sigma}_{rr}+ \bar{u}\partial_{\bar{r}}\bar{\sigma}_{rr}+\bar{w}\partial_{\bar{z}} \bar{\sigma}_{rr}-2(\partial_{\bar{r}}\bar{u})\bar{\sigma}_{rr}-2(\partial_{\bar{z}}\bar{u}) \bar{\sigma}_{rz}]=2\bar{\mu} \partial_{\bar{r}}\bar{u}, \end{gather}

(2.2b) \begin{gather} \begin{array}{@{}l@{}} \bar{\sigma}_{rz}+\bar{\lambda}[\partial_{\bar{t}}\bar{\sigma}_{rz}+\bar{u}\partial_{\bar{r}} \bar{\sigma}_{rz}+\bar{w}\partial_{\bar{z}}\bar{\sigma}_{rz}-(\partial_{\bar{r}}\bar{w}) \bar{\sigma}_{rr}-(\partial_{\bar{r}}\bar{u}+\partial_{\bar{z}}\bar{w}) \bar{\sigma}_{rz}-(\partial_{\bar{z}}\bar{u})\bar{\sigma}_{zz}] \\ \quad =\bar{\mu}(\partial_{\bar{z}}\bar{u}+\partial_{\bar{r}}\bar{w}),\end{array} \end{gather}

(2.2c) \begin{gather} \bar{\sigma}_{zz}+\bar{\lambda}[\partial_{\bar{t}}\bar{\sigma}_{zz}+\bar{u}\partial_{\bar{r}} \bar{\sigma}_{zz}+\bar{w}\partial_{\bar{z}}\bar{\sigma}_{zz}-2(\partial_{\bar{r}}\bar{w}) \bar{\sigma}_{rz}-2(\partial_{\bar{z}}\bar{w})\bar{\sigma}_{zz}]=2\bar{\mu}\partial_{\bar{z}}\bar{w}, \end{gather}

(2.2d) \begin{gather} \bar{\sigma}_{\theta\theta}+\bar{\lambda}[\partial_{\bar{t}} \bar{\sigma}_{\theta\theta}+\bar{u}\partial_{\bar{r}}\bar{\sigma}_{\theta\theta}+ \bar{w}\partial_{\bar{z}}\bar{\sigma}_{\theta\theta}]=\displaystyle\frac{2\bar{\mu} \bar{u}}{\bar{r}}, \end{gather}

where $\bar {\lambda }$ is the relaxation time, $\bar {\mu }$ is the dynamic viscosity and where the right-hand side of (2.2) are the components of the rate-of-strain tensor $\bar {\tau }$.

Boundary conditions at the tube wall, $\bar {r}=\bar {a}$, are no-slip,

(2.3)\begin{equation} \bar{u}=\bar{w}=0.\end{equation}

At the free surface, $\bar {r}=\bar {R}(\bar {z},\bar {t})$, the boundary conditions are given by: (i) continuity of tangential stress; (ii) jump in normal stress according to the Young–Laplace equation; and (iii) a kinematic boundary condition,

(2.4a)$$\begin{gather} \bar{\sigma}_{rz}+(\partial_{\bar{z}}\bar{R})(\bar{\sigma}_{rr}-\bar{\sigma}_{zz})- (\partial_{\bar{z}}\bar{R})^2\bar{\sigma}_{rz}=0, \end{gather}$$

(2.4b)$$\begin{gather} \begin{aligned}&-\bar{p}[1+(\partial_{\bar{z}}\bar{R})^2]+\bar{\sigma}_{rr}-2(\partial_{\bar{z}}\bar{R}) \bar{\sigma}_{rz}+(\partial_{\bar{z}}\bar{R})^2\bar{\sigma}_{zz}\\& =\bar{\gamma}[1+(\partial_{\bar{z}}\bar{R})^2]\left(\frac{1}{\bar{R} [1+(\partial_{\bar{z}}\bar{R})^2]^{1/2}}-\frac{\partial_{\bar{z}\bar{z}} \bar{R}}{[1+(\partial_{\bar{z}}\bar{R})^2]^{3/2}}\right),\end{aligned} \end{gather}$$

(2.4c)$$\begin{gather}\bar{u}=\partial_{\bar{t}}\bar{R}+\bar{w}\partial_{\bar{z}}\bar{R}, \end{gather}$$

where the background pressure has been set to zero for convenience.

Equations (2.1)–(2.4) may be non-dimensionalized by the scales

(2.5a{–}g) \begin{align} r=\frac{\bar{r}}{\bar{R}_0},\quad z=\frac{\bar{z}}{\bar{\varLambda}},\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{\varLambda}},\quad p=\frac{\bar{p}}{\bar{\rho} \bar{g}\bar{R}_0},\quad \sigma_{ij}=\frac{\bar{\sigma}_{ij}}{\bar{\rho} \bar{g}\bar{R}_0}, \end{align}

where $\bar {R}_0=\bar {a}-\bar {h}_0$ is the distance from the centre of the tube to the mean free surface of the film, $\bar {\varLambda }$ is a typical wavelength of a free-surface disturbance, and where $\bar {U}_0$ and $\bar {W}_0=\bar {\rho } \bar {g}\bar {R}_0^2/\bar {\mu }$ are velocity scales in the radial and axial directions, respectively. The model development that follows will make use of a long-wave assumption that exploits an assumed small ratio of length scales $\epsilon =\bar {R}_0/\bar {\varLambda }\ll 1$. The continuity equation then dictates that $\bar {U}_0=\epsilon \bar {W}_0$. This results in the following dimensionless governing equations:

(2.6a)$$\begin{gather} \epsilon^2 \widetilde{Re}(\partial_tu+u\partial_ru+w\partial_zu)={-}\partial_rp+ \frac{1}{r}\partial_r(r\sigma_{rr})+\epsilon\partial_z \sigma_{rz}{-\frac{\sigma_{\theta\theta}}{r}}, \end{gather}$$

(2.6b)$$\begin{gather}\epsilon \widetilde{Re}(\partial_tw+u\partial_rw+w\partial_zw)={-}\epsilon \partial_zp+\frac{1}{r}\partial_r(r\sigma_{rz})+\epsilon\partial_z\sigma_{zz}+1, \end{gather}$$

(2.6c)$$\begin{gather}\frac{1}{r}\partial_r(ru)+\partial_zw=0, \end{gather}$$

where $\widetilde {Re}=\bar {\rho } \bar {R}_0\bar {W}_0/\bar {\mu }$ is the Reynolds number. In dimensionless form, the constitutive equations are

(2.7a) \begin{gather} \sigma_{rr}+\epsilon \widetilde{We} [\partial_t\sigma_{rr}+u\partial_r\sigma_{rr}+w\partial_z \sigma_{rr}-2(\partial_ru)\sigma_{rr}-2\epsilon (\partial_zu)\sigma_{rz}]=2\epsilon\partial_ru, \end{gather}

(2.7b) \begin{gather} \begin{array}{@{}l@{}} \sigma_{rz}+ \widetilde{We}[\epsilon\partial_t\sigma_{rz}+\epsilon u\partial_r\sigma_{rz}+ \epsilon w\partial_z\sigma_{rz}-(\partial_rw)\sigma_{rr}-\epsilon(\partial_ru+\partial_zw) \sigma_{rz}-\epsilon^2(\partial_zu)\sigma_{zz}] \\ \quad =\epsilon^2\partial_zu+\partial_rw, \end{array} \end{gather}

(2.7c) $$\begin{gather} \sigma_{zz}+ \widetilde{We}[\epsilon\partial_t\sigma_{zz}+\epsilon u\partial_r\sigma_{zz}+ \epsilon w\partial_z\sigma_{zz}-2(\partial_rw)\sigma_{rz}-2\epsilon (\partial_zw)\sigma_{zz}]=2\epsilon\partial_zw, \end{gather}$$

(2.7d) $$\begin{gather}\sigma_{\theta\theta}+\epsilon \widetilde{We}[\partial_t\sigma_{\theta\theta}+ u\partial_r\sigma_{\theta\theta}+ w\partial_z\sigma_{\theta\theta}]=\frac{2\epsilon u}{r}, \end{gather}$$

where $\widetilde {We}=\bar {\lambda }\bar {W}_0/\bar {R}_0$ is the Weissenberg number, the value of which specifies the ratio of elastic forces to viscous forces in the problem. Boundary conditions at the tube wall, $r=a$, where $a=\bar {a}/\bar {R}_0$, are

(2.8)\begin{equation} u=w=0,\end{equation}

and at the free surface, $r=R(z,t)$,

(2.9a)$$\begin{gather} [1-\epsilon^2(\partial_zR)^2]\sigma_{rz}+\epsilon (\partial_zR)(\sigma_{rr}-\sigma_{zz})=0, \end{gather}$$

(2.9b)$$\begin{gather} -p[1+\epsilon^2(\partial_zR)^2]+\sigma_{rr}-2\epsilon (\partial_zR)\sigma_{rz}+ \epsilon^2(\partial_zR)^2\sigma_{zz}\nonumber\\ =\frac{1}{\widetilde{Bo}}[1+\epsilon^2(\partial_zR)^2] \left(\frac{1}{R[1+\epsilon^2(\partial_zR)^2]^{1/2}}- \frac{\epsilon^2\partial_{zz}R}{[1+\epsilon^2(\partial_zR)^2]^{3/2}}\right),\end{gather}$$

(2.9c)$$\begin{gather} u=\partial_tR+w\partial_zR, \end{gather}$$

where $\widetilde {Bo}=\bar {\rho } \bar {g}\bar {R}_0^2/\bar {\gamma }$ is the Bond number, whose value specifies the ratio of gravity forces to capillary forces. We thus have a total of five dimensionless parameters governing the flow: film thickness parameter $a$, Reynolds number $\widetilde {Re}$, Bond number $\widetilde {Bo}$, Weissenberg number $\widetilde {We}$ and aspect ratio $\epsilon$.

2.2. Leading-order equations and solution

Integrating the continuity equation (2.6c) across the fluid layer and applying the boundary conditions (2.8) and (2.9c) results in an evolution equation for the free surface $R(z,t)$,

(2.10)\begin{equation} \partial_tR=\frac{1}{R}\partial_z\int_R^awr\,{\rm d} r. \end{equation}

All that is needed to close the model is an expression for the axial velocity $w$. To find an approximation of $w$, a regular perturbation expansion in $\epsilon$ is used, with each variable expressed as

(2.11a)$$\begin{gather} w=w_0+\epsilon w_1+O(\epsilon^2),\quad p=p_0+\epsilon p_1+O(\epsilon^2), \end{gather}$$

(2.11b)$$\begin{gather}u=u_0+\epsilon u_1+O(\epsilon^2),\quad \sigma_{ij}=\sigma_{ij,0}+\epsilon \sigma_{ij,1}+O(\epsilon^2). \end{gather}$$

Substituting (2.11) into the governing equations (2.6) and truncating at leading order in $\epsilon$ results in

(2.12a{–}c)\begin{equation} \partial_rp_0=\frac{1}{r}\partial_r(r\sigma_{rr,0}){-\frac{\sigma_{\theta\theta,0}}{r}},\quad \frac{1}{r}\partial_r(r\sigma_{rz,0})={-}1,\quad \frac{1}{r}\partial_r(ru_0)+\partial_zw_0=0, \end{equation}

and the constitutive equations of the UCM model are given by

(2.13a)$$\begin{gather} \sigma_{rr,0}=0, \end{gather}$$

(2.13b)$$\begin{gather}\sigma_{rz,0}-\widetilde{We}(\partial_r w_0)\sigma_{rr,0}=\partial_rw_0, \end{gather}$$

(2.13c)$$\begin{gather}\sigma_{zz,0}-2\widetilde{We}(\partial_r w_0)\sigma_{rz,0}=0, \end{gather}$$

(2.13d)$$\begin{gather}\sigma_{\theta\theta,0}=0. \end{gather}$$

The no-slip boundary condition at $r=a$ dictates $u_0=w_0=0,$ and at $r=R(z,t)$, the boundary conditions are

(2.14) \begin{equation} \sigma_{rz}=0,\quad -p_0+\sigma_{rr,0}=\frac{1}{\widetilde{Bo}}\left(\frac{1}{R}-\epsilon^2R_{zz}\right), \end{equation}

where we have retained one term of higher order in $\epsilon$. While a full discussion of the retention of this term is omitted here, this term is commonly included in long-wave asymptotic models and has been shown to be the first-appearing term in the expansion that prevents shock formation; it also provides the correct cut-off wavenumber for linearly unstable modes and produces excellent agreement with experiments for Newtonian films.

The leading-order solution to (2.12a{–}c)–(2.14) is given by

(2.15) \begin{equation} \left.\begin{gathered} u_0=\frac{R(\partial_zR)}{4r}\left(r^2-a^2-2r^2\log\frac{r}{a}\right),\quad \sigma_{rz,0}=\frac{1}{2r}(R^2-r^2), \\ w_0=\frac{1}{4}\left(a^2-r^2+2R^2\log\frac{r}{a}\right),\quad \sigma_{zz,0}=\frac{\widetilde{We}}{2r^2}(R^2-r^2)^2,\\ p_0={-}\frac{1}{\widetilde{Bo}}\left(\frac{1}{R}-\epsilon^2\partial_{zz}R\right),\quad \sigma_{rr,0}=\sigma_{\theta\theta,0}=0. \end{gathered}\right\} \end{equation}

Note that the presence of elasticity does not modify the leading order base flow $(u,w)$, but does modify the $\sigma _{zz,0}$ component of the stress tensor. Equation (2.15) results in a leading-order model

(2.16)\begin{equation} {\partial_tR=\frac{1}{R}\partial_z\int_R^aw_0r\,{\rm d} r=\frac{1}{2} \left(R^2-a^2-2R^2\log\frac{R}{a}\right)\partial_zR.} \end{equation}

2.3. First-order equations and long-wave model

Next, the base flow may be used to find the first-order axial velocity $w_1$. The case of low-Reynolds-number flow, i.e. $\widetilde {Re}=O(\epsilon )$, is considered first. In this case, the inertial terms appear at $O(\epsilon ^2)$ and the governing equations at first-order in $\epsilon$ are given by

(2.17a)$$\begin{gather} \partial_rp_1=\frac{1}{r}\partial_r(r\sigma_{rr,1})+ \partial_z\sigma_{rz,0}{-\frac{\sigma_{\theta\theta,1}}{r}}, \end{gather}$$

(2.17b)$$\begin{gather}{\frac{1}{r}\partial_r(r\sigma_{rz,1})}=\partial_zp_0-\partial_z\sigma_{zz,0}, \end{gather}$$

(2.17c)$$\begin{gather}\frac{1}{r}\partial_r(ru_1)+\partial_zw_1=0. \end{gather}$$

The UCM model constitutive equations are

(2.18a) \begin{gather} \sigma_{rr,1}=2\partial_ru_0,\end{gather}

(2.18b) \begin{gather}\begin{array}{@{}l@{}} \sigma_{rz,1}+\widetilde{We}[\partial_t\sigma_{rz,0}+u_0 \partial_r\sigma_{rz,0}+w_0\partial_z\sigma_{rz,0} \\ \quad -(\partial_rw_0)\sigma_{rr,1}-(\partial_rw_1) \sigma_{rr,0}-(\partial_ru_0+\partial_zw_0)\sigma_{rz,0}]=\partial_rw_1, \end{array} \end{gather}

(2.18c) \begin{gather} \begin{array}{@{}l@{}} \sigma_{zz,1}+\widetilde{We}[\partial_t\sigma_{zz,0}+u_0\partial_r \sigma_{zz,0}+w_0\partial_z\sigma_{zz,0}\\ \quad -2(\partial_rw_0)\sigma_{rz,1}-2(\partial_rw_1) \sigma_{rz,0}-2(\partial_zw_0)\sigma_{zz,0}]=2\partial_zw_0,\end{array} \end{gather}

(2.18d) \begin{gather} \sigma_{\theta\theta,1}=\frac{2u_0}{r}. \end{gather}

Once again, no-slip boundary conditions apply at the wall $r=a$, while at the free surface $r=R(z,t)$, the first-order boundary conditions are

(2.19a)$$\begin{gather} \sigma_{rz,1}+(\partial_zR)(\sigma_{rr,0}-\sigma_{zz,0})=0, \end{gather}$$

(2.19b)$$\begin{gather}-p_1+\sigma_{rr,1}-2(\partial_zR)\sigma_{rz,0}=0. \end{gather}$$

While each equation and condition has been given through first-order, only (2.17b), (2.18a), (2.18b) and (2.19a) are needed to derive the first-order axial velocity $w_1$. In solving for $w_1$, note that (2.18b) contains $\partial _t\sigma _{rz,0}=RR_t/r$; to evaluate this contribution from the leading order stress, we use the leading order model (2.16) to substitute for $\partial _tR$. This approach is identical to that taken by e.g. Camassa et al. (2014) to incorporate inertial effects. Note that this approach may not be applicable to settings without a base flow, as considered by Halpern et al. (2010), for example. This results in the following expression for $w_1$:

(2.20)\begin{align} {w_1}&=\frac{1}{4\widetilde{Bo}}\left(r^2-a^2-2R^2\log\frac{r}{a}\right) \left(\frac{\partial_zR}{R^2}+\partial_{zzz}R\right) \nonumber\\ &\quad +\frac{\widetilde{We}R(\partial_zR)}{8} \left[r^2-a^2-R^2+\frac{a^2R^2}{r^2}+2(a^2-R^2)\log\frac{r}{a} \vphantom{\left(\log\frac{R}{a}+\log\frac{R}{r}\right)}\right. \nonumber\\ &\quad \left. +\;4R^2\log\frac{r}{a}\left(\log\frac{R}{a}+\log\frac{R}{r}\right)\right]. \end{align}

Using (2.15) and substituting $w_0+\epsilon w_1$ into (2.10), the model equation is

(2.21)\begin{equation} \partial_tR=\frac{1}{16R}{\partial_z}\left[f_1(R;a)\left(1-\frac{a^2}{Bo} \left(\frac{\partial_zR}{R^2}+\partial_{zzz}R\right)\right)-\frac{{We}R(\partial_zR)}{2a}f_2(R;a)\right], \end{equation}

where

(2.22a)$$\begin{gather} f_1(R;a)=a^4-4a^2R^2+3R^4-4R^4\log\frac{R}{a}, \end{gather}$$

(2.22b)$$\begin{gather}f_2(R;a)=3a^4-3R^4-4R^4\log\frac{R}{a}+16a^2R^2\log\frac{R}{a}+8R^4\left(\log\frac{R}{a}\right)^2, \end{gather}$$

and where we have dropped the $\epsilon$'s from the final model equation. This corresponds to returning to the original aspect ratio in the dimensionless problem; while the $\epsilon$'s are no longer in the model equation, it must be pointed out that the strict validity of the model derivation still relies on the assumption $\epsilon \ll 1$ is satisfied. Note that (2.21) has been expressed in terms of a rescaled Bond number ${Bo}=a^2\widetilde {Bo}$ and Weissenberg number ${We}=a\widetilde {We}$; these rescaled numbers ${Bo}=\bar {\rho }\bar {g}\bar {a}^2/\bar {\gamma }$ and ${We}=\bar {\lambda }\bar {\rho }\bar {g}\bar {a}/\bar {\mu }$ have the advantage of being independent of film thickness and depend only on parameters density, viscosity, tube radius, surface tension and relaxation time. It is interesting to note that the same model equation may be derived using the Oldroyd-B constitutive equations in place of the UCM model, so that the dynamical outcomes produced by either constitutive model are identical. Note that $f_i(R;a)>0$ (for $i=1,2$) since $0< R< a$.

In the case where ${Re}=O(1)$, the inertial terms appear at first-order; the details of the derivation in the case of a Newtonian film can be found from Camassa et al. (2014) and are omitted here as the resulting terms are identical (up to choice of scales):

(2.23)\begin{align} \partial_tR &= \frac{1}{16R}{\partial_z}\left[f_1(R;a)\left(1-\frac{a^2}{Bo} \left(\frac{\partial_zR}{R^2}+\partial_{zzz}R\right)\right) \right.\nonumber\\ &\quad \left.-\frac{{We}R(\partial_zR)}{2a}f_2(R;a)-\frac{Re}{a^3}f_3(R;a)\partial_zR\right], \end{align}

where

(2.24)\begin{align} f_3(R;a)&=\tfrac{59}{48}R^7-\tfrac{15}{16}a^2R^5-\tfrac{9}{16}a^4R^3+\tfrac{13}{48}a^6R- \tfrac{17}{4}a^2R^5\log (R/a)+\tfrac{7}{4}a^4R^3\log (R/a) \nonumber\\ &\quad -\tfrac{5}{2}R^7(\log(R/a))^2+\tfrac{5}{2}a^2R^5(\log (R/a))^2+2R^7(\log(R/a))^3, \end{align}

and where ${Re}=a^3\widetilde {Re}=\bar {\rho }^2\bar {g}\bar {a}^3/\bar {\mu }^2$ is a rescaled Reynolds number. As the focus here is on highly viscous films, model results will be primarily presented for (2.21), but comments will be given on (2.23) as well. It is again the case that $f_3(R;a)>0$ for $0< R< a$.

In the thin-film limit, the nonlinearities in the model simplify considerably. A film thickness $h$, scaled so that $h=0$ at the wall and $h=1$ at the mean free surface, may be defined by

(2.25)\begin{equation} h=\frac{a-R}{a-1}. \end{equation}

Substituting $R=a-\beta h$, with $\beta =a-1$, into (2.21) and expanding about $\beta =0$ results in

(2.26)\begin{equation} \partial_th+\frac{1}{3}{\partial_z}\left[h^3\left(1+\frac{1}{Bo^*}(\partial_zh+\partial_{zzz}h)\right)+ {We}^* h^4(\partial_zh)\right]=0, \end{equation}

where ${Bo}^*={Bo}/\beta$, ${We}^*=\beta ^2{We}$ and where time has been rescaled by $\beta ^2$. Equation (2.26) is similar to that of Frenkel (1992) but with viscoelasticity included, and is also similar to that of Kang & Chen (1995) for flow down an inclined plane but includes the effects of azimuthal curvature. Note that the viscoelastic term is proportional to $h^4$, while surface tension and advection terms contain $h^3$. In that respect – having a higher-order term that will be shown to drive instability growth – it bears some resemblance to the equation derived by Benney (1966) for film flow down an inclined plane, in which instability growth is provided by an inertial term proportional to $h^6$. For certain parameter values, unbounded solutions to this Benney equation can arise, which has no physical relevance for the falling-film-down-plane problem; this issue can be addressed through use of the WRIBL modelling approach (see, e.g. Ruyer-Quil & Manneville 1998; Dietze & Ruyer-Quil 2015). It would be interesting to conduct a detailed study of (2.26) for similar features, especially since unbounded solutions are not necessarily physically irrelevant given the potential for plug formation in the tube problem considered here, but nonlinear results will be presented here strictly for the model (2.21).

2.4. Parameter values

Before presenting model results, a brief discussion is given concerning the range of parameter values for Bo, ${Re}$, We and $a$ that are in mind here. As the model is derived to consider highly viscous flows, we briefly consider some experimental set-ups for which this model could apply.

In the Newtonian experiments of Camassa et al. (2014), a highly viscous silicone oil is used as the solvent, with density, dynamic viscosity and surface tension given approximately by

(2.27a{–}c)\begin{equation} \bar{\rho}=0.97\,\textrm{g}\,\textrm{cm}^{{-}3},\quad \bar{\mu}=125\,\textrm{P},\quad \bar{\sigma}=21.5\,\textrm{dyn}\,\textrm{cm}^{{-}1}, \end{equation}

respectively. The inner radius $\bar {a}$ of the tube in the experiments ranged from 0.5 cm down to $0.17$ cm; the thickness of the film $\bar {h}_0$ in viscous gravity-driven experiments varied, but most experiments had film thicknesses between 0.1 and 0.2 cm. Using these values results in the following estimated ranges for the dimensionless parameters in the model for a Newtonian fluid,

(2.28a{–}c) \begin{equation} {Bo}=1.3-11.1,\quad {Re}=3\times10^{{-}4}-7\times10^{{-}3},\quad a=1.25-2.5. \end{equation}

As ${Re}<0.01$ for all these cases, the focus will primarily be on the inertialess model (2.21), with values of Bond number primarily lying in the range $0.5<{Bo}<15$ and tube radius parameter $1< a<3$, occasionally showing results for other values of Bo or $a$ to illustrate a point about the model.

For viscoelastic films, Shaqfeh, Larson & Frederickson (1989) (see also Kang & Chen 1995) found that the Oldroyd-B model is capable of capturing the behaviour of polymeric thin films so long as ${We}<2.5$. Based on this criteria, and given that the model equation here is identical to that obtained using the Oldroyd-B equations, we will largely restrict results to ${We}<2.5$.

The addition of polymers to a Newtonian fluid may also modify the surface tension; e.g. an increase in ${We}$ may in some cases be accompanied by a corresponding increase in Bo. This effect was seen, e.g. in experiments conducted by Olander (2020), in which the falling film consisted of polybutene (PB) and various small concentrations of polyisobutylene (PIB); the PIB concentrations were small enough to result in negligible values of ${We}$ but still lower the value of $Bo$ through reduction of the surface tension. Some attention will thus be given to the increase in Bo required to offset a given increase in ${We}$ to produce an equivalent dynamical outcome.

One of the motivations for the experimental values chosen by Camassa et al. (2014) was the study of mucus lining human airways; the viscosities, surface tension and tube radii (and hence ${Bo}$) were chosen to be representative of upper airways. The rheological properties of human mucus can vary widely with patient, disease and airway generation, with elasticity decreasing with airway generation and increasing in the presence of airway diseases like cystic fibrosis (Lai et al. 2009).

3. Linear stability analysis

It is well known that the free surface of a film coating the interior or exterior of a tube is unstable to long-wave disturbances due to the Plateau–Rayleigh instability. The impact of viscoelasticity on this instability is explored next from both a temporal and spatiotemporal viewpoint, with the latter allowing classification of instabilities as convective or absolute.

3.1. Temporal linear stability analysis

To study the linear stability of a flat-interface solution $R(z,t)=1$, let the free surface be perturbed by a superposition of small-amplitude Fourier modes,

(3.1)\begin{equation} R=1+\hat{R}\exp({{\rm i}(kz-\omega t)}),\end{equation}

where $k$ is a spatial wavenumber and where it is assumed that $\hat {R}\ll 1$; while exploring the temporal stability of the free surface, the wavenumber $k$ will be taken to be real-valued. Substituting (3.1) into (2.21) results in the dispersion relation

(3.2)\begin{equation} \omega=\frac{k}{2}(a^2-1-2\log a)+\frac{{\rm i}a^2}{16{Bo}}({k^2}-k^4)f_1(1;a)+ \frac{{\rm i}{We}\,k^2}{32a}f_2(1;a). \end{equation}

The real part of (3.2), $\textrm {Re}[\omega ]/k$, determines the phase speed of disturbances, while the imaginary part, $\textrm {Im}[\omega ]$, dictates the growth of their amplitude. As with other thin-film models, the dispersion relation is similar to that of the Kuramoto–Sivashinsky (KS) equation in the dependence on wavenumber $k$. For reference, the dispersion relation for the model with inertial terms (2.23) and for the inertialess thin-film model (2.26) are also included:

(3.3) \begin{align} \omega=\frac{k}{2}(a^2-1-2\log a)+\frac{{\rm i}a^2}{16{Bo}}({k^2}-k^4)f_1(1;a)+ \frac{{\rm i}{We}\,k^2}{32a}f_2(1;a)+\frac{{\rm i}{Re} k^2}{16a^3}f_3(1;a) \end{align}

and

(3.4)\begin{equation} \omega=\frac{1}{3}\left(k+\frac{{\rm i}\beta }{Bo}(k^2-k^4)+{\rm i}{We}\beta^2k^2\right), \end{equation}

respectively.

Equation (3.2) has been well studied for Newtonian films (${We}=0$); similar to other settings where the Plateau–Rayleigh instability is the dominant one, the film is unstable to long-wave disturbances with cut-off wavenumber $k_c=1$ and the wavenumber of maximum growth rate, $k_{max}=1/\sqrt {2}$, set by a balance of destabilization of the free surface due to its azimuthal curvature and stabilization due to its axial curvature, both arising due to surface tension. The phase speed of disturbances is set by the balance of gravity and viscous forces. In the linear regime, the dimensionless growth rates of this model are determined by the $k^2$ and $k^4$ terms arising due to the Plateau–Rayleigh instability, with the growth rates increasing monotonically with $a$ and with $1/{Bo}$. The inclusion of base flow due to gravity in this first-order model does not impact these dimensionless growth rates as the hyperbolic terms in (2.21) associated with the effects of gravity arise solely in the real part of the dispersion relation, similar to many other first-order film flow models with base flow. This model dispersion relation has been shown to agree well with that of the Stokes equations (Camassa et al. 2014).

In the case where ${We}>0$, the presence of elasticity introduces a second mechanism of free-surface destabilization. Figure 2(a) shows the growth rates $\textrm {Im}[\omega (k)]$ (3.2) for $a=1.4$, ${Bo}=1$ and a variety of ${We}$ values; the growth rate increases with increasing ${We}$ for all wavenumbers $k$. The $k^2$ terms of (3.2) are identical to those of Zhou et al. (2014) in their long-wave analysis without surfactant. Both the wavenumber $k_{max}$ of maximum growth rate and the cutoff wavenumber $k_c$ increase as well, according to

(3.5a,b) \begin{equation} k_c=\sqrt{1+\frac{{We}\,{Bo}\, f_2(1;a)}{2a^3f_1(1;a)}},\quad k_{max}=k_c/\sqrt{2}, \end{equation}

so that the product ${Bo}\,{We}=\bar {\lambda }\bar {\rho }^2\bar {g}^2\bar {a}^3/\bar {\mu }\bar {\gamma }$ determines the change in $k_{max}$ and $k_c$ for fixed $a$. Figure 2(b) shows this dependence of $k_c$ and $k_{max}$ on ${We}$ for $a=1.4$. As the value of Bo increases, changing ${We}$ results in a larger change in $k_{max}$.

Figure 2. (a) Growth rates of (2.21) for $a=1.4$, ${Bo}=2$ and various values of We. Here, $\times$ symbols denote the maximum growth rate $GR_{max}$ and ${\cdot }$ symbols denote cutoff wavenumber $k_c$. (b) $k_{max}$ and $k_c$ for $a=1.4$ and various values of Bo and We. (c) Same as panel (a) but with $a=1.05$. (d) $f_2(1;a)/[a^3f_1(1;a)]$; maximum value obtained at $a\approx 1.375$. Open circles denote $a$ values used in panels (a,c).

As with the Newtonian film, the growth rates increase with $a$ and $1/{Bo}$, but the relative impact of increasing ${We}$ depends non-monotonically on film thickness parameter $a$. Figure 2(c) shows the growth rates for a much thinner film with $a=1.05$. As may be expected, the growth rates are smaller than for $a=1.4$, but it is also true that the growth rate curves are closer together than in figure 2(a), indicating that increasing ${We}$ has a smaller relative impact on the growth rates and on $k_{max}$ when $a=1.05$ than when $a=1.4$. This may be anticipated from (3.5a,b) or from the thin-film dispersion relation (3.4) in which the term proportional to ${We}$ occurs at higher order in $\beta$ than the other terms.

For moderate values of $a$, the impact of $a$ on $k_{max}$ must be assessed by examining $f_2(1;a)/[a^3f(1;a)]$ in (3.5a,b); this quantity is plotted in figure 2(d). This value reaches a maximum at $a\approx 1.375$; i.e. this value of $a$ produces the largest impact on $k_{max}$ and largest relative impact on the maximum growth rate $GR_{max}$ for small ${We}$ as determined by $[(\partial GR_{max}/\partial We)/GR_{max}]|_{We=0}=f_2(1;a)Bo/[a^3f_1(1;a)]$. We also note that for thick films (large $a$), the relative impact of viscoelasticity on $k_{max}$ and growth rates is again small, as $f_2(1;a)/[a^3f_1(1;a)]$ behaves like $1/a^3$ for large $a$.

In summary, linear instability is driven by both surface tension and viscoelasticity. While growth rates increase monotonically with $a$ and $1/{Bo}$, the relative impact of viscoelasticity on growth rates is highest for moderately thick films.

Note that while the base flow $w_0$ is unchanged by the presence of viscoelasticity, the normal stress is positive, $\sigma _{zz,0}>0$. This stress term plays a crucial role in determining the sign of the first-order velocity $w_1$ through (2.17b); $w_1$ may be expressed as the sum of three parts:

(3.6)\begin{equation} {w_1=w_{ST}+w_{NS}+w_{UCM},} \end{equation}

where

(3.7a) $$\begin{gather} {w_{1,ST}}=\frac{1}{4\widetilde{Bo}}\left(r^2-a^2-2R^2\log\frac{r}{a}\right) \left(\frac{\partial_zR}{R^2}+\partial_{zzz}R\right), \end{gather}$$

(3.7b) $$\begin{gather}{w_{1,NS}}=\frac{\widetilde{We}RR_z}{2}\left(r^2-a^2+2R^2\log\frac{a}{r} \left(1+\log\frac{a}{R}+\log\frac{r}{R}\right)\right), \end{gather}$$

(3.7c) \begin{gather}\begin{array}[b]{@{}l@{}} \displaystyle{w_{1,UCM}}=\frac{\widetilde{We}RR_z}{8}\left(3a^2-3r^2-R^2+\frac{a^2R^2}{r^2}+6R^2 \log\frac{r}{a} \right.\\ \displaystyle\phantom{{w_{1,UCM}}=}\left.+\;2a^2\log\frac{r}{a}-4R^2\log\frac{r}{a}\left(\log\frac{R}{a}+\log\frac{R}{r}\right)\right), \end{array} \end{gather}

where ST denotes surface tension, NS denotes normal stress corresponding to $\sigma _{zz,0}$ and UCM refers to the contribution of the terms in brackets in (2.18b). For $R< r< a$, note that if $R_z>0$, then $w_{1,NS}<0$, $w_{1,UCM}>0$ and $w_{1,NS}+w_{1,UCM}<0$. As a result, the normal stresses (specifically the gradient $\partial _z\sigma _{zz,0}$) generated by the base flow are responsible for the instability growth seen in the current model.

3.2. Absolute and convective instability

Linear stability analysis may also be conducted from a spatio-temporal viewpoint, which allows classification of instabilities as absolute or convective. In experiments, convectively unstable films are those in which visible instability growth is only evident downstream, far from the inlet, while absolutely unstable films exhibit visible instability growth right at or near the inlet.

Deissler, Oron & Lee (1991) first classified instabilities for a falling Newtonian film on the exterior of a tube by applying the methods developed and explored by Bers (1983), Briggs (1964) and Huerre & Monkewitz (1990) among others. They identified a critical velocity separating absolute instability from convective instability for a two-dimensional KS equation, with the critical velocity dependent on the azimuthal mode number. Duprat et al. (2007) applied these techniques to a strongly nonlinear model similar to (2.21) with ${We}=0$ for the exterior case. First finding a critical value for the thin-film parameter ${Bo}^*$ and then transforming the thin-film dispersion relation to the long-wave dispersion relation through a substitution, Duprat et al. (2007) identified a condition on Bond number ${Bo}^*$ and film thickness parameter $a$ that must be met for films to be absolutely unstable. Camassa et al. (2014) later applied this approach to identify a threshold thickness for Newtonian films inside a tube; films thicker than this threshold were shown to be absolutely unstable, while thinner films were convectively unstable.

We briefly state this threshold condition for the Newtonian case. The dispersion relation of the Newtonian case (denoted by subscripts of $N$),

(3.8)\begin{equation} \omega_N=\frac{k_N}{2}(a^2-1-2\log a)+\frac{{\rm i}a^2}{16{{Bo}_N}}({k_N^2}-k_N^4)f_1(1;a), \end{equation}

exhibits absolute instability if the criteria

(3.9)\begin{equation} {\frac{2a^2f_1(1;a)}{{Bo}_N(a^2-1-2\log a)}>8({-}17+7\sqrt{7})^{1/2}} \end{equation}

is met, where ${Bo}_N$ denotes the Bond number in the Newtonian case. This condition matches that found by Camassa et al. (2014) (adjusted for choice of scalings), and arises from considering branches of solutions to (3.8) found for complex wavenumber $k$ and identifying parameter values at which distinct branches coalesce (see, e.g. Huerre & Monkewitz (1990), Duprat et al. (2007) and Camassa et al. (2014) for further discussion and a thorough derivation of this criteria).

How does the presence of elasticity affect this condition? The dispersion relation of the Newtonian case (3.8) may be transformed into (3.1) through the substitution

(3.10a,b) \begin{align} {(k_N,\omega_N)=\left(1+\frac{{We}\, {Bo} f_2(1;a)}{2a^3f_1(1;a)}\right)^{{-}1/2}(k,\omega),\quad {Bo}_N=\left(1+\frac{{We}\, {Bo} f_2(1;a)}{2a^3f_1(1;a)}\right)^{{-}3/2}{Bo}.} \end{align}

Substituting (3.10a,b) into (3.9) results in the following condition for absolute instability:

(3.11)\begin{equation} \frac{2{Bo}^{1/2}f_1(1;a)}{a(a^2-1-2\log a)}\left(\frac{a^2}{Bo}+ \frac{{We} f_2(1;a)}{2af_1(1;a)}\right)^{3/2}>8({-}17+7\sqrt{7})^{1/2}. \end{equation}

Figure 3 shows the $a$-Bo boundary between absolute and convective instability for a variety of ${We}$ values. For ${We}=0$, the boundary between absolute and convective instability monotonically increases with increasing Bo. For large values of Bo (or $a$), the boundary is governed by ${Bo}\approx a^4/[4(-17+7\sqrt {7})^{1/2}]$, which may be found by using $f_1(1;a)=a^4+O(a^2)$ and $f_2(1;a)=3a^4+O(a^2)$ in (3.11).

Figure 3. (a) Boundary in Bo-$a$ space between absolute and convective instability for a variety of ${We}$ values. (b) Same as panel (a) but for a larger range of Bo values. Blue dots denote maximum values of $a$ obtained along each curve.

Increasing ${We}$ results in a greater region of absolute instability, with the effect most pronounced for larger Bo values. In contrast to the ${We}=0$ case, the boundary for ${We}>0$ exhibits non-monotonic behaviour, reaching some maximum attainable value of $a_c$ for a given fixed positive ${We}$. This maximum value, denoted by blue dots in figure 3(b), occurs when ${Bo}=4a^3f_1(1;a)/[{We} f_2(1;a)]$. Past this value of Bo, the boundary decreases monotonically as ${Bo}\rightarrow \infty$, the scaling of which may be found in the thin-film limit by substituting $a=1+\beta$ and $Bo=b\beta ^{\alpha }$ into (3.11) with $\beta \ll 1$. This results in

(3.12)\begin{equation} \frac{16\beta^{1-\alpha}}{3b}(1+{We} b \beta^{\alpha+1})^{3/2}>8({-}17+7\sqrt{7})^{1/2}. \end{equation}

There are two consistent scalings possible, corresponding to one of the two terms inside parentheses providing the dominant balance with the right-hand side of (3.12). If the second term dominates, then $\alpha =-5$ and ${Bo}={9(-17+7\sqrt {7})}/{4{We}^3\beta ^5}$, corresponding to the border between absolute and convective instability for large Bo. The other scaling, with $\alpha =1$ and ${Bo}=2\beta /[3(-17+7\sqrt {7})^{1/2}]$, corresponds to the border for small Bo and is independent of ${We}$.

4. Nonlinear results

In this section, the nonlinear equation (2.21) is solved numerically. Travelling wave solutions are found in § 4.2.

4.1. Evolution equation

A finite difference method was used to solve (2.21) numerically. Spatial derivatives were calculated using fourth-order approximations; time-stepping was done with an explicit second-order predictor-corrector method. An initial condition was prescribed by perturbing a constant free surface with a superposition of small-amplitude Fourier modes. In a typical simulation, 40 modes were used to perturb the free surface, though a variety of values were tried. The number of initial modes played a role in the transient dynamics as instabilities initially grew, but once the evolution had settled into a quasi-steady state, there was – in almost all cases – no clear signature of the initial conditions. The one exception to this is for films with a mean thickness very near the critical thickness – in which plug formation would only be triggered by wave mergers and not linear wave growth – where a change in initial conditions could result in a change of the final state (plug or no plug) in the finite-time simulations conducted. This was also noted in a study of Newtonian films with wall slip; see Schwitzerlett et al. (2023) for numerical tests and further discussion of this point.

At $z=0$, inlet boundary conditions were prescribed by populating several ghost nodes upstream of the inlet using a periodic extension of the initial condition moving at the linear phase speed of the free-surface disturbances. At $z=L$, outlet boundary conditions were prescribed by populating several ghost nodes downstream of the inlet using a reflection of the final nodes within the domain. Other outlet boundary conditions were also tried; the exact prescription of these ghost nodes downstream only had a significant impact on the final few nodes near $z=L$, which will be largely ignored in the results. Simulations were also conducted using periodic boundary conditions. Since the dynamical outcomes of the simulations did not change when the boundary conditions were varied, these results are omitted, and results using inlet/outlet boundary conditions will be presented here.

Two dynamical outcomes are evident in simulation snapshots shown in figure 4(a–c). The baseline case in figure 4(a) shows a snapshot of the film in the tube using ${Bo}=10$, ${We}=0$ and $a=1.8$. This snapshot is taken near $t=500$, long after the dynamics have settled into a quasi-steady state. The gently perturbed film entering the domain at the inlet (left) is clearly convectively unstable, with disturbances growing as the film falls down the tube. Between $z=100$ and $z=150$, this growth saturates, with the remainder of the tube coated with a wavy film consisting of small-amplitude travelling-wave-like ripples.

Figure 4. (a–c) Snapshots of the film evolution for $a=1.8$ and for various Bo and We values. In panels (b,c), the snapshot is taken just prior to the formation of the first plug. (d) Maximum and minimum values of $R$ over the domain in each simulation as a function of $t$.

A second outcome is shown in figure 4(b) with identical parameter values to figure 4)(a) except for ${Bo}=5$. The film is still convectively unstable, but visible growth is seen closer to the inlet, e.g. near $z\approx 50$; the free-surface waves also have larger amplitude than with ${Bo}=10$. However, the simulation is halted near $t=150$ as one of the waves – near $z=220$ – begins to undergo accelerated growth, with the free surface at the wave crest approaching $R=0$ in finite time. Figure 4(d) shows the evolution of $\max R$ and $\min R$ over the domain as a function of time. As with other studies of asymptotic models, $R$ approaching zero may be taken as the model's prediction that plug formation is about to occur. No attempt is made here to continue running the model forward past this time, as the nonlinearities in (2.21) include logarithms and inverse powers of $R$, preventing the model from being run forward as it is, though there has been interesting recent work in this direction. Most notably, Dietze et al. (2020) modified a WRIBL model for film flow inside a tube to prevent $R$ from reaching a value of zero during plug formation; this was achieved through inclusion of a novel pressure term that allowed for arbitrarily narrow filaments of core fluid (air) to form where the model exhibited plug formation. Here, the focus is on whether plugs may be expected to form or not, and their propagation after formation is left for future work.

Figure 4(a,b) indicates that for $a=1.8$, there is some critical Bond number $5< {Bo}_c<10$ such that for ${Bo}>{Bo}_c$, plugs do not form; while for ${Bo}<{Bo}_c$, plugs do form. This is consistent with previous work (Camassa et al. 2014; Dietze et al. 2020; Ogrosky 2021) showing that stronger base flow inhibits plug formation in Newtonian films. The impact of viscoelasticity is examined next.

Figure 4(c) shows a snapshot of the film with ${Bo}=10$, ${We}=0.9$ and $a=1.8$; note that all three panels of figure 4 have the same value of film thickness parameter $a$. Once again, the film is convectively unstable, though instability growth is visible even closer to the inlet, consistent with the higher linear growth rates found in § 3. The free-surface waves with $We>0$ appear to have shorter typical wavelength than those with $We=0$; e.g. in the simulation snapshots shown in figure 4, the average wavelength from $z=100$ to $z=300$ was approximately 9.5 in figure 4(b) and 7.7 in figure 4(c). This shorter average wavelength may be expected given the impact of viscoelasticity on wavelength shown in (3.5a,b). Similar to figure 4(b), the simulation was halted (near $t=95$) as one wave – near $z=150$ – underwent accelerated amplitude growth, resulting in $R\rightarrow 0$ indicating plug formation.

The results above indicate that either increasing ${We}$ or decreasing Bo past some critical value while holding the film thickness constant results in plug formation. Said another way, for fixed Bo and ${We}$, there is a critical value of $a$ such that if $a$ exceeds this value, plug formation may be expected to occur; increasing ${We}$ or decreasing Bo leads to a reduction in this critical value of $a$.

4.2. Travelling wave solutions

Since the free surface in figure 4 evolves into a travelling wave train, travelling wave solutions are sought next. Substituting $Q(Z)=Q(z-ct)=R(z,t)$ into (2.21) and integrating once results in a third-order ordinary differential equation,

(4.1)\begin{equation} 8cQ^2+K+\left[1-\frac{a^2}{Bo}\left(\frac{Q'}{Q^2}+Q'''\right)\right]f_1(Q;a)-\frac{{We}\,QQ'}{2a}f_2(Q;a)=0, \end{equation}

with speed $c$ and constant of integration $K$. The approach used to find solutions to this ODE follows the approach outlined by Camassa et al. (2016) and is briefly outlined here. First, fixed point solutions are found for the three-dimensional dynamical system associated with (4.1); these correspond to constant solutions $Q=Q_0$. If the wave speed $c$ is varied, a Hopf bifurcation in the solution family is identified using XPP/AUTO; XPP is a numerical tool for simulating and analysing dynamical systems (Ermentrout 2002) and also provides a front end to the numerical continuation and bifurcation package AUTO (Doedel et al. 2008). This bifurcation is a zero-Hopf bifurcation, which can be somewhat difficult to find numerically; to avoid this issue, a small viscosity-like term $\eta Q''$ was added to (4.1) which makes continuation onto a family of periodic solutions easier. Once periodic solutions have been found, the viscosity parameter $\eta$ may be taken to zero. The mean value of $Q^2$ is ensured to be equal to one by imposing an integral condition. As the focus here is on waves with profiles similar to those seen in transient solutions like figure 4, only solution families with such profiles are presented here; other families, e.g. were considered by Zhou et al. (2016) for an integral boundary layer model.

Figure 5(a) shows several families of travelling wave solutions with period $4{\rm \pi}$ for a variety of values of Bo and ${We}$; the period $4{\rm \pi}$ was chosen as it is approximately the period of waves for ${We}=0$ seen in figure 4(a) especially near the bottom of the tube, e.g. from $z=150$ to $z=250$. Each curve shows the wave amplitude of solutions within a family as a function of film thickness parameter. Before discussing the impact of viscoelasticity on these curves, we first note that each family has a turning point (denoted by a red dot) at some critical value $a_c$ of the thickness parameter $a$. For ${We}=0$, this critical thickness $a_c$ has previously been shown to be a reliable proxy for the critical thickness required for plugs to form (Camassa et al. 2014, 2016; Ding et al. 2019; Dietze et al. 2020; Camassa et al. 2021). The critical thickness associated with the turning point solution depends on the value of Bo, with higher values of Bo corresponding to higher critical thickness, consistent with the results of § 4.1 and with the results of Dietze et al. (2020) showing that while the base flow does not contribute to the linear growth/decay of small disturbances, it plays a nonlinearly stabilizing role. Note that the horizontal axis shows $(a-1)/a$, with values close to 0 representing thin films, and values close to 1 representing thick films (since in the limit $a\rightarrow 1$ (thin films), $(a-1)/a\rightarrow 0$, while in the limit $a\rightarrow \infty$ (thick film), $(a-1)/a\rightarrow 1$). For each thickness $(a-1)/a$ less than the turning-point thickness, there are two travelling wave solutions found, one lying on a ‘lower’ branch with relatively small amplitude and one lying on an ‘upper’ branch with relatively large amplitude. Lower branch solutions are shown for several values of ${We}$ in figure 5(c); these lower branch solutions correspond well to the waves seen in numerical simulations like those of figure 4(a–c). Upper branch solutions have a qualitatively similar shape, but are omitted here as waves with such large amplitude are not seen in numerical simulations (but see examples of these waves in, e.g. Camassa et al. 2014 for a Newtonian film).

Figure 5. (a) Families of travelling wave solutions with ${Bo}=1$ and ${Bo}=10$ for four values of We. Turning points indicated by red dots. (b) Travelling wave solutions with ${Bo}=10$ corresponding to red dots in panel (a). (c) Travelling wave solutions with ${Bo}=10$ and $a=1.5$ corresponding to blue dots in panel (a).

How does the presence of viscoelasticity impact the critical thickness associated with the turning point of solution families? Figure 5(a) shows that as ${We}$ increases, this turning-point thickness – which serves as a proxy for the critical thickness required for plugs to form – decreases, again consistent with the results of § 4.1 showing that viscoelasticity promotes plug formation. The impact is most pronounced for larger values of Bo. Figure 5(b) shows the solution corresponding to the turning point denoted in red in panel  (a) for ${Bo}=10$; as ${We}$ increases, the support of these waves decreases significantly. While some of this change in support could be attributed to the different film thickness of each wave, viscoelasticity also plays a role; to demonstrate this, figure 5(c) shows travelling wave solutions (corresponding to blue dots in figure 5c) for fixed $a$ and ${Bo}$ and various ${We}$. Despite having uniform thickness, the change in support persists, most noticeable when one observes the markedly different distance from wave crest to the trough (capillary ripple) that precedes the crest.

How does the turning point in travelling wave solution families depend on the chosen period size? Figure 6(a) shows the dependence of $a_c$ on period for ${We}=0$ and ${We}=1.5$. For both curves, there appears to be a value of $a_c$ which the curves approach asymptotically in the limit of large period, while for small period, the value of $a_c$ increases rather sharply. The period $4{\rm \pi}$ is denoted along both curves by a (blue or red) dot.

Figure 6. (a) Dependence of $a_c$ on period size for $We=0$ and $We=1.5$. (b) Families of travelling wave solutions with ${Bo}=10$ for $We=0$ and $We=1.5$; dashed grey and black lines show period $4{\rm \pi}$ solutions (corresponding to the dots in panel a), solid grey line shows period ${\approx }8.5$ solutions (corresponding to red $\times$ in panel a). Thick dots denote the maximum value of $h_{max}$ during numerical simulations of (2.21); see text for details. Vertical (red and blue) dashed lines indicate thickness at which plugs formed in nonlinear simulations.

Given that figure 4(a–c) indicates that elasticity produces waves with shorter wavelength, it may be that a more accurate estimate of the critical thickness would be provided by a turning point found for shorter-wavelength travelling wave solutions for ${We}>0$. Figure 6(b) shows the ${Bo}=10$ travelling wave solution families for ${We}=0$ and ${We}=1.5$ for period $4{\rm \pi}$ (dashed lines), as well as the solution family for ${We}=1.5$ and period ${\approx }8.5$ (solid line). This shorter wavelength produces a turning point at a larger value of $a_c$ than the $4{\rm \pi}$ domain waves, as expected.

To assess which period length better predicts the critical thickness in model solutions, an additional test was conducted. Model solutions were found using an initial condition with a thickness well below the critical thickness. Once the free surface had settled into a quasi-steady state, the film thickness was increased by a small amount. This process was repeated until the simulation indicated plug formation through $R\rightarrow 0$ in finite time. The average value of $\min _z R$ was calculated over a period once a quasi-steady state had been achieved for each thickness; the resulting data points are plotted in figure 6(b) for ${We}=0$ and ${We}=1.5$, with the thickness resulting in plug formation indicated by a vertical dashed line. It appears that for both the Newtonian and non-Newtonian films, the $4{\rm \pi}$ solutions provide a reasonable estimate of the critical thickness.

How does $a_c$ depend on Bo for various values of ${We}$? Figure 7 shows the dependence of the critical thickness parameter $a_c$ on Bo and ${We}$; panel  (a) shows this dependence with various fixed Bo while in panel  (b), ${We}$ is held fixed. For ${We}=0$, the critical thickness increases arbitrarily with increasing Bo. Ogrosky (2021) noted that for Newtonian flow (${We}=0$), there appeared to be a scaling law between Bo and $a_c$ in the limit of strong base flow (weak surface tension), i.e. $Bo\rightarrow \infty$; an empirical fit was found with a 1/5 scaling. Here, we note that the scaling appears to lie somewhere between 1/5 and 1/6 in figure 7(b). A brief justification (relying partially on numerical observations) for the existence of a scaling law is given in Appendix A, where the relationship

(4.2)\begin{equation} {Bo}\propto \frac{a_c^6}{\log a_c} \end{equation}

is identified for $a_c\rightarrow \infty$. Thus, in the absence of elasticity, a film may be arbitrarily thick and still not form plugs, provided the base flow is strong enough relative to the effects of surface tension.

Figure 7. (a) Critical thickness $a_c$ as a function of We for various fixed values of Bo. (b) $a_c$ versus ${Bo}$ for various fixed values of ${We}$. Grey dashed lines show slopes $m=1/5$ and $m=1/6$; black dashed line shows the approximate critical thickness for vanishing Bo.

For small values of ${Bo}$, the presence of elasticity only slightly modifies the value of $a_c$, consistent with the results of figure 5(a). The curves for various fixed ${We}$ nearly overlap for small ${Bo}$ and approach a value corresponding to a film thickness of 12 % of the tube radius (denoted with a grey horizontal dashed line). This percentage has been identified in previous studies with thin-film models as an approximate estimate of the thickness required for plugs to form (Gauglitz & Radke 1988) in a surface-tension dominated limit.

In contrast to the Newtonian case, however, there is now a maximum value of $a_c$ that is attainable for some ${Bo}$ at a given fixed ${We}>0$; past this ${Bo}$ value, the critical thickness $a_c$ decreases monotonically. Some intuition for this result may be found by revisiting the temporal linear stability analysis. Recall that the maximum growth rate is given by

(4.3)\begin{equation} {GR_{max}=\frac{a^2f_1(1;a)}{16{Bo}}(k_{max}^2-k_{max}^4)+\frac{{We} f_2(1;a)}{32a}k_{max}^2}, \end{equation}

where $k_{max}$ is given in (3.5a,b). Substituting (3.5a,b) into (4.3), taking a derivative with respect to ${Bo}$ and setting equal to zero, a local minimum in $GR_{max}$ as a function of ${Bo}$ may be found at ${Bo}=2a^3f_1(1;a)/{We} f_2(1;a)$ (or alternately, when $k_{max}=1$). This is depicted in figure 8, in which growth rate curves for ${We}=0.8$ and $a=1.4$ are shown for various ${Bo}$ values. The dependence of $GR_{max}$ on ${Bo}$ (through $k_{max}$ in (3.5a,b)) is shown by the dashed line. For these fixed values of $a$ and ${We}$, there is some value of ${Bo}$ for which the growth rate is minimized over all ${Bo}$. While a minimum in the (linear) growth rate cannot be expected to exactly correspond to a maximum in critical thickness, it seems reasonable to conjecture that the existence of such a minimum may indicate a value of ${Bo}$ for which linear instabilities may more readily saturate in the nonlinear regime, suppressing plug formation and allowing for a relatively large value of $a_c$.

Figure 8. Growth rates of (2.21) for $a=1.4$, ${We}=0.8$ and various values of Bo. $\times$ symbols denote the maximum growth rate $GR_{max}$, and ${\cdot }$ symbols denote cutoff wavenumber $k_c$. Black dashed line indicates dependence of $GR_{max}$ on ${Bo}$ through $k_{max}$ defined in (3.5a,b).

It is important to note that the addition of polymers to a previously Newtonian film not only introduces elasticity (modifying ${We}$) but can also change the surface tension of the film (modifying Bo). Thus, if polymers reduce the film's surface tension, it is possible that their presence could actually inhibit plug formation, despite the plug-promoting role of elasticity. It may be important, therefore, to quantify the dependence of $a_c$ on both Bo and ${We}$, particularly for ${We}\ll 1$. Similar to figure 7, figure 9(a) shows curves containing turning point solutions but now for various fixed values of $a_c$. Each of these curves reaches a maximum value of ${We}$, consistent with the discussion regarding figure 7. The blue dashed lines are tangent lines to each curve at ${We}=0$, and indicate the slope or $(\partial Bo/\partial {We})|_{{{We}=0}}$ for isolines of $a_c$. Those slopes are also shown in figure 9(b) by the three blue dots, with the black curve showing the dependence of this slope on $Bo$. This dependence indicates that the slope $(\partial Bo/\partial {We})|_{{We}=0}$ scales with $Bo^2$.

Figure 9. (a) Critical thickness $a_c$ isolines as a function of We and Bo. Also shown are the maximum value of We for fixed $a_c$ (red dots) and the tangent line to each curve at $We=0$ (blue dashed lines). (b) Derivative $\partial Bo_c/\partial {We}_c|_{{We}=0,a=a_c}$ plotted as a function of Bo (see text for interpretation of this term).

For example, if polymers are introduced to a Newtonian film with $a_c=1.5$, the impact on plug formation is as follows: if the polymers increase Bo by an amount greater than ${We}\times (\partial Bo/\partial {We})|_{{We}=0}\approx 1.0 {We}$, plug formation will be suppressed. If, however, Bo increases by less than this, or is unchanged or decreases, then plug formation will be promoted by the addition of polymers.

Before concluding, we summarize the model results on plug formation and instability classification (absolute versus convective) in figure 10. The combinations of Bo and $a$ values that produce each outcome are shown in figure 10(a) for several values of ${We}$ and small to moderate Bo; figure 10(b) shows this for higher values of Bo, highlighting the dramatic effect an increase in ${We}$ has for large values of Bo. Note also that, for some small Bo (e.g. ${Bo}\approx 0.125$ for ${We}=0$), the AI/CI boundary and the plug/no-plug boundary cross, indicating the possibility of absolutely unstable films that do not form plugs for ${Bo}\lesssim 0.125$.

Figure 10. (a) Boundary separating plug formation from no plug formation as determined by turning points in travelling wave families (red lines), and the boundary separating absolutely unstable films from convectively unstable films (grey lines), shown for three values of ${We}$. (b) Same as panel (a) but for larger values of Bo.

5. Conclusion

We have developed a long-wave asymptotic model for a viscoelastic falling film inside a tube. An upper convected Maxwell model was used to incorporate viscoelasticity; the resulting model equation is identical to that obtained using an Oldroyd-B model (up to first-order). Linear stability analysis allowed the impact of elasticity on the cutoff wavenumber and wavenumber of maximum growth rate to be determined analytically. Similarly, the impact of elasticity on the classification of instabilities as absolute or convective was determined, with elasticity promoting absolute instability. The impact of elasticity was more pronounced as the Bond number increased. In contrast to the Newtonian case, for ${We}>0$, the AI/CI boundary exhibits non-monotonic behaviour with increasing Bo.

Numerical solutions to the model equation show that elasticity promotes plug formation. Travelling wave solution families were found, and the turning points in those families may be used as a proxy for the critical thickness required for plugs to form. By continuation of these turning points, the impact of elasticity on plug formation is easily demonstrated throughout parameter space. In contrast to the Newtonian version of the model (in which plug formation can be suppressed in a film of any thickness, provided the base flow is strong enough), any amount of elasticity creates a maximum critical thickness, past which plugs form regardless of the base flow strength.

The addition of polymers to a Newtonian film may not only introduce elasticity but also modify the surface tension. The trade-off between increasing We (which promotes instability growth and plug formation) and increasing Bo (which inhibits instability growth and plug formation) was explored.

It would be interesting to test the impact of other constitutive equations on these results, particularly on the large impact that increasing ${We}$ has on $a_c$ for large values of ${Bo}$; this is left for future work.

One potential application of this information relates to the airway surface liquid lining human airways. It is well established that this highly viscous fluid is viscoelastic, and the information provided here may be important in accounting for plug formation and rupture in human airway models, specifically in upper airways where gravity may be expected to play a role in the lining's evolution.