The effect of viscoelasticity in an airway closure model

Abstract

The closure of a human lung airway is modelled as a pipe coated internally with a liquid that takes into account the viscoelastic properties of mucus. For a thick-enough coating, the Plateau–Rayleigh instability blocks the airway by the creation of a liquid plug, and the preclosure phase is dominated by the Newtonian behaviour of the liquid. Our previous study with a Newtonian-liquid model demonstrated that the bifrontal plug growth consequent to airway closure induces a high level of stress and stress gradients on the airway wall, which is large enough to damage the epithelial cells, causing sublethal or lethal responses. In this study, we explore the effect of the viscoelastic properties of mucus by means of the Oldroyd-B and FENE-CR model. Viscoelasticity is shown to be very relevant in the postcoalescence process, introducing a second peak of the wall shear stresses. This second peak is related to an elastic instability due to the presence of the polymeric extra stresses. For high-enough Weissenberg and Laplace numbers, this second shear stress peak is as severe as the first one. Consequently, a second lethal or sublethal response of the epithelial cells is induced.

1. Introduction

Bronchi and bronchioles are human airways which depart from the windpipe as a network of tubular branches contained in the lungs. Each bifurcation or generation of the airway implies a reduction of the tubular cross-section, reducing the airway diameter down to microscopic scales at which the bronchioles connect to microscopic sacs of air termed alveoli. In total, an adult lung consists of approximately 23 generations with the trachea being the zeroth generation and the terminal bronchioles the last one. This bifurcation process is well described, for the first 15 generations, by a power law which characterizes the airway radius: $a_n=a_0 2^{-n/3}$, where $a_0$ is the trachea's radius (for adults it normally ranges in $a_0\in [0.8, 1]$ cm) and $a_n$ the radius of the airway at generation $n$ (Weibel & Gomez 1962).

On the inside, the airways are lined with an annular liquid film whose thickness is normally between 2 % and 4 % of the airway diameter. In the first 15–16 generations, the liquid film is a bilayer consisting of mucus on the inner side and serous, or periciliary, on the outer side (Widdicombe et al. 1997), whereas from the 17th generation on, the liquid film becomes a single-layered waterish fluid. It is well known that the mucus is a non-Newtonian fluid which exhibits viscoelastic properties (Yeates, Crystal & West 1990), possesses a yield stress and shows shear thinning characteristics (Basser, McMahon & Griffith 1989; Quraishi, Jones & Mason 1998). In some pathological conditions, the liquid film thickness increases substantially and may even reach 40 % of the airway radius (Sackner & Kim 1987). This abnormal film thickness exceeds the critical threshold established by Gauglitz & Radke (1988), after which an infinitesimally small varicose perturbation of the liquid–air interface leads to a Plateau–Rayleigh instability forming a liquid plug which occludes the airway. This phenomenon, known as airway closure, may lead to the collapse of the airway (Macklem, Proctor & Hogg 1970; Greaves, Hildebrandt & Hoppin 1986), hence to the blockage of gas exchange between trachea and distal airways.

For these reasons, modelling airway closure must take into account pathologies which induce accumulation of liquid due to infections or edema, hypersecretion of mucus and surfactant deficiency or dysfunction. Typical lung diseases are pneumonia (Gunther et al. 1996), asthma (Veen et al. 2000), bronchiolitis (Dargaville, South & Mcdougall 1996), chronic obstructive pulmonary disease (known as COPD; Guerin et al. (1997)), cystic fibrosis (Griese et al. 2004) and acute respiratory distress syndrome (known as ARDS; Viola et al. (2019)). For literature reviews of respiratory airway closure, liquid plug propagation and rupture, we refer to Grotberg (1994, 2001, 2011). Moreover, Grotberg (2019) has recently pointed out that mechanical stresses and strains induced by airway closure and reopening can cause, themselves, lung disease or injury. Further investigation along this line are considered in our study.

When a pathology induces mucus hypersecretion, the importance of non-Newtonian behaviours of the bilayer liquid lining the airway increases and viscoelasticity, yield stress and shear-dependent viscosity may affect the whole closure phenomenon (Halpern, Fujioka & Grotberg 2010), or prevent reopening once an airway is completely occluded by a liquid plug. This has important implications on the respiratory system and might even have lethal consequences (Synek et al. 1994).

Owing to the importance of airway closure, several investigations have tried to identify and study the various causes which might induce the occlusion of respiratory ways, including the capillary instability of the lining liquid film and the mechanical instability of the tube walls (see e.g. Halpern & Grotberg 1992). The compliance of airway walls has been extensively studied by Heil, Hazel & Smith (2008), who took into account cross-sectional instabilities of the respiratory ways and demonstrated the importance of non-axisymmetric mechanical and hydrodynamic instabilities, which act in favour of the airway occlusion. Other investigations focused on the role of surfactants in the process, and, correspondingly, on the negative impact caused by surfactant deficiencies on the airway closure. Experimental and numerical studies on this topic have been carried out by Cassidy (1999) and Halpern et al. (2008), respectively. Thereafter, Halpern et al. (2010) investigated the effect of viscoelasticity on the airway closure, by modelling the bilayer thin film using the Oldroyd-B model. They focused on the effect of the Weissenberg number on the closure time and the maximum shear stress. Several of the aforementioned studies employ a thin-film approximation valid under the assumptions of the lubrication theory. Even though such assumptions are valid prior to coalescence, when the liquid film is thin compared with the airway radius, the lubrication theory provides a very good estimate of the precoalescence phase and of the closure time. However, several important insights in the airway closure phenomenon can still be disclosed by reproducing the whole closure phenomenon, including the formation of a stationary plug by bifrontal plug growth that has been recently shown to be potentially a major cause of epithelial cell damage in the entire air closure process (Romanò et al. 2019).

When the airway closure occurs, a liquid plug of mucus and serous occludes the respiratory way and gets advected by inhaling air during the respiratory phase. Propagation of the liquid plug leads to rupture of liquid meniscus if the leading edge film is thinner than the trailing edge one (airway reopening, see e.g. Hassan et al. (2011) and Muradoglu et al. (2019)), or if the plug reaches an airway bifurcation. This normally causes a pressure wave which can be auscultated by a stethoscope (Piirila & Sovijarvi 1995). The airway reopening has been extensively investigated due to the damaging effects of the very high stress levels at the airway wall, which are induced by the plug propagation and the successive plug rupture. As pointed out by Bilek, Dee & Gaver (2003), Kay et al. (2004) and Huh et al. (2007), the mechanical stress on the epithelial cells due to high wall shear stress, and axial wall stress gradients may cause serious injuries to the epithelial cells, and even leading to lethal consequences. The connection between plug rupture and high wall stress levels has been numerically confirmed by Fujioka & Grotberg (2004), Fujioka & Grotberg (2005), Fujioka, Takayama & Grotberg (2008) and Muradoglu et al. (2019), who investigated the airway reopening in a rigid pipe and a rigid channel. Further confirmations have been provided by the in vivo experiments by Muscedere et al. (1994) and Taskar et al. (1997), who considered animal models and excised lungs to prove that a sequence of plug ruptures severely damages the airway tissues. Additional considerations are reserved to the effect of non-uniform surface tension and yield stress on plug rupture (Ghadiali & Gaver 2008; Muradoglu et al. 2019; Hu, Romanò & Grotberg 2020), and compliance of airway walls, which enhance the level of stress on the epithelial cells (Zheng et al. 2009).

