2. Problem formulation

Consider a core–annular flow of two immiscible fluids with an incompressible Newtonian fluid in the core region and an incompressible shear-thinning viscoelastic liquid layer in the annular region flowing through a tube of radius $R$. The core fluid possesses a constant dynamic viscosity $\mu _1^*$, density $\rho _1$ and occupies region $r^* \in [0,HR]$, where $0 \leq H \leq 1$. Here, the superscript ‘$*$’ signifies a dimensional quantity, which will be scaled in the following discussion. The ASL has shear-rate-dependent dynamic viscosity $\mu ^*_2$ and relaxation time $\lambda ^*_2$, constant shear modulus $G$, density $\rho _2$ and occupies region $r^* \in [HR, R]$. The schematic of the flow under consideration is shown in figure 1 in terms of the lengths scaled by the radius of the tube.

Figure 1.

Schematic of the model core–annular flow geometry in dimensionless coordinates. Waves at the air–ASL interface signify the development of interfacial instability.

To include the shear-thinning nature of the ASL, the dimensional ASL viscosity ($\mu ^*_2$) and relaxation time ($\lambda ^*$) are assumed to be functions of the second invariant ($\dot \gamma ^{*}$) of the strain-rate tensor $\dot {\boldsymbol {\gamma }}^{*(2)}=(\boldsymbol {\nabla }\boldsymbol {v}^{*(2)})+(\boldsymbol {\nabla }\boldsymbol {v}^{*(2)})^{\rm T}$. Following earlier studies (Wilson, Renardy & Renardy Reference Wilson, Renardy and Renardy1999; Patne Reference Patne2021), we employ power-law model

(2.1)$$\begin{gather} \mu^*_2(\dot \gamma^{*}) = K \dot \gamma^{*(n-1)}, \end{gather}$$

(2.2)$$\begin{gather}\lambda^*_2(\dot \gamma^{*}) =\frac{\mu^*_2(\dot \gamma{*})}{G_f}= \frac{K}{G} \dot \gamma^{*(n-1)}, \end{gather}$$

(2.3)$$\begin{gather}\dot \gamma^*=\sqrt{\tfrac{1}{2} \dot \gamma_{ij}^{*(2)} \dot \gamma_{ji}^{*(2)}}. \end{gather}$$

Here, $K$ is the consistency parameter in the power-law model and $n$ is the power-law index.

The velocities and lengths are scaled by $V_c$ and $R$, respectively, where $V_c$ is the velocity of the core fluid at the centre of the airway. The viscosity and relaxation time are scaled by viscosity scale $\mu _{s}=K ( {V_c}/{R} )^{(n-1)}$ and $\mu _s/G$, respectively. The pressure ($p^{(i)}$) and stress components ($\boldsymbol {\tau }^{(i)}$) are non-dimensionalised by $\mu _s V_c/R$. Thus, the scaled ASL viscosity and relaxation time are

(2.4a)\begin{equation} \mu_2(\dot \gamma)= \lambda(\dot \gamma)= \dot \gamma^{(n-1)}. \end{equation}

Let $\boldsymbol {v}^{(i)}=(v_r^{(i)},v_\theta ^{(i)},v_z^{(i)})$ be the dimensionless velocity fields in fluid where $i=1$ for the core fluid and $i=2$ for the annular fluid. The scaled continuity and Cauchy momentum equations are

(2.4b)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}^{(i)}=0, \end{gather}$$

(2.4c)$$\begin{gather}\rho^{(i)} \,Re \left[\frac{\partial \boldsymbol{v}^{(i)}}{\partial t}+(\boldsymbol{v}^{(i)}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{v}^{(i)} \right]=-\boldsymbol{\nabla} p^{(i)} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau}^{(i)}, \end{gather}$$

where $Re=\rho _2 V_c R/\mu _s$ is the Reynolds number and $\boldsymbol {\nabla }$ is the gradient operator. The parameter $\rho ^{(i)}$ is $\rho ^{(1)} = \rho _r$ for the core fluid and $\rho ^{(2)}=1$ for the annular fluid. The core fluid is a Newtonian fluid; thus we assume $\boldsymbol {\tau }^{(1)} = \mu _r \dot {\boldsymbol {\gamma }}^{(1)}$, where $\mu _r = \mu _1/\mu _s$ and $\dot {\boldsymbol {\gamma }}^{(1)}=(\boldsymbol {\nabla }\boldsymbol {v}^{(1)})+(\boldsymbol {\nabla }\boldsymbol {v}^{(1)})^{\rm T}$ is the strain-rate tensor. For the annular viscoelastic fluid, we employ the WM model (Wilson & Rallison Reference Wilson and Rallison1999; Patne Reference Patne2021)

(2.4d)\begin{equation} \boldsymbol{\tau}^{(2)}+W \lambda(\dot \gamma) \left[ \frac{\partial \boldsymbol{\tau}^{(2)}}{\partial t} + \boldsymbol{v}^{(2)} \boldsymbol{\cdot} \boldsymbol{\tau}^{(2)}- (\boldsymbol{\nabla}\boldsymbol{v}^{(2)})^{\rm T} \boldsymbol{\cdot}\boldsymbol{\tau}^{(2)}-\boldsymbol{\tau}^{(2)} \boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{v}^{(2)}) \right]= \mu_2 (\dot \gamma) \dot{\boldsymbol{\gamma}}^{(2)},\end{equation}

where $W=\mu _s V_c/ (G R)$ is the Weissenberg number.

At the air–ASL interface ($r=H$), the above equations are supplemented with the kinematic boundary condition, velocity continuity and tangential and normal stress balance conditions. No-slip and impermeability conditions are applied at $r=1$, i.e. at the ASL–airway muscle interface.

Note that the oscillatory airflow due to breathing and the branching network of the bifurcating airways will certainly modify the predictions of the present study. These factors have been neglected in the flow considered here in order to simplify the extremely complicated flows in the airways and make the flow amenable to theoretical treatment. Considering an oscillatory airflow necessitates linear Floquet stability analysis, which becomes a cumbersome mathematical exercise when applied to investigate the spatio-temporal stability analysis. The present study aims to provide a theoretical basis for the experimental observations using, admittedly, a simplified model by focusing on a single airway.

For a laminar flow of a Newtonian fluid, the length of the pipe ($L$) necessary to achieve a fully developed flow is (White Reference White2003; Incropera et al. Reference Incropera, DeWitt, Bergman and Lavine2007)

(2.5a)\begin{equation} L \approx 0.06 Re_D D, \end{equation}

where $Re_D = {\rho U D}/{\mu }$ is the Reynolds number based on the diameter $D$ of the pipe and uniform inlet velocity $U$. For the higher-generation airways, $Re_D \sim O(10^{-3}-10)$ implies that a fully developed flow exists in the airways for $L > 0.6 D$. In the Stokes flow ($Re_D \rightarrow 0$) limit, the development length is comparable to the airway diameter, implying the existence of a fully developed flow for $L > D$ (Badr et al. Reference Badr, Dennis, Bates and Smith1985; Stocker & Duck Reference Stocker and Duck1995; Abrahams, Davis & Smith Reference Abrahams, Davis and Smith2008). The airway length is typically six times larger than its radius in higher airway generations (Crystal Reference Crystal1997). Thus, assuming that the above calculations apply to the flow under consideration, a fully developed flow exists for the most part of the airways for the range of $Re_D$ considered here. For a steady-state, axisymmetric, fully developed base state, the quantities are

(2.5b)$$\begin{gather} \bar v_z^{(1)} = \frac{(r^2-H^2)}{4\mu_r} \frac{\partial \bar p}{\partial z} + \frac{n}{2^{1/n} (1+n)} \left[ \left(-\frac{\partial \bar p}{\partial z} \right)^{1/n} - H \left(-H\frac{\partial \bar p}{\partial z} \right)^{1/n}\right], \end{gather}$$

(2.5c)$$\begin{gather}\bar v_z^{(2)} = \frac{n}{2^{1/n} (1+n)} \left[ \left(-\frac{\partial \bar p}{\partial z} \right)^{1/n} - r \left(-r\frac{\partial \bar p}{\partial z} \right)^{1/n}\right], \end{gather}$$

(2.5d)$$\begin{gather}\frac{\partial \bar p}{\partial z} =-\frac{4 \mu_r}{\mu_r + H^2 (1 - \mu_r)}, \end{gather}$$

(2.5e)$$\begin{gather}\bar \tau_{rz}^{(2)}= \bar \mu_2(\bar{\dot{\gamma}}) \frac{{\rm d} \bar v_z^{(2)}}{{\rm d}r}; \quad \bar \tau_{zz}^{(2)}=2W \bar \lambda(\bar{\dot{\gamma}})\bar \tau_{rz}^{(2)} \frac{{\rm d} \bar v_z^{(2)} }{{\rm d}r}, \end{gather}$$

where an overbar signifies a base-state quantity and $\bar \mu _2(\bar {\dot {\gamma }})= \bar \lambda (\bar {\dot {\gamma }})= ( -{{\rm d} \bar v_z^{(2)}}/{{\rm d}r} )^{(n-1)}$.

2.1. Linearised perturbation equations

