1. Introduction
The human oral area is lined with saliva (thickness 
$\sim$1–10 
$\mathrm {\mu }$m) (Collins & Dawes Reference Collins and Dawes1987; Zussman, Yarin & Nagler Reference Zussman, Yarin and Nagler2007) while the airways, from the trachea through the bronchioles to the alveoli, are lined with mucus of varying thicknesses from 
$\sim$0.1 
$\mathrm {\mu }$m in the alveoli to 
$\sim$10 
$\mathrm {\mu }$m in the trachea (Widdicombe & Widdicombe Reference Widdicombe and Widdicombe1995; Widdicombe Reference Widdicombe1997; Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). The saliva lines the oral area, including the tongue, possessing a shear modulus of 
${\sim }10^{3}$ Pa (Cheng et al. Reference Cheng, Gandevia, Green, Sinkus and Bilston2011; Nakamori et al. Reference Nakamori2020) while the mucus lines the mucosa and submucosa layers of trachea, having a shear modulus of 
${\sim }10^{4}$ Pa (Wang, Mesquida & Lee Reference Wang, Mesquida and Lee2011) and is thus deformable. The present study aims to understand the role of the air, shear-thinning viscoelastic saliva and mucus, and deformable muscles lined by saliva or mucus in introducing elastohydrodynamic instabilities. These instabilities could play a role in airway closure or in the formation of fluid particles (i.e. droplets and aerosols) expelled from a person during breathing, talking, coughing and sneezing (Tang et al. Reference Tang, Nicolle, Pantelic, Koh, Wang and Amin2012, Reference Tang, Nicolle, Pantelic, Koh, Wang and Amin2013; Theriault et al. Reference Theriault, Eddy, Borgaonkar, Babar and Manos2018; Mittal, Ni & Seo Reference Mittal, Ni and Seo2020), which may lead to the spread of respiratory diseases.
1.1. Flows of non-Newtonian liquids
The saliva and mucus exhibit strong viscoelastic properties owing to their constituents (Hwang, Litt & Forsman Reference Hwang, Litt and Forsman1969; Sims & Horne Reference Sims and Horne1997; Stokes & Davies Reference Stokes and Davies2007; Zussman et al. Reference Zussman, Yarin and Nagler2007; Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Hamed & Fiegel Reference Hamed and Fiegel2014; Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). The viscoelasticity of the mucus is caused by the mucins and non-mucin proteins secreted by the mucous glands (Wu & Carlson Reference Wu and Carlson1991; Moriarty & Grotberg Reference Moriarty and Grotberg1999; Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). These characteristics help the mucus to protect the airways from abrasion of the particles and pathogens carried by the air inhaled during breathing or exhaled during coughing and sneezing (Stokes & Davies Reference Stokes and Davies2007; Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Hamed & Fiegel Reference Hamed and Fiegel2014). The thickness of the mucus can substantially increase or mucus properties could be drastically altered by respiratory diseases such as the common cold and cystic fibrosis (Shah et al. Reference Shah, Scott, Knight, Marriott, Ranasinha and Hodson1996). The elastic nature of the saliva and mucus could substantially affect the exhibited instabilities. Thus, it is essential to first review the instabilities exhibited by flows of viscoelastic liquids.
The study of the onset of transition from laminar to turbulent flow in viscoelastic liquids has received renewed interest in part due to the unambiguous demonstration of ‘early transition’ in the flow of polymer solutions through tubes at Reynolds numbers (
$Re \sim 800$) much lower than 
$Re \sim 2000$ at which the Newtonian transition is typically observed in experiments (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Chandra, Shankar & Das Reference Chandra, Shankar and Das2018; Choueiri, Lopez & Hof Reference Choueiri, Lopez and Hof2018). In addition, there have been several reports of instabilities in the flow of highly concentrated polymer solutions in tubes and channels, and the Reynolds number at which this transition is observed is again significantly lower than 
$Re \sim 2000$ (Bodiguel et al. Reference Bodiguel, Beaumont, Machado, Martinie, Kellay and Colin2015; Poole Reference Poole2016; Picaut et al. Reference Picaut, Ronsin, Caroli and Baumberger2017; Wen et al. Reference Wen, Poole, Willis and Dennis2017; Chandra et al. Reference Chandra, Shankar and Das2018). The linear stability of the plane Couette flow of an upper convected Maxwell (UCM) liquid was first studied by Gorodtsov & Leonov (Reference Gorodtsov and Leonov1967), who predicted the existence of two stable discrete modes in the creeping-flow limit. Renardy & Renardy (Reference Renardy and Renardy1986) studied the same flow, but for non-zero Reynolds number, and concluded that the flow is linearly stable for arbitrary 
$Re$.
Viscoelastic liquid flows sheared by the air can be encountered in coating processes (Kistler & Schweizer Reference Kistler and Schweizer1997). The air that flows past the coated film during the drying process could lead to an instability which may result in the formation of undesirable patterns. The formation of these patterns on the surface of a coating film degrades the quality of the final product and plays a crucial role in coating technology (Weinstein & Ruschak Reference Weinstein and Ruschak2004). Thus, the present study is also relevant to processes where air cooling is employed to dry the coating.
Besides viscoelasticity, saliva and mucus also exhibit shear thinning that depends on the concentration of the mucins and other constituent proteins (Lafforgue et al. Reference Lafforgue, Seyssiecq-Guarante, Poncet and Favier2017). As the concentration of the mucin in the mucus increases, the shear-thinning nature of the mucus increases. If the shear dependence of the viscosity is modelled using a power-law model, Lafforgue et al. (Reference Lafforgue, Seyssiecq-Guarante, Poncet and Favier2017) showed that the power-law index (
$n$) characterising the shear-thinning nature of the mucus is given by 
$n=1-0.346 C$ where 
$C$ (in wt %) is the concentration of the mucin. The power-law index for a healthy individual is 
$n \sim 0.1\text {--}0.5$ (Lai et al. Reference Lai, Wang, Wirtz and Hanes2009), thus it is a highly shear-thinning fluid. Nouar & Frigaard (Reference Nouar and Frigaard2009) studied the combined plane Couette–Poiseuille flow of a Carreau fluid in a rigid channel to understand the role of the wall velocity and shear-thinning behaviour on the stability. Their analysis predicted that the wall velocity and thus Couette component has a stabilising effect on plane Couette–Poiseuille flow and that increasing the shear-thinning behaviour stabilises the long-wavelength modes. Furthermore, they concluded that the plane Couette flow of an inelastic Carreau fluid in a rigid channel is linearly stable.
Additionally, at the air–mucus interface, a pulmonary surfactant secreted by the alveolar type II cells (Grotberg Reference Grotberg2001) is also present, which leads to the Marangoni stresses. The Marangoni stresses caused by the secreted surfactant play a major role in regulating the air–mucus interface tension. Thus, the absence of the surfactant at the air–mucus interface could lead to respiratory distress syndrome due to a high air–mucus interfacial tension, which in turn may lead to the inhibition of normal breathing (Grotberg Reference Grotberg2001). Blyth & Pozrikidis (Reference Blyth and Pozrikidis2004, Reference Blyth and Pozrikidis2006), Bassom, Blyth & Papageorgiou (Reference Bassom, Blyth and Papageorgiou2012) and Halpern & Grotberg (Reference Halpern and Grotberg1992) predicted a strong effect of the surfactant along the fluid–fluid interface on the stability of flows where the surfactant effect arises due to the Marangoni stresses caused by the inhomogeneity of the surface tension due to its dependence on the concentration of the surfactant.
1.2. Flows past a deformable-solid layer
The surfaces lined by the saliva and mucus are deformable owing to their shear modulus of 
$O(10^{3}\text {--}10^{4})$ Pa (Cheng et al. Reference Cheng, Gandevia, Green, Sinkus and Bilston2011; Wang et al. Reference Wang, Mesquida and Lee2011; Nakamori et al. Reference Nakamori2020). The deformability of the wall plays a major role in fluid flows past a deformable-solid layer. Thus, next, we review the instabilities exhibited by the flows past a deformable-solid layer.
The stability of plane Poiseuille flow through a channel with deformable walls was studied by Hains & Price (Reference Hains and Price1962), Gajjar & Sibanda (Reference Gajjar and Sibanda1996) and Davies & Carpenter (Reference Davies and Carpenter1997). They used the plate-and-spring model for the deformable walls, assuming that the solid is made up of plates connected by Hookean springs. However, the continuum models, such as linear viscoelastic, neo-Hookean models, describe the dynamics of the solid better than the plate-and-spring model (Malvern Reference Malvern1969; Holzapfel Reference Holzapfel2000). Thus, Kumaran, Fredrickson & Pincus (Reference Kumaran, Fredrickson and Pincus1994) used the linear viscoelastic model for the solid to study the linear stability of the plane Couette flow past a deformable solid in the creeping-flow limit (
$Re=0$) and predicted the presence of the ‘viscous instability.’ The viscous instability arises as a result of the shear work done by the flowing fluid on the deformable solid (Kumaran et al. Reference Kumaran, Fredrickson and Pincus1994). The experiments of Kumaran & Muralikrishnan (Reference Kumaran and Muralikrishnan2000) confirmed the existence of the viscous instability predicted by Kumaran et al. (Reference Kumaran, Fredrickson and Pincus1994).
To account for the finite base-state deformations of the solid, Gkanis & Kumar (Reference Gkanis and Kumar2003) used the frame-invariant neo-Hookean model for the solid and predicted the ‘short-wave instability’ in addition to the viscous instability in plane Couette flow past a deformable solid. Neelamegam, Giribabu & Shankar (Reference Neelamegam, Giribabu and Shankar2014) studied the stability of the plane Couette flow past a two-layered gel and showed an excellent agreement between theory (using the neo-Hookean model) and experiments, thereby proving the usefulness of the neo-Hookean model to describe the dynamics of the deformable solid. Gkanis & Kumar (Reference Gkanis and Kumar2005), Gaurav & Shankar (Reference Gaurav and Shankar2010) and Patne & Shankar (Reference Patne and Shankar2017) subsequently studied the plane Poiseuille flow through a channel lined with neo-Hookean solid and predicted the existence of the short-wave instability in the creeping-flow limit. The analysis of Gaurav & Shankar (Reference Gaurav and Shankar2010) also predicted a new class of unstable modes at low 
$Re$ and two classes of modes, viz. wall modes, and inviscid modes at high 
$Re$. Experiments on plane Poiseuille flow past a deformable solid by Verma & Kumaran (Reference Verma and Kumaran2013) found good agreement with theory provided that the solid is modelled using the neo-Hookean model and the change in the shape of the deformable wall due to a decrease in pressure along the length of the channel and its effect on the base-state velocity profile are considered.
The stability of the plane Couette flow of viscoelastic liquids past a deformable-solid layer has been studied by Shankar & Kumar (Reference Shankar and Kumar2004) and Chokshi & Kumaran (Reference Chokshi and Kumaran2008a,Reference Chokshi and Kumaranb). Shankar & Kumar (Reference Shankar and Kumar2004) predicted the destabilisation of the stable modes predicted by Gorodtsov & Leonov (Reference Gorodtsov and Leonov1967) in the creeping-flow limit due to the deformability of the solid. A subsequent study by Joshi & Shankar (Reference Joshi and Shankar2019) showed the existence of a new mode of instability that results from the elastic nature of the viscoelastic liquid and solid. Thus, Shankar & Kumar (Reference Shankar and Kumar2004) and Joshi & Shankar (Reference Joshi and Shankar2019) showed the importance of considering the deformability of the solid in viscoelastic liquid flows past a deformable-solid layer.
The effect of the shear thinning on the stability of the plane Couette flow past a deformable solid in the creeping-flow limit was first studied by Roberts & Kumar (Reference Roberts and Kumar2006). Their analysis predicted a non-monotonic effect of the variation in the power-law index on the viscous and short-wave modes of instability. Their work was extended to arbitrary Reynolds number by Giribabu & Shankar (Reference Giribabu and Shankar2017) and Tanmay, Patne & Shankar (Reference Tanmay, Patne and Shankar2018).
1.3. Liquid layer sheared by the air
In the oral area and airways, the saliva and mucus layers are sheared by the air, which becomes the driving force to set in the instabilities. The stability of a Newtonian liquid layer sheared by the air was first studied by Miles (Reference Miles1960). He considered the shear stress exerted by the flowing air and subsequent shear flow generated in the liquid. His study predicted critical Reynolds number 
$Re_c \sim 203$ for the existence of the instability. Later, Smith & Davis (Reference Smith and Davis1982) pointed out a ‘missing term’ in the normal stress continuity boundary condition of Miles (Reference Miles1960). Upon inclusion of the missing term, Smith & Davis (Reference Smith and Davis1982) predicted 
$Re_c \sim 34$ in the absence of the air–liquid interfacial tension. Furthermore, they observed that the air–liquid interfacial tension leads to a strong stabilisation.
1.4. Present study
As explained above, mucus and saliva exhibit a strong shear-thinning viscoelastic behaviour. Thus, a fluid model which considers both of their properties is necessary to successfully capture the dynamics of these fluids sheared by the air. Furthermore, in the literature, the power-law index and relaxation time values for mucus are available (see Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Lafforgue et al. Reference Lafforgue, Seyssiecq-Guarante, Poncet and Favier2017), which indicate that a model must have a power-law variation of the viscosity. These requirements are fulfilled by the phenomenological White & Metzner (Reference White and Metzner1963) (WM) model. The WM model introduces shear thinning or thickening and strain softening or hardening by assuming dependences of the viscosity, shear modulus and relaxation time on the second invariant of the strain-rate tensor. Thus, data reported in the literature for the rheological parameters of the WM model such as viscosity, relaxation time and shear modulus, can be readily fit into the power-law model, thereby making it easy to relate the experimental observations to theoretical predictions. The linear stability analysis of pressure-driven channel flow of a WM fluid was first studied by Wilson & Rallison (Reference Wilson and Rallison1999), Wilson & Loridan (Reference Wilson and Loridan2015) and Castillo & Wilson (Reference Castillo and Wilson2017). Their analysis predicted the existence of a new shear-thinning elastic instability.
Halpern & Grotberg (Reference Halpern and Grotberg1992) analysed the stability of a Newtonian liquid lining the inner surface of an elastic tube. The role of the pulmonary surfactant on the fluid-elastic instabilities in a Newtonian liquid-lined elastic tube was studied by Halpern & Grotberg (Reference Halpern and Grotberg1993), where they predicted a delay in the airway closure due to the surfactant. It must be noted that Halpern & Grotberg (Reference Halpern and Grotberg1992, Reference Halpern and Grotberg1993) utilised the plate-and-spring model to describe the dynamics of the deformable wall. Moriarty & Grotberg (Reference Moriarty and Grotberg1999), motivated by the understanding of the instabilities related to mucus clearance, studied the stability of a bilayer (formed by the elastic sheet of the mucin and serous layers) sheared by the air. They modelled the elastic sheet formed by the mucin layer using the Kelvin–Voigt model while the serous layer lying below the mucus layer was modelled as a viscous liquid. Their analysis predicted a minimum air speed necessary for the existence of the instabilities in the airways. Heil & White (Reference Heil and White2002), by using nonlinear shell theory for the deformable wall, studied the airway closure due to elastic instabilities for a Newtonian fluid. A thin-film model for the airway closure due to the surface-tension-driven instabilities was later developed by White & Heil (Reference White and Heil2005). Halpern, Fujioka & Grotberg (Reference Halpern, Fujioka and Grotberg2010) studied the stability of a viscoelastic film coating a rigid tube and flowing under the action of gravity. Their analysis predicted a destabilising influence of viscoelasticity. Zhou et al. (Reference Zhou, Peng, Zhang and Zhuge2014) extended the study of Halpern et al. (Reference Halpern, Fujioka and Grotberg2010) to include the effect of imposed shear at the free surface and the presence of an insoluble surfactant. Their analysis predicted a destabilising effect of the surfactant.
The study of Moriarty & Grotberg (Reference Moriarty and Grotberg1999) suffers from some important drawbacks as follows. The differentiation of the mucus into the elastic sheet formed by the mucins at the air–liquid interface and serous layers assumes a sharp interface between the mucin surface layer and serous layer which is unphysical since the mucins and other solutes are in fact dissolved in the serous solution (Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). Furthermore, Halpern & Grotberg (Reference Halpern and Grotberg1992, Reference Halpern and Grotberg1993), Moriarty & Grotberg (Reference Moriarty and Grotberg1999), Heil & White (Reference Heil and White2002), White & Heil (Reference White and Heil2005) and Halpern et al. (Reference Halpern, Fujioka and Grotberg2010) neglect the shear-thinning viscoelastic nature of the mucus, saliva and presence of the surfactant at the air–mucus interface. A consistent physical model requires treating the mucin and serous layers as a single layer with viscoelastic and shear-thinning properties possessing a finite air–liquid interfacial tension and surfactant at the free surface. Furthermore, Moriarty & Grotberg (Reference Moriarty and Grotberg1999), Halpern et al. (Reference Halpern, Fujioka and Grotberg2010) and Zhou et al. (Reference Zhou, Peng, Zhang and Zhuge2014) neglected the dynamics of the deformable mucosa and submucosa layers lined by the mucus. However, as shown in the present study, the dynamics of these deformable layers can introduce a whole new mode of instability, thereby showing their importance.
Studies that take into consideration the deformable nature of the walls of the airways are due to Halpern & Grotberg (Reference Halpern and Grotberg1992, Reference Halpern and Grotberg1993), Heil & White (Reference Heil and White2002) and White & Heil (Reference White and Heil2005). They employ a plate-and-spring model or shell theory to describe the dynamics of the deformable wall. However, continuum models, such as linear viscoelastic, neo-Hookean models, offer a better description of the deformable wall dynamics instead of lumped parameter models, such as the plate-and-spring model or shell theory (Malvern Reference Malvern1969; Holzapfel Reference Holzapfel2000). Furthermore, the plate-and-spring model allows only vertical oscillations of the deformable wall. However, as noted by Thaokar & Kumaran (Reference Thaokar and Kumaran2002), tangential motion of deformable media is necessary for the existence of a class of instabilities exhibited by the flows past deformable surfaces. The airflow in the airways and oral area is turbulent (Moriarty & Grotberg Reference Moriarty and Grotberg1999), however, Halpern & Grotberg (Reference Halpern and Grotberg1992, Reference Halpern and Grotberg1993), Moriarty & Grotberg (Reference Moriarty and Grotberg1999), Heil & White (Reference Heil and White2002), White & Heil (Reference White and Heil2005) and Halpern et al. (Reference Halpern, Fujioka and Grotberg2010) assume a laminar flow of air.
The three-layer system consisting of the air, saliva or mucus, and underlying muscles requires a consistent model that will take into account the dynamics and properties of the three layers, as shown in figure 1. The present study removes the shortcomings of the previous studies by using the WM and neo-Hookean models to describe the dynamics of the mucus and saliva and underlying muscles, respectively, while the turbulent airflow is accounted for by the shear stress exerted at the air–liquid interface. Additionally, the effect of the surfactant and its impact on the stability analysis is also considered.
Figure 1. Schematic of the flow geometry in dimensionless coordinates. The liquid flows as a consequence of the driving force imparted by the air blowing past the liquid layer at 
$y=1$. The flowing liquid exerts stress on the deformable solid to cause deformation, which leads to the coupling. The boundary at 
$y=-H$ is assumed to be perfectly bonded to a rigid and impermeable surface, a reasonable assumption for the trachea due to the presence of the relatively rigid outer musculocartilaginous and adventitia layers (Wang et al. Reference Wang, Mesquida and Lee2011).
The airflow through airways could be better described in a cylindrical flow geometry due to the tubular structure of the airways. However, we assume a planar flow in the present study for the following reasons. The airflow through airways is similar to the sliding Couette flow with a deformable wall studied by Patne & Shankar (Reference Patne and Shankar2020). The sliding Couette flow with a deformable wall consists of an outer stationary and deformable cylinder and a steadily moving inner cylinder. The stability analysis of Patne & Shankar (Reference Patne and Shankar2020) showed that if the parameter 
$\xi =r_i/r_o$ approaches unity, then the predicted results approach the results for the plane Couette flow, i.e. the effect of the curvature can be neglected. Here, 
$r_i$ is the radius of the inner cylinder and 
$r_o$ is the sum of the radius of the inner cylinder and fluid layer thickness. Their analysis also predicted that the axisymmetric mode is most unstable for the sliding Couette flow. The flow geometry in the present study is essentially the same as the sliding Couette flow except that the inner cylinder should be replaced by the flowing air. The thickness of the mucus is considerably smaller than the air domain, implying 
$\xi \rightarrow 1$, thereby justifying the planar flow assumption. It must be noted that the planar flow assumption neglects the possible stabilising/destabilising mechanism caused by the surface tension at the air–liquid interface in a core–annular geometry which may quantitatively alter the growth rate of the instabilities predicted in the present study.
The rest of the paper is arranged as follows. The perturbation governing equations and boundary conditions are derived in § 2. The numerical methodology employed here to carry out the linear stability analysis is presented in § 3. The results obtained from the stability analysis and the mechanism of the instabilities predicted here are presented in § 4. Finally, the major conclusions obtained in the present study are presented in § 5.
2. Governing equations
Consider an incompressible viscoelastic liquid of a constant shear modulus 
$G_f$ is flowing past a deformable layer having shear modulus 
$G_s$, as shown in figure 1. For the processes considered here, the same density 
$\rho$ of the liquid and deformable solid is assumed. The driving force for the liquid is provided by the air flowing past the liquid at 
$y=1$. The resulting flow in the liquid then induces a deformation field in the deformable solid. The deformation leads to a coupling between the flowing liquid and deformable solid which may result in the elastohydrodynamic instabilities (Kumaran Reference Kumaran1995; Kumaran & Muralikrishnan Reference Kumaran and Muralikrishnan2000). The Cartesian reference frame chosen here is such that the liquid flows in the 
$x-z$ plane while the 
$y$-axis is normal to the flowing liquid.
The dimensional fluid viscosity (
$\mu ^{*}_p$) and relaxation time (
$\lambda ^{*}_p$) are assumed to be functions of the second invariant (
$\dot \gamma ^{*}$) of the strain-rate tensor (
$\dot {\boldsymbol {\gamma }}^{\boldsymbol {*}}$; refer to (2.5f)). For simplicity, we use the power-law model for the polymer viscosity (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1977; Wilson & Rallison Reference Wilson and Rallison1999), which relates the viscosity to the shear rate as
$$\begin{gather} \mu^{*}_p(\dot \gamma^{*}) = K \dot \gamma^{*(n-1)}, \end{gather}$$
$$\begin{gather}\dot \gamma^{*}=\sqrt{\tfrac{1}{2} \dot \gamma_{ij}^{*} \dot \gamma_{ji}^{*}}. \end{gather}$$
Here, 
$K$ denotes the consistency parameter in the power-law model and 
$n$ is the power-law index. The superscript ‘
$*$’ in the above equation and in the subsequent discussion denotes dimensional quantities while subscript ‘
$p$’ denotes a polymer quantity. By using the above definition of viscosity, the relaxation time (
$\lambda ^{*}_p$) for the WM fluid becomes
\begin{equation} \lambda^{*}_p(\dot \gamma^{*}) =\frac{\mu^{*}_p(\dot \gamma^{*})}{G_f}= \frac{K}{G_f} \dot \gamma^{*(n-1)} . \end{equation}
At the air–liquid interface the surfactant is also present. Following Blyth & Pozrikidis (Reference Blyth and Pozrikidis2004, Reference Blyth and Pozrikidis2006) and Bassom et al. (Reference Bassom, Blyth and Papageorgiou2012), we assume a linear dependence of the surface tension (
$\varSigma ^{*}$) on the surfactant concentration (
$\zeta ^{*}$) such that 
$\varSigma ^{*} = \varSigma _0^{*} - \gamma ^{*} (\zeta ^{*} - \zeta _0^{*})$, where 
$\zeta _0^{*}$ is the reference concentration of the surfactant and 
$\gamma ^{*} = -{\textrm {d} \varSigma ^{*}}/{\textrm {d} \zeta ^{*}}$.
Before proceeding further, we present the non-dimensionalisation scheme for the problem under consideration. The power-law model does not have a natural viscosity scale, thus we define the viscosity scale as 
$\mu ^{*}_{sc}=K ( {V_m}/{R} )^{(n-1)}$, where 
$V_m$ and 
$R$ are respectively the maximum velocity and height of the liquid layer. The maximum velocity of the liquid is 
$V_m=\tau _a R/\mu$, where 
$\tau _a$ is the shear stress exerted by the air at the air–liquid interface. We non-dimensionalise lengths, velocities, viscosities and pressures and stresses, respectively, by 
$R$, 
$V_m$, 
$\mu ^{*}_t=\mu ^{*}_s+\mu ^{*}_{sc}$ and 
$\mu ^{*}_t V_m/R$, where 
$\mu ^{*}_t$ and 
$\mu ^{*}_s$ are the total (solution) and solvent viscosities. To relate the perturbed state equations with the Oldroyd-B equations for a non-shear-thinning viscoelastic fluid, the viscosity of the polymer is non-dimensionalised by 
$\mu ^{*}_{sc}$. Hence, the dimensionless viscosity for the power-law model (2.1) becomes
\begin{equation} \mu_p(\dot \gamma)= \dot \gamma^{(n-1)}. \end{equation}
The liquid layer extends from 
$y=0$ to 
$y=1$ while the deformable layer extends from 
$y=-H$ to 
$y=0$. Thus, liquid–deformable-solid form an interface at 
$y=0$.
Let the velocity field in the liquid be 
$\boldsymbol {v}=(v_x,v_y,v_z)$. The dimensionless continuity and Cauchy momentum equations are
$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}=0, \end{gather}$$
$$\begin{gather}Re \left(\frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{v} \right) ={-}\boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau }, \end{gather}$$
$$\begin{gather}\boldsymbol{\tau }=\boldsymbol{\tau}^{\boldsymbol{s}}+\boldsymbol{\tau}^{\boldsymbol{p}}, \end{gather}$$
where, 
$\boldsymbol {v}=(v_1,v_2, v_3)$ is the velocity field, 
$v_i$ is the component of the velocity in the 
$i$th direction, the Reynolds number 
$Re=\rho V_m R/\mu ^{*}_t$, 
$p$ is the pressure, 
$\boldsymbol {\tau }$, 
$\boldsymbol {\tau ^{s} }$, 
$\boldsymbol {\tau ^{p} }$ are respectively the total, solvent and polymer stress tensors and 
$\nabla$ is the gradient operator. To relate the stress tensor with the shear-rate tensor, we use following constitutive relations:
$$\begin{gather} \boldsymbol{\tau^{s} }=\beta \dot{\boldsymbol{\gamma}}, \end{gather}$$
$$\begin{gather}\boldsymbol{\tau}^{\boldsymbol{p}}+W \lambda(\dot \gamma) \left( \frac{\partial \boldsymbol{\tau}^{\boldsymbol{p}}}{\partial t} + \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\tau}^{\boldsymbol{p}}- (\boldsymbol{\nabla}\boldsymbol{v})^{T} \boldsymbol{\cdot}\boldsymbol{\tau}^{\boldsymbol{p}}-\boldsymbol{\tau}^{\boldsymbol{p}} \boldsymbol{\cdot} (\boldsymbol{\nabla}\boldsymbol{v}) \right)=(1-\beta) \mu_p(\dot \gamma) \dot{\boldsymbol{\gamma}}, \end{gather}$$
$$\begin{gather}\dot{\boldsymbol{\gamma}}=(\boldsymbol{\nabla}\boldsymbol{v})+(\boldsymbol{\nabla}\boldsymbol{v})^{T}. \end{gather}$$
Here, 
$\beta =\mu ^{*}_s/\mu ^{*}_t$. The shear-rate-dependent dimensionless relaxation time in the above equation is 
$\lambda (\dot \gamma )={\mu _p(\dot \gamma )}/{G}$, which is non-dimensionalised by using scale 
$\lambda ^{*}_{sc}=\mu ^{*}_{sc}/G$. Also, 
$W=\lambda ^{*}_{sc} V_m/ R$ is the Weissenberg number. The constitutive relation for 
$\boldsymbol {\tau }^{\boldsymbol {p}}$ is given by the WM model, which differs from the Oldroyd-B model due to the dependence of the polymer viscosity on the shear rate.
In between the mucus and the deformable muscle layers, another liquid layer of much less viscous fluid, the periciliary layer (PCL), is also present. As discussed in § 4, the existence of the PCL will not affect the qualitative picture of the instabilities predicted here.
For the solid, consider a representative particle with position vector 
$\boldsymbol {X}=(X_{1},X_{2},X_{3})$ in the undeformed solid. Assuming the solid to be incompressible and letting liquid flow past it, the representative particle will now assume position vector, 
$\boldsymbol {x}=(x_{1},x_{2},x_{3})$. The current and undeformed position vectors are related by,
\begin{equation} \boldsymbol{x}(\boldsymbol{X})=\boldsymbol{X}+\boldsymbol{u}(\boldsymbol{X}), \end{equation}
where 
$\boldsymbol {u}(\boldsymbol {X})$ is the Lagrangian displacement in the solid. Hence, the deformation gradient is, 
$\boldsymbol {F}={\partial \boldsymbol {x}}/{\partial \boldsymbol {X}}$. The incompressibility condition is given by
\begin{equation} \det{(\boldsymbol{F})}=1. \end{equation}
The dimensionless Cauchy stress (
$\boldsymbol {\sigma }$) for a purely elastic neo-Hookean solid is (Macosko Reference Macosko1994; Holzapfel Reference Holzapfel2000; Patne, Giribabu & Shankar Reference Patne, Giribabu and Shankar2017),
\begin{equation} \boldsymbol{\sigma}={-}p_{g} \boldsymbol{I}+\frac{1}{\varGamma} \boldsymbol{F} \boldsymbol{\cdot} \boldsymbol{F}^{T}, \end{equation}
where 
$p_{g}$ is the pressure field in the solid and 
$\varGamma =\mu V_m/(G R)$ is the dimensionless maximum velocity of the base flow. The parameter 
$\varGamma$ is also measure of the deformation in the deformable solid. The non-dimensionalised momentum balance equation (Malvern Reference Malvern1969; Holzapfel Reference Holzapfel2000) is
$$\begin{gather} Re \frac{\partial^{2} \boldsymbol{u}}{\partial t^{2}}=\nabla_{\bar{\boldsymbol{x}}} \boldsymbol{\cdot} \boldsymbol{P}, \end{gather}$$
$$\begin{gather}\boldsymbol{P}=\boldsymbol{F}^{{-}1} \boldsymbol{\cdot} \boldsymbol{\sigma}, \end{gather}$$
where 
$\boldsymbol {P}$ is the first Piola–Kirchoff stress tensor. At the rigid-solid–deformable-solid interface (
$\bar y=-H$) there will be no perpendicular or tangential motion as it is assumed that the neo-Hookean solid is perfectly bonded to the rigid solid. Hence, the boundary condition at this interface becomes 
$\boldsymbol {u}=0$. At the liquid–solid interface, continuity conditions are imposed for the velocity and stresses. Similarly, at the air–liquid interface, the boundary conditions are the kinematic boundary condition and continuity of the tangential and normal stresses. The presence of the pulmonary surfactant at the air–liquid interface does not affect the base state.
Here, we have considered an axisymmetric deformation of the airways which allows us to consider the assumption of a planar flow geometry due to the small thickness of the saliva and mucus. However, the airways may undergo non-axisymmetric deformation which is not considered in the present study. Such a deformation in a cylindrical flow geometry with a neo-Hookean model for the mucosa and submucosa layers leads to a base-state deformation which is a function of the angular coordinate 
$\theta$. The stability analysis of such a flow is mathematically a cumbersome exercise and will necessitate global stability analysis. Thus, for the sake of simplicity in the present study, we have not considered the non-axisymmetric nature of airway deformation.
The trachea wall is lined with the cells containing small hair-like projections called cilia, with tips reaching the bottom of the mucus layer (Bustamante-Marin & Ostrowski Reference Bustamante-Marin and Ostrowski2017). The cilia layer helps in the entrapment and ejection of the pathogens and food particles that have strayed into the trachea by the mucociliary clearance mechanism. The present study assumes the cilia layer is also a part of the deformable-solid layer and neglects phenomena involving the active role of the cilia such as the mucociliary clearance mechanism owing to the related mathematical intricacies.
2.1. Base state
The fully developed, steady-state base-state velocity profile for the sheared flow is (Miles Reference Miles1960; Smith & Davis Reference Smith and Davis1982)
\begin{equation} \bar v_x=y; \quad \bar v_y=0; \quad \bar v_z=0; \quad \bar p=0. \end{equation}In the above equation and henceforth, an overbar signifies a base-state quantity. By using the above velocity profiles and constitutive equations, the base-state stresses are given by
\begin{equation} \bar \tau^{p}_{xy}=(1-\beta)\bar \mu_p(\bar{\dot{\gamma}}) \frac{\textrm{d} \bar v_x}{\textrm{d} y}; \quad \bar \tau^{p}_{xx}=2W\bar \lambda(\bar{\dot{\gamma}}) \bar \tau^{p}_{xy} \frac{\textrm{d} \bar v_x}{\textrm{d} y}; \quad \bar \tau^{p}_{yy}=0. \end{equation}
Here, 
$\bar \mu _p(\bar {\dot {\gamma }})$ and 
$\bar \lambda (\bar {\dot {\gamma }})$ are the dimensionless shear-rate-dependent viscosity and relaxation time of the polymer in the base state.
For the deformable layer, in the base state, the motion is, 
$\bar {\boldsymbol {x}}(\boldsymbol {X})=\boldsymbol {X}+\bar {\boldsymbol {u}}(\boldsymbol {X})$, which leads to the base-state deformation gradient 
$\bar {\boldsymbol {F}}={\partial \bar {\boldsymbol {x}}}/{\partial \boldsymbol {X}}$. Substitution of 
$\bar {\boldsymbol {F}}$ in (2.7a)–(2.7d), and following the procedure outlined by Patne et al. (Reference Patne, Giribabu and Shankar2017), the steady base-state deformation in the solid in terms of the pre-stressed state coordinates is
\begin{equation} \bar u_{x}= \varGamma (\bar y+H); \quad \bar p_{g} =0. \end{equation}2.2. Linearised perturbation equations
Next, we impose two-dimensional infinitesimally small perturbations on the base state (2.8). It must be noted that Squire's theorem, in general, is not applicable for the present three-layer system. However, for a Newtonian liquid layer sheared by the air past a rigid solid, Squire's theorem is applicable (Smith & Davis Reference Smith and Davis1982). For the planar flows past a deformable solid from the analysis of Patne & Shankar (Reference Patne and Shankar2019), two-dimensional perturbations are more unstable than the corresponding three-dimensional perturbations. For planar flows of an Oldroyd-B fluid past a rigid surface, Squire's theorem is applicable (Bistagnino et al. Reference Bistagnino, Boffetta, Celani, Mazzino, Puliafito and Vergassola2007), i.e. two-dimensional disturbances are more unstable than the corresponding three-dimensional disturbances. A similar demonstration is not possible for the WM fluid sheared by air, but we nonetheless assume, in the interests of simplicity, that two-dimensional perturbations to be more unstable than three-dimensional ones, and we restrict our attention to two-dimensional perturbations. For the perturbed state, the motion is described by 
$\boldsymbol {x}(\bar {\boldsymbol {x}})=\bar {\boldsymbol {x}}+\boldsymbol {u}'(\bar {\boldsymbol {x}},t)$. Here, 
$\boldsymbol {u}'(\bar {\boldsymbol {x}},t)$ is the Lagrangian displacement of the particle from base state and prime indicates that it is a perturbation quantity (henceforth, perturbation quantities will be indicated by a prime). The deformation gradient becomes 
$\boldsymbol {F}={\partial \boldsymbol {x}}/{\partial \boldsymbol {X}}$. The incompressibility condition for perturbed state after using base-state incompressibility condition becomes 
$\det {(\boldsymbol {F'})}=1$. Further, we assume two-dimensional disturbances with normal modes of the form
\begin{equation} f'(\boldsymbol{x},t)= \tilde{f}(y) \, \textrm{e}^{\textrm{i}(kx-\omega t)}, \end{equation}
where 
$f'(\boldsymbol {x},t)$ is the perturbation to a quantity and 
$\tilde {f}(y)$ is the corresponding eigenfunction. Here, 
$k$ is the (real-valued) wavenumber, while 
$\omega =\omega _r+\textrm {i} \omega _i$ is the complex frequency of the disturbances. The flow is linearly unstable if at least one eigenvalue satisfies the condition 
$\omega _i>0$.
Substituting (2.9) in the linearised governing equations, we obtain the following linearised continuity and momentum equations
$$\begin{gather} \textrm{i}k \tilde{v}_x+D\tilde{v}_y=0, \end{gather}$$
$$\begin{gather}Re[\textrm{i}(k\bar v_x-\omega)\tilde{v}_x+ D\bar v_x \tilde{v}_y]={-}\textrm{i}k\tilde p+\beta(D^{2}-k^{2})\tilde{v}_x+\textrm{i}k\tilde{\tau}^{p}_{xx}+D\tilde{\tau}^{p}_{xy}, \end{gather}$$
$$\begin{gather}\textrm{i}\,Re(k\bar v_x-\omega)\tilde{v}_y={-}D\tilde p+\beta(D^{2}-k^{2})\tilde{v}_y+\textrm{i}k\tilde{\tau}^{p}_{xy}+D\tilde{\tau}^{p}_{yy}, \end{gather}$$
where 
$D=\textrm {d}/{\textrm {d} y}$. The constitutive equation (2.5e) gives
$$\begin{gather} \hspace{-5pc}W\bar
\lambda(\bar{\dot{\gamma}}) [\textrm{i}(k\bar v_x-\omega)
\tilde{\tau}^{p}_{xx}+\tilde{v}_y D\bar
\tau^{p}_{xx}-2\textrm{i}k \bar \tau^{p}_{xx} \tilde{v}_x-2
D\bar v_x \tilde{\tau}^{p}_{xy}-2 \bar \tau^{p}_{xy} D
\tilde{v}_x]\nonumber\\ -\underline{2 W \bar
\tau^{p}_{xy} D \bar v_x \left(\frac{\textrm{d}
\lambda}{\textrm{d} \dot \gamma}\right)_{\dot
\gamma=\bar{\dot{\gamma}}}
[\textrm{i}k\tilde{v}_y+D\tilde{v}_x]}+\tilde{\tau}^{p}_{xx}=2\textrm{i}k(1-\beta)\bar
\mu_p(\bar{\dot{\gamma}})\tilde{v}_x,
\\ \hspace{-3pc}W\bar \lambda(\bar{\dot{\gamma}})
[\textrm{i}(k\bar v_x-\omega)
\tilde{\tau}^{p}_{xy}+\tilde{v}_y D\bar \tau^{p}_{xy}
-\textrm{i}k \bar \tau^{p}_{xy} \tilde{v}_x-D\bar v_x
\tilde{\tau}^{p}_{yy}-\textrm{i}k \bar \tau^{p}_{xx}
\tilde{v}_y-\bar \tau^{p}_{xy} D\tilde{v}_y]\nonumber\\
\hspace{-3pc}+\tilde{\tau}^{p}_{xy}=(1-\beta) \underline{
\left[\bar \mu_p(\bar{\dot{\gamma}})
+\left(\frac{\textrm{d} \mu_p}{\textrm{d} \dot
\gamma}\right)_{\dot \gamma=\bar{\dot{\gamma}}}
\bar{\dot{\gamma}} \right]
(\textrm{i}k\tilde{v}_y+D\tilde{v}_x)},\end{gather}$$
$$\begin{gather}W\bar \lambda(\bar{\dot{\gamma}}) [\textrm{i}(k\bar v_x-\omega) \tilde{\tau}^{p}_{yy}-2\textrm{i}k \bar \tau^{p}_{xy} \tilde{v}_y]+\tilde{\tau}^{p}_{yy}=2(1-\beta)\bar \mu_p(\bar{\dot{\gamma}})D\tilde{v}_y, \end{gather}$$
where 
$({\textrm {d} \lambda }/{\textrm {d} \dot \gamma })_{\dot \gamma =\bar {\dot {\gamma }}}$ and 
$({\textrm {d} \mu _p}/{\textrm {d} \dot \gamma })_{\dot \gamma =\bar {\dot {\gamma }}}$ are the derivatives of the polymer viscosity and relaxation time with respect to 
$\dot \gamma$ evaluated at the base state. The underlined terms in (2.10d)–(2.10e) are a consequence of the WM model. Additionally, the WM model modifies all the terms of Oldroyd-B equations due to multiplication by 
$\bar \lambda (\bar {\dot {\gamma }})$ for the elastic terms and 
$\bar \mu _p(\bar {\dot {\gamma }})$ for the viscous terms. The explicit expressions for 
$\bar \mu _p, \bar \lambda (\bar {\dot {\gamma }}),({\textrm {d} \mu _p}/{\textrm {d} \dot \gamma })_{\dot \gamma =\bar {\dot {\gamma }}}$ and 
$({\textrm {d} \lambda }/{\textrm {d} \dot \gamma })_{\dot \gamma =\bar {\dot {\gamma }}}$ for a WM fluid are
\begin{equation} \bar \lambda(\bar{\dot{\gamma}})=\bar \mu_p(\bar{\dot{\gamma}})=\bar{\dot{\gamma}}^{(n-1)}; \quad \left(\frac{\textrm{d} \lambda}{\textrm{d} \dot \gamma}\right)_{\dot \gamma=\bar{\dot{\gamma}}}=\left(\frac{\textrm{d} \mu_p}{\textrm{d} \dot \gamma}\right)_{\dot \gamma=\bar{\dot{\gamma}}}=(n-1)\bar{\dot{\gamma}}^{(n-2)}. \end{equation}Similarly for the solid, the incompressibility condition is
\begin{equation} \textrm{i}k \tilde{u}_{x}+D \tilde{u}_{y}=0, \end{equation}and momentum balance equations are
$$\begin{gather} -Re \varGamma \omega^{2} \tilde{u}_x={-}\textrm{i}k \varGamma \tilde{p}_{g}+(D^{2}+2\textrm{i}k \varGamma D-k^{2} (1+\varGamma^{2})) \tilde{u}_{x}, \end{gather}$$
$$\begin{gather}-Re \varGamma \omega^{2} \tilde{u}_y ={-}\varGamma D \tilde{p}_{g}+(D^{2}+2\textrm{i}k \varGamma D-k^{2} (1+\varGamma^{2})) \tilde{u}_{y}, \end{gather}$$
where 
$D={\textrm {d}}/{\textrm {d} \bar y}$.
The above equations (2.10) and (2.11) are subject to the following boundary conditions. At the air–liquid interface, let the free surface be specified by 
$y=1+h(x,t)$, where 
$h(x,t)$ is the displacement of the surface from mean position 
$y=1$. At the air–liquid interface, the dimensionless equation governing the transport of the surfactant is (Halpern & Frenkel Reference Halpern and Frenkel2003)
\begin{equation} \frac{ \partial}{\partial t} (\sqrt{1+h_x^{2}} \,\zeta)+ \frac{ \partial }{\partial t} (\sqrt{1+h_x^{2}} \, \zeta \, v_x) = \frac{1}{Pe} \frac{\partial}{\partial x} \left( \frac{1}{\sqrt{1+h_x^{2}}} \frac{\partial \zeta}{\partial x} \right), \end{equation}
where 
$h_x = {\partial h}/{\partial x}$ and 
$Pe = V_m R/D_s$ is the Péclet number with 
$D_s$ as the surface diffusivity of the surfactant. We assume an insoluble pulmonary surfactant for which 
$D_s/(V_m R) \ll 1$, thus the diffusion term can be neglected. At 
$y=1$, the boundary conditions are the kinematic boundary condition, the surfactant evolution equation and the continuity of the tangential and normal stresses (Smith & Davis Reference Smith and Davis1982)
$$\begin{gather} \textrm{i}(k\bar v_x-\omega) \tilde h = \tilde{v}_y, \end{gather}$$
$$\begin{gather}\textrm{i}(k\bar v_x-\omega) \tilde \zeta ={-}k \tilde{v}_x, \end{gather}$$
$$\begin{gather}\tilde{\tau}^{p}_{xy}+\beta (\textrm{i}k\tilde{v}_{y}+D\tilde{v}_{x})-2\textrm{i}k(1-\beta)W \tilde h=\textrm{i}k Ma \tilde \zeta, \end{gather}$$
$$\begin{gather}-\tilde{p} + \tilde{\tau}^{p}_{yy} + 2\beta D\tilde{v}_{y} - 2\textrm{i}k \tilde h={-} T k^{2} \tilde h, \end{gather}$$
where 
$h(x,t)=\tilde h \,\textrm {e}^{\textrm {i}(kx-\omega t)}, \zeta (x,t)=\tilde {\zeta } \,\textrm {e}^{\textrm {i}(kx-\omega t)}$, 
$Ma={\gamma ^{*} \zeta _0^{*}}/{\mu _{t}^{*} V_m}$ is the Marangoni number and 
$T=T^{*}/(\mu V_m)$ is the dimensionless interfacial tension with 
$T^{*}$ as the dimensional interfacial tension. At 
$y=0$, i.e. at the liquid–deformable-solid interface, velocity and stress continuity give
$$\begin{gather} \tilde{v}_{y}={-}\textrm{i} \omega \tilde{u}_{y}, \end{gather}$$
$$\begin{gather}\tilde{v}_{x}+ \tilde{u}_{y}={-}\textrm{i} \omega \tilde{u}_{x}, \end{gather}$$
$$\begin{gather}\frac{1}{\varGamma} D \tilde{u}_{x}+ \frac{\textrm{i}k}{\varGamma} \tilde{u}_{y}= \tilde{\tau}^{p}_{xy}+\beta (\textrm{i}k\tilde{v}_{y}+D\tilde{v}_{x}) -2\textrm{i}k(1-\beta)W \tilde{u}_{y}, \end{gather}$$
$$\begin{gather}-\tilde{p}_{g}+\frac{2}{\varGamma} (\textrm{i}k \varGamma+D) \tilde{u}_{y} ={-}\tilde{p} + \tilde{\tau}^{p}_{yy} + 2\beta D\tilde{v}_{y}. \end{gather}$$
At 
$y=-H$, as deformable solid is perfectly bonded to the rigid solid, the displacement of the deformable solid will be zero i.e.
\begin{equation} \tilde{u}_{x}=0, \quad \tilde{u}_{y}=0. \end{equation} The present study aims to probe the instabilities induced by the airflow in the airways and oral area and the role of the air, saliva, mucus and muscles lined by the saliva or mucus in causing such instabilities. Thus, the parameter values used in probing the instabilities will be restricted to the physically relevant range. The dimensional and dimensionless parameters relevant in the processes under consideration are listed in tables 1 and 2. From table 2, it is clear that the mode of instability predicted by Smith & Davis (Reference Smith and Davis1982), which becomes unstable for 
$Re>34.2$, will not be relevant for the processes under consideration.
Table 1. Order of magnitude of dimensional parameters, thickness of the viscoelastic liquid (
$R$), thickness of the deformable muscle layer (
$HR$), relaxation time of the viscoelastic liquid (
$\lambda$), air speed past the liquid (
$V_a$) and the shear modulus of the deformable layer (
$G$) for various parts of the airways and oral area relevant during breathing, talking, coughing and sneezing actions. For the oral area, the data are from Collins & Dawes (Reference Collins and Dawes1987), Cheng et al. (Reference Cheng, Gandevia, Green, Sinkus and Bilston2011), Nakamori et al. (Reference Nakamori2020), Stokes & Davies (Reference Stokes and Davies2007), Zussman et al. (Reference Zussman, Yarin and Nagler2007) and Geoghegan et al. (Reference Geoghegan, Laffra, Hoogendorp, Taylor and Jermy2017), for the trachea from Sims & Horne (Reference Sims and Horne1997), Lawrence et al. (Reference Lawrence, Branson, Oliva and Rubinowitz2014), Hwang et al. (Reference Hwang, Litt and Forsman1969), Wang et al. (Reference Wang, Mesquida and Lee2011), Tang et al. (Reference Tang, Nicolle, Pantelic, Koh, Wang and Amin2012), Chen et al. (Reference Chen, Zhong, Luo, Deng, Hu and Song2019) and Shields & Jeffery (Reference Shields and Jeffery1987) and for bronchi and bronchioles from Hwang et al. (Reference Hwang, Litt and Forsman1969), Wang et al. (Reference Wang, Mesquida and Lee2011), Theriault et al. (Reference Theriault, Eddy, Borgaonkar, Babar and Manos2018) and Lai et al. (Reference Lai, Wang, Wirtz and Hanes2009). Other properties such as density and viscosity of the liquids are 
${\sim }10^{3}$ kg m
$^{-3}$ and 
${\sim }10^{2}\text {--}10^{-3}$ Pa s, respectively and air–liquid interfacial tension is 
${\sim } 10^{-1}\text {--}10^{-2}$ N m
$^{-1}$ (Hwang et al. Reference Hwang, Litt and Forsman1969; Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Hamed & Fiegel Reference Hamed and Fiegel2014; Chen et al. Reference Chen, Zhong, Luo, Deng, Hu and Song2019). For the oral area, the viscoelastic liquid is saliva while for the airways the viscoelastic liquid is the mucus. For the trachea, the parameters are for the mucosa and submucosa layers since the shear modulus of the musculocartilaginous layer and adventitia is high enough to treat them as a rigid solid (Wang et al. Reference Wang, Mesquida and Lee2011). Using 
$V_a$, the maximum velocity 
$V_m$ is 
$10^{-2}\text {--}10^{-6}$ m s
$^{-1}$. The power-law index for the mucus of a healthy individual is 
$n \sim 0.1\text {--}0.5$ (Lai et al. Reference Lai, Wang, Wirtz and Hanes2009).
Table 2. Order of magnitude of the dimensionless parameters using the dimensional parameters listed in table 1. The range of Reynolds number is 
$10^{-2}\text {--}10^{-5}$ while the range of 
$T$ is 
$10^{3}\text {--}10^{-1}$.
3. Numerical methodology
3.1. For 
$Re=0,n=1$ and $\beta =0$
In the creeping-flow limit (
$Re = 0$), for a Newtonian fluid (
$W=0$ and 
$n=1$), the perturbation equations for the liquid (2.10) can be analytically solved to give
\begin{equation} \tilde{v}_y = \textrm{e}^{{-}k y} (a_1 + y a_2) + \textrm{e}^{k y} (a_3 + y a_4), \end{equation}
where 
$a_1,a_2,a_3$ and 
$a_4$ are constants. For a UCM fluid (
$n=1$ and 
$\beta =0$), the reduced governing equation is (Gorodtsov & Leonov Reference Gorodtsov and Leonov1967)
\begin{equation} (y_n^{2} D^{2} - 2y_n D-k^{2}y_n^{2}+2)(D^{2}+2\textrm{i}kD-2k^{2}W^{2}-k^{2})\tilde{v}_y =0, \end{equation}
where 
$y_n=y-\omega /k-\textrm {i}/(kW)$. The solution of the above differential equation is
\begin{align} \tilde{v}_y &= \frac{ a_1(\omega-ky_n)}{2 k^{3} W (\textrm{i} + W)} \exp(\omega + \textrm{i}/W - k y_n ) + \frac{a_2 (\omega - k y_n)}{4 k^{4} W (-\textrm{i} + W)} \exp(-\omega - \textrm{i}/W + k y_n ) \nonumber\\ &\quad + a_3 \exp\left[\frac{-(\textrm{i} k W + \sqrt{k^{2} (1 + W^{2})}) (-\textrm{i} - \omega W + k W y_n )}{k W} \right] \nonumber\\ &\quad + a_4 \exp\left[\frac{(-\textrm{i} k W + \sqrt{k^{2} (1 + W^{2})}) (-\textrm{i} - \omega W + k W y_n )}{k W} \right]. \end{align}
Similarly, for the solid, the perturbation governing equations (2.11) can be solved to obtain the solution for 
$\tilde {u}_y$ as
\begin{equation} \tilde{u}_y = \textrm{e}^{{-}k y} [c_1 + \textrm{e}^{-\textrm{i} \varGamma k y} (\textrm{e}^{2 k y} c_2 + c_3) + \textrm{e}^{2 k y} c_4], \end{equation}
where 
$c_1,c_2,c_3$ and 
$c_4$ are constants.
The dispersion relation can be obtained as follows. Substitute the solutions (3.3) and (3.4) in boundary conditions (2.12). The corresponding expressions for the liquid quantities 
$\tilde {v}_x,\tilde {p}$ and solid quantities 
$\tilde {u}_x,\tilde {p}_g$ can be obtained from the respective continuity and 
$x$-component of the momentum equations. Upon substitution, the result is eight equations using which a matrix can be formed as
\begin{equation} \boldsymbol{M}\boldsymbol{a}=0, \end{equation}
where 
$\boldsymbol {M}$ and 
$\boldsymbol {a}$ are the matrix formed of the coefficients of the eight equations and a vector containing the eight constants 
$a_1,\ a_2,\ a_3,\ a_4,\ c_1,\ c_2,\ c_3$ and 
$c_4$, respectively. The determinant of the matrix 
$\boldsymbol {M}$ then gives the required dispersion relation. The general dispersion relation derived above is unwieldy but can be reduced by specifying the value of the parameters and then solving the resulting equation for 
$\omega$. Thus, the stability analysis reduces to finding an eigenvalue 
$\omega$ with 
$\omega _i >0$ as the parameters 
$k,H,\varGamma$ and 
$W$ are varied.
3.2. Arbitrary $Re,n$ and $\beta$
To carry out the linear stability analysis of (2.10) and (2.11) subject to boundary conditions (2.12) at an arbitrary 
$Re$, the pseudo-spectral method is employed, in which the eigenfunctions corresponding to each dynamic field are expanded into series of the Chebyshev polynomials as
\begin{equation} \tilde{f}(y)=\sum_{m=0}^{m=N} a_m {\mathcal{T}}_m (y), \end{equation}
where 
${\mathcal {T}}_m(y)$ are Chebyshev polynomials of degree 
$m$ and 
$N$ 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
\begin{equation} \boldsymbol{A}\boldsymbol{e} + \omega \boldsymbol{B}\boldsymbol{e} + \omega^{2} \boldsymbol{C}\boldsymbol{e}=0, \end{equation}
where 
$\boldsymbol {A}$, 
$\boldsymbol {B}$ and 
$\boldsymbol {C}$ are matrices obtained from the discretisation procedure and 
$\boldsymbol {e}$ is the vector containing the coefficients of all series expansions.
Further details of the discretisation of the governing equations and boundary conditions, and of the construction of the matrices 
$\boldsymbol {A}$, 
$\boldsymbol {B}$ and 
$\boldsymbol {C}$ can be found in the standard procedure described by Trefethen (Reference Trefethen2000) and Schmid & Henningson (Reference Schmid and Henningson2001). We use the polyeig MATLAB routine to solve the constructed, generalised eigenvalue problem given by (3.7). To filter out the spurious modes from the genuine, numerically computed spectrum of the problem, the latter is determined for 
$N$ and 
$N+2$ collocation points, and the eigenvalues are compared with an a priori specified tolerance, e.g. 
$10^{-4}$. The genuine eigenvalues are 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 sixth significant digit, the same number of collocation points is used to determine the critical parameters of the system. As 
$W$ increases, the number of collocation points to capture the unstable mode also increases. Thus, in the subsequent analysis, we restrict the analysis to 
$W<1$.
4. Results and discussion
Before proceeding further, we validate the numerical methodology utilised here in two ways. First, for arbitrary 
$Re$, the curve predicted by Smith & Davis (Reference Smith and Davis1982) for vanishing air–liquid interfacial tension is reproduced and compared with the digitally extracted data points in figure 2. The critical Reynolds number (
$Re_c$) for the unstable Newtonian mode for 
$T=0$ is predicted to be 
$34.2$, in excellent agreement with Smith & Davis (Reference Smith and Davis1982). Second, the unstable mode predicted by using the pseudo-spectral method in creeping flow is compared with the one obtained using the analytical method outlined in § 3.1, thereby validating the pseudo-spectral code in two ways.
Figure 2. The neutral stability curve (
$\omega _i=0$) at 
$W=0$, 
$n=1$, 
$Ma=0$ and 
$T=0$ corresponding to the curve for vanishing air–liquid interfacial tension in figure 3 of Smith & Davis (Reference Smith and Davis1982). The data points from Smith & Davis (Reference Smith and Davis1982) are digitally extracted. The excellent agreement between the present numerical approach and data points of Smith & Davis (Reference Smith and Davis1982) validates the former.
To understand the individual roles of the liquid and deformable solid in triggering the instability for a viscoelastic liquid sheared by air and flowing past a deformable layer, four different cases are considered here. First, a Newtonian liquid sheared by the air and flowing past a rigid solid, henceforth abbreviated as ‘NR’, where N stands for Newtonian and ‘R’ stands for the rigid. This case has been studied in detail by Miles (Reference Miles1960) and Smith & Davis (Reference Smith and Davis1982). Their analyses reported the existence of an instability for 
$Re>34.2$. Thus, the NR case will not be studied further but the related results will be frequently referred to from Miles (Reference Miles1960) and Smith & Davis (Reference Smith and Davis1982). The instability mode predicted by Smith & Davis (Reference Smith and Davis1982) will be henceforth referred to as the ‘Newtonian mode.’ Second, a shear-thinning viscoelastic liquid sheared by the air and flowing past a rigid solid, henceforth abbreviated as ‘VR’, will be considered where ‘V’ stands for viscoelastic. The VR case will shed light on the role of the elasticity of the liquid in triggering the instability. To understand the effect of the solid elasticity, the third case considers a Newtonian liquid sheared by the air and flowing past a deformable solid, hereafter abbreviated as ‘ND’, where ‘D’ stands for deformable. Finally, the fourth case is a shear-thinning viscoelastic liquid sheared by the air and flowing past a deformable solid, hereafter abbreviated as ‘VD’. This case will show the combined effect of the elastic nature of the liquid and solid. The effect of the pulmonary surfactant on the instabilities predicted in the present study is discussed in § 4.4.
4.1. Viscoelastic liquid sheared by the air and flowing past a rigid solid (VR)
To understand the role of the elasticity of the liquid in introducing instability, first, we study a shear-thinning viscoelastic liquid layer sheared by the air and flowing past a rigid solid, i.e. 
$G \rightarrow \infty$ or 
$\varGamma \rightarrow 0$. For the VR case, the boundary conditions at 
$y=0$ are the no slip and impermeability, i.e. 
$\tilde {v}_x =0$ and 
$\tilde {v}_y=0$. For ease of the discussion and presentation, the results have been divided into two sections dealing with stability in the creeping-flow limit and non-zero 
$Re$.
4.1.1. Creeping flow
A typical spectrum of the VR in the creeping flow is shown in figure 3(a). The figure shows the existence of a downstream travelling unstable mode (
$\omega _r>0$). From the analysis of Smith & Davis (Reference Smith and Davis1982), the NR case is linearly stable in the creeping flow. However, from figure 3(a), there is an unstable mode for a viscoelastic liquid layer, hereafter termed the ‘liquid elastic mode’. This clearly shows that the elasticity of the liquid leads to the predicted instability. Another discrete mode in figure 3(a) is 
$\omega =k$, which remains stable for the parameters considered here. Thus, its variation with the other parameters will not be explored in the following discussion. Using the analytical solutions (3.3) and (3.4), the shooting code employing the vpasolve MATLAB function gives 
$\omega = 0.433075 + 0.023960i$, in excellent agreement with the eigenvalue obtained using the pseudo-spectral method thereby validating the latter. This further affirms the genuine nature of the liquid elastic mode.
Figure 3. Panel (a) shows the spectrum of the flow at 
$Re=0$, 
$k=0.5$, 
$T=1$, 
$n=1$, 
$Ma=0$, 
$\beta =0$ and 
$W=0.3$ for the VR case obtained using the pseudo-spectral method. A downstream travelling unstable mode (
$\omega = 0.433075 + 0.023960i$) exists in the spectrum since 
$\omega _i>0$. Such an unstable mode is absent for the Newtonian liquid layer sheared by the air (Smith & Davis Reference Smith and Davis1982), thereby showing that the liquid elasticity is responsible for the destabilisation of this mode. Besides the unstable mode, a neutrally stable mode with 
$\omega =k$ is also present. The overlap of the unstable mode for 
$N=50$ and 
$N=75$ confirms the genuine nature of the mode. The continuous spectrum line at 
$\omega _i=-3.3333$ is made up of multiple stable modes. The number of modes in the continuous spectrum depends on the number of collocation points utilised in obtaining the spectrum. The same unstable mode has been also predicted by the shooting code. Panel (b) illustrates the movement of the continuous spectrum due to variation in 
$W$. For 
$W=0.3$ and 
$W=0.1$, the continuous spectrum is located at 
$\omega _i\sim -3.3333$ and 
$\omega _i \sim -10$, respectively, which confirms the rule for the continuous spectrum as 
$\omega _i=-1/W$.
The continuous spectrum in figure 3(a) corresponds to 
$\omega _i \sim -3.3333$, which is an equivalent of the high-frequency Gordotsov–Leonov line for the plane Couette flow reported in the literature (Renardy & Renardy Reference Renardy and Renardy1986; Joshi & Shankar Reference Joshi and Shankar2019). It must be noted that the plane Couette flow of a viscoelastic liquid exhibits a high-frequency Gordotsov–Leonov line at 
$\omega _i \sim -1/(2W)$ (Joshi & Shankar Reference Joshi and Shankar2019). However, in the present case from figure 3(b), we observe that the continuous spectrum line obeys the rule 
$\omega _i \sim -1/W$. Therefore, a viscoelastic liquid layer sheared by the air not only introduces a new mode of instability but also affects the continuous spectrum compared with the plane Couette flow of a viscoelastic liquid.
To further understand the new mode of instability, the variation of 
$\omega _i$ with 
$k$ for select values of 
$W,n,\beta$ and 
$T$ is shown in figure 4. For 
$T=0$, Smith & Davis (Reference Smith and Davis1982) predicted that, for the Newtonian mode to become unstable, 
$Re>34.2$, thus 
$Re_c=34.2$, thereby showing that a minimum 
$Re$ is needed for the existence of the Newtonian mode. However, from figure 4, the flow is unstable at an arbitrary 
$W$, implying an absence of the minimum 
$W$ needed for the instability to exist. Thus, a viscoelastic liquid layer sheared by the air and flowing past a rigid solid is unconditionally unstable. However, from figure 4(a), the parameter 
$W$ does affect the growth rate of the perturbations. Also, at low 
$k$, the growth rate shows characteristic scaling 
$\omega _i \sim Wk^{2}$. According to the analysis of Smith & Davis (Reference Smith and Davis1982), 
$Re_c$ for a Newtonian liquid rapidly increases with an increasing 
$T$. However, from figure 4(b), the unstable mode still exists, albeit with a lesser growth rate and at low values of 
$k$.
Figure 4. Variation in the growth rate 
$\omega _i$ with the disturbance wavenumber 
$k$ in the creeping-flow limit at 
$Ma=0$. Panel (a) shows that the liquid layer is unconditionally unstable since a minimum 
$W$ is unnecessary for the instability to exist. The Weissenberg number only affects the growth rate of the disturbances. The curve for 
$W=0.1$ is incomplete due to the elimination of the negative values of 
$\omega _i$ to use a log scale. Panel (b) shows the effect of the surface tension on the dispersion curves. An increasing surface tension suppresses the instability at moderate and high wavenumbers but fails to stabilise at low wavenumbers, thereby showing that the flow remains unconditionally unstable. Panel (c) illustrates the destabilising effect of the shear-thinning nature of the fluid. Panel (d) shows the decrease in the growth rate of the perturbations as the contribution of the solvent viscosity signified by 
$\beta$ is increased. Parameters are: (a) 
$n=1$, 
$\beta =0$ and 
$T=0$; (b) 
$n=1$, 
$\beta =0$ and 
$W=0.1$; (c) 
$W=0.1$, 
$\beta =0$ and 
$T=1$; (d) 
$n=1$, 
$T=1$ and 
$W=0.1$.
The effect of the shear thinning on the growth rate is illustrated in figure 4(c). A decreasing power-law index signifies increasing shear thinning. From figure 4(c), as the shear-thinning nature of the fluid increases the growth rate of the liquid elastic mode increases. The increasing contribution of the solvent to the solution, i.e. 
$\beta$ leads to a decrease in the growth rate of the liquid elastic mode, as shown in figure 4(d). This is to be expected since the NR case is linearly stable in the creeping-flow limit and increasing 
$\beta$ implies an increase in the Newtonian solvent contribution.
The mechanism of the liquid elastic mode can be understood as follows. For simplicity, consider the case of a UCM fluid sheared by air which corresponds to 
$n=1$ and 
$\beta =0$. For the existence of an instability, a term is necessary which can lead to an energy transfer from the base-state quantities to the perturbation quantities. From (2.10) for 
$n=1$ and 
$\beta =0$, there are several terms including the convective terms 
$\textrm {i}kWy\tilde {\tau }_{xx}$, 
$\textrm {i}kWy\tilde {\tau }_{xy}$ and 
$\textrm {i}kWy\tilde {\tau }_{yy}$ and other interaction terms 
$-4\textrm {i}kW^{2} \tilde {v}_x$, 
$-2\textrm {i}kW^{2} \tilde {v}_y$, 
$-2WD\tilde {v}_x$, 
$-2W\tilde {\tau }_{xy}$ in the constitutive equations that satisfy this criterion. However, similar terms are also present for the plane Couette flow of UCM liquid yet it is linearly stable in the creeping-flow limit (Gorodtsov & Leonov Reference Gorodtsov and Leonov1967; Renardy & Renardy Reference Renardy and Renardy1986; Joshi & Shankar Reference Joshi and Shankar2019). This indicates that the terms responsible for the instability exhibited by VR may not be the terms present in the linearised constitutive equation. The other terms that represent the interaction between the base-state and perturbed state quantities are in the boundary conditions at 
$y=1$, i.e. 
$(-2\textrm {i}k\tilde h)$, which is also present for the Newtonian liquid, and 
$(-2\textrm {i}kW\tilde h)$, which is only present for the viscoelastic liquid. Since the NR case is stable in the creeping flow, the instability predicted in the present study could be due to the term 
$-2\textrm {i}kW\tilde h$.
To verify the above hypothesis, eigenvalues are obtained with/without the term 
$(-2\textrm {i}kW\tilde h)$ in the boundary condition (2.12d), as shown in table 3. Table 3 shows that the above hypothesis is indeed correct and the term responsible for the existence of an unstable liquid elastic mode is 
$(-2\textrm {i}kW\tilde h)$. The term 
$(-2\textrm {i}kW\tilde h)$ in the basic form is 
$(\textrm {i}k (\bar \tau _{yy}-\bar \tau _{xx}) \tilde h)$, where the quantity 
$(\bar \tau _{yy}-\bar \tau _{xx})$ is the first normal stress difference across the air–liquid interface in the base state. Thus, the term 
$(-2\textrm {i}kW\tilde h)$ originates from the first normal stress difference, a feature unique to a viscoelastic liquid, thereby showing that the liquid elastic mode is driven by the first normal stress difference exhibited by the viscoelastic liquid. The plots of the normalised perturbations for the liquid elastic mode are shown in figure 5. Since the energy for the perturbations is supplied by the term 
$(-2\textrm {i}kW\tilde h)$ at the air–liquid interface, the perturbations exhibit the highest variation at 
$y=1$.
Figure 5. The normalised perturbations, for the liquid at 
$W=0.3$, 
$\beta =0$, 
$n=1$, 
$T=1$, 
$k=0.5$, 
$Ma=0$ and 
$Re=0$ for the eigenvalue 
$\omega =0.433075 + 0.023960$. Here, 
$v'_x=Re[ \tilde {v}_x \,\textrm {e}^{\textrm {i}kx} ]$ and 
$v'_y=Re[ \tilde {v}'_y \,\textrm {e}^{\textrm {i}kx} ]$. For convenience, the axes have been normalised to the interval 
$[0,1]$. The length of the domain in the 
$x$-direction is equal to a wavelength (
$2{\rm \pi} /k$) of the perturbations. The plots shows that the perturbations exhibit maximum variation near the air–liquid interface, thereby indicating the destabilisation caused by the term 
$(-2\textrm {i}kW\tilde h)$. The perturbations vanish at 
$y=0$ as needed by the no-slip and impermeability conditions at 
$y=0$; (a) 
$v'_x$ and (b) 
$v'_y$.
Table 3. The most unstable eigenvalue (
$\omega =\omega _r +\textrm {i} \omega _i$) for the VR case with/without the 
$(-2\textrm {i}kW\tilde h)$ term in the tangential stress boundary condition (2.12d) at 
$n=1$, 
$\beta =0$, 
$Ma=0$ and 
$T=1$. The presence of the stable eigenvalues in the second column of the table shows that the term responsible for the predicted instability is 
$(-2\textrm {i}kW\tilde h)$. The term 
$(-2\textrm {i}kW\tilde h)$ originates from the first normal stress difference across the air–liquid interface. Thus, the liquid elastic mode is triggered by the first normal stress difference exhibited by a viscoelastic liquid.
4.1.2. A simple model for liquid elastic mode
The above discussion shows that the driving force for the liquid elastic mode is the first normal stress difference exhibited by a viscoelastic solid across the air–liquid interface. In the following analysis, we derive a simple model which will incorporate the basic driving force for the liquid elastic mode based on the long-wave analysis. Thus, we scale 
$x=X/\epsilon$ and 
$t=\tau /\epsilon$ and expand the fluid quantities as
$$\begin{gather} v_x = \bar v_x + \epsilon v'_x +\cdots, \end{gather}$$
$$\begin{gather}v_y = \epsilon v'_y +\cdots, \end{gather}$$
$$\begin{gather}p = p'+\cdots, \end{gather}$$
where 
$\epsilon \ll 1$. Further, we assume creeping flow (
$Re=0$) and UCM fluid in the absence of the surfactant which implies 
$n=1,\beta =0$ and 
$Ma=0$. The above expansion closely follows the linear stability analysis derivation discussed in § 2. Next, we substitute the above expansions in (2.5). At 
$O(1)$, recalling 
$\bar p=0$ for the present problem, we obtain the base-state quantities specified in (2.8) and the air–liquid interface is flat. At 
$O(\epsilon )$ we obtain the perturbed state equations which can be further simplified due to the long-wave approximation
$$\begin{gather} \frac{\partial^{2} v'_x}{\partial y^{2}}-\frac{\partial p'}{\partial X}=0, \end{gather}$$
$$\begin{gather}\frac{\partial p'}{\partial y}=0. \end{gather}$$
The boundary conditions at 
$O(\epsilon )$ are
$$\begin{gather} \text{at} \ y=0,\quad v'_x=0, \end{gather}$$
$$\begin{gather}\text{at} \ y=1+\epsilon h(X,\tau),\quad \frac{\partial v'_x}{\partial y} - \underline{2\hat W \frac{\partial h}{\partial X}}=0; \quad p'={-}\hat T \frac{\partial^{2} h}{\partial X^{2}}, \end{gather}$$
where 
$\hat W=\epsilon W$ and 
$\hat T=\epsilon ^{3} T$ and the underlined term is due to the first normal stress difference exhibited by the viscoelastic fluid. Thus, 
$v'_x$ at 
$O(\epsilon )$ is
\begin{equation} v'_x = \frac{1}{2} y \left[ 4\hat W \frac{\partial h}{\partial X} -\hat T (y-2) \frac{\partial^{3} h}{\partial X^{3}} \right]. \end{equation}
The kinematic boundary condition at 
$O(\epsilon )$ is
\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial h}{\partial X} + \frac{\partial}{\partial X} \int_0^{1} v'_x \,\textrm{d} y=0. \end{equation}Substituting (4.8) in the above equation yields
\begin{equation} \frac{\partial h}{\partial \tau} + \frac{\partial h}{\partial X} + \frac{\partial}{\partial X} \left[ \frac{1}{3} \left(3\hat W\frac{\partial h}{\partial X} +\hat T \frac{\partial^{3} h}{\partial X^{3}} \right) \right]=0. \end{equation}Rescaling the above equation yields
\begin{equation} \frac{\partial h}{\partial t} + \frac{\partial h}{\partial x} + \frac{\partial}{\partial x} \left[ \frac{1}{3} \left(3 W\frac{\partial h}{\partial x} + T \frac{\partial^{3} h}{\partial x^{3}} \right) \right]=0, \end{equation}
where 
$h$ is still a small deviation of the air–liquid interface from the base-state equation 
$y=1$. Equation (4.11) is strictly applicable in the limit of small perturbations.
To carry out the linear stability analysis, we substitute 
$h=\tilde h \,\textrm {e}^{\textrm {i}(kx-\omega t)}$ and simplify to yield
\begin{equation} \omega = k - \tfrac{1}{3} \textrm{i} k^{4} T + \textrm{i} k^{2} W. \end{equation}
A comparison between the dispersion curves predicted by using the simple model and the numerical approach is shown in figure 6. The above model correctly predicts 
$\omega _i \sim Wk^{2}$ at low 
$k$, as shown in figure 4(a), thereby capturing the importance of the first normal stress difference exhibited by the viscoelastic fluid at the air–liquid interface. Also, the simple model succeeds in capturing the stabilisation due to the increasing air–liquid interfacial tension shown in figure 4(b). However, for 
$k>0.1$, both the models disagree, which clearly points out the stabilising influence of the other terms at high 
$k$, a feature missing in the simple model.
Figure 6. A comparison of the growth rate predicted by the proposed simple model and numerical prediction for the liquid elastic mode at 
$n=1$, 
$W=0.1$, 
$T=0.1$ and 
$\beta =0$ in the creeping-flow limit. Panel (a) shows the dispersion curves predicted by both methods for the whole range of 
$k$. Panel (b) illustrates the excellent agreement between two methods for 
$k<0.1$. (a) Complete curves and (b) magnified low 
$k$ region.
It must be noted that, here, we have assumed a planar flow geometry due to the small thickness of the saliva and mucus layers which removes the curvature effects through the air–liquid interfacial tension. From the analysis of Zhou et al. (Reference Zhou, Peng, Zhang and Zhuge2014), the surface tension plays a destabilising role on the long-wave (
$k<1$) instabilities for a viscoelastic film coating a tube and flowing under the action of gravity and imposed shear at the free surface. Thus, given the presence of the long-wave instabilities in the present study if we consider a flow coating a deformable tube and sheared by the airflow then the surface tension may result in the growth rate enhancement of the liquid elastic mode.
4.1.3. Non-zero $Re$
The effect of an increasing inertia, i.e. 
$Re$, on the growth rate of the perturbations is shown in figure 7. For the NR case, Smith & Davis (Reference Smith and Davis1982) predicted a destabilising effect of increasing inertia on the Newtonian mode. However, from figure 7, an increasing inertia has a stabilising effect on the liquid elastic mode. The stabilising effect of increasing inertia further affirms the purely elastic nature of the liquid elastic mode.
Figure 7. The effect of variation in 
$Re$ on the growth rate of the perturbations for the liquid elastic mode at 
$n=1$, 
$\beta =0$, 
$Ma=0$, 
$T=1$ and 
$W=0.1$. The figure shows the stabilising effect of an increasing inertia on the elastic mode.
4.2. Newtonian liquid sheared by the air and flowing past a deformable solid (ND)
To understand the role of the deformable-solid layer, next, we consider the ND case for which 
$n=1$ and 
$\lambda \rightarrow 0$ or 
$W\rightarrow 0$. Similar to the VR case, the results are divided into two sections dealing with the results in the creeping-flow limit and at finite 
$Re$.
4.2.1. Creeping flow
For the ND case, there are three eigenvalues in the creeping-flow limit, as shown in table 4. The unstable mode is a downstream travelling mode while the remaining two stable modes are upstream travelling. As per the stability analysis of Smith & Davis (Reference Smith and Davis1982), the NR case is linearly stable in the creeping-flow limit. Thus, the instability predicted here is a pure consequence of the elasticity of the deformable layer. Henceforth the mode of instability introduced by the elastic nature of the solid will be referred to as the ‘solid elastic mode.’
Table 4. Typical eigenspectrum in the creeping-flow limit at 
$k=0.8$, 
$H=10$, 
$\varGamma =0.5$ and 
$T=1$ for the ND case by using the analytical solution. The other two modes always remain stable. The unstable mode predicted by the pseudo-spectral method is 
$0.607433 + 0.052304i$, which shows an excellent agreement with the eigenvalue predicted by the analytical solution thereby validating the former.
To further understand the role of the deformable layer in introducing the solid elastic mode, dispersion curves are plotted for select values of 
$\varGamma$ and 
$T$ in figure 8. A plane Couette flow of a Newtonian liquid past a deformable layer in the creeping-flow limit exhibits a finite-wave instability referred to as the ‘viscous instability’ in the literature (Kumaran et al. Reference Kumaran, Fredrickson and Pincus1994; Kumaran & Muralikrishnan Reference Kumaran and Muralikrishnan2000; Gkanis & Kumar Reference Gkanis and Kumar2003; Patne & Shankar Reference Patne and Shankar2017; Patne et al. Reference Patne, Giribabu and Shankar2017; Joshi & Shankar Reference Joshi and Shankar2019). For the existence of the viscous instability, a minimum value of 
$\varGamma$ is necessary, thus a finite critical 
$\varGamma$, i.e. 
$\varGamma _c$ exists. However, the dispersion curves of figure 8(a) illustrate that, for an arbitrary 
$\varGamma$, the flow is unstable, thereby ruling out the existence of 
$\varGamma _c$. Therefore, a Newtonian liquid sheared past a neo-Hookean solid is unconditionally unstable. The effect of the air–liquid interfacial tension 
$T$ on the solid elastic mode is shown in figure 8(b). An increasing 
$T$ has a stabilising effect on the high-wavenumber part of the solid elastic mode but fails to stabilise the low-wavenumber part of the solid elastic mode. From figure 8(b), the growth rate increases with an increasing 
$H$. However, for 
$H>10$ the growth rate is relatively insensitive to variation in 
$H$.
Figure 8. Variation in the growth rate 
$\omega _i$ with the disturbance wavenumber 
$k$ in the creeping-flow limit for the ND case. Panel (a) shows the dispersion plots for 
$Re=0$, 
$Ma=0$ and 
$T=0$. The plots show that, similar to the viscoelastic liquid layer sheared past a rigid solid, the Newtonian liquid layer sheared past a deformable layer is unconditionally unstable since a minimum 
$\varGamma$ is not needed for the instability to exist. The parameter 
$\varGamma$ only affects the growth rate of the disturbances. Panel (b) shows the effect of the air–liquid interfacial tension on the dispersion curves for 
$\varGamma =0.5$. An increasing surface tension suppresses the instability at moderate and high wavenumbers but fails to stabilise at low wavenumbers, thereby showing that the flow remains unconditionally unstable. Panel (c) shows the effect of variation in the dimensionless thickness of the deformable solid on the growth rate of the solid elastic mode; (a) 
$H=10$ and 
$T=0$, (b) 
$H=10$ and 
$\varGamma =0.5$ and (c) 
$T=1$ and 
$\varGamma =0.5$.
The mechanism of the solid elastic mode can be understood as follows. Kumaran et al. (Reference Kumaran, Fredrickson and Pincus1994), Kumaran & Muralikrishnan (Reference Kumaran and Muralikrishnan2000) and Gkanis & Kumar (Reference Gkanis and Kumar2003) showed that the viscous instability is triggered due to the shear work done by the flowing liquid on the deformable layer. This shear work is transferred via the tangential velocity continuity interface condition (2.12g) at the liquid–deformable-solid interface at 
$y=0$. The particular term responsible for the transfer of the shear work in (2.12g) is 
$(\tilde {u}_y)$, which in basic form is 
$(D\bar v_x \tilde {u}_y)$. However, in the present study for the ND case, another term 
$(-2\textrm {i}k\tilde h)$ at 
$y=1$ in the normal stress continuity boundary condition (2.12e) is equally important, as shown in table 5. In basic form, the term 
$(-2\textrm {i}k\tilde h)$ can be written as 
$(-2\textrm {i}k \bar \tau _{xy} \tilde h)$. But at 
$y=1$, in the base state, the shear stress in the liquid is equal to the shear stress exerted by the air thus, in dimensionless terms 
$\bar \tau _{xy} = \tau _a R/(\mu V_m)$, modifying the term to 
$(-2\textrm {i}k \tau _a R/(\mu V_m) \tilde h)$. Therefore, the energy for the solid elastic mode is provided by the shear stress exerted by the air on the liquid. To conclude, the energy for the destabilisation of the solid elastic mode is provided by the shear stresses exerted by the air at the air–liquid interface and by the liquid at the liquid–deformable-solid interface.
Table 5. The most unstable eigenvalue for the ND case with/without 
$(-2\textrm {i}k\tilde h)$ term in the normal stress boundary condition (2.12e) at 
$H=10$, 
$Ma=0$ and 
$T=1$. From the second column of the table the solid elastic mode is stable in the absence of 
$(-2\textrm {i}k\tilde h)$, thereby showing the importance of the same in triggering the solid elastic mode. The shear stress exerted at the air–liquid interface is responsible for the term 
$(-2\textrm {i}k\tilde h)$. A similar vital role is also played by the term 
$(\tilde {u}_y)$ in the tangential velocity continuity condition (2.12g) at the liquid–deformable-solid interface (
$y=0$). Therefore, the solid elastic mode is combined effect of the shear stresses exerted by the air at the air–liquid interface and by the liquid at the liquid–deformable-solid interface.
The plots of the normalised velocity perturbations 
$v'_x$ and 
$v'_y$ in the liquid and displacement perturbations 
$u'_x$ and 
$u'_y$ in the deformable solid are shown in figure 9. Since the destabilisation mechanism of the solid elastic mode is the shear stresses exerted by the air at the air–liquid interface and by the liquid at the liquid–deformable-solid interface, the velocity perturbations exhibit maximum variation near both the air–liquid and liquid–deformable-solid interfaces. However, the displacement perturbations exhibit highest variation near the liquid–deformable-solid interface due to the shear stresses exerted by the liquid on the deformable solid while they vanish at 
$y=-H$, which corresponds to 
$y=0$ in figures 9(c) and 9(d), as required by the boundary condition (2.12j).
Figure 9. The normalised perturbations for the liquid and deformable-solid layer for the ND case at 
$\varGamma =0.5$, 
$T=1$, 
$k=0.8$, 
$H=10$, 
$Ma=0$ and 
$Re=0$ for the solid elastic mode 
$\omega =0.433075 + 0.023960$. Here, 
$u'_x=Re[ \tilde {u}_x \,\textrm {e}^{\textrm {i}kx} ]$ and 
$u'_y=Re[ \tilde {u}'_y \,\textrm {e}^{\textrm {i}kx} ]$. For convenience, the axes have been normalised to the interval 
$[0,1]$. The length of the domain in the 
$x$-direction is equal to a wavelength (
$2{\rm \pi} /k$) of the perturbations. The velocity perturbations (panels a and b) show maximum variation both near the air–liquid (
$y=1$) and liquid–deformable-solid (
$y=0$) interface, thereby indicating the destabilisation caused by the terms 
$(-2\textrm {i}k\tilde h)$ and 
$(\tilde {u}_y)$. For the solid displacement perturbations (panels c and d), the perturbations exhibit maximum at the liquid–deformable-solid (
$y=1$) interface and vanish at the rigid-solid–deformable-solid (
$y=0$) interface; (a) 
$v'_x$, (b) 
$v'_y$, (c) 
$u'_x$ and (d) 
$u'_y$.
It must be noted that, similar to the first normal stress difference exhibited by the viscoelastic fluid at the air–liquid interface, the neo-Hookean solid also exhibits a first normal stress difference across the liquid–deformable-solid interface. As discussed in § 4.1, the first normal stress difference exhibited by the viscoelastic fluid across the air–liquid interface gives rise to the liquid elastic mode. Similarly, Gkanis & Kumar (Reference Gkanis and Kumar2003) predicted that the first normal stress difference exhibited by the neo-Hookean solid also exhibits a new mode of instability, termed the ‘short-wave instability’, for 
$\varGamma >1$. However, for the processes under consideration here, 
$\varGamma <1$, thereby ruling out the presence of the short-wave instability and thus the role of the first normal stress difference exhibited by the neo-Hookean solid at the liquid–deformable-solid interface.
4.2.2. Non-zero $Re$
The effect of the inertia on the solid elastic mode is shown in figure 10. Similar to the liquid elastic mode, increasing inertia has a stabilising effect on the solid elastic mode, thereby confirming its purely elastic nature. The viscous mode of instability predicted by Kumaran et al. (Reference Kumaran, Fredrickson and Pincus1994), Kumaran & Muralikrishnan (Reference Kumaran and Muralikrishnan2000) and Chokshi & Kumaran (Reference Chokshi and Kumaran2007) for the plane Couette flow past a deformable solid in the creeping-flow limit can be destabilised by increasing inertia. For 
$Re>10$, the viscous mode of instability exhibits a characteristic scaling 
$\varGamma _c \sim Re^{-1}$ (Chokshi & Kumaran Reference Chokshi and Kumaran2008a; Giribabu & Shankar Reference Giribabu and Shankar2017; Tanmay et al. Reference Tanmay, Patne and Shankar2018). However, from figure 10, the inertia has a stabilising effect on the solid elastic mode, which further sets apart the viscous mode of instability from the solid elastic mode predicted here. The stabilising effect of the increasing inertia also implies the purely elastic nature of the solid elastic mode similar to the liquid elastic mode.
Figure 10. The effect of variation in 
$Re$ on the growth rate of the solid elastic mode at 
$H=10$, 
$T=1$, 
$Ma=0,$ and 
$\varGamma =0.5$. The figure shows the stabilising effect of an increasing inertia on the solid elastic mode.
4.3. Viscoelastic liquid sheared by the air and flowing past a deformable solid (VD)
In §§ 4.1 and 4.2, the individual influence of the solid and liquid elasticity was shown to manifest in the form of two elastic modes, viz. the liquid elastic mode due to the elasticity of the viscoelastic liquid layer in the absence of the deformable solid and the solid elastic mode due to the elasticity of the solid for which we considered a Newtonian liquid sheared by the air and flowing past a deformable solid. Both modes of instability were found to be unconditionally unstable. In this section, we consider the VD case, which involves the elastic liquid and solid relevant to the processes under consideration. Similar to the VR and ND cases, here, we divide the discussion into two sections dealing with creeping-flow and finite 
$Re$ results.
4.3.1. Creeping flow
From §§ 4.1 and 4.2, the liquid and solid elastic modes are unconditionally unstable. For the VD case, the elastic characteristic of the solid and liquid can either lead to two different unstable modes representing the elasticity of the liquid and solid, or the two modes could reinforce each other, i.e. undergo a resonance to result in a more destructive mode of instability. Figure 11 shows that, for the VD case, the resonance pathway is preferred, thereby leading to a single mode resulting from the resonance of the liquid and solid elastic modes thus the name ‘resonance mode.’ For the resonance mode, both the elastic modes reinforce each other, which leads to a higher growth rate. The resonance between both modes for the whole range of 
$k$ and select values of the other parameters is shown in figure 12.
Figure 11. The solid elastic mode (
$\varGamma =0.2,W=0$), liquid elastic mode (
$\varGamma =0,W=0.2$) and resonance mode (
$\varGamma =0.2,W=0.2$) at 
$n=1$, 
$\beta =0$, 
$Re=0$, 
$k=0.5$, 
$Ma=0$, 
$T=1$ and 
$H=10$ obtained using the pseudo-spectral method. Here, the resonance mode refers to the unstable mode for the VD case. The spectrum shows that the resonance mode exhibits a higher growth rate compared with the solid and liquid elastic modes as a result of the resonance between the solid and liquid elastic modes, thus the name resonance mode. To obtain the liquid elastic mode, 
$\varGamma =0$ and for the solid elastic mode 
$W=0$.
Figure 12. Variation in the growth rate 
$\omega _i$ with the disturbance wavenumber 
$k$ in the creeping-flow limit for the VD case at 
$Ma=0$. The figure illustrates an increase in growth rate due to the resonance between the liquid and solid elastic modes. Panel (a) illustrates the growth rate increase as 
$W$ increases, thereby showing the proof of the resonance between the liquid and solid elastic modes over whole 
$k$. A similar effect is also exhibited if 
$W$ is fixed and 
$\varGamma$ is increased. The stabilising effect of the air–liquid interfacial tension on the resonance mode is illustrated in panel (b). The destabilising effect of increasing 
$H$ on the resonance mode is shown in panel (c). Panel (d) illustrates destabilising effect of the shear thinning on the resonance mode. The stabilising effect of increasing 
$\beta$ is shown in panel (e). The figure illustrates that the resonance mode exhibits features of both the liquid and solid elastic modes, thereby affirming its origin: (a) 
$H=10$, 
$\varGamma =0.2$, 
$\beta =0$, 
$n=1$ and 
$T=1$; (b) 
$H=10$, 
$W=0.2$, 
$\beta =0$, 
$n=1$ and 
$\varGamma =0.2$; (c) 
$\varGamma =0.2$, 
$\beta =0$, 
$n=1$, 
$W=0.2$ and 
$T=1$; (d) 
$T=1$, 
$W=0.2$, 
$\beta =0$, 
$H=10$ and 
$\varGamma =0.2$; (e) 
$T=1$, 
$W=0.2$, 
$H=10$, 
$n=1$ and 
$\varGamma =0.2$.
The resonance mode has combined characteristics of both the modes, thus it is unconditionally unstable. Additionally, the growth rate of the resonance mode will increase with increasing 
$\varGamma ,W$ and 
$H$ or decreasing 
$n$ while its growth rate will decrease with increasing 
$T$ and 
$\beta$, as shown in figure 12. As 
$\beta \rightarrow 1$, the dispersion curve approaches the corresponding curve for the purely solid elastic mode explored in § 4.2. The physical mechanism responsible for the resonance mode is a combination of the mechanisms responsible for the destabilisation of the liquid and solid elastic modes. Thus, the first normal stress difference exhibited by the viscoelastic liquid across the air–liquid interface and the shear stresses exerted by the air at the air–liquid interface and by the liquid on the deformable solid at the liquid–deformable-solid interface are responsible for the existence of the resonance mode.
In addition to the first normal stress difference across the air–liquid interface, there is also the first normal stress difference across the liquid–deformable-solid interface due to the viscoelastic nature of the liquid. From § 4.1, the first normal stress difference across the air–liquid interface leads to the liquid elastic mode. However, the first normal stress difference across the liquid–deformable-solid interface due to the viscoelastic fluid does not introduce a new mode of instability, instead, it only affects the growth rate of the unstable modes and thus does not play a major role in determining the stability of the flow.
Please note that in the present study we have neglected the existence of the PCL, which may lead to some modifications in the predictions obtained in the present study regarding the solid elastic mode due to its proximity to the deformable muscle layers. However, it will not affect the liquid elastic mode which originates at the air–liquid interface. Additionally, the solid elastic mode exists even for a Newtonian fluid layer subjected to airflow induced shear, thus the solid elastic mode will exist even in the presence of the PCL. Thus, the PCL layer can affect the growth rate predicted for the solid elastic and resonance modes but may not affect the liquid elastic mode.
4.3.2. Non-zero $Re$
Since the resonance mode is a result of the resonance between the liquid and solid elastic modes, the inertia is expected to lead to the stabilisation. The effect of increasing inertia on the resonance mode as illustrated in figure 13 is indeed stabilising, which further affirms its purely elastic nature.
4.4. Effect of pulmonary surfactant
The surface tension at the mucus–air interface in the mucus could be lowered due to the presence of the pulmonary surfactant. The surfactant is secreted by the alveolar type II cells (Grotberg Reference Grotberg2001; Grotberg & Jensen Reference Grotberg and Jensen2004), which lower the air–mucus interfacial tension by introducing Marangoni stresses. In this section, we try to understand the effect of the pulmonary surfactant on the stability of the flow. The destabilising effect of increasing 
$Ma$ on the resonance mode is illustrated in figure 14. It must be noted that an increasing Marangoni stresses at the air–mucus interface leads to the corresponding decrease in the surface tension 
$T$ (Grotberg & Jensen Reference Grotberg and Jensen2004), which in turn fuels the growth of the perturbations. Thus, the presence of the pulmonary surfactant in the airways leads to an increase in the severity of the unstable perturbations.
5. Conclusions
Saliva and mucus are known to have a strong shear-thinning viscoelastic nature. Furthermore, these viscoelastic liquids are supported by a deformable-solid layer such as the tongue, mucosa and submucosa layers. Thus, a consistent fluid mechanical model for this three-layer system needs to consider the dynamics of the turbulent air, viscoelastic saliva or mucus and deformable muscles. This three-layer system is modelled here as a shear-thinning viscoelastic liquid layer sheared by the air and flowing past a deformable solid. The dynamics of the liquid and solid is described using the White & Metzner (Reference White and Metzner1963) and neo-Hookean models, respectively.
To understand the role of the elasticity of the liquid and solid, the stability analysis has been divided into four cases as follows. The first case is the Newtonian liquid (
$n=1$ and 
$W\rightarrow 0$) sheared by the air and flowing past a rigid solid (NR) known to be unstable for Reynolds numbers 
$Re>34.2$ in the absence of the air–liquid interfacial tension (Smith & Davis Reference Smith and Davis1982). The second case is the shear-thinning viscoelastic liquid sheared by the air and flowing past a rigid solid (VR) which highlights the shear thinning and elasticity of the liquid. The stability analysis carried out here reveals that the VR case is unconditionally unstable, i.e. does not possess a critical Weissenberg number 
$W=\lambda _{sc} V_m/R$ or 
$Re$ even in the presence of the air–liquid interfacial tension 
$T$. The elastic mode arising due to the elasticity of the liquid for the VR case is termed the ‘liquid elastic mode.’ The growth rate of the liquid elastic mode increases with increasing 
$W$ and decreasing 
$T$ and 
$n$ (power-law index). The destabilisation mechanism for the liquid elastic mode is found to be the first normal stress difference across the air–liquid interface, a unique feature of viscoelastic liquids.
The third case is the Newtonian liquid (
$n=1$ and 
$W\rightarrow 0$) sheared by the air and flowing past a deformable solid (ND) which sheds light on the elasticity of the solid. The stability analysis carried out here reveals that the ND case is also unconditionally unstable, thereby showing that the elastic characteristics of the liquid and deformable solid give rise to the unconditionally unstable modes. The elastic mode caused by the deformable solid is termed the ‘solid elastic mode.’ The growth rate of the solid elastic mode increases with increasing dimensionless speed of the air 
$\varGamma$ and decreasing 
$T$. The mechanism responsible for the destabilisation of the solid elastic mode is the shear stresses exerted by the air on the liquid at the air–liquid interface and by the liquid on the deformable solid at the liquid–deformable-solid interface.
The fourth case considers a viscoelastic liquid sheared by the air and flowing past a deformable solid layer (VD) which shows the combined effect of the liquid and solid elasticity. The analysis shows that the liquid and solid elastic modes undergo resonance to result in the ‘resonance mode’ of instability. The resonance mode exhibits a much higher growth rate than the liquid and solid elastic modes and exhibits characteristics of both liquid and solid elastic modes. Additionally, the effect of the pulmonary surfactant on the growth rate of the resonance mode has been also investigated. The analysis predicts a growth rate enhancing effect of Marangoni stresses arising due to the pulmonary surfactant due to lowering the air–mucus interfacial tension. To conclude, the present analysis shows that the shearing action of the air induces instabilities in the airways and oral area. These instabilities originate as a result of the elastic nature of the fluids lining the oral and airways organs such as saliva and mucus and the muscle layers supporting the liquid layers.
Furthermore, we hypothesise that the predicted instability may lead to the fluid particle generation in the airways and oral area. The fluid particle size exhaled during breathing, talking and coughing is in the range of 
$10^{-5}\text {--}10^{-7}$ m with a preponderance of particles with diameters less than 
$10^{-6}$ m (Papineni & Rosenthal Reference Papineni and Rosenthal1997). From figure 12(b), the resonance mode achieves maximum growth rate for 
$k_m\sim 0.2\text {--}1$ depending on the air–liquid interfacial tension 
$T$. This implies that the dimensional wavelength of the most unstable disturbances is 
$\lambda _m \sim 10 {\rm \pi}R - 2 {\rm \pi}R$. From table 1, 
$R\sim 10^{-5}\text {--}10^{-8}$ m which yields 
$\lambda _m \sim 10^{-4}\text {--}10^{-7}$ m, which is in good correlation with the exhaled fluid particle size distribution measured by Papineni & Rosenthal (Reference Papineni and Rosenthal1997).
To further understand the fluid particle formation due to the elastic instabilities, it would be useful to carry out a nonlinear stability analysis which could shed light on the process through which the elastic instabilities lead to the formation of the fluid particles. The exhaled fluid particles could also be a result of the turbulent bursts caused during coughing and sneezing, thus the present analysis needs to be extended to include such turbulent bursts. Considering the thin layer of the mucus and saliva, a thin-film analysis could be carried out which could give a better idea of the thin-film rupture that results in fluid particle formation.
Figure 13. The effect of increasing inertia in the growth rate on the resonance mode at 
$\varGamma =0.2$, 
$n=1$, 
$\beta =0$, 
$H=10$, 
$W=0.2$, 
$Ma=0$ and 
$T=1$. Similar to the liquid and solid elastic modes, the resonance mode is stabilised by increasing inertia, implying the purely elastic nature of the resonance mode.
Figure 14. The destabilising effect of the pulmonary surfactant on the growth rate of the resonance mode at 
$\varGamma =0.2$, 
$n=0.5$, 
$\beta =0$, 
$H=10$, 
$W=0.2$ and 
$T=1$. The lowering of the air–mucus interfacial tension due to Marangoni stresses results in the growth rate enhancement of the resonance mode.
Acknowledgements
The author would like to thank the anonymous reviewers for their careful reading of the manuscript and their many insightful comments and suggestions which helped in improving the manuscript.
Declaration of interest
The author reports no conflict of interest.