Even though a rich literature can be found regarding airway reopening, not many studies have focused on airway closure and the stress levels associated with it. Recently, our experimental and numerical studies (Bian et al. 2010; Tai et al. 2011; Romanò et al. 2019) have addressed the airway closure as a potentially damaging phenomenon for the airways, showing that the level of stress due to the formation of a liquid plug in rigid microscopic pipes may be comparable to the injuring threshold indicated by Bilek et al. (2003) and Huh et al. (2007) for the mechanical stress on epithelial cells. In particular, the recent computational study by Romanò et al. (2019) has clearly demonstrated that the highest peak of the wall stresses is reached a few instants after the coalescence. They have thus showed that, although it is often ignored in many theoretical and computational studies, the postcoalescence phase is the most important phase in the entire airway closure process since it may induce from 300 % to 600 % larger mechanical stresses compared with the precoalescence values.

In our previous study (Romanò et al. 2019), we modelled the annular film coating the airway by means of a Newtonian fluid. The density and dynamic viscosity were assumed to be the weighted mean of the physical properties of mucus and periciliary liquid that form the multilayer liquid coating the airways. In Romanò et al. (2019), Tai et al. (2011) and Bian et al. (2010) we demonstrate that in a first phase of the instability, the radial velocity dominates forming a liquid bulge that grows because of the Plateau–Rayleigh instability mechanism up to forming a liquid plug that occludes the airway. Immediately after coalescence, the velocity profile inside the liquid phase is dominated by the axial velocity, resembling two receding air fingers. The sharp topological change leading to the liquid plug, as well as the growth of the symmetric liquid plug, deserve the highest attention in a Newtonian process since they induce the highest wall stresses. To stress such a peculiar quick phase of the airway closure, Romanò et al. (2019) termed it bifrontal plug growth.

The effect of initial film thickness, liquid dynamic viscosity and surface tension has been investigated by Romanò et al. (2019). Upon an increase of the initial film thickness $\varepsilon$, the liquid plug formation speeds up (see also Bian et al. 2010; Tai et al. 2011). No remarkable changes with $\varepsilon$ have been observed in terms of tangential and normal wall stresses, and shear stress gradients, whereas Romanò et al. (2019) demonstrated for the first time that the initial film thickness strongly influences the wall pressure gradient. In our previous study, we also showed that the dependence of the closure time $t_{c}$ on the surface tension $\sigma$ is slightly sublinear, as well as the dependence of $t_{c}$ on $\mu _{L}^{-2}$, where $\mu _{L}$ denotes the dynamic viscosity of the liquid phase. In terms of non-dimensional wall stresses, the variations observed by changing $\sigma$ and $\mu _{L}$ are not very significant.

The simulations by Romanò et al. (2019) model the liquid film as a Newtonian fluid. However, since the mucus is a non-Newtonian fluid with viscoelastic and viscoplastic properties, the present study aims to refine the physical model of Romanò et al. (2019) by including the effects of viscoelasticity on the airway closure flow. Major qualitative and quantitative differences are expected between the viscoelastic airway closure and the Newtonian one, hence comparing the viscoelastic model with the Newtonian model will help us in recognizing the net effect of viscoelasticity.

A particular emphasis of this study is placed on the effects of viscoelasticity on the mechanical stresses over the airway wall where the epithelial cells are located. For this purpose, the liquid film is modelled as a viscoelastic fluid using the Oldroyd-B model. A similar study has been carried out by Halpern et al. (2010), but they employed a lubrication approximation for investigating only the precoalescence dynamics. In this paper we simulate the whole process, including precoalescence and postcoalescence phases for biologically relevant conditions. Simulations are also performed using the FENE-CR model to examine the effects of polymer extensibility.

The rest of the paper is organized as follows. The problem formulation and the corresponding numerical discretization are reported in §§ 2 and 3, respectively. Section 4 presents the simulation results and compares the Newtonian airway closure of Romanò et al. (2019) with the non-Newtonian numerical results of this study, highlighting the phases when viscoelastic behaviours dominate the wall stresses. Finally, § 5 summarizes our results and draws the conclusions of our study, while the appendix A reports a validation of the numerical solver and gathers the considerations in terms of the polymer extensibility parameter used in the FENE-CR model.

2. Problem formulation

We model the airway closure employing a cylindrical rigid tube of radius $a$ and length $L$, lined with a viscoelastic non-Newtonian liquid film of initial average thickness $h$, total dynamic viscosity $\mu _{L}$ and density $\rho _{L}$. Inside the rigid pipe, the core gas has constant dynamic viscosity $\mu _{G}$ and density $\rho _{G}$ and is surrounded by the thin viscoelastic film as depicted in figure 1. The surface tension $\sigma$, acting at the interface between the liquid phase and the gas phase, is assumed to be constant.

Figure 1. Schematic of the geometry of the airway model: the rigid tube has radius $a$, length $L$ and is coated by a liquid film (light blue) of average thickness $h$ and surrounded by a gas core. The interface is initially located at a distance $R_{I}$ from the axis of the pipe.

The initial radial location of the interface is subject to a varicose perturbation with a first Fourier mode. The amplitude of the perturbation is 10 %, hence the radius of the interface is initialized as

(2.1)\begin{equation} r = R_{I} = a - h [1-0.1 \times \cos(2 {\rm \pi}z /L)], \end{equation}

where $z$ and $r$ denote the axial and radial coordinates, respectively. To non-dimensionalize the governing equations, we make use of a capillary scaling, i.e. length, time, pressure and velocity are scaled with $a$, $\mu _{L}a/\sigma$, $\sigma /a$, $\sigma /\mu _{L}$, respectively. Assuming that the flow is incompressible, the resulting single-field dimensionless equations (Tryggvason, Scardovelli & Zaleski 2011; Popinet 2018) read

(2.2a)\begin{gather} \partial_t \tilde{\varrho}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\tilde{\varrho}=0, \end{gather}

(2.2b)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}

(2.2c)\begin{gather}{La} \tilde{\varrho}\left(\partial_t \boldsymbol{u} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}\right)={-} \boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot} \tilde{\boldsymbol{\tau}} + \chi \boldsymbol{n} \delta_{s}, \end{gather}

where $p$ denotes the pressure, $\boldsymbol {u}=(u_r,u_\phi ,u_z)$ is the velocity vector, $\chi =\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {n}$ is the total local curvature of the interface, $\boldsymbol {n}$ refers to the outward unit normal at the interface. The surface Dirac $\delta$-function $\delta _{s}$ equals zero everywhere, but at the two-phase interface, $\tilde {\varrho }$ is the variable density field required by the single-field approach to include the effect of the gas-to-liquid density ratio. In the liquid phase $\tilde {\varrho }=1$, while $\tilde {\varrho }=\varrho$ in the gas phase. In (2.2), $\tilde {\boldsymbol {\tau }}$ denotes the non-dimensional deviatoric component of the stress tensor. In the gas phase $\tilde {\boldsymbol {\tau }}=\boldsymbol {\tau }_{G}=\mu (\boldsymbol {\nabla } \boldsymbol {u} + \boldsymbol {\nabla }^T\boldsymbol {u})$, where $\boldsymbol {\tau }_{G}$ is the deviatoric part of the Newtonian stress tensor with $\mu$ being the gas-to-liquid dynamic viscosity ratio. In the liquid phase, $\tilde {\boldsymbol {\tau }}$ consists of the Newtonian and the viscoelastic extra stresses. Following Halpern et al. (2010), we model $\tilde {\boldsymbol {\tau }}$ in the liquid using the Oldroyd-B model in order to investigate the effect of the viscoelastic properties of the liquid film. This non-Newtonian model assumes that the liquid film behaves as a diluted polymer (solute) immersed in a Newtonian liquid (solvent). The polymer is normally modelled as a linear elastic dumb-bell and, based on kinetic theory, a constitutive equation for the polymeric stresses is derived. Consistently, $\tilde {\boldsymbol {\tau }}$ is assumed to be a superposition of the solvent and the solute deviatoric stresses as follows:

(2.3)\begin{equation} \tilde{\boldsymbol{\tau}} = \tilde{\boldsymbol{\tau}}_{L} = \mu_{S}\left(\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla}^T\boldsymbol{u}\right) + \boldsymbol{S}, \end{equation}

where $\tilde {\boldsymbol {\tau }}_{L}$ is the deviatoric component of the stress tensor in the liquid phase, $\mu _{S}$ is the solvent-to-total dynamic viscosity ratio in the liquid film and $\boldsymbol {S}$ is the deviatoric stress tensor due to the polymeric part of the viscoelastic fluid. The polymeric tensor $\boldsymbol {S}$ is governed by the constitutive equation

(2.4)\begin{equation} We \left[\partial_t\boldsymbol{S} + \left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \right) \boldsymbol{S} - \left(\boldsymbol{\nabla}\boldsymbol{u} \right)\boldsymbol{S} - \boldsymbol{S}\left(\boldsymbol{\nabla}^{T}\boldsymbol{u} \right)\right] + \boldsymbol{S} = \mu_{P}\left(\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla}^T\boldsymbol{u}\right), \end{equation}

where $\mu _{P}=1-\mu _{S}$ is the solute-to-total dynamic viscosity ratio in the liquid film. Several independent non-dimensional groups arise from the momentum equation: the Laplace number $La$, the Weissenberg number $We$, the gas-to-liquid density ratio $\varrho$ and the gas-to-liquid dynamic viscosity ratio $\mu$, and the solvent-to-total dynamic viscosity ratio in the liquid film $\mu _{S}$. Besides, two additional aspect ratios are required to completely define the problem parameters: the length-to-radius aspect ratio $\lambda$ and the normalized average film thickness $\varepsilon$. The non-dimensional parameters can be summarized as follows:

(2.5a–g)\begin{equation} \left.\begin{gathered} La = \frac{\rho_{L}\sigma a}{\mu_{L}^2}, \quad We = \frac{\varLambda\sigma}{a\mu_{L}}, \quad\varrho = \frac{\rho_{G}}{\rho_{L}}, \quad \mu = \frac{\mu_{G}}{\mu_{L}}, \quad \mu_{S} = \frac{\mu_{L,S}}{\mu_{L}},\\ \lambda = \frac{L}{a}, \quad \varepsilon = \frac{h}{a}, \end{gathered}\right\} \end{equation}

where $\varLambda$ is a characteristic relaxation time, $\mu _{L,S}$ is the solvent dynamic viscosity and $\mu _{L}=\mu _{L,S}+\mu _{L,P}$ is the total viscosity with $\mu _{L,P}$ being the polymeric viscosity.

Since the typical Weissenberg numbers involved in pulmonary flows are very high, the strong elasticity of the non-Newtonian fluid makes the constitutive equations very stiff and may produce regions of high stresses and fine flow structures which are known to induce numerical instabilities. To overcome this difficulty, we first recast (2.4) by introducing the conformation tensor defined as $\boldsymbol {A}=\boldsymbol {S}+\mu _{P}We^{-1}\boldsymbol {I}$, where $\boldsymbol {I}$ is the identity matrix. Then following Fattal & Kupferman (2004) and Fattal & Kupferman (2005), and introducing $\boldsymbol {\varPsi }=\log (\boldsymbol {A})$, the log-conformation representation of the viscoelastic constitutive equation can be written as

(2.6)\begin{equation} We \left[\partial_t\boldsymbol{\varPsi} + \left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} \right)\boldsymbol{\varPsi} - \left(\boldsymbol{\varOmega} \boldsymbol{\varPsi} - \boldsymbol{\varPsi}\boldsymbol{\varOmega}\right) - 2\boldsymbol{B}\right] = \left(\textrm{e}^{-\boldsymbol{\varPsi}}-\boldsymbol{I}\right), \end{equation}

where $\boldsymbol {\varOmega }$, $\boldsymbol {N}$ and $\boldsymbol {B}$ are two antisymmetric and one symmetric tensors, respectively, resulting from the decomposition $\boldsymbol {\nabla }\boldsymbol {u}=\boldsymbol {\varOmega } + \boldsymbol {B} + \boldsymbol {N}\boldsymbol {A}^{-1}$. Equation (2.6) is solved numerically together with the flow equations and the conformation tensor is then obtained using the inverse transformation as $\boldsymbol {A} = \textrm {e}^{\boldsymbol {\varPsi }}$.

The mathematical problem (2.2) is closed by enforcing periodic boundary conditions in axial direction, and no-slip and no-penetration conditions along the wall

(2.7a,b)\begin{equation} \boldsymbol{u}\left(z=0 \right) = \boldsymbol{u}\left(z=\lambda\right), \quad \boldsymbol{u}\left(r=1 \right) = \boldsymbol{0}. \end{equation}

3. Numerical method

The airway closure is modelled as an axisymmetric phenomenon, enforcing $\partial _\phi =0$, $u_\phi =0$. The governing equations (2.2) are integrated in time by means of a fractional step method consisting of a second-order pressure-correction projection scheme. The viscous term is treated implicitly, whereas the nonlinear convective term is discretized explicitly by means of the Bell–Collela–Glaz advection scheme (Bell, Colella & Glaz 1989). The spatial discretization is carried out employing a second-order finite volume method on a staggered grid which solves for the pressure at the cell centre and for the velocity at the cell faces.

The two-phase flow is dealt with using the volume of fluid (VOF) method, which captures the liquid–gas interface by means of a fraction field $f(r,z,t)$. The fraction field represents the volume fraction of liquid in a computational cell, i.e. it equals 1 if the cell contains only liquid, 0 if the cell is filled only with gas and $f(r,z,t) \in (0,1)$ if the cell is filled partially by liquid and partially by gas. The condition $f(r,z,t) \in (0,1)$ identifies the computational cells which contain the liquid–gas interface. Building on the fraction field, the variable fluid properties in the single-field approach are expressed as

(3.1a)\begin{gather} \tilde{\varrho}\left[f\left(r,z,t\right)\right] = f + (1-f)\varrho, \end{gather}

(3.1b)\begin{gather}\tilde{\boldsymbol{\tau}}\left[f\left(r,z,t\right)\right] = f{\boldsymbol{\tau}}_{L} + (1-f){\boldsymbol{\tau}}_{G}, \end{gather}

where the physical properties of the interface are dealt with by the same approach of Popinet (2003) and Popinet (2008). An additional equation is required to complete the description of the fraction field dynamics, modelled, for immiscible fluids, as simple advection of $f$ due to the flow velocity $\boldsymbol {u}$,

(3.2)\begin{equation} \partial_t f + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\,f \boldsymbol{u}\right)=0. \end{equation}

Equation (3.2) is discretized by employing a staggered approach in time making use of a second-order scheme as done by Popinet (2008). The liquid–gas interface is discretized using a piecewise-linear geometrical VOF method, which solves (3.1) representing the interface as a straight line within each computational cell. The interface is then advected taking into account its local unit normal (mixed young centred method, see Aulisa, Manservisi & Scardovelli (2006)) and the local indicator function $f$. Consistently with the projection method, (2.6) is also discretized by the Bell–Collela–Glaz method (Bell et al. 1989).