We impose infinitesimal perturbations on the base state (2.5) and substitute in the governing equations (2.4). The resulting equations are then linearised around the base state (2.5) to obtain the linearised perturbation equations. Following Guillot et al. (Reference Guillot, Colin, Utada and Ajdari2007) and Camassa et al. (Reference Camassa, Ogrosky and Olander2014), here, we consider axisymmetric perturbations. In the linearised perturbation equations, we substitute

(2.6)\begin{equation} {\mathcal{F}}'(\boldsymbol{x},t)= \tilde{{\mathcal{F}}}(r) \, {\rm e}^{{\rm i}(kz -\omega t)}, \end{equation}

where ${\mathcal {F}}'(\boldsymbol {x},t)$ and $\tilde {{\mathcal {F}}}(r)$ are the perturbation to a quantity and the corresponding eigenfunction, respectively. Here, $k = k_r+{\rm i} k_i$ is the wavenumber, and $\omega =\omega _r+{\rm i} \omega _i$ is the complex frequency of the perturbations. The base state (2.5) is temporally unstable if at least one eigenvalue satisfies the condition $\omega _i>0$ for $k_i=0$. For the flow to be absolutely unstable, the requirement will be explained in § 7.

After substituting (2.6) in the linearised governing equations, we obtain

(2.7a)$$\begin{gather} D \tilde v_{r}^{(i)} +\frac{\tilde v_{r}^{(i)} }{r}+{\rm i}k \tilde v_{z}^{(i)} =0, \end{gather}$$

(2.7b)$$\begin{gather}\rho^{(i)}\,Re [ {\rm i}\tilde v_{z}^{(i)} (k\bar v_{z}^{(i)} - \omega)+D \bar v_{z}^{(i)} \tilde v_{r}^{(i)}]=-{\rm i}k\tilde p^{(i)} +{\rm i}k\tilde \tau_{zz}^{(i)} +D\tilde \tau_{rz}^{(i)} + \frac{1}{r} \tilde \tau_{rz}^{(i)}, \end{gather}$$

(2.7c)$$\begin{gather}{\rm i} \rho^{(i)} \,Re \,\tilde v_{r}^{(i)} (k\bar v_{z}^{(i)} -\omega) =-D\tilde p^{(i)} +{\rm i}k\tilde \tau_{rz}^{(i)} +D\tilde \tau_{rr}^{(i)} + \frac{1}{r} (\tilde \tau_{rz}^{(i)} - \tilde \tau_{\theta \theta}^{(i)} ), \end{gather}$$

where $D={\rm d}/{\rm d}r$. For the core fluid, the stress tensor components are

(2.7d)$$\begin{gather} \tilde \tau_{zz}^{(1)} = 2 {\rm i} k \mu_r \tilde v_{z}^{(1)}; \quad \tilde \tau_{rz}^{(1)} = \mu_r( D \tilde v_{z}^{(1)} + {\rm i} k \tilde v_{r}^{(1)}), \end{gather}$$

(2.7e)$$\begin{gather}\tilde \tau_{rr}^{(1)} = 2 \mu_r D \tilde v_{r}^{(1)}; \quad \tilde \tau_{\theta \theta}^{(1)} = \frac{2 \mu_r}{r} \tilde v_{r}^{(1)}. \end{gather}$$

The constitutive equation (2.4d) gives

(2.7f)\begin{gather} W\bar \lambda(\bar{\dot{\gamma}}) [({\rm i}k\bar v_z^{(2)}-\omega) \tilde \tau_{zz}^{(2)}+\tilde v_r^{(2)} D\bar \tau_{zz}^{(2)} -2{\rm i}k \bar \tau_{zz}^{(2)} \tilde v_z^{(2)} -2 D\bar v_z^{(2)} \tilde \tau_{rz}^{(2)} -2 \bar \tau_{rz}^{(2)} D \tilde v_z^{(2)}]\nonumber\\ \quad -\,2 W \bar \tau_{rz}^{(2)} D \bar v_z^{(2)} \left(\frac{{\rm d} \lambda}{{\rm d} \dot \gamma}\right)_{\dot \gamma=\bar{\dot{\gamma}}} [{\rm i}k\tilde v_r^{(2)}+D\tilde v_z^{(2)}]+\tilde \tau_{zz}^{(2)} =2{\rm i}k\bar \mu_2(\bar{\dot{\gamma}}) \tilde v_z^{(2)}, \end{gather}

(2.7g)\begin{gather} W\bar \lambda(\bar{\dot{\gamma}}) [{\rm i}k(\bar v_z^{(2)}-c) \tilde \tau_{rz}^{(2)}+\tilde v_r^{(2)} D\bar \tau_{rz}^{(2)} -{\rm i}k \bar \tau_{rz}^{(2)} \tilde v_z^{(2)}-D\bar v_z^{(2)} \tilde \tau_{rr}^{(2)} -{\rm i}k \bar \tau_{zz}^{(2)} \tilde v_r^{(2)} -\bar \tau_{rz}^{(2)} D\tilde v_r^{(2)}]\nonumber\\ \quad +\,\tilde \tau_{rz}^{(2)}=n ({\rm i}k\tilde v_r^{(2)}+D\tilde v_z^{(2)}), \end{gather}

(2.7h)$$\begin{gather} W\bar \lambda(\bar{\dot{\gamma}}) [{\rm i}k(\bar v_z^{(2)}-c) \tilde \tau_{rr}^{(2)}-2{\rm i}k \bar \tau_{rz}^{(2)} \tilde v_r^{(2)}]+\tilde \tau_{rr}^{(2)} = 2\bar \mu_2(\bar{\dot{\gamma}}) D\tilde v_r^{(2)}, \end{gather}$$

(2.7i)$$\begin{gather}W\bar \lambda(\bar{\dot{\gamma}}) [{\rm i}k(\bar v_z^{(2)}-c) \tilde \tau_{\theta \theta}^{(2)}]+\tilde \tau_{\theta \theta}^{(2)}=2\bar \mu_2(\bar{\dot{\gamma}}) \frac{\tilde v_r^{(2)}}{r}. \end{gather}$$

To solve the above (2.7), we utilise the following boundary conditions. At the centre of the tube, for the axisymmetric disturbances

(2.8a)\begin{equation} u^{(1)}=0 \quad \text{and} \quad \partial_r v^{(1)}=0. \end{equation}

At the air–ASL interface, let the interface be specified by $r=H+h(z,t)$. Here, $h(z,t)$ is the imposed perturbation on the interface. At $r=H$, the relevant boundary conditions include the kinematic boundary condition and the continuity of the velocity and stress components (Wei Reference Wei2005a; Wei & Rumschitzki Reference Wei and Rumschitzki2005; Bassom, Blyth & Papageorgiou Reference Bassom, Blyth and Papageorgiou2012; Gaurav & Shankar Reference Gaurav and Shankar2013)

(2.8b)$$\begin{gather} {\rm i}(k\bar v_z^{(1)}-\omega)\tilde h = \tilde v_r^{(1)}, \end{gather}$$

(2.8c)$$\begin{gather}\tilde v_r^{(1)} = \tilde v_r^{(2)}, \end{gather}$$

(2.8d)$$\begin{gather}\tilde v_z^{(1)} + D\bar v_z^{(1)} \tilde h = \tilde v_z^{(2)} + D\bar v_z^{(2)} \tilde h, \end{gather}$$

(2.8e)$$\begin{gather}\tilde \tau_{rz}^{(1)}=-2{\rm i}kW (D\bar v_z^{(2)})^2 \tilde h + \tilde \tau_{rz}^{(2)}, \end{gather}$$

(2.8f)$$\begin{gather}-\tilde p^{(1)} + \tilde \tau_{rr}^{(1)} +\tilde p^{(2)} - \tilde \tau_{rr}^{(2)}= \frac{1}{Ca} \left( \frac{1}{H^2}-k^2\right) \tilde h, \end{gather}$$

where $h(z,t)=\tilde h \,{\rm e}^{{\rm i} (kz-\omega t)}$ and $Ca=\mu _s V_c/\varSigma _0^*$ is the capillary number. At $r=1$, i.e. at the ASL–muscle interface, we assume no slip and impermeability of the airway wall

(2.8g)\begin{equation} \tilde v_{r}=0, \quad \tilde v_{z}=0. \end{equation}

The present study aims to understand the ASL layer dynamics in the higher-generation airways. Thus, first, we estimate the relevant parameters. In these airways, the thickness of the ASL layer in a healthy lung is ${\sim } O(10^{-6}-10^{-7})$ m while the airway lumen radius is ${\sim } O(10^{-3}-10^{-4})$ m; thus $H \sim 1$ (Widdicombe & Widdicombe Reference Widdicombe and Widdicombe1995; Widdicombe Reference Widdicombe1997; Moriarty & Grotberg Reference Moriarty and Grotberg1999; Montaudon et al. Reference Montaudon, Desbarats, Berger, De Dietrich, Marthan and Laurent2007; Furlow & Mathisen Reference Furlow and Mathisen2018; Karamaoun et al. Reference Karamaoun, Sobac, Mauroy, Van Muylem and Haut2018). The respiratory diseases could substantially increase the thickness of the ASL layer, thus $H$ could assume values much less than unity (Shah et al. Reference Shah, Scott, Knight, Marriott, Ranasinha and Hodson1996; Moriarty & Grotberg Reference Moriarty and Grotberg1999; Heil et al. Reference Heil, Hazel and Smith2008). The power-law index ($n$) for the ASL is 0.1–0.5 (Lai et al. Reference Lai, Wang, Wirtz and Hanes2009).

The density and viscosity of the ASL and air–ASL interfacial tension are ${\sim }10^{3}$ kg m$^{-3}$, ${\sim } 10^{-1}-10^{-2}$ N m$^{-1}$ and ${\sim } 10^{2}-10^{-2}$ Pa${\cdot }$s, respectively (Hwang, Litt & Forsman Reference Hwang, Litt and Forsman1969; Li, Xu & Kumacheva Reference Li, Xu and Kumacheva2000; Schatz & Neitzel Reference Schatz and Neitzel2001; Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Hamed & Fiegel Reference Hamed and Fiegel2014; Hill et al. Reference Hill2014; Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). The maximum velocity of air in the higher-generation airways is $V_c \sim 10^{-3}-1$ m s$^{-1}$ (Pedley Reference Pedley1977). The velocity decreases as the air travels to the higher generation of airways. This leads to a laminar airflow in the higher-generation airways (Lin et al. Reference Lin, Tawhai, McLennan and Hoffman2007; Kleinstreuer & Zhang Reference Kleinstreuer and Zhang2009, Reference Kleinstreuer and Zhang2010; Xia et al. Reference Xia, Tawhai, Hoffman and Lin2010). The order of magnitude of the dimensionless parameters is given in table 1.

Table 1.

Order of magnitude of the dimensionless parameters using the dimensional parameters listed in the text. Note that the viscosity and density of the air are ${\sim } 10^{-5}$ Pa${\cdot }$s and $\sim$1 kg m$^{-3}$.

3. Numerical methodology

To solve the eigenvalue problem (2.7) subjected to the boundary conditions (2.8), we employ the pseudo-spectral method. Each eigenfunctions are expanded in terms of the Chebyshev polynomials

(3.1)\begin{equation} \tilde{f}(y)=\sum_{m=0}^{m=M} a_m {\mathcal{T}}_m (y), \end{equation}

where ${\mathcal {T}}_m(y)$ are Chebyshev polynomials of degree $m$ and $M$ is the highest degree of the polynomial in the series expansion or, equivalently, the number of collocation points. The series coefficients $a_m$ are the unknowns to be solved for. The generalised eigenvalue problem is constructed in the form

(3.2)\begin{equation} \boldsymbol{Ae} + \omega \boldsymbol{Be} =0, \end{equation}

where $\boldsymbol {A}$ and $\boldsymbol {B}$ are the discretised matrices and $\boldsymbol {e}$ is the vector of $a_m$.

The standard procedure for the discretisation procedure of the matrices can be found in Trefethen (Reference Trefethen2000) and Schmid & Henningson (Reference Schmid and Henningson2001). To solve the so-formed eigenvalue problem, we employ the eig MATLAB routine. Along with the genuine eigenvalues, the pseudo-spectral code also predicts numerically spurious eigenvalues. To filter out these spurious eigenvalues, we compare the numerically computed spectrum for $N$ and $N+2$ collocation points with a priori specified tolerance, e.g. $10^{-6}$. The genuine eigenvalues are further verified by increasing the number of collocation points by $25$ and monitoring the variation of the obtained eigenvalues. Whenever the eigenvalue does not change up to a prescribed precision, e.g. to the fifth significant digit, the same number of collocation points is used to determine the critical parameters of the system. The above-described numerical methodology has been previously utilised to solve similar problems by Patne (Reference Patne2021, Reference Patne2024), Patne & Shankar (Reference Patne and Shankar2020) and Patne & Chandarana (Reference Patne and Chandarana2023). To validate the numerical methodology described above, the results predicted by the numerical method will be compared with the analytically obtained solution in § 6.

4. General linear stability analysis

From previous studies of Chen (Reference Chen1991), Hu & Patankar (Reference Hu and Patankar1995), Gaurav & Shankar (Reference Gaurav and Shankar2013), Wei (Reference Wei2005a,Reference Weib), Bassom et al. (Reference Bassom, Blyth and Papageorgiou2012), Joseph et al. (Reference Joseph, Bai, Chen and Renardy1997) and Wei & Rumschitzki (Reference Wei and Rumschitzki2005), the present flow geometry can exhibit the Yih mode due to viscosity stratification, the elastic mode due to the viscoelastic nature of the ASL and the capillary mode due to the interfacial tension. The existence of these modes is demonstrated in figure 2. From the analysis of Yih (Reference Yih1967), the Yih or viscosity stratification mode exhibits a non-zero growth rate for a finite $Re$, which can be positive depending on the viscosity ratio and interface position. From figure 2(a) for $\mu _r<1$ and $H>0.8$, the growth rate is positive. The existence of the capillary mode is demonstrated in figure 2(b). In agreement with the previous studies (Hu & Patankar Reference Hu and Patankar1995; Joseph et al. Reference Joseph, Bai, Chen and Renardy1997; Wei Reference Wei2005a,Reference Weib; Wei & Rumschitzki Reference Wei and Rumschitzki2005; Bassom et al. Reference Bassom, Blyth and Papageorgiou2012; Gaurav & Shankar Reference Gaurav and Shankar2013), the capillary mode exhibits an increasing growth rate for a decreasing capillary number. The elastic mode shown in figure 2(c) arises due to the first normal stress difference across the air–ASL interface (Chen Reference Chen1991). From figure 2, the three instability modes relevant to the present study are long wave. Thus, we proceed with deriving the evolution equation to carry out the stability of the flow.

Figure 2.

Variation in the growth rate $\omega _i$ with the disturbance wavenumber $k$ at $\mu _r=0.1$, $H=0.8$, $n=1$ and $\rho _r=0.001$. Note that $\omega _i >0$ implies an unstable mode. Panel (a) shows the destabilising effect of an increasing Reynolds number $Re$ for $W=0$ and $Ca=0.1$, thereby demonstrating the existence of the Yih mode. Panel (b) demonstrates the capillary mode for $W=0$ and $Re=0$. Panel (c) demonstrates the elastic mode for $Re=0$ and $Ca=0.1$. A finite value of $Ca$ is considered in panels (a) and (c) to remove the unphysical arbitrarily high wavenumber instability.

5. Evolution equation

To derive the long-wave evolution equation describing the dynamics of the liquid layers for the disturbances, let $\lambda$ be the wavelength of the long-wave disturbances such that $R=\epsilon \lambda$, where $\epsilon \ll 1$. Now, following the procedure described in Oron, Davis & Bankoff (Reference Oron, Davis and Bankoff1997), Samanta (Reference Samanta2013), Patne, Agnon & Oron (Reference Patne, Agnon and Oron2021) and Patne & Chandarana (Reference Patne and Chandarana2023), we scale the $z$-axis and time as: $z \to Z \epsilon ^{-1}$ and $t \to \tau \epsilon ^{-1}$. While the fluid dynamical quantities are expanded as

(5.1a)$$\begin{gather} v^{(i)}_z = u_0^{(i)}+ \epsilon u_1^{(i)}+\cdots, \end{gather}$$

(5.1b)$$\begin{gather}v^{(i)}_r = \epsilon v_0^{(i)}+ \epsilon^2 v_1^{(i)}+\cdots, \end{gather}$$

(5.1c)$$\begin{gather}p^{(i)} = \frac{1}{\epsilon} \varPi_0^{(i)}+ \varPi_1^{(i)}+\cdots. \end{gather}$$

The elastic mode arises due to the first normal stress difference across the fluid–fluid interface, while the other elastic terms in the linearised constitutive equations (2.7f)–(2.7i) involving $W$ provide a negligible quantitative correction, as shown in the table 2 in the absence of shear thinning. Due to this, we have neglected elastic terms in the linearised constitutive equations (2.7f)–(2.7i). Thus, we incorporate the elastic effect only through the air–ASL interface tangential stress boundary condition. This assumption considerably simplifies the analysis while retaining the elastic characteristic of the ASL. Note that the analysis demonstrated here is restricted to $n=1$ for the ASL since, for $n \neq 1$, the expressions become unwieldy to show here. For $O(\epsilon )$ calculations, MATHEMATICA fails to solve the equations for $n \neq 1$ due to an extremely high number of terms. Thus, the effect of the shear thinning on the inertial terms will not be shown here.

Table 2.

The most unstable (least stable) numerically computed eigenvalue with/without elastic terms in the linearised constitutive equations (2.7f)–(2.7i) for $Re=0$, $H=0.8$, $W=10$, $\mu _r=0.001$, $n=1$ and $Ca=0.1$. A comparison of the last two columns indicates the negligible role of the elastic terms in the linearised constitutive equations (2.7f)–(2.7i).

At the leading order $O(\epsilon ^0)$, after substitution of the expansions (5.1) into the governing equations, (2.4) reduce to

(5.2a,b)\begin{equation} {-} \partial_Z \varPi^{(i)}_0 + \frac{\mu^{(i)}}{r} \partial_r ( r \partial_r u^{(i)}_0 )=0,\quad \partial_r \varPi^{(i)}_0 = 0. \end{equation}

Equations (5.2) are subjected to the following boundary conditions. At the centre of the tube ($r=0$) and tube wall ($r=1$), the boundary conditions are

(5.3a)$$\begin{gather} \text{at} \ r=0, \quad u^{(1)}_0 \ \text{is finite and} \ \partial_r v^{(1)}_0=0, \end{gather}$$