The last discretization step is operated on the surface-tension term on the right-hand side of (2.2c). Numerical issues are normally associated with the discretization of $\chi \boldsymbol {n}\delta _{s}$, which might induce parasitic currents (see Brackbill, Kothe & Zemach 1992; Popinet & Zaleski 1999) either in front-tracking or in front-capturing methods. Combining a balanced-force approach and a height-function estimator, Popinet (2008) showed that numerical parasitic currents can be avoided, and second-order convergence rates are achieved even for challenging benchmarks. In the followings, the approach of Popinet (2008) is adopted.

All the simulations presented hereinafter are computed on a Cartesian grid designed for capturing the wall stresses and their gradients for biologically relevant parameters in airway closure problems. Romanò et al. (2019) performed a grid convergence study for the Newtonian airway closure model and showed that grid-independent solutions are obtained using a grid resolution containing $512\times 86$ cells in the axial and radial directions, respectively. Although not shown here due to space considerations, a similar grid convergence study has been performed and the same grid resolution is found to be sufficient to obtain grid-independent solutions for all the simulations reported in the present paper. The numerical implementation of all the discretization schemes in space and time is implemented in the free software package called Basilisk (Popinet 2014; http://basilisk.fr) which has been extensively tested for a large number of benchmark cases and successfully used for a wide range of multiphase flows of practical interest (Popinet 2014). In addition, Romanò et al. (2019) have validated the method against the experimental measurements of Bian et al. (2010) and the numerical results obtained using the front-tracking method of Muradoglu et al. (2019). The log-conformation approach applied to the Oldroyd-B model is provided as an embedded library in Basilisk and has been validated against the results of Figueiredo et al. (2016) for a viscoelastic axisymmetric droplet impacting and spreading on a plane wall. An additional validation is reported in the appendices, comparing the results of Basilisk with the corresponding simulations carried out with the code of Izbassarov & Muradoglu (2015). As discussed in the appendices, the results obtained using Basilisk and the front-tracking method are in very good agreement providing a further verification for the accuracy of the present simulations.

4. Results and discussion

The effect of the viscoelastic properties of the liquid film is studied by varying the two independent non-dimensional groups characterizing the viscoelasticity, namely the Weissenberg number $We$ and the concentration of solvent which controls the relative dynamic viscosities $\mu _{S}$ and $\mu _{p}=1-\mu _{S}$. Based on the works of Guo & Kanso (2017), Ross (1971), Smith, Gaffney & Blake (2006) and Mitran (2007), the viscoelastic relaxation time of mucus varies in the range of $\varLambda _{M} \in [0.1, 10]$ s in healthy conditions while Lauga (2007) and Gilboa & Silberberg (1976) suggest that $\varLambda _{M}$ normally grows in diseased conditions (see, e.g. Hwang, Litt & Forsman (1969), Gilboa & Silberberg (1976), who measured the relation time in the range $\varLambda _{M} \in [30, 100]$ s). Hence, the range of interest for the relaxation time of the mucus layer is $\varLambda _{M} \in [0.1, 100]$ s. The value of the solver-to-total viscosity ratio $\mu _{S}$ is usually assumed as an investigation parameter (see, e.g. Halpern et al. 2010) and we consider it in the range $\mu _{S}\in [0.25, 0.9]$. The serous layer is usually assumed to be Newtonian, but Guo & Kanso (2017), Tarran et al. (2001), Button et al. (2012) and Boucher (2004) pointed out that its weak viscoelasticity may play a non-negligible role. In the present study, we use a homogenization approach in which the liquid layer is assumed to be a homogeneously viscoelastic fluid ignoring details of the highly viscoelastic mucus and the weakly viscoelastic serous layers, and consider a single Weissenberg number in the range $We\in [5, 1000]$.

The choice of the other non-dimensional groups refers to biologically relevant parameters for the Plateau–Rayleigh instability causing the airway closure at the ninth-to-10th branching generations in the lungs. Considering an adult lung airway, the bronchioles’ radius at ninth-to-10th generation is approximately $a \in [0.05, 0.1]$ cm (see e.g. Crystal 1997) and a normal length-to-radius ratio of the airway is $\lambda =6$. These parameters are used to quantify the wall stresses in dimensional terms. Following Romanò et al. (2019), the density of the multilayer liquid made out of mucus and serous is not far from that of water, i.e. $\rho _{L} = 1\ \textrm {g}\ \textrm {cm}^{-3}$, hence a representative gas-to-liquid density ratio is $\varrho = 10^{-3}$. Moreover, we assume that the serous layer is very watery hence its dynamic viscosity is similar to that of water, i.e. ${\approx }0.01$ poise. In addition, we consider that the dynamic viscosity of mucus may vary over several orders of magnitude, ranging from 10 to 10000 times that of water (see e.g. Lai et al. 2009) depending on physiological functions, pathological conditions and age. For these reasons Romanò et al. (2019) suggested to deal with the bilayer liquid film as a homogeneous medium with dynamic viscosity in the range $\mu _{L}\in [0.126, 1.26]$ poise. Thereafter, considering that the dynamic viscosity of air at $37\,^\circ \textrm {C}$ is $\mu _{G} = 1.89 \times 10^{-4}$ poise, these estimates lead to $\mu \in [0.15, 1.5]\times 10^{-3}$. Romanò et al. (2019) showed that a change of the dynamic viscosity ratio over the two-orders-of-magnitude interval $\mu \in [0.15, 1.5]\times 10^{-3}$ has a weak effect on the wall stresses and increasing the liquid viscosity plays a significant role only in slowing down the airway closure process. Hence, since our study focuses on the effect of viscoelasticity on the wall stresses, we do not vary the gas-to-liquid dynamic viscosity ratio, and keep it fixed at $1.5 \times 10^{-4}$ for all the results presented in this paper. On the other hand, an important role in the bifrontal plug growth is played by the Laplace number (Romanò et al. 2019). The same range of $La$ used by Romanò et al. (2019) is also employed in the present study, i.e. $La\in \lbrace 2, 20, 200\rbrace$. Finally, the initial non-dimensional average film thickness is set to $\varepsilon =0.25$ for all the simulations presented in this paper.

4.1. Comparison of the Newtonian and viscoelastic liquid airway models

The comparison between the viscoelastic (solid lines) and the Newtonian (dashed lines) airway closure models is presented in figure 2 in terms of the excursion of the wall shear stress as a function of time, $\Delta \tau _{w} (t) = \max _z \tau (r=1) - \min _z \tau (r=1)$. All the results in figure 2 are obtained for $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$ and three different Laplace numbers, i.e. $La = 200$ (black), 20 (red) and 2 (blue). For the Oldroyd-B model, we use $\mu _{S}=0.5$ and $We=100$.

Figure 2. Excursion of the tangential stress distribution along the wall, $\Delta \tau _{w} (t) = \max _z \tau (r=1) - \min _z \tau (r=1)$ for $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$. Three Laplace numbers are considered: $La = 200$ (black), 20 (red) and 2 (blue). Two constitutive models for the liquid phase, i.e. the Newtonian (dashed lines) and the Oldroyd-B model (solid lines, $\mu _{S}=0.5$ and $We=100$), are compared.

For both the constitutive models, an increase of the Laplace number increases the peak of the shear stress level, that rises from $\max (\Delta \tau _{w})\approx 0.4$ for $La = 2$ to $\max (\Delta \tau _{w})\approx 0.6$ for $La = 200$. As will be further discussed when presenting the results for different Weissenberg numbers, the peak of the shear stress level at the wall is almost insensitive to the viscoelastic properties of the fluid. This is well understood considering the fact that the initial peak observed in the tangential stress results from the immediate reaction of the fluid to the bifrontal plug growth. Hence, the sharp peak of the curves is due to the Newtonian part of the liquid and extra stresses (whenever present, i.e. for $We\neq 0$), play a very minor role right after the coalescence. A more detailed analysis of the wall stresses for the Oldroyd-B model is presented later, including the effect of $We$ and $\mu _{S}$ on the Newtonian part of the tangential stress, the viscoelastic extra stresses and on the pressure.

Another effect of the viscoelasticity on the airway closure is to reduce the closure time by approximately 20 %. This is consistent with the theoretical prediction of Halpern et al. (2010). However, the viscoelastic properties of mucus play their most significant role sometime after the coalescence event. Figure 2 depicts the qualitative difference between the Newtonian and viscoelastic cases: the tangential stresses in a Newtonian liquid relax to approximately $\Delta \tau _{w} \approx 0.1$ after the very quick transient due to the bifrontal plug growth. On the other hand, when viscoelastic effects are taken into account, the elastic reaction of the polymeric phase generates a second increase of the wall shear stress excursion, that can be of the same order of magnitude of the first Newtonian peak. This qualitatively and quantitatively very different postcoalescence dynamics has never been investigated before and certainly deserves a special attention. To better understand and characterize the viscoelastic effect, in the following sections we present a typical scenario of viscoelastic plug formation and we carry out a parametric study over $La$, $We$ and $\mu _{S}$.

4.2. Analysis of viscoelastic effects in a typical airway closure scenario

The surface tension between liquid and gas reduces the pressure proportionally to the cross-sectional curvature ($R_{I}^{-1}$). This is a destabilizing effect because a perturbed interface will have local pressure minima at the lowest $R_{I}$, with a consequent drainage of liquid from the near-wall regions (highest $R_{I}$) towards the bulge tip (lowest $R_{I}$) up to liquid plug formation or coating film pinch-off. On the other hand, the in-plane curvature ($\partial ^2_z R_{I}$) is a stabilizing effect that tends to smear out the interface perturbations because it creates pressure minima at the highest $R_{I}$, hence draining fluid from the bulge to the coating film. These two competing effects are at the basis of the Plateau–Rayleigh instability that may occlude the airways by formation of a liquid plug. Increasing the surface tension ($La$) favours the destabilizing effect leading to a quicker airway closure. Decreasing the liquid film thickness ($\varepsilon$) slows down the liquid plug formation as it becomes increasingly difficult to drain thinner and thinner films. Along the same line, also increasing the liquid viscosity slows down the closure process, as extensively proved in the literature (Bian et al. 2010; Tai et al. 2011; Dietze & Ruyer-Quil 2015; Romanò et al. 2019).

All these considerations can be inferred by studying a two-phase Newtonian system that does not include any viscoelastic effect in the liquid phase. Since we are interested in the viscoelastic effects, the formation of a liquid plug is here first investigated for $La = 2$, $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$, $\mu _{S} = 0.25$ and $We=500$. Under these physiologically relevant conditions, we expect to observe significant viscoelastic effects.

The distribution of pressure and in-plane extra-stress components, $p$ and $S_{rz}$, $S_{rr}$ and $S_{zz}$, is plotted in figure 3 at five time instants, two of them before and three after the airway closure, occurring at $t=t_c\in [199, 200]$. Owing to the symmetries of the problem, figure 3 depicts the four fields for only one quarter of the domain. The pressure distribution in the liquid is depicted at the top-left quadrant and is symmetric with respect to the vertical solid line in figure 3. Even though we consider strong viscoelastic effects, $p$ is qualitatively and quantitatively very similar to what observed for Newtonian fluids, as also reported in Halpern et al. (2010). For moderate deformations of the liquid film (see $t=150\approx 3 t_c/4$ in figure 3), the pressure is almost independent of the radial coordinate, consistently with the leading-order approximation of the lubrication theory for thin films (see Halpern et al. 2010). Owing to the bulge growth and the film drainage, the interface deformation increases, the liquid film configuration deviates further from the hypothesis of the thin-film approximation, and the pressure develops a radial dependence whose minimum is observed at the bulge tip (see $t=199 \lesssim t_c$ in figure 3). The quick bifrontal plug growth observed for $t > t_c$ in figure 3 leads to a capillary wave for the thin film near the wall, whose correlation with the wall pressure gradient has already been discussed in Romanò et al. (2019). The pressure $p$ increases by one order of magnitude within the quick topological change of the domain occupied by the liquid phase.

Figure 3. Time evolution of the plug formation by the Plateau–Rayleigh instability and bifrontal plug growth. The contour plots of $p$ (top left), $S_{rz}$ (top right), $S_{rr}$ (bottom left) and $S_{zz}$ (bottom right) are shown at two times before the closure $t = 150$ and $t = 199\lesssim t_c$ and at three times after the closure $t = 200 \gtrsim t_c$, $t = 360$ and $t = 450$. The simulation parameters are $La = 2$, $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$, $We = 500$ and $\mu _{S} = 0.25$.

The in-plane diagonal extra stresses $S_{rr}$ and $S_{zz}$ are depicted at the bottom-left and the bottom-right quadrant, respectively, and are symmetric with respect to the vertical solid line in figure 3. The extra stress $S_{rr}$ presents the highest values near the bulge for the whole preclosure phase. The highest extra stresses $S_{rr}$ remain localized near the axis also during the postcoalescence phase, concentrating more and more towards the symmetry line denoted by the solid black line. The order of magnitude of $S_{rr}$ does not change over the whole airway closure process. Almost the opposite trend is observed for $S_{zz}$, that increases by approximately one order of magnitude during the liquid-plug formation. Except for a quick dynamics right after the closure, figure 3 shows that the highest values of $S_{zz}$ are always localized near the airway wall, where the thin film connects to the bulge (for $t<t_c$) or to the liquid plug (for $t>t_c$). The in-plane extra-diagonal extra stress $S_{rz}$ is shown at the top-right quadrant and is antisymmetric with respect to the vertical solid line in figure 3. Also for $S_{rz}$ we observe an increase of one order of magnitude during the airway closure and, similarly to $S_{zz}$, the highest values of $S_{rz}$ are always localized near the airway wall. The top-right quadrants of figure 3 also demonstrate the correlation of $S_{rz}$ with the highest in-plane curvature of the interface. We anticipate that this feature is key to the understanding of the walls-stress dynamics.

Beside the characterization of pressure and extra-stress fields, the major focus of our study is on the effect of viscoelasticity on the wall stresses since it has direct implications for the epithelial cells living on the airway walls. Keeping the set of Newtonian parameters constant, i.e. $La = 2$, $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$, as well as the solvent-to-total viscosity ratio $\mu _{S} = 0.25$, we vary the relaxation time of the viscoelastic liquid, hence its Weissenberg number $We$. Figure 4(b) depicts the minimum radial location of the interface $R_{min}$ (solid line), and the excursion of the normal stress distribution along the wall, $\Delta p_{w} (t) = \max _z p(r=1) - \min _z p(r=1)$ (dashed line), for $We=\lbrace 5, 10, 50, 100, 500, 1000\rbrace$. Increasing the Weissenberg number, the airway closure speeds up significantly because the extra stresses play a destabilization role (see figure 3). The excursion in wall normal stresses, however, is only weakly affected by $We$, that reduces the pressure level of less than 10 %. This is consistent with the pressure field depicted at the top-left quadrant of figure 3, that agrees qualitatively and quantitatively with the pressure distribution for Newtonian liquids reported in Romanò et al. (2019). For these reasons, the following parametric study will rather focus only on the wall tangential stresses, that are significantly influenced by the viscoelastic properties of the liquid phase.

Figure 4. (a) Excursion of the tangential stress distribution along the wall, $\Delta \tau _{w} (t) = \max _z \tau (r=1)$ $- \min _z \tau (r=1)$ (solid line), consisting of the Newtonian ($\Delta \tau _{w}^{N}$, dashed line) and the extra-stress component ($\Delta S_{w}$, dashed–dotted line). (b) Excursion of the normal stress distribution along the wall, $\Delta p_{w} (t) = \max _z p(r=1) - \min _z p(r=1)$ (dashed line), and minimum of the radial location of the interface, $R_I$ (solid line). The Weissenberg number is varied between 5 and 1000, whereas the other simulation parameters are fixed: $La = 2$, $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$ and $\mu _{S} = 0.25$.

As observed when comparing the wall tangential stress for a Newtonian and a viscoelastic liquid film (see figure 2), the amplitude of the first peak of $\Delta \tau _{w}$, right after $t=t_c$, is a weak function of the viscoelastic non-dimensional groups $We$ and $\mu _{S}$. On the other hand, a remarkable qualitative difference between Newtonian and viscoelastic cases is observed in the long-term dynamics after the coalescence. To better understand it, figure 4(a) depicts the tangential stress excursion at the wall, $\Delta \tau _{w}$ (solid line), given by the sum of its two components, i.e. the Newtonian stress due to the solvent $\Delta \tau _{w}^{N}$ (dashed–dotted line) and the extra stress due to the polymers $\Delta S_{w}$ (dashed line). Figure 4 demonstrates that the first peak of the tangential stresses at the wall is due to the effect of the Newtonian dynamics that instantaneously reacts to the bifrontal plug growth. A significant viscoelastic contribution is observed for the highest Weissenberg numbers, long enough after the airway closure. The delay of the extra stresses increases with the Weissenberg number, that represents a non-dimensional measure of the relaxation time of the polymer. Moreover, for high-enough $We$, a peculiar dynamics is observed in $\Delta S_{w}$. The wall extra stress undergoes several oscillations whose amplitude is comparable to the first Newtonian peak of $\Delta \tau _{w}$. The greatest contribution to $\Delta S_{w}$ is given, in fact, by $\max _z S_{rz}(r=1) - \min _z S_{rz}(r=1)$. To understand the origin of such oscillations, we should consider the stretching of the polymeric chains along curved streamlines, where stretching is known to induce a purely elastic instability of type hoop stress (Pakdel & McKinley 1996; Groisman & Steinberg 1998). This non-Newtonian instability has been extensively investigated in some recent studies (see, e.g. Lin-Gibson et al. 2004; Steinberg & Groisman 1998; Graham 2003) which also focused on the definition of a parameter that can characterize the dominance of the elastic over the Newtonian instabilities (Groisman & Steinberg 2000). Recent experimental evidence demonstrated, however, that a combination of both inertial and elastic effects may significantly influence the onset of the elastic instability leading, de facto, to an elasto-inertial instability. In flows with curved streamlines such instability induces pronounced flow oscillations whose onset is governed by a modified elasticity number $El_{m}$. In fact, Casanellas et al. (2016) and Kim et al. (2017) showed that the non-dimensional group of interest that governs the elasto-inertial instability links the relaxation time $\varLambda$ to the viscous time scale $a^2/\nu$ and it scales like the square root of the polymer-to-solvent dynamic viscosity ratio. In our system, the interface curvature after closure induces the streamline curvature needed to stretch the polymeric chains that leads to the instability and, in our scaling, the modified elasticity number $El_{m}$ corresponds to $El_{m}=We\sqrt {1-\mu _{S}}/La=\varLambda \nu \sqrt {1-\mu _{S}}/a^2$. If $El_{m}$ is too large or too small, the elasto-inertial instability does not set in (see, e.g. blue symbols in figure 6 of Kim et al. (2017)) and purely elastic or inertial instabilities may still occur. On the other hand, if $El_{m}$ falls in the unstable regime, the elasto-inertial instability generates oscillations corresponding to the ones we also report in our study in terms of $\Delta S_{w}$. A detailed characterization of the elasto-inertial stability region is reported in § 4.5.

4.3. Effect of the viscoelastic parameters: $We$ and $\mu _{S}$

The effects of the two viscoelastic parameters are investigated by comparing the results obtained for three polymer concentrations and six Weissenberg numbers. The Weissenberg number is varied in the range $5 \leqslant We\leqslant 1000$ for the three different solvent-to-total dynamic viscosity ratios of $\mu _s = 0.25$, $0.5$ and 0.9.

In figure 5, excursion of the tangential stress distribution along the wall, $\Delta \tau _{w}$ is depicted for $La=2$, $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$. For $\mu _{S} = 0.9$ (panel (c)), upon an increase of two orders of magnitude of the Weissenberg number, no significant change of the wall tangential stress $\Delta \tau _{w} (t)$ is observed, other than a relatively weak speed up of the liquid plug formation. The extra stresses, plotted in dashed lines, are approximately one order of magnitude smaller than the Newtonian stresses (dashed–dotted lines). As a result, the total shear stress very well resembles the typical scenario observed for the Newtonian airway closure. This is well understood because, for $\mu _{S} = 0.9$, the polymeric concentration is low and viscoelastic phenomena are not very relevant. It is, however, remarkable that $\Delta S_{w}$ experiences a qualitative change upon an increase of $We$, passing from a time-decaying trend for low Weissenberg numbers, to a growing trend for high-enough Weissenberg numbers. We stress that the postcoalescence growth of $\Delta S_{w}$ for $\mu _{S} = 0.9$ has features consistent with the dominant elasto-inertial instability observed for $\mu _{S} = 0.25$.

Figure 5. Excursion of the tangential stress distribution along the wall, $\Delta \tau _{w} (t) = \max _z \tau (r=1) - \min _z \tau (r=1)$ (solid line), consisting of the Newtonian ($\Delta \tau _{w}^{N}$, dashed line) and the extra-stress component ($\Delta S_{w}$, dashed–dotted line). The Weissenberg number is varied between 5 and 1000, while the other parameters are fixed at $La = 2$, $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$. In each panel, the values of $\mu _{S}$ are indicated: $\mu _{S}=0.25$ (a), $\mu _{S}=0.5$ (b) and $\mu _{S}=0.9$ (c).

For $La = 2$ and $\mu _{S} = 0.5$, the viscoelastic component of the total stress starts playing an important role in the postcoalescence phase. Even if $\Delta S_{w}$ is important, the postcoalescence stress excursion does not differ qualitatively from the Newtonian case for $We \leqslant 10$. This is well understood considering the fact that, for low-enough Weissenberg numbers, the extra stresses have a Newtonian-like behaviour in the Oldroyd-B model, i.e. $\boldsymbol {S} \approx \mu _{P}(\boldsymbol {\nabla }\boldsymbol {u} + \boldsymbol {\nabla }^T\boldsymbol {u})$. Moreover, the total shear stress relaxation after the quick bifrontal plug growth distributes, almost equally, the stresses between $\tau _{w}^{N}$ and $S_{w}$. The equal distribution of the total wall tangential stresses $\tau _{w}$ between $\tau _{w}^{N}$ and $S_{w}$ is a direct consequence that, for $\mu _{S} = 0.5$, the solvent dynamic viscosity and the corresponding polymeric counterpart $\mu _{S}$ are equal, i.e. $\mu _{S}=\mu _{P}=0.5$.

For $La = 2$, $\mu _{S} = 0.5$ and $We \geqslant 50$ we note the viscoelastic footprint induced by the elastic instability. In fact, for large-enough Weissenberg numbers, the upper-convective derivative of the extra stresses gains increasing importance and overcomes the part of $\boldsymbol {S}$ proportional to the strain tensor. This elastic component of the total tangential stress at the wall is observed in $\Delta \tau _{w}$ for $t>t_{e}=1.6 t_{c}$, where $t_{c}$ denotes the coalescence time of the liquid plug (the first time at which $R_{min}=0$) occurring right before the Newtonian peak of $\Delta \tau _{w}$, and $t_{e}$ denotes the time at which the elastic oscillations due to the instability become significant, i.e. it is responsible of at least 10 % of the tangential wall stress excursions.

4.4. Effect of the Laplace number

The effect of the Laplace number is investigated by varying $La$ within $2\leqslant La \leqslant 200$, for the Weissenberg number in the range $5 \leqslant We \leqslant 1000$ and for three different solvent-to-total liquid dynamic viscosity ratios, i.e. $\mu _{S} = 0.25$, 0.5 and 0.9. In all the cases, upon an increase of the Laplace number, the peak value of $\Delta \tau _{w}$ right after the coalescence grows by approximately 50 % rising from ${\approx } 0.4$ for $La = 2$ up to ${\approx } 0.6$ for $La = 200$. This is demonstrated in figure 6, where three Weissenberg numbers among the six considered are depicted. The first-peak growth is not significantly affected by the viscoelastic parameters $We$ and $\mu _{S}$ because the shear stress peak right after the coalescence is a Newtonian peak, as also discussed in the previous sections. The same trend of the first peak of $\Delta \tau _{w}$ as function of $La$ is, in fact, also observed for the Newtonian airway closure, see figure 2.

Figure 6. Excursion of the tangential stress distribution along the wall, $\Delta \tau _{w} (t) = \max _z \tau (r=1) - \min _z \tau (r=1)$ for $La = 200$ (solid line), $La = 20$ (dashed line) and $La = 2$ (dashed–dotted line). Three Weissenberg numbers are depicted, i.e. 10 (blue), 100 (cyan) and 1000 (magenta). Three $\mu _{S}$ are considered: $\mu _{S}=0.25$ (a), $\mu _{S}=0.5$ (b) and $\mu _{S}=0.9$ (c). The fixed parameters of the simulations are $\mu =1.5 \times 10^{-4}$, $\varrho = 10^{-3}$, $\epsilon = 0.25$, $\lambda = 6$.

A second effect of the Laplace number is observed during the late postcoalescence phase. Increasing the Laplace number corresponds to an increase of importance of the surface tension over the viscous effects. Since the surface tension is proportional to the shear stress that drives the elastic instability, higher Laplace numbers correspond to stronger excitation forces of the viscoelastic medium, hence larger elastic oscillation amplitudes, hence larger viscoelastic-induced wall shear stress excursions. We stress that such an interplay between the Newtonian inertial effects (proportional to $La$) and the viscoelastic inertial effects (proportional to $We$) is a nonlinear component of the elastic instability we observe, and it is not taken into account by the classic linear theory of the hoop stress, developed in the context of creeping flows. As a result, for large-enough Laplace numbers ($La\geqslant 20$, solid and dashed lines in figure 6), large-enough polymeric concentration ($\mu _{S} \leqslant 0.5$, figure 6a,b) and large-enough Weissenberg numbers ($We\geqslant 100$, cyan and magenta lines in figure 6), the viscoelastic effects become very relevant in the wall shear stress budget, up to surpassing the Newtonian peaks observed right after the coalescence. The high Weissenberg numbers required for significant viscoelastic effects explain why the significant elastic oscillations are observed only in a late postcoalescence phase. Indeed, the higher the Weissenberg number, the larger the relaxation time of the viscoelastic fluid. Consequently, the elastic part of the extra stresses reacts with a certain delay to the flow shear rate, and this delay grows proportionally to $We$ in the Oldroyd-B model. A weak effect of the Laplace number on the viscoelastic delay is also observed, since higher $La$ increases the excitation force applied on the viscoelastic medium, resulting in a quicker elastic response, i.e. $La\uparrow$ induces $t_{e}/t_{c}\downarrow$.

4.5. Stability diagram

The postcoalescence stability properties of the viscoelastic flow are characterized by considering the excursion of the wall extra stress $\Delta S_{w}$. The growth rate of the elasto-inertial instability is quantified by fitting $\Delta S_{w}(t)$ using a local exponential fitting function $\Delta S_{w}(t)\approx C_0 + C_1\exp (\sigma t)$. A least squares fit is performed in combination with a pattern recognition algorithm that finds the best match between $\Delta S_{w}$ and the fitting function over the largest time interval that grants a minimum confidence level of 90 %. Three examples are given in the right-hand panels of figure 7 for $\mu _{S}=0.25$, $La = 200$ and $We=10$ (bottom), $We=100$ (middle), and $We=1000$ (top). The black lines show the results of our numerical simulations $\Delta S_{w}$, while the thick lines depict the local exponential fit, either for stable ($\sigma <0$, red) or unstable ($\sigma >0$, blue) conditions. The exponential growth rates $\sigma$ are, thereafter, quantified and plotted against $El_{m}$. An example is depicted in the inset of the left-hand panel of figure 7 for $La = 200$, $\mu _{S}\in [0.25, 0.9]$ and $We\in [5, 2000]$. The neutral conditions (violet markers) $El_{m,n}$ are computed by finding the zeros of $\sigma (El_{m})$ at constant $La$. The neutral stability curve is, thereafter, approximated by spline interpolation of the three neutral stability points (dashed violet line). As shown in figure 7, both, inertial ($La$) and elastic ($We \sqrt {1-\mu _{S}}$) effects have an impact on the neutral stability curve, and employing the modified elasticity number $El_{m}=We \sqrt {1-\mu _{S}}/La$ provides an effective parameter for characterizing the onset of the elasto-inertial instability.

Figure 7. Stability diagram for the postcoalescence elasto-inertial instability. The left-hand panel represents the stable (red) and unstable (blue) conditions for $\mu _{S}=0.25$ (circles), $\mu _{S}=0.5$ (squares) and $\mu _{S}=0.9$ (triangles). The violet markers denote the neutral conditions extrapolated in terms of the modified elastic parameter $El_{m}$. An example of such an extrapolation is reported in the inset for $La = 200$. The dashed violet line denotes the neutral stability curve approximated by spline interpolation of the three neutral stability points (violet markers). The three right-hand panels show the local exponential fit $\Delta S_{w} (t)\approx C_0 + C_1\exp (\sigma t)$ for $\mu _{S}=0.25$, $La = 200$ and three Weissenberg numbers, i.e. $We=10$ (bottom), $We=100$ (middle) and $We=1000$ (top). The black lines denote the extra-stress excursion and the thick lines denote the local exponential fit, either for stable ($\sigma <0$, red) or unstable ($\sigma >0$, blue) conditions.

We stress that the instability-induced flow oscillations are also observed when changing the constitutive model employed for the extra-stress tensor. This is demonstrated by the viscoelastic models employed in Kim et al. (2017) and Casanellas et al. (2016) and further confirmed in the appendices where simulations are performed using the FENE-CR model of Chilcott & Rallison (1988). The results show that the elasto-inertial instability also occurs in the FENE-CR model for physically relevant polymer extensibility values. We therefore conclude that the instability is not an artefact of the viscoelastic model in use. Moreover, we anticipate that limiting the extensibility of polymers, the elasto-inertial instability gets suppressed, which further confirms the importance of stretching the polymeric chains to support the elastic nature of the instability mechanism (see appendices).

5. Conclusions

The airway closure in human lungs has been investigated by computationally studying the Plateau–Rayleigh instability of a thin liquid film lining a rigid tube. The parameters are selected based on the physiological values at the ninth or 10th branching generation of a typical adult lung, where the airway diameters are very small and surface tension dominates over the gravitational and viscous forces. The interfacial instability occurring in our millimetric tubes leads to the formation of a liquid plug that, within the model framework, represents the airway closure typically observed in small bronchioles of human lung.

We here extend our previous study (Romanò et al. 2019) taking into account the effect of mucus viscoelasticity by means of the Oldroyd-B model. The two new parameters included in the liquid phase model, i.e. the Weissenberg number and the solvent-to-total dynamic viscosity ratio, have been investigated for physiologically relevant ranges. This same modelling approach was used by Halpern et al. (2010), who employed the lubrication theory to investigate the capillary instability during the precoalescence phases. They observed a speed up of the airway closure and a weak change of the wall stresses upon the increase of $We$ and $\mu _{S}$. Their same trends are qualitatively and quantitatively reproduced in our study, demonstrating that the major effect of viscoelasticity during the precoalescence phase consists of speeding up the instability process.

The postcoalescence phase, simulated in our study for the first time, shows, however, that viscoelastic effects may be very significant after the formation of a liquid plug and after the quick bifrontal plug growth. In fact, if $We$, $La$ and $\mu _{P}=1-\mu _{S}$ are high enough, the shear rate and the curvature at the interface induce an elasto-inertial instability of the viscoelastic liquid coating the airway. This phenomenon was not observed in Newtonian models since it is a peculiar behaviour of viscoelastic fluids due to an unstable elastic response of the liquid induced by the stretching of the polymeric chains. As demonstrated by extending the viscoelastic model to the FENE-CR extra-stress constitutive equation (see appendix B), such an instability is also predicted by viscoelastic models other than the Oldroyd-B. A similar unstable phenomenon has recently been reported in the literature by Zhou et al. (2016), who employed the upper convective Maxwell model to perform a stability analysis of a viscoelastic fluid coating the walls of a flexible pipe in a gravitational field directed along the pipe axis. For small Deborah (Weissenberg) numbers, they find that the viscoelasticity of the fluid strengthens the wave dispersion of the steady travelling waves, which, in turn, weakens the capillary ripples. They also find that the large velocity gradients near the wave troughs stretch the polymeric molecules leading to thinner capillary ripples, and more pronounced local polymeric stretching. This same principle applies to our instability, where the stretching of the polymeric chains is mainly induced by the curvature of the streamlines after closure. This leads to a thin film with very strong shear and extra stresses (as for Zhou et al. (2016)) that builds the feedback mechanism leading to the elasto-inertial instability we report. Moreover, the same elasto-inertial instability has been experimentally reported by Kim et al. (2017) in a $90^\circ$ bent microchannel for which they characterize the stability limits in terms of the modified elastic number defined as $We\sqrt {1-\mu _{S}}/Re$, where $Re$ denotes the Reynolds number in their channel flow. In our case, the corresponding modified elastic number reads $El_{m}=We\sqrt {1-\mu _{S}}/La$. Consistently with Kim et al. (2017), we showed that the local growth rate of the elasto-inertial instability is well described in terms of $El_{m}$, and an approximation of the neutral stability boundary has been reported in figure 7. For physiologically relevant parameters, we further demonstrated that $El_{m}$ mainly affects the wall shear-stress excursions $\Delta \tau _{w}$, whereas the wall normal-stress excursions $\Delta p_{w}$ seem to be almost insensitive to it.

Our simulations provide significant insights and can be used to estimate the viscoelastic effect on the wall stresses, and, in turn, on the epithelial cells that cover the airway walls. As pointed out by Romanò et al. (2019), the Newtonian bifrontal plug growth can induce a lethal response of the epithelial cells since the wall stresses and their gradients may equal $\max _{t,z}(\vert \tau _{w}\vert )\approx 250\ \textrm {dyn}\ \textrm {cm}^{-2}$, $\max _{t,z}(\vert \partial _z p_{w}\vert )\approx 4.50\times 10^4\ \textrm {dyn} \ \textrm {cm}^{-3}$ and $\max _{t,z}(\vert \partial _z \tau _{w}\vert )\approx 8\times 10^3\ \textrm {dyn} \ \textrm {cm}^{-3}$, which are far beyond the damaging threshold experimentally established by Bilek et al. (2003) ($\max _{t,z}(\vert \tau _{w}\vert ) > 12.9\ \textrm {dyn} \ \textrm {cm}^{-2}$, $\max _{t,z}(\vert \partial _z \tau _{w}\vert ) > 2.1\times 10^3\ \textrm {dyn} \ \textrm {cm}^{-3}$ and $\max _{t,z}(\vert \partial _z p_{w}\vert ) > 3.21\times 10^4\ \textrm {dyn} \ \textrm {cm}^{-3}$) and Huh et al. (2007) ($\max _{t,z}(\vert \tau _{w}\vert ) > 98.58\ \textrm {dyn} \ \textrm {cm}^{-2}$). With the present study we demonstrate that the Newtonian peak of the wall stresses induced by the quick bifrontal plug growth is, in fact, almost insensitive to the viscoelastic parameters of the Oldroyd-B model, i.e. to the Weissenberg number $We$ and the solvent-to-total dynamic viscosity ratio $\mu _{S}$. This consideration further stresses that the bifrontal plug growth may severely damage the epithelial cells, regardless of the viscoelastic properties of mucus. On top of that, when the $We$, $La$, $We/La$ and $\mu _{P}=1-\mu _{S}$ are sufficiently high, the elastic instability in the viscoelastic liquid induces additional shear stresses that further contribute to a damaging response of the epithelial cells, inducing $\Delta S_{w}$ of the same order of magnitude observed during the Newtonian peak. We further remark that such results are robust also in terms of the airway aspect ratio. Indeed, preliminary simulations in which we double the periodic domain ($\lambda =12$) show no remarkable change of the instability mechanisms and of the amplitude of the stress peaks. Only a shorter closure time and, consequently, a time shift of the elasto-inertial instability is observed passing from $\lambda =6$ to $\lambda =12$.

Our study provided crucial insights on the effect of the viscoelastic properties of mucus. To focus on them, we decided to simplify the airway model by considering rigid walls. We, however, point out that the elastic deformation of the airway wall can significantly favour the capillary instability leading to airway closure (see Halpern & Grotberg 1992). Moreover, elastic airway walls can even result in a structural instability of the airway cross-section, which also favours the airway closure (see White & Heil 2005; Heil et al. 2008). We, moreover, speculate that including the compliance of airway walls may provide further insights on the postcoalescence phase by unravelling the eventual interplay between mechanical deformations and the elasto-inertial instability.

Future studies will include the effect of viscoplasticity and, finally, the interplay between viscous, elastic and plastic behaviours of the mucus. Adding one effect at a time will allow us to really highlight the net role of non-Newtonian behaviours of mucus and will help constructing reduced-order models. On top of that, we think that the interaction with deformable walls and surfactant is a relevant topic for all these models applied to airway closure, hence it will be addressed in our future works.