(5.3b)$$\begin{gather}\text{at} \ r=1, \quad u^{(2)}_0=0 \ \text{and} \ v^{(2)}_0=0. \end{gather}$$

At the fluid–fluid interface $y=\zeta (z,t)$, the continuity of the velocity continuity and stress balance gives

(5.3c)$$\begin{gather} u^{(1)}_0 = u^{(2)}_0, \end{gather}$$

(5.3d)$$\begin{gather}v^{(1)}_0 = v^{(2)}_0, \end{gather}$$

(5.3e)$$\begin{gather}\mu_r \partial_r u^{(1)}_0 - \partial_r u^{(2)}_0=-2\hat{W}_m \partial_Z \zeta, \end{gather}$$

(5.3f)$$\begin{gather}\varPi^{(1)}_0-\varPi^{(2)}_0=-\widehat{Ca}^{-1} \left(\frac{\zeta}{H^2} + \partial^2_Z \zeta \right), \end{gather}$$

where

(5.3g)\begin{equation} \hat{W}_m=\epsilon W (D \bar v_z^{(2)})^2 \quad \text{and} \quad \widehat{Ca}=\epsilon^{-3}Ca, \end{equation}

are the scaled Weissenberg and capillary numbers, respectively.

The solution of the governing equations (5.2) using boundary conditions (5.3) leads to

(5.4a)$$\begin{gather} u^{(1)}_0 = a_1 + a_2 \ln r + \frac{1}{2 \mu_r} \partial_Z\varPi^{(1)}_0, \end{gather}$$

(5.4b)$$\begin{gather}u^{(2)}_0 = a_3 + a_4 \ln r + \frac{1}{2} \partial_Z \varPi^{(2)}_0, \end{gather}$$

where $a_i$ are the integration constants

(5.4c)\begin{gather} a_1=\frac{a'_1}{4 Ca \, H^2 (-\mu_r + (-1 + \mu_r) \zeta^4)}, \end{gather}

(5.4d) \begin{gather} a'_{1} =-16 Ca \,H^2 \mu_r q - 8 \zeta \, Ca \, H^2 \mu_r \hat{W}_m (1 + \ln \zeta) \partial_Z \zeta \nonumber\\ \quad +\,8 \zeta^3 \,Ca \, H^2 (-1 + 2 \mu_r) \hat{W}_m \partial_Z \zeta + 8 \zeta^5 \,Ca \, H^2 (-1 + \mu_r) \hat{W}_m (-1 + \ln \zeta) \partial_Z \zeta \nonumber\\ \quad +\,(-3 + 4 \mu_r) \zeta^4 (\partial_Z \zeta + H^2 \partial^3_Z \zeta) + 2 (-1 + \mu_r) (-1 + \ln \zeta) \zeta^6 (\partial_Z \zeta + H^2 \partial^3_Z \zeta)\nonumber\\ \quad +\,(1 - 2 \mu_r - 2 \mu_r \ln \zeta) \zeta^2 \partial_Z \zeta + H^2 \zeta^2 (16 Ca\, (-1 + \mu_r) q + (1 - 2 \mu_r - 2 \mu_r \ln \zeta)\partial^3_Z \zeta, \end{gather}

(5.4e)\begin{gather} a_2=0, \end{gather}

(5.4f)\begin{gather} a_3 = \frac{\partial_Z \zeta}{4 Ca \, H^2} + \frac{\partial^3_Z \zeta}{4Ca} + \frac{\mu_r H^2 (-16 Ca \,q + (-1 + \zeta^2)^2 \partial^3_Z \zeta)}{4 Ca \, H^2 (-\mu_r + (-1 + \mu_r) \zeta^4)} \nonumber\\ \quad +\,\frac{\mu_r [(-1 + \zeta^2) (-1 + 8 Ca \, H^2 \hat{W}_m \zeta + \zeta^2) \partial_Z \zeta]}{4 Ca \, H^2 (-\mu_r + (-1 + \mu_r)\zeta^4)}, \end{gather}

(5.4g)\begin{gather} a_4 = \frac{\zeta^2 }{2 Ca \, H^2} (\partial_Z \zeta + H^2 \partial^3_Z \zeta) + 2 \hat{W}_m \zeta \partial_Z \zeta, \end{gather}

where $q$ is the volumetric flow rate obtained using the base-state velocity profiles (2.5).

The rescaled kinematic boundary condition at the interface up to $O(1)$ is

(5.5)\begin{equation} \partial_\tau \zeta =-\frac{1}{\zeta} \partial_Z \int_0^\zeta r u^{(1)}_0 \,{\rm d}r. \end{equation}

On substitution of the velocity profiles derived above into (5.5), the interface evolution equation can be obtained

(5.6)\begin{equation} {\mathcal{F}}_1 (\zeta,\partial_Z \zeta,\partial_Z^2 \zeta,\partial_Z^3 \zeta,\partial_Z^4 \zeta) + {\mathcal{F}}_2 (\zeta) \partial_t \zeta=0, \end{equation}

where ${\mathcal {F}}_1 (\zeta,\partial _Z \zeta,\partial _Z^2 \zeta,\partial _Z^3 \zeta,\partial _Z^4 \zeta )$ and ${\mathcal {F}}_2 (\zeta )$ are known functions.

To include the effect of the inertia, we proceed with $O(\epsilon )$ calculations. Adding the inertial effects leads to a better agreement between the long-wave model and numerically computed eigenvalues. At $O(\epsilon )$ the governing equations lead to

(5.7a)$$\begin{gather} \rho_r \,Re ( \partial_\tau u^{(i)}_0+ u^{(i)}_0 \partial_Z u^{(i)}_0 + v^{(i)}_0 \partial_r u^{(i)}_0 ) =-\partial_Z \varPi^{(i)}_1 + \frac{\mu^{(i)}}{r} \partial_r ( r \partial_r u^{(i)}_1 ), \end{gather}$$

(5.7b)$$\begin{gather}\partial_r \varPi^{(i)}_1 = 0. \end{gather}$$

The radial velocity $v^{(i)}_0$ can be derived using $u^{(i)}_0$ by using the continuity equation at $O(1)$. The boundary conditions at $O(\epsilon )$ are

(5.8a)$$\begin{gather} \text{at} \ r=0, \quad u^{(1)}_1 \ \text{is finite and} \ \partial_r v^{(1)}_1=0, \end{gather}$$

(5.8b)$$\begin{gather}\text{at} \ r=1, \quad u^{(2)}_1=0 \ \text{and} \ v^{(2)}_1=0. \end{gather}$$

At the fluid–fluid interface, the interface conditions are

(5.8c)$$\begin{gather} u^{(1)}_1 = u^{(2)}_1, \end{gather}$$

(5.8d)$$\begin{gather}v^{(1)}_1 = v^{(2)}_1, \end{gather}$$

(5.8e)$$\begin{gather}\mu_r \partial_r u^{(1)}_1 - \partial_r u^{(2)}_1=0, \end{gather}$$

(5.8f)$$\begin{gather}\varPi^{(1)}_1-\varPi^{(2)}_1=0. \end{gather}$$

The system of (5.7) is solved using the boundary conditions (5.8) to obtain the $O(\epsilon )$ contribution. At $O(\epsilon )$, equations turn out to be quite large; thus, we do not provide explicit expressions for the quantities. The rescaled kinematic boundary condition at the interface up to $O(\epsilon )$ is

(5.9)\begin{equation} \partial_\tau \zeta =-\frac{1}{\zeta} \partial_Z \int_0^\zeta r (u^{(1)}_0+ \epsilon u^{(1)}_1) \,{\rm d}r. \end{equation}

Using the above kinematic boundary condition, we obtain the evolution equation for the thickness of the layers. Before proceeding further, we scale out the $\epsilon$ parameter by rescaling the obtained equation. Thus, the parameters used below will be without hat subscript. The subsequent temporal and spatio-temporal stability analysis will be performed using the rescaled evolution equation.

6. Temporal stability analysis

As mentioned previously, due to the extreme density and viscosity stratification and the small thickness of the ASL in the higher-generation airways, the results from the previous studies dealing with the core–annular flow will not be readily relevant. Thus, first, we determine the relevance of the three modes of instability discussed in § 4 as follows. We linearise the evolution equation by substituting $\zeta =H + \zeta '$, where $\zeta ' \ll 1$. Next, we use normal modes $\zeta ' = \tilde \zeta \,{\rm e}^{{\rm i}(kz-\omega t)}$ in the resulting linearised equation to obtain the dispersion relation

(6.1)\begin{align} \omega&=k f_0(H,\mu_r,n)+{\rm i} k^2 \left[ f_1 (H,\mu_r,\rho_r,n) \, Re + f_2 (H,\mu_r,n) W + \frac{f_3 (H,\mu_r,n)}{Ca} \right]\nonumber\\ &\quad - {\rm i} H^2 k^4 \frac{f_3 (H,\mu_r,n)}{Ca}. \end{align}

Here, the volumetric flow rate $q$ has been eliminated by using relation $q=\int _0^H \bar v_z^{(1)} r \,{\rm d}r+\int _H^1 \bar v_z^{(2)} r\,{\rm d}r$. The functions $f_i$ are

(6.2a)$$\begin{gather} f_0 = \frac{(H^2-1) \mu_r}{H^4 (\mu_r -1) - \mu_r}, \quad f_1=\frac{- H^4 ( \mu_r -1) g_1}{g_2}, \quad f_2 = \frac{2 H^4 \mu_r^2 g_3}{g_4}, \end{gather}$$

(6.2b)$$\begin{gather}f_3 = \frac{H}{16} \left[ \frac{g_5}{H^4 ( \mu_r -1) - \mu_r} - 4 \ln{H} \right], \end{gather}$$

and

(6.2c) \begin{align} g_1 &=-(H^2 -1) \{\mu_r^2 (2 \mu_r^2 - H^{12} ( \mu_r -1) (-13 + 22 \mu_r) \nonumber\\ &\quad +H^2 \mu_r (-9 + 56 \mu_r) + H^{10} ( \mu_r -1) (-83 + 152 \mu_r) +\cdots\nonumber\\ &\quad H^4 (-24 + (217 - 314 \mu_r) \mu_r) +\cdots\nonumber\\ &\quad H^8 (-193 + 612 \mu_r - 434 \mu_r^2) +5 H^6 (27 + 4 \mu_r (-31 + 28 \mu_r))) +\cdots\nonumber\\ &\quad 3 H^2 (H^8 (-3 + 7 H^2) + H^4 (5 - 29 H^2 + 59 H^4 - 35 H^6) \mu_r +\cdots\nonumber\\ &\quad 2 H^2 (H^2 -1)^2 (6 - 33 H^2 + 35 H^4) \mu_r^2 -\cdots\nonumber\\ &\quad 6 (H^2 -1)^4 ( 11 H^2 -1) \mu_r^3 +24 (H^2 -1)^5 \mu_r^4) \rho_r \}\nonumber\\ &\quad +24 H^2 [H^4 ( \mu_r -1) - \mu_r] \{3 \mu_r^3 (\rho_r - 1 ) \cdots\nonumber\\ &\quad +6 H^2 \mu_r^2 (-1 + 2 \mu_r + \rho_r - 2 \mu_r \rho_r) -\cdots\nonumber\\ &\quad 6 H^6 ( \mu_r -1) \mu_r (-\mu_r - \rho_r + 2 \mu_r \rho_r) +\cdots\nonumber\\ &\quad H^8 ( \mu_r -1) (\rho_r + \mu_r (-\mu_r + 3 ( \mu_r -1) \rho_r)) +\cdots\nonumber\\ &\quad 2 H^4 \mu_r (\rho_r + \mu_r (6 - 7 \mu_r + 9 ( \mu_r -1) \rho_r))\} \ln{H}, \end{align}

(6.2d)$$\begin{gather} g_2 = 144 [H^4 ( \mu_r -1) - \mu_r]^3 [-H^2 ( \mu_r -1) + \mu_r]^2, \end{gather}$$

(6.2e)$$\begin{gather}g_3 =-H^2 + H^4 - 2 (H^2 -1)^2 \mu_r + 2 (H^4 ( \mu_r -1) - \mu_r) \ln{H}, \end{gather}$$

(6.2f)$$\begin{gather}g_4 = [H^4 ( \mu_r -1) - \mu_r] [-H^2 ( \mu_r -1) + \mu_r]^2, \end{gather}$$

(6.2g)$$\begin{gather}g_5 = (H^2 -1) [1 - 4 \mu_r + H^2 (4 \mu_r -3 )]. \end{gather}$$

The function $f_0 (H,\mu _r,n)$ is the wave speed of the perturbations. The functions $f_0 (H,\mu _r,n), f_1 (H,\mu _r,\rho _r,n)$ and $f_3 (H,\mu _r,n)$ are in agreement with the long-wave asymptotic analysis of Gaurav & Shankar (Reference Gaurav and Shankar2013) for $\rho _r=1$, thereby validating the present analysis.

For $n \neq 1$, the functions $f_0 (H,\mu _r,n), f_2 (H,\mu _r,\rho _r,n)$ and $f_3 (H,\mu _r,n)$ are

(6.3a)$$\begin{gather} f_0 = \frac{2^{-2 - 1/n} (H^2 -1)}{\mu_r - H^4 \mu_r + H^4 n} g'_0; \quad f_2 = \frac{2^{-2 - 1/n} H^{3 + 1/n} g'_3}{g'_4}, \end{gather}$$

(6.3b)$$\begin{gather}f_3 =-\frac{H}{16n} \left[ \frac{g'_5}{\mu_r - H^4 \mu_r + H^4 n} + 4 \ln{H} \right], \end{gather}$$

where

(6.3c)$$\begin{gather} g'_0 =-2^{1 + 1/n} \mu_r H (H^2 -1) \left(-H \frac{\partial \bar p}{\partial z} \right)^{1/n}\frac{\partial \bar p}{\partial z} [H^2 + \mu_r + H^4 (-1 - \mu_r + n)], \end{gather}$$

(6.3d)$$\begin{gather}g'_3 = \left(-\frac{\partial \bar p}{\partial z} \right)^{1/n+1} [ (H^2 -1) (2 (H^2 -1) \mu_r - H^2 n) +2 (\mu_r - H^4 \mu_r + H^4 n) \ln{H} ], \end{gather}$$

(6.3e)$$\begin{gather}g'_4 = n (\mu_r - H^4 \mu_r + H^4 n), \end{gather}$$

(6.3f)$$\begin{gather}g'_5 = (H^2 -1) [n - 4 \mu_r + H^2 (4 \mu_r -3n )], \end{gather}$$

(6.3g)$$\begin{gather}\frac{\partial \bar p}{\partial z} =-\frac{4 \mu_r}{\mu_r + H^2 (1 - \mu_r)}. \end{gather}$$

The above equations reduce to (6) for $n=1$. As stated above, owing to computational difficulty, the effect of shear thinning on the inertial term will not be explored. An examination of terms in (6.1) provides valuable insight into three modes of instability as follows. The first term, $f_1 (H,\mu _r, \rho _r,n) Re$, is clearly a function of the density and viscosity stratification and leads to the Yih mode. The second term, $f_2 (H,\mu _r,n) W$, arises due to the first normal stress difference across the air–ASL interface. Interestingly, this term is also a function of the viscosity stratification, as indicated by the presence of $\mu _r$ in $f_2 (H,\mu _r,n)$. The elastic instability mode predicted by Chen (Reference Chen1991) for the core–annular flow of viscoelastic fluids is due to $f_2 (H,\mu _r,n) W$. The third term, ${f_3 (H,\mu _r,n)}/{Ca}$, arises due to the air–ASL interfacial tension and is responsible for the capillary mode. Additionally, all three functions depend on the power-law index, underlining the importance of the ASL shear thinning on the instabilities.

To demonstrate the utility of the evolution equation, we compare the growth rate predicted by the numerical and long-wave analyses in figure 3. Excellent agreement between the results predicted by both approaches for $k<0.2$ certifies the applicability of the evolution equation. Since the analytical results have been validated by directly comparing the functions with Gaurav & Shankar (Reference Gaurav and Shankar2013), figure 3 also validates the present numerical methodology.

Figure 3.

The variation of the growth rate ($\omega _i$) with the wavenumber ($k$) for $Re=10$, $W=10$, $Ca=0.1$, $n=1$, $\mu _r=0.1$, $H=0.8$ and $\rho _r=1$. The figure demonstrates excellent agreement between the numerical and asymptotic analyses for $k<0.2$, thereby demonstrating the ability of the long-wave analysis. The semilog axis scale is used to magnify the low wavenumber region of the dispersion curve.

The variation of $f_1(H,\mu _r,\rho _r,n)$ with $H,\rho _r$ and $\mu _r$ is shown in figure 4. For the Yih mode to be unstable $f_1(H,\mu _r,\rho _r,n) >0$, since $Re>0$ in the present analysis. For $\mu _r<1$, for low values of $H$, $f_1(H,\mu _r,\rho _r,n)<0$ and the opposite case for high $H$, as shown in figure 4(a) in the absence of the density stratification, thus demonstrating the existence of the stable and unstable zones. From figure 4(b), the opposite scenario regarding $f_1(H,\mu _r,\rho _r,n)$ values is predicted for $\mu _r>1$.

Figure 4.

Variation of function $f_1(H,\mu _r,\rho _r,n)$ with the relative position of the liquid–liquid interface $H$ and the viscosity ratio $\mu _r$ for $n=1$. Note that $f_1(H,\mu _r,\rho _r,n)>0$ implies an unstable Yih mode. Panels (a) and (b) show the presence of stable and unstable zones for $\mu _r<1$ and $\mu _r>1$, respectively. Panel (c) demonstrates the effect of variation in $\rho _r$ on stable and unstable zones. Note the absence of stable zones for low $\rho _r$. Panel (d) shows the $f_1(H,\mu _r,\rho _r,n)$ values for the relevant parameters for the core–annular flow in the higher-generation airways.

For the present problem, the inner fluid layer is air while the outer layer is ASL, thus $\rho _r \sim 10^{-3}$. From figure 4(c), $f_1(H,\mu _r,\rho _r,n)>0$ for the whole range of $H$, implying the removal of the stable zone due to the density stratification. Depending on the diseases, the viscosity of the ASL also varies, and $H$ could differ substantially from unity (Moriarty & Grotberg Reference Moriarty and Grotberg1999; Hamed & Fiegel Reference Hamed and Fiegel2014). Thus, for these relevant parameters, $f_1(H,\mu _r,\rho _r,n)$ values are shown in figure 4(d), where $f_1(H,\mu _r,\rho _r,n)>0$ irrespective of parametric combination. This implies that the Yih mode will lead to an unstable flow in the higher-generation airways.

The variation of $f_2(H,\mu _r,n)$ and $f_3(H,\mu _r,n)$ with $H$ and $\mu _r$ is shown in figure 5. Note that $f_2(H,\mu _r,n)$ is the coefficient of $W$, thus representing elastic mode, while $f_3(H,\mu _r,n)$ is the coefficient of $1/Ca$, thus representing capillary mode. Similar to $f_1(H,\mu _r,\rho _r,n)$, $f_2(H,\mu _r,n)$ also exhibits negative values for low $H$, as shown in figure 5(a). The elastic and capillary modes are unstable if $f_2(H,\mu _r,n)>0$ and $f_3(H,\mu _r,n)>0$, respectively. For the parameters relevant to the higher-generation airways, from figure 5(b) $f_2(H,\mu _r,n)>0$, implying the destabilising role of the ASL elasticity in the base flow. From figure 5(c), the capillary mode is unstable for an arbitrary value of $H$. The value of $f_3(H,\mu _r,n)$ for the parameters relevant to the core–annular flow in the higher-generation airways is shown in figure 5(d). To conclude from figures 4 and 5, the viscosity and density stratification, ASL elasticity and capillarity lead to the destabilisation of the core–annular flow in the higher-generation airways.

Figure 5.

Variation of functions $f_2(H,\mu _r,n)$ and $f_3(H,\mu _r,n)$ with the relative position of the liquid–liquid interface $H$ and the viscosity ratio $\mu _r$ for $n=1$. Panel (a) demonstrates the existence of the stable and unstable zones for the elastic mode similar to $f_1(H,\mu _r,\rho _r,n)$. Panel (c) shows the existence of only an unstable zone for the capillary mode. Panels (b) and (d) show the values of the functions $f_2(H,\mu _r,n)$ and $f_3(H,\mu _r,n)$ relevant to the higher-generation airways.

The effect of shear thinning on the elastic mode is shown in figure 6(a). For $\mu _r<0.4$, an increasingly shear-thinning annular fluid has a strong stabilising effect on the elastic mode and the opposite for $\mu _r >0.4$. For the capillary mode, as the shear-thinning index decreases and the viscosity of the ASL increases, the value of the function $f_3(H,\mu _r,n)$ increases, indicating a destabilising effect. This aspect is shown in figure 6(b) for select values of the power-law index $n$.

Figure 6.

Variation of functions $f_2(H,\mu _r,n)$ (a) and $f_3(H,\mu _r,n)$ (b) with the relative viscosity $\mu _r$ for select values of power-law index $n$ at $H=0.9$. The figure demonstrates the stabilising effect of the shear-thinning nature of the ASL on the elastic mode and the opposite on the capillary mode.

7. Spatio-temporal stability analysis ($k_i \neq 0$)

In §§ 4 and 6, we assumed a complex frequency and real wavenumber, implying a temporal linear stability analysis. The flow is temporally unstable if the imaginary part of the frequency, $\omega _i >0$, for at least one eigenvalue. Thus, temporal stability analysis provides information only about the temporal growth of the perturbations. Similarly, one can assume real frequency and complex wavenumber to carry out the spatial stability analysis. The flow will be spatially unstable if the imaginary part of the wavenumber, $k_i<0$, for at least one eigenvalue. Thus, spatial stability provides information about the spatial growth of the perturbations. To analyse the growth of the perturbations in space and time, we must treat both frequency and wavenumber as complex numbers. The temporal and spatial stability analyses are straightforward and can be readily deduced from the eigenspectra. To analyse the spatio-temporal stability analysis, one can utilise methods developed by Briggs (Reference Briggs1964) and Kupfer et al. (Reference Kupfer, Bers and Ram1987).

The spatio-temporal stability analysis classifies temporal instability as an absolute or convective instability depending on the value of $\omega _i$ at the cusp point. An absolutely unstable mode implies the development of the corresponding mode in up- and down-stream directions while a convectively unstable mode signifies that the disturbances develop only in the downstream direction from the source of the disturbances (Briggs Reference Briggs1964). To demonstrate the absolute and convective instabilities, we assume an arbitrary dispersion relation $\omega ={\mathcal {F}}(k)$, where ${\mathcal {F}}(k)$ is a continuous and differentiable function of $k$. Next, we differentiate the dispersion relation with respect to $k$ to obtain the group velocity of the perturbations. To determine the existence of absolute or convective instability, the group velocity of the disturbances is equated to zero so that ${\partial \omega }/{\partial k}=0$ (Huerre & Monkewitz Reference Huerre and Monkewitz1990; Schmid & Henningson Reference Schmid and Henningson2001). The solution of this equation gives $k$ values at which saddle points exist (denoted by $k_0$) in the $k$ or $k_r - k_i$ plane. The calculated saddle point must also obey the causality principle for the existence of the absolute instability (Huerre & Monkewitz Reference Huerre and Monkewitz1990). For the first-order saddle point in the $k$-plane, a local Taylor series expansion about this point gives $(\omega -\omega _{0})\sim (k-k_{0})^{2}$, where $\omega _{0}$ is the cusp point in the $\omega$ or $\omega _r - \omega _i$ plane. This relationship implies that the mapping from the $k$-plane to the $\omega$-plane is characterised by angle doubling (i.e. phase doubling). The solution of the equation ${\partial \omega }/{\partial k}=0$ could lead to multiple solutions in terms of $k$. Among these solutions, if at least one saddle point satisfies the condition $\omega _{0i}>0$, then the flow is absolutely unstable else convectively unstable (Kupfer et al. Reference Kupfer, Bers and Ram1987; Huerre & Monkewitz Reference Huerre and Monkewitz1990).

Briggs (Reference Briggs1964) devised a method to determine the saddle point in the $k$-plane and the cusp point in the $\omega$-plane by using analytic continuation. This method requires systematic movement of the curves in both planes from the temporal stability curves. This necessitates a considerable computation, which could become unwieldy for realistic flows. Kupfer et al. (Reference Kupfer, Bers and Ram1987) suggested an alternative method in which only the cusp point formation in the $\omega$-plane is accomplished by varying $k_i$ value. If, for a finite $k_i$, the cusp point in the $\omega$-plane possesses $\omega _{i0}>0$ and follows causality, then the flow is absolutely unstable. Given the relative ease of the Kupfer et al. (Reference Kupfer, Bers and Ram1987) method, it will be utilised in the present study.

The Kupfer et al. (Reference Kupfer, Bers and Ram1987) procedure to determine the genuine cusp point is as follows. We begin with the temporal curve in the $\omega$-plane. This curve corresponds to $k_i=0$ and varying $k_r$. Next, a finite negative value is assumed for $k_i$, and $k_r$ varies similarly to the temporal curve. The same procedure is repeated until a cusp point forms in the $\omega$-plane. To determine whether the obtained cusp point follows causality, a straight line parallel to the $\omega _i$ axis is drawn from the cusp point until it crosses the temporal curve. If the line crosses the temporal curve an odd number of times, then the obtained cusp point follows causality and is termed a genuine cusp point. Otherwise, it is an evanescent cusp point (Yeo, Khoo & Zhao Reference Yeo, Khoo and Zhao1999, Reference Yeo, Khoo and Zhao2001; Patne & Shankar Reference Patne and Shankar2017; Patne & Ramon Reference Patne and Ramon2020; Patne & Chandarana Reference Patne and Chandarana2023).

To obtain the saddle point, we differentiate the dispersion relation (6.1) with respect to $k$ and solve the resulting equation for $k$ to obtain three roots

(7.1a)$$\begin{gather} k_0 =-\frac{{\rm i} F_1}{3^{1/3} \varDelta} + \frac{{\rm i} \varDelta}{2\, 3^{2/3} F_2}, \end{gather}$$

(7.1b)$$\begin{gather}k_0 = \frac{({\rm i} - \sqrt{3}) F_1}{2\, 3^{1/3}\varDelta} - \frac{({\rm i} + \sqrt{3}) \varDelta}{4\, 3^{2/3} F_2}, \end{gather}$$

(7.1c)$$\begin{gather}k_0=\frac{({\rm i} + \sqrt{3}) F_1}{2\, 3^{1/3}\varDelta} - \frac{({\rm i} - \sqrt{3}) \varDelta}{4\, 3^{2/3} F_2}, \end{gather}$$

where

(7.1d)$$\begin{gather} \varDelta = \left(9 f_0 F_2^2 + \sqrt{3} \sqrt{8 F_1^3 F_2^3 + 27 f_0^2 F_2^4}\right)^{1/3}, \end{gather}$$

(7.1e)$$\begin{gather}F_1 = f_1 (H,\mu_r,\rho_r,n) Re + f_2 (H,\mu_r,n) W + \frac{f_3 (H,\mu_r,n)}{Ca} + f_4 (H,\mu_r) Ma, \end{gather}$$

(7.1f)$$\begin{gather}F_2 = H^2 \frac{f_3 (H,\mu_r,n)}{Ca}. \end{gather}$$

The first root (7.1a) is purely imaginary since the factors $F_1, F_2$ and $\varDelta$, are real numbers, thus a spurious saddle point. The second root (7.1b) will be a complex number with a negative real part. The third root (7.1c) is also a complex number but with a positive real part. In § 4, we considered positive real values of $k$ while analysing the temporal stability analysis. Thus, the third root (7.1c) possessing a positive real part is the relevant root. Note that the second root (7.1b) will provide the same cusp point, albeit for negative values of $k_r$. Thus, the conclusions obtained using the third root (7.1c) will be applicable to the second root (7.1b) as well. To obtain the cusp point $\omega _0$, the third root (7.1c) is substituted into the dispersion relation (6.1). Next, we investigate the effect of capillarity, inertia, ASL elasticity and viscosity in introducing an absolute instability.

7.1. Effect of capillarity ($Ca$)

The cusp points for select values of $Ca$ and parameters relevant to the higher generation airways are shown in figure 7. Note that we will neglect the effect of inertia and ASL elasticity in this section to focus on the capillary effect. The cusp points are shown here to demonstrate the genuine nature of the cusp point as follows. To determine the genuine cusp point, a straight line is drawn from the cusp point parallel to the $\omega _i$ axis. From figure 7, for the values of $Ca$ considered here, the line crosses the temporal curve only once, i.e. an odd number of times. From the above discussion, this implies the cusp points shown in figure 7 are genuine.

Figure 7.

Cusp point formation in the $\omega$-plane at $\mu _r=0.001$, $H=0.95$, $Re=0$, $n=1$ and $W=0$. The locations of the saddle and cusp points in the $k$ and $\omega$ planes for all panels are as follows: (a) $k_0= 0.7444 - 0.006477 i$, $\omega _0= 0.00008908 + 0.001279 i$; (b) $k_0= 0.7523 - 0.06299 i$, $\omega _0= 0.00008940 + 0.0001241 i$; (c) $k_0= 0.9565 - 0.3468 i$, $\omega _0= 0.0001012 - 0.00001278 i$; (d) $k_0= 1.754 - 0.9167 i$, $\omega _0= 0.0001648 - 0.00007712 i$. Note that, from dispersion relation (6.1), the flow is temporally unstable for a finite value of $Ca$.

The role of capillary instability in causing absolute instability has been previously demonstrated by Guillot et al. (Reference Guillot, Colin, Utada and Ajdari2007, Reference Guillot, Colin and Ajdari2008) and Utada et al. (Reference Utada, Fernandez-Nieves, Stone and Weitz2007). This does not imply an absolutely unstable flow for an arbitrary value of $Ca$, as shown in figure 7. For low $Ca$, the flow is absolutely unstable and convectively unstable for large values of $Ca$. Thus, there exists a critical value of the capillary number ($Ca_c$) for given values of $H$ and $\mu _r$ beyond which the flow is convectively unstable. For the parameters specified in figure 7, $Ca_c=0.6585$. From the definition of $Ca$, this implies a minimum ASL viscosity and/or core velocity is necessary to avoid absolute instability. Drawing an analogy with the experiments of Guillot et al. (Reference Guillot, Colin, Utada and Ajdari2007) and Utada et al. (Reference Utada, Fernandez-Nieves, Stone and Weitz2007), the absolute instability predicted here could lead to the formation of air bubbles and ASL plugs in the airways. These ASL plugs could function similarly to the ASL bridges or plugs formed due to film collapse to result in airway closure. Thus, the information about the $Ca_c$ could be utilised to predict the airway closure triggered by the absolute instability. Pulmonary diseases lead to an increase in the viscosity and thickness of the ASL (Hwang et al. Reference Hwang, Litt and Forsman1969; Shah et al. Reference Shah, Scott, Knight, Marriott, Ranasinha and Hodson1996; Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). In the core–annular flow under consideration, the dimensionless thickness of the ASL layer is $1-H$. Thus, we determine the variation of $Ca_c$ with $1-H$ and $\mu _r$ in figure 8.

Figure 8.

The variation of the critical capillary number $Ca_c$ with the dimensionless thickness of the ASL layer $1-H$ for $Re=0$, $n=1$ and $W=0$. The curves exhibit scaling $Ca_c \sim (1-H)^2$.

For a healthy lung, in the higher-generation airways, the thickness of the ASL layer is ${\sim } O(10^{-6}-10^{-7})$ m while the airway lumen radius is ${\sim } O(10^{-3}-10^{-4})$ m; thus $H \sim 0.999$ (Widdicombe & Widdicombe Reference Widdicombe and Widdicombe1995; Widdicombe Reference Widdicombe1997; Moriarty & Grotberg Reference Moriarty and Grotberg1999; Montaudon et al. Reference Montaudon, Desbarats, Berger, De Dietrich, Marthan and Laurent2007; Furlow & Mathisen Reference Furlow and Mathisen2018; Karamaoun et al. Reference Karamaoun, Sobac, Mauroy, Van Muylem and Haut2018). This corresponds to $1-H = 0.001$ in figure 8. For $H \sim 0.999$, the calculations show that the flow is absolutely unstable provided that $Ca< 0.00027$. To achieve such a low $Ca={\mu _s V_c}/{\varSigma }$, the velocity of the air ($V_c$) must drop significantly since the ASL viscosity and air–ASL interfacial tension are constant. Such a severe drop in $V_c$ could occur towards the end of expiration in the breathing cycle (Pedley Reference Pedley1977). Thus, even in a healthy lung, absolute instability can exist, which could lead to airway closure in the airways during the end phase of the expiration in agreement with the previous studies (Burger & Macklem Reference Burger and Macklem1968; Heil et al. Reference Heil, Hazel and Smith2008). These liquid blockages normally rupture in the beginning phase of inspiration. As inspiration begins, the $Ca$ value will increase and flow once again becomes convectively unstable, thereby prohibiting the airway closure.

In a lung suffering from diseases (e.g. cystic fibrosis and pulmonary oedema), there could be an increase in the thickness and viscosity of the ASL mainly due to the alteration in the properties and increasing thickness of the mucus layer (Hwang et al. Reference Hwang, Litt and Forsman1969; Shah et al. Reference Shah, Scott, Knight, Marriott, Ranasinha and Hodson1996; Heil et al. Reference Heil, Hazel and Smith2008; Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). Previous studies have indicated a higher chance of airway closure in diseased lungs (Heil et al. Reference Heil, Hazel and Smith2008; Halpern et al. Reference Halpern, Fujioka and Grotberg2010; Erken et al. Reference Erken, Fazla, Muradoglu, Izbassarov, Romanò and Grotberg2023). An increase in the ASL thickness and viscosity leads to the exacerbation of the absolute instability in two ways. An increasing ASL thickness implies a high value of $1-H$. From figure 8, as the ASL thickness increases, larger values of $Ca$ become accessible for the flow to exhibit absolute instability. Simultaneously, an increasing viscosity of the ASL implies a decreasing $\mu _r$. From figure 8, a decreasing $\mu _r$ shifts the absolute instability boundary to higher $Ca$. Thus, flow can exhibit absolute instability even at a higher $Ca$. To conclude, a pulmonary disease affecting the ASL layer enhances the possibility of an absolutely unstable flow in the higher-generation airways, thereby increasing the possibility of airway closure, in excellent agreement with the experimental and clinical observations (Mall Reference Mall2016).

Also, from figure 8, if one were to plot $\mu _r Ca_c$ instead of $Ca_c$, then the three curves would overlap, indicating the scaling $Ca_c \sim 1/\mu _r$. Thus, the critical capillary number for the flow to become convectively unstable scales with the dimensionless thickness of the ASL and viscosity ratio as $Ca_c \sim (1-H)^2/\mu _r$.

7.2. Effect of inertia ($Re$)

Due to considerable air velocity in a certain generation of airways (Pedley Reference Pedley1977), there are appreciable inertial effects, as indicated by high values of $Re$ in table 1. Thus, in this section, we analyse the effect of a finite $Re$ on the existence of absolute instability. From figure 7(c), the flow is convectively unstable for $Ca=1$ for $Re=0$. If we increase $Re$ to $100$ while keeping other parameters the same as in figure 7(c), then the cusp point changes to $\omega _0 = 0.0001172 + 9.821 \times 10^{-6} i$, thus there is an absolutely unstable flow. This indicates the supporting role of the inertial effects in introducing absolute instability.

To further explore the effect of inertia, the variation of $Ca_c$ with $1-H$ is shown in figure 9 for $Re=100$. The inertial effects curtail the convectively unstable region for a higher thickness of the ASL. The airway cross-sectional area becomes narrower for a diseased lung with a high ASL layer thickness. Now, to maintain the appropriate gas exchange balance, the same volume of gas must flow through a lesser cross-section. This can result in an increase in the air velocity associated with rapid breathing. This increases the inertial effects in the airways since $Re={\rho _2 V_c R}/{\mu _s}$ and increase in the capillary number since $Ca= {\mu _s V_c}/{\varSigma }$. From figure 9, if the capillary number increases faster than the Reynolds number and the ASL thickness is not appreciable, then the flow could become convectively unstable, thereby avoiding airway closure due to absolute instability. However, once the thickness of the ASL crosses a threshold value, the effect of the increased inertia comes into play, and the flow becomes absolutely unstable for an arbitrarily high value of $Ca$, thus leading to airway closure.

Figure 9.

The variation of the critical capillary number ($Ca_c$) with the dimensionless thickness of the ASL ($1-H$) for $Re=100$, $W=0$, $\mu _r=0.001$, $n=1$ and $\rho _r=0.001$. From the location of the turning point of the curve, the flow is absolutely unstable for arbitrarily high $Ca$ for $1-H > 0.054$ or $H<0.946$.

For a healthy lung, during exercise, an increased air velocity due to rapid breathing will increase $Re$ and $Ca$. However, this will decrease the airway closure possibility due to an increase in $Ca_c$ since the thickness $1-H \sim 0.001$ for which the inertial effects do not affect the absolute instability (see figure 9). To conclude, increasing air velocity could increase the risk of airway closure due to inertial effects in diseased lungs while having the opposite effect on healthy lungs.

7.3. Effect of shear thinning ($n$)

The ASL exhibits a strong shear thinning with power-law index $n \sim 0.1\unicode{x2013}0.5$ (Lai et al. Reference Lai, Wang, Wirtz and Hanes2009). In this section, we will analyse the effect of ASL shear thinning on the absolute instability parametric regions. To obtain a flavour of the modifications brought about by the shear thinning, consider a parameter set $Ca=10$, $H=0.95$, $\mu _r =0.001$, $Re=0$, $W=0$ and $n=1$ for which the cusp point is $\omega _0 = 0.000101166 - 0.0000127818 i$. Now, if we introduce shear thinning by varying $n=0.3$, the cusp point shifts to $\omega _0 = 0.000113344 + 0.0000272691 i$. Thus, the shear thinning increases the threshold $Ca$ at which the flow becomes convectively unstable. The effect of the power-law index on the absolute instability region is shown in figure 10. The shear thinning enables a higher $Ca$ range accessible for the absolute instability. Note that pulmonary diseases (e.g. cystic fibrosis) can exacerbate the shear-thinning nature of the ASL owing to the increased concentration of the mucins in the mucus. Thus, in addition to the increased viscosity and thickness, pulmonary diseases increase the probability of absolute instability induced airway closure through enhanced shear thinning.

Figure 10.

The variation of the critical capillary number $Ca_c$ with the dimensionless thickness of the ASL layer $1-H$ for $Re=0$, $\mu _r=0.001$ and $W=0$. The figure shows that ASL shear thinning enlarges the absolutely unstable region.

7.4. Effect of elasticity ($W$)

The mucus layer in the ASL exhibits strong viscoelastic properties (Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Hill et al. Reference Hill2014). Thus, in this section, we explore the effect of variation in the ASL elasticity (i.e. $W$) on the absolute instability. Table 3 shows the variation in the cusp point location due to variation in the Weissenberg number $W$. Even as $W$ changes by two orders, there is a negligible change in the $\omega _{i0}$ value, indicating a negligible effect of the ASL elasticity on the absolute instability in the higher-generation airways. Note that a pulmonary disease also increases ASL elasticity along with the viscosity and thickness (Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Hill et al. Reference Hill2014). An increasing ASL elasticity does not majorly affect the absolute instability, while from § 7.1, increasing ASL viscosity and thickness have a strong effect on the existence of absolute instability. Additionally, for $n<1$, the elastic mode exhibits a further decrease in the growth rate for the parameters relevant to the higher-generation airways, implying a decreasing role of the elastic mode in introducing absolute instability. Thus, pulmonary diseases primarily affect airway closure probability through the increase in the viscosity, shear thinning and thickness of the ASL.

Table 3.

The location of the cusp point ($\omega _0$) for select values of the Weissenberg number ($W$). The other parameters are $Ca=0.1$, $Re=0$, $\mu _r=0.001$, $n=1$ and $H=0.95$. The table demonstrates the weak effect of the variation in ASL elasticity on absolute instability.

8. Summary

In the present study, we analyse the temporal and spatio-temporal stability of flow in the higher-generation airways. Unlike turbulent airflow in the proximal or lower-generation airways, the airflow in the higher-generation airways is laminar. Thus, we treat this as a core–annular flow with air as the core fluid and ASL as the annular fluid. Owing to the shear-thinning viscoelastic properties of the ASL, we utilise the WM model for the annular fluid while the core fluid is treated as a Newtonian fluid. Additionally, we consider density and viscosity stratified flow due to the large difference in the densities and viscosities of the air and ASL.

The general linear temporal stability analysis predicts three unstable long-wave modes, in agreement with the previous studies. First, the Yih mode arises due to the viscosity stratification. Second, the elastic mode is due to the first normal stress difference across the air–ASL interface. Third, the capillary mode is due to the air–ASL interfacial tension. The temporal stability analysis reveals that, for the parameters relevant to the higher-generation airways, the three modes mentioned above will contribute to the destabilisation of the ASL. The shear thinning destabilises the capillary mode and stabilises the elastic mode.

The experiments of Guillot et al. (Reference Guillot, Colin, Utada and Ajdari2007, Reference Guillot, Colin and Ajdari2008) and Utada et al. (Reference Utada, Fernandez-Nieves, Stone and Weitz2007) on core–annular flows showed the existence of a continuous core-fluid jet, core-fluid drops and annular liquid slug formation depending on the flow parameters. A lubrication-theory-based model of Guillot et al. (Reference Guillot, Colin, Utada and Ajdari2007) concluded that the formation of the core-fluid drops and annular liquid slug is a consequence of an absolute instability. Thus, if the flow in the higher-generation airways exhibits absolute instability, then it could lead to the formation of air bubbles and ASL slugs. The presence of air bubbles or ASL slugs implies the presence of liquid bridges similar to those originating from the film collapse described in Heil et al. (Reference Heil, Hazel and Smith2008), or ASL plugs predicted by Romano et al. (Reference Romano, Fujioka, Muradoglu and Grotberg2019) that are responsible for airway closure. Thus, an absolutely unstable flow, if it exists in the higher-generation airways, may lead to airway closure. However, the existence of an absolutely unstable flow in these airways cannot be immediately concluded from the previous studies due to the extreme parametric range relevant to the higher-generation airways that remains unexplored in the previous analyses due to extreme density and viscosity stratification, the effect of inertial forces and the shear-thinning viscoelastic nature of the ASL. The present study is aimed at addressing this important missing gap.

In the spirit of the lubrication-theory-based model of Guillot et al. (Reference Guillot, Colin, Utada and Ajdari2007) and long-wave equation of Camassa et al. (Reference Camassa, Ogrosky and Olander2014) owing to the long-wave nature of the unstable modes, here, we derive an evolution equation for the long-wave perturbations and analyse the temporal and spatio-temporal stability using the resulting dispersion relation. To probe for absolute instability, we utilise a combination of the analytical calculations and the Kupfer et al. (Reference Kupfer, Bers and Ram1987) method.

The present analysis reveals the possibility of an absolutely unstable flow in the higher-generation airways in certain parametric regimes. This leads to the existence of the critical capillary number $Ca_c$ such that, for $Ca < Ca_c$, the flow will be absolutely unstable for vanishing Reynolds number $Re$. The parameter $Ca_c$ exhibits the scaling of $Ca_c \sim (1-H)^2/\mu _r$, where $1-H$ is the dimensionless thickness of the ASL, and $\mu _r$ is the ratio of the air viscosity to the ASL viscosity.

In agreement with the clinical observations of Burger & Macklem (Reference Burger and Macklem1968) and suggestions of Heil et al. (Reference Heil, Hazel and Smith2008), for a healthy lung, where the ASL thickness is $H \sim 0.999$, our analysis predicts absolute instability only at the end of expiration during a breathing cycle owing to the drop in $Ca$ value triggered by the low air velocity. In diseased lungs, the ASL thickness and viscosity can drastically increase. From the scaling $Ca_c \sim (1-H)^2/\mu _r$, one can immediately conclude that this will enhance the possibility of absolutely unstable flow and thus airway closure, in agreement with the experimental studies (Mall Reference Mall2016).

Due to the substantial air velocity in the earlier generations of the distal airways, we also analyse the effect of the Reynolds number on the absolute instability. An increasing $Re$ curtails the convectively unstable region in $(1-H) - Ca_c$ parametric space, thereby increasing the possibility of airway closure. A patient with a diseased lung can exhibit a higher $Re$ due to rapid breathing. This could further enhance the occurrence of airway closure. Due to the small thickness of the ASL in healthy lungs, an increasing air velocity (e.g. due to exercise) reduces the airway closure probability. The ASL also exhibits a strong shear thinning with the power-law index $n \sim 0.1\unicode{x2013}0.5$. The analysis shows that an increasing shear thinning leads to a bigger absolute instability region, thereby enhancing the probability of airway closure. The ASL elasticity has a negligible effect on the absolute instability.

The above conclusions must be considered in light of the following shortcomings of the present study. Due to breathing, an oscillatory flow exists in the airways, which could modify the occurrence of absolute instability. Additionally, repeated airway bifurcation can affect the convective to absolute transition. These factors have been neglected in the flow model considered here in order to simplify the complicated flow conditions and make the flow amenable to theoretical treatment. The present study could be extended by considering these modifications. Also, in the present analysis, the instability has not been tracked until the airway closure by considering the nonlinear development of the perturbations. Rather, the conclusions of the present study are based on the experimental outcomes of Guillot et al. (Reference Guillot, Colin, Utada and Ajdari2007). We have treated the ASL as a single layer by considering the average properties of the mucus and serous layers. One could extend the present study by considering a three-layer flow of air, mucus and serous layers. Careful experiments could validate the predictions of the present analysis.