1. Introduction
The human lungs form a branching network with circular-like cross-sections, where each bifurcation results in progressively narrower tubes. This branching hierarchy starts with the trachea and proceeds through the bronchi, bronchioles, respiratory bronchioles and finally the alveoli (Grotberg Reference Grotberg2011). The airway diameters of the first 15 generations satisfy the relation
$a_i=a_02^{-i/3}$
, where
$a_0$
denotes the trachea’s radius which typically falls within the range
$a_0 \in [0.8,1]$
cm for adults (Weibel & Gomez Reference Weibel and Gomez1962), and
$a_i$
represents the radius of the airway at generation
$i$
. The liquid film that coats the airways in the first 16 generations is composed of two distinct layers: a mucus layer and a serous layer, the latter in contact with the airway walls (Widdicombe et al. Reference Widdicombe, Bastacky, Wu and Lee1997). Typically, the thickness of this film constitutes around
$2\,\%{-}4\,\%$
of the airway diameter, though, in diseased states, it may thicken to as much as
$20\,\%$
(Codd et al. Reference Codd, Lambert, Alley and Pack1994; Yager et al. Reference Yager, Cloutier, Feldman, Bastacky, Drazen and Kamm1994). The role of the fluid film on the airway surface is to help remove pollutants, protect against viruses and bacteria and keep the epithelial cells moist (Grotberg Reference Grotberg1994).
The formation, propagation and rupture of the plugs in respiratory airways have significant physiological consequences (Burger Jr & Macklem Reference Burger and Macklem1968; Hughes, Rosenzweig & Kivitz Reference Hughes, Rosenzweig and Kivitz1970; Dargaville, South & McDougall Reference Dargaville, South and McDougall1996; Cassidy et al. Reference Cassidy, Halpern, Ressler and Grotberg1999; Grotberg Reference Grotberg2011). It is well established that a thicker liquid lining has the potential to occlude the airways and to reduce gas exchange for millimetres or even smaller airway diameters. The corresponding airway closure arises from the formation of a liquid plug, driven by the capillary instability of Plateau–Rayleigh (Gauglitz & Radke Reference Gauglitz and Radke1988; Grotberg & Jensen Reference Grotberg and Jensen2004; Shemilt et al. Reference Shemilt, Thompson, Horsley, Whitfield and Jensen2025). Once a liquid plug is formed and the airway is occluded, the affected distal airways are prevented from exchanging gas between the upper and the lower generations. Airway reopening is the process by which a pressure difference caused by breathing or coughing pushes the movement of the liquid plug along the airway until it ruptures (Grotberg & Jensen Reference Grotberg and Jensen2004; Ryans et al. Reference Ryans, Fujioka, Halpern and Gaver2016; Grotberg Reference Grotberg2019; Elias-Kirma et al. Reference Elias-Kirma, Artzy-Schnirman, Sabatan, Dabush, Waisman and Sznitman2021). Large mechanical stresses resulting from the corresponding interfacial dynamics can damage the pulmonary epithelium (see Huh et al. Reference Huh, Fujioka, Tung, Futai, Paine, Grotberg and Takayama2007; Tavana et al. Reference Tavana, Zamankhan, Christensen, Grotberg and Takayama2011; Fujioka et al. Reference Fujioka, Halpern, Ryans and Gaver2016; Grotberg Reference Grotberg2019; Romanò et al. Reference Romanò, Fujioka, Muradoglu and Grotberg2019, Reference Romanò, Muradoglu, Fujioka and Grotberg2021, Reference Romanò, Muradoglu and Grotberg2022; Dietze Reference Dietze2024). Moreover, the complex rheological response of the mucus layer for airway reopening, can lead to an effective delay or even prevent the reopening of airways (Zamankhan et al. Reference Zamankhan, Helenbrook, Takayama and Grotberg2012; Grotberg Reference Grotberg2019; Bahrani et al. Reference Bahrani, Hamidouche, Moazzen, Seck, Duc, Muradoglu, Grotberg and Romanò2022), with profound respiratory implications that may prove fatal in some cases (Synek et al. Reference Synek, Beasley, Goulding, Holloway and Holgate1994).
The serous layer primarily exhibits Newtonian or weakly viscoelastic behaviours, while mucus is frequently described as a markedly non-Newtonian liquid (Girod et al. Reference Girod, Zahm, Plotkowski, Beck and Puchelle1992; Cone Reference Cone2009; Lai et al. Reference Lai, Wang, Wirtz and Hanes2009; Lavalle et al. Reference Lavalle, Mergui, Grenier and Dietze2021). The mucus layer is a gel-like fluid, which contains mucins secreted by caliciform cells and glands. This gel consists of
$95\,\%$
water and contains
$1\,\%{-}2\,\%$
proteins,
$1\,\%$
glycoproteins (called mucins) and
$1\,\%$
lipids (Vasquez & Forest Reference Vasquez and Forest2014; Spagnolie Reference Spagnolie2015; Lafforgue et al. Reference Lafforgue, Seyssiecq, Poncet and Favier2018). The rheological properties of mucus are intrinsically linked to its solid and polymeric concentration, which is governed by the body’s physiological state. For instance, while the solid content of healthy pulmonary mucus typically hovers around
$2\,\%$
wt, it rises to
$8\,\%$
in patients with cystic fibrosis (CF) (Hill et al. Reference Hill2014). Corresponding non-Newtonian model parameters are reported in Erken et al. (Reference Erken, Fazla, Muradoglu, Izbassarov, Romanò and Grotberg2023) who considered the Saramito–Herschel–Bulkley (Saramito–HB) constitutive law (Saramito Reference Saramito2009) and fitted the rheological properties of sputum across asthma, chronic obstructive pulmonary disease (COPD), CF and healthy, based on the experimental results of Patarin et al. (Reference Patarin, Ghiringhelli, Darsy, Obamba, Bochu, Camara, Quétant, Cracowski, Cracowskir and Robert de Saint Vincent2020). Their findings revealed that the sputum associated with CF and COPD are distinguished by their significant elastic characteristics and yield stress, indicating that such pathological conditions tend to return to their equilibrium state after deformation and considerably enhance flow resistance. In contrast, the sputum of asthma patients displays a moderate elasticity and yield stress, and the mucus of healthy individuals exhibits the weakest traces of rheological complexity among all. Similar considerations are reported also in the other literature sources (Vasquez & Forest Reference Vasquez and Forest2014; Lafforgue et al. Reference Lafforgue, Seyssiecq, Poncet and Favier2018; Patarin et al. Reference Patarin, Ghiringhelli, Darsy, Obamba, Bochu, Camara, Quétant, Cracowski, Cracowskir and Robert de Saint Vincent2020).
The viscoelastic properties, yield stress and shear-dependent viscosity play a crucial role in forming liquid plugs. Shemilt et al. (Reference Shemilt, Horsley, Jensen, Thompson and Whitfield2022) investigate the instability of the viscoplastic Bingham fluid layer in an airway closure model using numerical solutions of the long-wave evolution equations. It is found that when the layer has finite thickness, the critical thickness for the formation of a liquid plug increases with the capillary Bingham number. Halpern, Fujioka & Grotberg (Reference Halpern, Fujioka and Grotberg2010) carried out a corresponding stability analysis based on a thin film asymptotic approach, while Romanò et al. (Reference Romanò, Muradoglu, Fujioka and Grotberg2021) simulated the whole closure process modelling the effects of mucus viscoelasticity by using the Oldroyd-B and FENE-CR (finite extendable nonlinear elastic – Chilcott and Rallison) constitutive laws. Both these studies reported an increased growth rate of the capillary instability that ultimately leads to airway closure. Moreover, Romanò et al. (Reference Romanò, Muradoglu, Fujioka and Grotberg2021) identified an elastoinertial instability produced during the postcoalescence phase for high elastic numbers. Erken et al. (Reference Erken, Fazla, Muradoglu, Izbassarov, Romanò and Grotberg2023) analysed airway closure effects using an elastoviscoplastic (EVP) liquid layer, modelled with the Saramito–HB model. Their results indicate that the liquid layer remains largely unyielded before closure, suggesting that elasticity and solvent viscosity are dominant factors. The elastic properties are crucial in determining closure time and occurrence of liquid plug formation. Additionally, wall stresses in COPD and CF do not relax smoothly post-closure, maintaining high levels similar to peak closure values, which is consistent with the discoveries of Romanò et al. (Reference Romanò, Muradoglu, Fujioka and Grotberg2021). This has been further investigated by Fazla et al. (Reference Fazla, Erken, Izbassarov, Romanò, Grotberg and Muradoglu2024), who included the effect of isotropic kinematic hardening in the constitutive law. They identified a few pathological conditions for which the energy dissipation before yielding represented by the kinematic hardening can ease airway closure to occur.
In contrast to the plug formation, the rupture of a liquid plug does not rely on instability, but rather results from the deposit of a trailing film thicker than the leading coating film. Under these conditions the liquid plug keeps being drained all along its propagation until it ruptures. Bilek et al. (Reference Bilek, Dee and Gaver2003) experimentally explored the impact of airway reopening on lung epithelial cell damage through a model involving a semi-infinite bubble travelling through a narrow, fluid-filled channel lined with pulmonary epithelial cells. Their findings revealed that cell damage is intensified as reopening velocity is decreased, though the presence of pulmonary surfactant acts as a protective agent against injury. Additionally, the steep pressure gradient near the edge of the leading air finger was identified as the primary cause of cellular damage. The study of Kay et al. (Reference Kay, Bilek, Dee and Gaver2004) further indicated that repeated cycles of reopening and closure resulted in cellular injury, even under conditions that would not cause significant damage during a single reopening event. Huh et al. (Reference Huh, Fujioka, Tung, Futai, Paine, Grotberg and Takayama2007) delved into the mechanical trauma inflicted upon primary human small airway epithelial cells (SAEC) by the propagation and rupture of liquid plugs. They lined a microfluidic channel with SAEC and a thin liquid film. When subjected to physiological fluid mechanical stresses akin to those encountered during surfactant-deficient airway reopening, the SAEC experienced considerable cellular damage. This damage intensified with the increasing frequency of plug propagation and rupture. Viola et al. (Reference Viola2024) generalized these findings to semicircular cross-sections and investigated the role of the liquid viscosity. Fujioka, Takayama & Grotberg (Reference Fujioka, Takayama and Grotberg2008) studied the dynamics of plug propagation and discussed the effects of film thickness, initial liquid plug length, driving pressure and surface tension on the behaviour of the liquid plug and the wall stress prior to rupture. Hassan et al. (Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011) numerically simulated the topological changes of the plug interface during the liquid plug rupture. They confirmed the prior-to-rupture results of Fujioka et al. (Reference Fujioka, Takayama and Grotberg2008) and reported that the peak wall stresses are intensified as the initial precursor film thickness diminishes. Following rupture, the peak mechanical stresses diminish as the blockage dissipates and the flow fully develops. It appears that a reduction in the pressure drop and an increase in the Laplace number contribute to a delay in rupture, as fluids with significant inertial forces take longer to respond to minor pressure drops. Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019) performed a computational investigation of surfactant-laden liquid plug propagation and rupture in the ninth-to-tenth generation airway, employing a front-tracking method. The numerical approach was initially validated for surfactant-free scenarios, showing strong agreement with previous simulations by Fujioka et al. (Reference Fujioka, Takayama and Grotberg2008) and Hassan et al. (Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011). The study revealed that surfactant markedly decreases pressure and shear stress while delaying rupture.
As observed for the airway closure, the airway reopening dynamics are significantly affected by the non-Newtonian characteristics of the airway surface liquid. The process of pulmonary airway reopening in two-dimensional channels of Bingham fluid plugs was numerically investigated by Zamankhan et al. (Reference Zamankhan, Helenbrook, Takayama and Grotberg2012), and the influence of the capillary number and Bingham number was discussed, emphasizing that the Bingham number slows down the propagation of the plug. Using lubrication theory, Jalaal & Balmforth (Reference Jalaal and Balmforth2016) analysed the formation of a viscoplastic film on the pipe wall by long bubbles moving under the effect of ambient fluid flow and compared it with simulation results, and found that the presence of yield stresses may affect the flow behaviour of the end of the bubbles and thicken the liquid film. Hu, Romanò & Grotberg (Reference Hu, Romanò and Grotberg2020) numerically investigated the effects of surface tension and yield stress on the rupture of a mucus plug in a rectangular cross-section channel. The yield stress was found to slow down the rupture and had no significant effect on the wall shear stress. Bahrani et al. (Reference Bahrani, Hamidouche, Moazzen, Seck, Duc, Muradoglu, Grotberg and Romanò2022) studied the rupture and propagation of mucus plugs by injecting synthetic mucus into prewetted capillaries and found that thickening of the fluid film due to non-Newtonian effects accelerates plug rupture, while yield stresses impede rupture. Despite its significance, a detailed numerical investigation of airway reopening under complex rheological conditions is still missing.
In this study, we account for the elastic, yield-stress and shear-thinning characteristics of mucus using the Saramito–HB model. The remainder of the paper is structured as follows. Sections 2 and 3 detail the problem formulation and numerical method, respectively. Section 4 presents the results and discussion of this study, organizing them into three main parts: (i) at first the viscoelastic results are discussed; (ii) thereafter we present and discuss the viscoplastic effects; (iii) finally, we combine viscoelastic and viscoplastic effects to investigate the interplay between them and draw conclusion on the robustness standalone mechanisms identified in (i) and (ii). A verification of the robustness of such an analysis with respect to the presence of surfactant on the air–mucus interface closes the results section. Finally, the conclusions of the study are given in § 5.
2. Problem formulation
The airway reopening model employed in this study is shown in figure 1. It consists of a rigid cylindrical tube characterized by radius
$a$
and length
$L$
. The leading and trailing air fingers are surrounded by the liquid film of thickness
$h$
, with reference dynamic viscosity
$\mu _L$
(evaluated at a reference shear rate) and density
$\rho _L$
. At
$t=0$
, a liquid plug of length
$l_p$
is centred at
$z=z_p$
. The gas phase depicted in white has constant dynamic viscosity
$\mu _G$
and density
$\rho _G$
. Additionally, there is a surface tension
$\sigma$
, initially assumed constant (no surfactant), acting at the gas–liquid interface. The flow is presumed to be incompressible and exhibits axial symmetry. The gas is Newtonian, whereas the liquid is non-Newtonian.
Figure 1. Schematic of the plug propagation in a human airway (a) and the corresponding EVP model with insoluble interfacial surfactant (b).
2.1. Governing equations
The governing equations are non-dimensionalized by using a capillary scaling, i.e. length, time, pressure and velocity are scaled with
$a$
,
$\mu _La/\sigma$
,
$\sigma /a$
,
$\sigma /\mu _L$
, respectively. The dimensionless form of the Navier–Stokes and continuity equations yields (Popinet Reference Popinet2018; Romanò et al. Reference Romanò, Muradoglu, Fujioka and Grotberg2021)
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} & = 0 , \end{align}
\begin{align} {La}\tilde {\rho }(\partial _t \boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{\nabla }\boldsymbol{u}) &=-\boldsymbol{\nabla }p+\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{\tilde {\tau }}+\chi \boldsymbol{n}\delta_{s} , \end{align}
where
${La}$
is the Laplace number,
$\boldsymbol{u}=(u_r,0,u_z)$
is the velocity vector field,
$p$
is the pressure field,
$t$
is the time. The total local curvature of the interface is
$\chi =\boldsymbol{\nabla }\boldsymbol{\cdot } \boldsymbol{n}$
, where
$\boldsymbol{n}$
refers to the outward unit normal at the interface. The surface Dirac delta function
$\delta_{s}$
equals zero everywhere except at the gas–liquid interface,
$\tilde {\rho }$
is the variable density field required by the one-field approach to include the effect of the gas-to-liquid density ratio
$\rho =\rho _G/\rho _L$
. In the liquid phase
$\tilde {\rho } = 1$
, while
$\tilde {\rho } =\rho$
in the gas phase.
In (2.1),
$\boldsymbol{\tilde {\tau }}$
is the non-dimensional deviatoric part 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
$\mu =\mu _G/\mu _L$
is the gas-to-liquid dynamic viscosity ratio. The liquid is modelled using the Saramito–HB model (Saramito Reference Saramito2009), a combination of the Oldroyd-B and HB models. This model effectively captures the EVP character of the mucus, i.e. viscoelasticity, viscoplasticity and shear-thinning. The liquid film is assumed to behave as a diluted polymer (solute) immersed in a Newtonian liquid space (solvent). Then
$\tilde {\boldsymbol{\tau }}$
in the liquid phase is expressed as a superposition of both solvent viscous stress and solute extra stresses
$\boldsymbol{S}$
, as follows:
\begin{equation} \tilde {\boldsymbol{\tau }} =\boldsymbol{\tau }_L=\mu_{S}(\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla} ^T \boldsymbol{u})+\boldsymbol{S}, \end{equation}
where
$\mu_{S}=\mu _{\textit{L,S}}/\mu _L$
is the reference solvent-to-total dynamic viscosity ratio in the liquid film and
$\mu _{\textit{L,S}}$
is the dynamic viscosity of the solvent. The Oldroyd-B model is used to describe the viscoelastic character of the mucus where a typical polymer chain is modelled as a linear elastic dumbbell. The primary limitation of the Oldroyd-B model lies in its assumption of infinite polymer extensibility (
$\mathcal{L}$
), i.e.
$\mathcal{L} \rightarrow \infty$
. Romanò et al. (Reference Romanò, Muradoglu, Fujioka and Grotberg2021) used a finite extension nonlinear elastic model, in particular FENE-CR, to verify the robustness of the Oldroyd-B model to elastic stresses build-up mechanisms by testing
$\mathcal{L}\in \mbox{[2, 50]}$
, which is within physiologically relevant regimes (
$\mathcal{L}\lt 150$
, see Wagner, Bourouiba & McKinley Reference Wagner, Bourouiba and McKinley2015). The results by Romanò et al. (Reference Romanò, Muradoglu, Fujioka and Grotberg2021) demonstrated that, for
$\mathcal{L}\gt 30$
, the deviations in outcomes between FENE-CR and Oldroyd-B were negligible for viscoelastic-induced instabilities reported in airway closure. We therefore selected the Oldroyd-B model to carry out our investigation of the viscoelastic effects in this study.
By applying kinetic theory, a constitutive equation can be derived to characterize the stresses (Fattal & Kupferman Reference Fattal and Kupferman2005; Hao & Pan Reference Hao and Pan2007) as
\begin{equation} \textit{Wi} \big[ \partial _t \boldsymbol{S}+(\boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{\nabla })\boldsymbol{S} -(\boldsymbol{\nabla }\boldsymbol{u})\boldsymbol{S}-\boldsymbol{S}(\boldsymbol{\nabla} ^T \boldsymbol{u}) \big] +\mu_{P}^{-1}\boldsymbol{S}= \left (1-\mu_{S}\right )(\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla} ^T \boldsymbol{u}), \end{equation}
where
$\mu_{P}$
is the polymeric dynamic viscosity ratio, defined as the ratio between the local polymeric viscosity and the reference polymeric viscosity
$\mu _{\textit{L,P}}$
, which equals 1 for the Oldroyd-B model and varies as a field in Saramito–HB model. The reference total viscosity is computed by
$\mu _L=\mu _{\textit{L,S}}+\mu _{\textit{L,P}}$
.
In order to investigate the effects of viscoplasticity and shear-thinning, the Saramito–HB polymeric viscosity ratio
$\mu_{P}$
is defined by (Saramito Reference Saramito2007, Reference Saramito2009)
\begin{equation} \mu_{P}= \cfrac {1}{1-\mu_{S}}\, \max \left [ 0,\frac {K_c {\left | \boldsymbol{S}^{d} \right |}^n}{\left | \boldsymbol{S}^{d} \right | - Bi} \right ]^\frac {1}{n}, \end{equation}
where
$n$
is power-law index, and the magnitude of
$\boldsymbol{S}^{d}$
can be written as
$\left | \boldsymbol{S}^{d} \right |=\sqrt {S_{\textit{ij}}^{d}S_{\textit{ij}}^{d}/2}$
, where
$d$
denotes the deviatoric part of the tensor, i.e.
$\boldsymbol{S}^{d} = \boldsymbol{S} - tr(\boldsymbol{S})\ \boldsymbol{I}/3$
. At
$t=0$
,
$\mu_{P}$
is set to
$\mu_{P}(t=0) = 1$
. The dimensionless consistency parameter
$K_c$
is computed using
$K_c=K \sigma ^{n-1}a^{-n+1}\mu _L^{-n}$
, where
$K$
is consistency parameter given by
$\mu _{\textit{L,P}}=\tau _y a/u_c+K(a/u_c)^{n-1}$
, with
$u_c=\mu _L/\sigma$
standing for the characteristic velocity, and
$\tau _y$
being the yield stress.
Several independent non-dimensional groups arising from the momentum (2.1b
): the Laplace number
${La}$
, the Weissenberg number
$\textit{Wi}$
, the Bingham number
$\textit{Bi}$
, the gas-to-liquid density ratio
$\rho$
and the gas-to-liquid dynamic viscosity ratio
$\mu$
, and the reference 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
$\epsilon$
. The non-dimensional parameters can be summarized, as follows:
\begin{eqnarray} \begin{array}{cc} {La}=\cfrac {\rho _L \sigma a}{\mu ^2_L}, \quad \textit{Wi}=\cfrac { \sigma (1-\mu_{S})}{a G}, \quad \textit{Bi}=\cfrac {\tau _y a}{\sigma }, \quad \mu_{S}=\cfrac {\mu _{\textit{L,S}}}{\mu _L}, \\ \rho =\cfrac {\rho _G}{\rho _L} , \quad \mu =\cfrac {\mu _{G}}{\mu _L}, \quad \lambda =\cfrac {L}{a}, \quad \epsilon =\cfrac {h}{a}, \quad L_p =\cfrac {l_p}{a}, \quad \\[16pt] \Delta p=\cfrac {a\left [p(z=0) - p(z=\lambda )\right ]}{\sigma }, \end{array} \end{eqnarray}
where
$G$
is the elastic modulus, then the reference relaxation time is calculated as
$\varLambda _r=\mu _{\textit{L,P}}/G$
.
The pronounced elasticity of the non-Newtonian fluid renders the numerical solution of the constitutive equation difficult to integrate in time (Fattal & Kupferman Reference Fattal and Kupferman2004, Reference Fattal and Kupferman2005), due to the elevated Weissenberg numbers characterizing pulmonary flows. The flow is supposed, in fact, to embed regions of high stress and fine features, which are known to induce numerical instabilities. This potential instability imposes practical constraints on the values of the relaxation time of the EVP fluid
$\varLambda$
. To surmount this challenge, (2.3) is recast using the log-conformation representation by Fattal & Kupferman (Reference Fattal and Kupferman2004, Reference Fattal and Kupferman2005).
Finally, to mimic inhaling conditions, a constant pressure difference
$\Delta p = p(z=0) - p(z=\lambda )\gt 0$
between the left and right ends of the airway is enforced to drive the liquid plug moving from left to right. At
$t=0$
, the left and right film thicknesses are fixed to
$\epsilon$
, and the capillary pressure jump across the liquid–gas interface is included as
$(1-\epsilon )^{-1}$
. The mathematical problem (2.1) is closed by enforcing: a piecewise-constant Dirichlet and homogeneous Neumann conditions for pressure and velocity, respectively, at the inlet and outlet, i.e. at
$z=0$
and
$z=\lambda$
. No-slip and no-penetration conditions are enforced along the wall, i.e. at
$r=1$
, while axisymmetric conditions are set at
$r=0$
.
2.2. Physiological parameters
This study aims to model the propagation and rupture of liquid plugs in adult human lungs, considering the EVP properties of mucus. All parameters are selected to be compatible with the physiological ranges. We consider healthy and pathological conditions, and extend the non-dimensional parametric ranges slightly beyond the reference values reported in the literature in order to test the robustness of our results. For an airway characterized by a sufficiently thick coating and a narrow radius, the Plateau–Rayleigh instability generates a liquid plug that obstructs the airway, hindering distal gas exchange. In accordance with Guidotti (Reference Guidotti1997), we adopt an airway radius of
$a\in [0.65,\ 1.2]$
mm as representative of the typical dimensions found in the eighth to the tenth generation of adult lung airways where surface tension dominates and gravity is negligible (Romanò et al. Reference Romanò, Muradoglu and Grotberg2022). The dimensionless initial plug length is set to
$L_p(t=0)=1$
, while the centre of the plug is initially set at
$z_p(t=0)=2.95a$
and the length-to-radius ratio
$\lambda$
of the airway is extended to 16 in order to reduce the influence of inflow and outflow conditions. The initial dimensionless coating liquid film thickness
$\epsilon$
spans 0.05–0.1 for both the leading and trailing films, as indicated by Hassan et al. (Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011) and Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019). Following Romanò et al. (Reference Romanò, Fujioka, Muradoglu and Grotberg2019), the multilayer liquid’s density comprising mucus and serous components, closely matches that of water (
$\rho _L=1$
$\text{g cm}^{-3}$
), resulting in a gas-to-liquid density ratio of
$\rho =10^{-3}$
. The serous layer exhibits a dynamic viscosity similar to water, around 0.01 poise. Meanwhile, the dynamic viscosity of mucus can vary widely, spanning orders of magnitude. To simplify, Tai et al. (Reference Tai, Bian, Halpern, Zheng, Filoche and Grotberg2011) and Romanò et al. (Reference Romanò, Fujioka, Muradoglu and Grotberg2019) propose treating the bilayer liquid film as a homogeneous medium, with dynamic viscosity (
$ \mu _L$
) ranging from 0.126 to 1.26 poise. We here specialize this approach taking into account the healthy and pathological mucus measurements of Patarin et al. (Reference Patarin, Ghiringhelli, Darsy, Obamba, Bochu, Camara, Quétant, Cracowski, Cracowskir and Robert de Saint Vincent2020) and the corresponding data postprocessing by Erken et al. (Reference Erken, Fazla, Muradoglu, Izbassarov, Romanò and Grotberg2023). By using the HB fit of Erken et al. (Reference Erken, Fazla, Muradoglu, Izbassarov, Romanò and Grotberg2023), we identified
$\mu _L\in [0.01,\ 0.02]$
Pa s for healthy conditions and
$\mu _L\in [0.02,\ 0.2]$
Pa s for pathological conditions that include asthma, COPD and CF. Additionally, dynamic viscosity of air at
$37\,^{\circ}$
C (
$ \mu _G = 1.89 \times 10^{-4}$
poise) is considered. Given the study’s emphasis on elastoviscoplasticity effect on wall stresses, we maintain a constant characteristic gas-to-liquid dynamic viscosity ratio of
$\mu =1.5 \times 10^{-3}$
throughout our analysis. For a detailed analysis of shear stresses with
$\mu$
, we refer to Romanò et al. (Reference Romanò, Fujioka, Muradoglu and Grotberg2019). The surface tension
$\sigma$
depends mainly on the amount of surfactant and ranges in
$\sigma \in [0.02,\ 0.03]$
N m−1 in healthy condition (Clements, Brown & Johnson Reference Clements, Brown and Johnson1958; Moriarty & Grotberg Reference Moriarty and Grotberg1999) and in
$\sigma \in [0.04,\ 0.06]$
N m−1 for surfactant deficient conditions (Schürch et al. Reference Schürch, Gehr, Im Hof, Geiser and Green1990; Moriarty & Grotberg Reference Moriarty and Grotberg1999; Bian et al. Reference Bian, Tai, Halpern, Zheng and Grotberg2010). The Laplace number emerges as a critical determinant in plug propagation and rupture (Fujioka et al. Reference Fujioka, Takayama and Grotberg2008; Hassan et al. Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011; Muradoglu et al. Reference Muradoglu, Romanò, Fujioka and Grotberg2019). In the current investigation utilizes the range
$ {La} \in \left \{ 50, 100, 200, 300, 500 \right \}$
. Based on the measurements of Patarin et al. (Reference Patarin, Ghiringhelli, Darsy, Obamba, Bochu, Camara, Quétant, Cracowski, Cracowskir and Robert de Saint Vincent2020), fitted by Erken et al. (Reference Erken, Fazla, Muradoglu, Izbassarov, Romanò and Grotberg2023) and Fazla et al. (Reference Fazla, Erken, Izbassarov, Romanò, Grotberg and Muradoglu2024), the elastic modulus
$G$
is approximately 0.1 Pa in healthy conditions and
$\approx [0.1,\ 2.8]$
Pa in disease conditions. In this investigation, we explore effects of
$G$
ranging from 0 to 0.5, which covers both healthy and diseased conditions. Consistently, the fits of Erken et al. (Reference Erken, Fazla, Muradoglu, Izbassarov, Romanò and Grotberg2023), Fazla et al. (Reference Fazla, Erken, Izbassarov, Romanò, Grotberg and Muradoglu2024) provide
$\tau _y\approx 0.5$
Pa in healthy conditions and
$\approx [1,\ 12]$
Pa for pathological conditions. Finally, they also identify the power-law index
$n$
to range in
$n\in [0.5,\ 0.7]$
for healthy conditions and in
$n\in [0.3,\ 0.9]$
for pathological conditions. We therefore consider ranges of
$\tau _y\in [0,\ 3.75]$
and
$n\in [0.3,\ 1]$
in order to cover most of the physiologically relevant regimes. The range of bulk and interfacial dimensional and non-dimensional parameters covered in this study is summarized in table 1. As for the solvent-to-total characteristic viscosity ratio,
$\mu_{S}$
is a key investigation parameter as well (Halpern et al. Reference Halpern, Fujioka and Grotberg2010), and we consider an extensive range which varies from 0.25 to 0.9 in order to study the impact of polymeric concentration, yet assuming that
$\mu_{S}$
does not vary during the course of a simulation.
Table 1. Ranges of dimensional and non-dimensional parameters used in this study, compared with the reference values for healthy and pathological mucus conditions. Based on the values of
$G$
,
$\mu _L$
and
$\mu_{S}$
, the corresponding range for the polymer relaxation time investigated in this study is computed by
$\varLambda _r = (1-\mu_{S})\mu _L/G \in [0.002,\ 2]$
s.
3. Numerical method
The airway reopening phenomenon is modelled in an axisymmetric fashion, with the enforcement of
$\partial _{\phi }=0$
and
$u_{\phi }=0$
. Integrating the governing equation (2.2) over time is achieved through a fractional step method employing a second-order pressure-correction projection scheme. Implicit treatment is applied to the viscous term, while explicit discretization is employed for the nonlinear convective term using the Bell–Collela–Glaz advection scheme (Bell, Colella & Glaz Reference Bell, Colella and Glaz1989).
The spatial discretization is conducted using 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. For handling the two-phase flow, the volume of fluid (VOF) method is utilized, capturing the liquid–gas interface via a fraction field
$f(r,z,t)$
. This fraction field represents the volume fraction of liquid within a computational cell: 1 for solely liquid, 0 for solely gas and
$f(r,z,t)\in (0,1)$
for cells containing both liquid and gas. The condition
$f(r,z,t)\in (0,1)$
identifies computational cells housing the liquid–gas interface. Leveraging the fraction field, the variable fluid properties in the one-field approach (Popinet Reference Popinet2003, Reference Popinet2009) are expressed as
\begin{align} \tilde {\rho }[f(r,z,t)] &=f+(1-f)\rho, \end{align}
\begin{align} \tilde {\boldsymbol{\tau }}[f(r,z,t)] &=f{\boldsymbol{\tau }}_L+(1-f){\boldsymbol{\tau }}_G .\end{align}
To complete the description of the dynamics of the fraction field for immiscible fluids, modelled as the simple advection of
$f$
due to the flow velocity
$\boldsymbol{u}$
, an additional equation is needed,
\begin{equation} \partial _tf+\boldsymbol{\nabla }\boldsymbol{\cdot }(f\boldsymbol{u})=0. \end{equation}
The discretization of (3.2) follows a staggered time approach, employing the second-order scheme developed by Popinet (Reference Popinet2009). A piecewise-linear geometrical VOF method is utilized to reconstruct the liquid–gas interface, representing the interface as a straight line within each computational cell by solving (3.1). Advection of the interface incorporates the local unit normal using the mixed Young-centred method (Aulisa, Manservisi & Scardovelli Reference Aulisa, Manservisi and Scardovelli2006) and the local indicator function
$f$
. In alignment with the projection method, (2.5) is discretized using the Bell–Collela–Glaz method (Bell et al. Reference Bell, Colella and Glaz1989). Upon use of the one-field approach, surface tension numerical issues typically arise from the discretization of
$\chi \boldsymbol{n}\delta_{s}$
, potentially causing parasitic currents, as discussed Brackbill, Kothe & Zemach (Reference Brackbill, Kothe and Zemach1992) and Popinet & Zaleski (Reference Popinet and Zaleski1999), in either front-tracking or front-capturing methods. Popinet (Reference Popinet2009) demonstrated that a combination of a balanced-force approach and a height-function estimator can significantly reduce these numerical artefacts, enabling the attainment of second-order convergence rates even in demanding benchmarks. Subsequently, the methodology proposed by Popinet (Reference Popinet2009) is employed in this study.
In the following sections, the simulations demonstrate the utilization of biologically relevant parameters in airway reopening problems. A Cartesian grid serves as the computational mesh, customized to capture key quantities of interest in our study, including pressure and stress. As the focus of our study is on the liquid plug dynamics and the wall stresses, we optimize the numerical resources by employing a tailored adaptive mesh refinement, which is a technique that dynamically adjusts the grid resolution based on user-defined criteria, allowing for efficient simulations that focus computational resources where they are most needed. Figure 2 shows an example of the time evolution of the mesh for one of our simulations at three instants in time. The spatial and temporal discretization schemes are numerically implemented using the open-source software Basilisk (Popinet & COllaborators Reference Popinet2013–2025). This same solver has been tested for similar multiphase flows. In particular, in our previous study (Romanò et al. Reference Romanò, Fujioka, Muradoglu and Grotberg2019) we validated it against experiments for an airway closure model, while in Viola et al. (Reference Viola2024) we reproduced the pressure excursion and the plug velocity of microfluidic experiments modelling airway reopening. Moreover, the computational framework of viscoelasticity has undergone validation in modelling airway closure, as detailed in our prior work (Romanò et al. Reference Romanò, Muradoglu, Fujioka and Grotberg2021) and we present a further validation for our airway reopening problem in Appendix A.
Figure 2. Mesh evolution over time. The black and red lines represent the grid and the gas–liquid interface, respectively:
$l_{max}=11$
,
$l_{min}=9$
,
$\Delta p=1$
,
$L_p(t=0)=1$
;
${La}=100$
;
$\textit{Wi}=100$
;
$\mu_{S}=0.5$
;
$\lambda = 16$
;
$\epsilon =0.05$
.
4. Results and discussion
Considering the complex rheological characteristics of the liquid plug, this study employs the Saramito–HB model to build upon the rationalization of the effects of viscoelasticity and elastoviscoplasticity. The viscoelasticity is explored in § 4.1 by fixing
$\textit{Bi}=0$
and
$n=1$
, which turns the Saramito–HB model into the Oldroyd-B model. Investigating the viscoelastic liquid plug involves varying two critical dimensionless parameters: the Weissenberg number (
$\textit{Wi}$
) and the characteristic solvent concentration, which determine the relative dynamic viscosity ratio,
$\mu_{S}$
. In addition, the Laplace number
${La}$
is varied to explore the relevance of capillary effects on the non-Newtonian features. Then in § 4.2,
$\textit{Bi}$
and
$n$
are varied to investigate the effects of viscoplasticity when
$\textit{Wi}={O}(1)$
. In the final results and discussion section, § 4.3, we investigate the regime in which viscoelasticity and viscoplasticity are both significant.
4.1. Viscoelastic airway reopening
4.1.1. Comparison of the Newtonian and viscoelastic airway reopening models
In this section, we compare a viscoelastic and Newtonian airway reopening focusing on the wall pressure excursion over time. Figure 3(a) compares the viscoelastic (solid lines) and Newtonian (dashed lines)
$\Delta p_w=\max p_w(r=1)-\min p_w(r=1)$
for
$\mu =1.5 \times 10^{-3}, \rho =10^{-3}, \epsilon =0.05$
and three Laplace numbers:
$ La = 100$
(blue); 200 (green); 300 (red). For the Oldroyd-B model, parameters
$ \mu_{S} = 0.5$
and
$ We = 100$
are employed for figures 3 and 4. In figure 3(a), the wall pressure excursion peak consistently rises upon an increase of the Laplace number for the Newtonian model, i.e.
$\max (\Delta p_w)$
increases form around
$\max (\Delta p_w)=4.2$
for
${La}=100$
to
$\max (\Delta p_w)=5.6$
for
${La}=300$
in Newtonian condition. The peak of
$\Delta p_w$
normally occurs right after plug rupture, which is induced by the suddenly changing topology of the liquid film at the front meniscus because the plug bounces back to the wall (Hassan et al. Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011). Figure 3(b) shows the time evolution of the wall pressure (
$p_w$
) distribution and the gas–liquid interface corresponding to the red solid line in figure 3(a), which illustrates well the relationship between the peak of
$\Delta p_w$
and the liquid film topology. As shown in figure 3(b) at
$t=\left \{ 30, 60, 150 \right \}$
, the
$p_w$
highest gradient always occurs close to the front meniscus where the curvature changes sharply. Right after the plug rupture (
$t=212$
), the
$p_w$
distribution at the front plug is more intense due to the plug recoil towards the wall. This moment also corresponds to the
$\Delta p_w$
peak in figure 3(b). As the liquid film stops moving and the curvature of the front interface becomes smooth (
$t=280$
),
$p_w$
also becomes flat corresponding to
$\Delta p_w$
decaying at figure 3(a). Additionally, as shown in figure 4, the decrease of
$\epsilon$
leads to an increase in
$\Delta p_w$
consistently across both the constitutive models as expected by the scaling predicted by lubrication theory (Hori Reference Hori2006; Halpern et al. Reference Halpern, Fujioka and Grotberg2010).
Figure 3. (a) Comparison of the time evolution of the wall pressure difference for
${La}= \{100, 200, 300\}$
when
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
,
$\Delta p=1$
and
$\epsilon =0.05$
. Two constitutive models are used for the liquid phase: Newtonian (
$\textit{Wi}=0, \mu_{S}=0$
) and Oldroyd-B model (
$\textit{Wi}=100, \mu_{S}=0.5$
) (b) Time evolution of pressure distribution for
${La} = 300$
and Oldroyd-B model.
Figure 4. Comparison of the time evolution of the maximum wall pressure difference for
$\epsilon = \{0.05, 0.07, 0.10\}$
when
${La}=100$
,
$\mu =1.5 \times 10^{-3}$
and
$\rho =10^{-3}$
. Two constitutive models are used for the liquid phase: Newtonian (
$\textit{Wi}=0, \mu_{S}=0$
) and Oldroyd-B model (
$\textit{Wi}=100, \mu_{S}=0.5$
).
Notably, the plug will not rupture in all cases where
$\epsilon =0.10$
, which is attributed to a thick leading film slowing down the net mass loss of the plug, as reported in Fujioka et al. (Reference Fujioka, Takayama and Grotberg2008), Hassan et al. (Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011) and Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019). For all ruptured cases, the effect of the viscoelasticity speeds up the plug rupture approximately
$20\,\%$
compared with Newtonian fluids. As mentioned above, the viscoelastic properties of the fluid stem from the presence of dissolved polymers in the solvent, which exhibits most significant impact during plug propagation. Figures 3(a) and 4 illustrate the qualitative and quantitative disparity between Newtonian and viscoelastic conditions: in the Newtonian liquid, wall pressure experiences a plateau under all conditions. Conversely, when considering viscoelastic effects, the elastic response of the polymeric phase leads to an amplification of wall pressure excursion during propagation, which can reach a magnitude comparable to the rupture peak. This feature of the viscoelastic plug propagation is here reported for the first time and it deserves to be investigated further. In the subsequent sections, we elucidate a typical viscoelastic plug propagation and rupture scenario and conduct a parametric investigation across
${La}$
,
$\textit{Wi}$
and
$\mu_{S}$
to comprehensively understand and characterize the viscoelastic effect.
4.1.2. Analysis of viscoelastic effects for propagation and rupture of liquid plug
The airway reopening is driven by the pressure difference between the front and the rear air fingers. Increasing
${La}$
causes delayed rupture in capillary time units due to increased inertia. Decreasing the initial liquid film thickness
$\epsilon$
elevates wall pressure and shear stress as expected by the thin film theory. Moreover, augmenting the inlet-to-outlet pressure difference accelerates plug rupture and enhances wall stress. The above considerations have been extensively demonstrated using the Newtonian two-phase model (Bilek et al. Reference Bilek, Dee and Gaver2003; Kay et al. Reference Kay, Bilek, Dee and Gaver2004; Hassan et al. Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011; Muradoglu et al. Reference Muradoglu, Romanò, Fujioka and Grotberg2019). Given our emphasis on viscoelastic effects, we first focus on two cases: a low-(
$\textit{Wi}=10$
) and a high-(
$\textit{Wi}=100$
) Weissenberg number case in order to observe qualitatively distinct viscoelastic effects. The other parameters are kept as
$\mu_{S}=0.5$
,
${La}=300$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
and
$\rho =10^{-3}$
.
Figure 5. Comparison of plug propagation and airway reopening time dynamics. Here (a) and (b) correspond to
$\textit{Wi}=10$
and
$\textit{Wi}=100$
, respectively, and the contour plots of pressure and extra stress components (
$S_{\textit{zz}}$
,
$S_{\textit{zr}}$
and
$S_{\textit{yy}}$
) are shown at three times before the rupture and at one time after the rupture. The other non-dimensional groups are
${La}=300$
,
$\textit{Bi}=0$
,
$n=1$
,
$\mu_{S}=0.5$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
Figure 5 illustrates the pressure distribution
$p$
and the in-plane extra-stress components,
$ S_{rz}$
,
$ S_{rr}$
and
$ S_{\textit{zz}}$
at four times, three of them preceding and one following airway reopening for
$\textit{Wi}=10$
(figure 5
a) and
$\textit{Wi}=100$
(figure 5
b). The white area in the picture represents air. The fields from top to bottom are
$p$
,
$S_{\textit{zz}}$
,
$S_{rz}$
and
$S_{rr}$
, while the red rectangle marks the zoomed-in area used to identify significant qualitative difference. During the plug propagation, the tail meniscus thickens inducing a continuous mass loss from the plug bulk. Simultaneously, the plug squeezes the front meniscus, resulting in a minimum liquid film thickness where the maximum wall pressure drops and wall shear stress increases. Despite the presence of strong viscoelastic effects, the pressure
$p$
and interface evolutions exhibit similarities to those observed in Newtonian fluids (see Hassan et al. Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011; Muradoglu et al. Reference Muradoglu, Romanò, Fujioka and Grotberg2019).
For
$\textit{Wi}=10$
, the polymeric stress in the trailing liquid film is almost uniformly distributed. This results in a relatively constant liquid film thickness without significant fluctuations. The qualitative evolution of the interface near the plug shows a reduced amplitude of changes and lower extreme values compared with
$\textit{Wi}=100$
. In contrast, for
$\textit{Wi}=100$
, in the early status of plug propagation
$t=\left \{80,130 \right \}$
, the maximum values of
$S_{\textit{zz}}$
and
$S_{rr}$
always appear in a specific narrow area of the trailing meniscus and at the airway wall, respectively, and their size and distribution area increase with time. The minimum value of
$S_{rz}$
appears in the interface of trailing meniscus. Note that the maximum value of
$S_{rz}$
moves from the plug centre area to the front meniscus gas–liquid interface during the propagation process (from
$t=80$
to
$t=200$
). Up until about the rupture
$(t=200)$
, both
$S_{\textit{zz}}$
and
$S_{rr}$
increase, with
$S_{rz}$
experiencing an order of magnitude growth. The qualitative difference between
$S_{rr}$
and
$S_{\textit{zz}}$
is that
$S_{rr}$
develops a high extra stress area near the front and trailing meniscus interface while the stress concentration area of
$S_{\textit{zz}}$
is always close to the trailing meniscus and the wall. After rupture
$(t=240)$
, the three extra stresses decay rapidly. It is noteworthy that an enlarged view of the trailing liquid film exhibits a distinct concentration distribution of
$S_{rz}$
and
$S_{\textit{zz}}$
within the film. As the plug propagates, this distribution intensifies and then diminishes following rupture. This behaviour is primarily due to the significant stretching of the polymer at the interface and the wall, where the polymer responds to the accumulated elastic stress most intensely. This also explains the absence of stress concentration areas in the front liquid film where the fluid is compressed. Such polymer stress distribution also causes liquid film thickness fluctuations, as shown at
$t=200$
in the case of
$\textit{Wi}=100$
. For the comparison, the case with
$\textit{Wi}=10$
exhibits a similar qualitative evolution of the interface near the plug, while the trailing film thickness gradients and peaks are reduced. The qualitative difference in the stress peak between
$\textit{Wi}=10$
and
$\textit{Wi}=100$
predominantly occurs in the trailing liquid film rather than in the plug. Hence, the observed fluctuations in polymer stress and liquid film thickness in case of
$\textit{Wi}=100$
indicate further points towards an elastically driven dynamics during plug propagation, which will be explored in more detail in § 4.1.3.
4.1.3. Critical conditions for the kinematic stretching in the trailing film
In order to quantitatively characterize the differences reported in figure 5, we have to focus on the thin-film dynamics. During plug propagation, we consider the excursion of the interface extra-stress (
$\boldsymbol{S}_i$
, where the subscript
$i$
denotes the interface) in three components (
$S_{i,rz}$
,
$S_{i,rr}$
and
$S_{i,zz}$
). Figure 6 shows the two aforementioned cases to compare the low-
$\textit{Wi}$
(
$\textit{Wi}=10$
; figure 6
a) and the high-
$\textit{Wi}$
(
$\textit{Wi}=100$
; figure 6
b) cases at five instants in time. The black line in the figure represents the gas–liquid interface using the left-hand and bottom axes to represent the
$z$
and
$r$
coordinates, respectively. The three polymer stress components at the interface are shown as blue, yellow and magenta lines, using the right-hand axis to represent the extra stress value. As shown in figure 6(a), the distribution of extra stress consistently correlates with the interface curvature at all times, which increases smoothly near the liquid plug and approaches zero far from the plug. This is considered to be typical interfacial polymer stress behaviour for subcritical flow. For comparison, figure 6(b) illustrate the high-
$\textit{Wi}$
case. At the early stage of plug propagation (
$t=30$
), the three polymer stress components are induced by the interface curvature, and their respective peaks always appear at the meniscus shoulders which is qualitatively consistent with the low-
$\textit{Wi}$
case. However, as the plug propagates (
$t=60$
), the extra stress of the trailing liquid film significantly exceeds the extra stress of the leading film. This can be well understood when comparing our plug dynamics with Bothe et al. (Reference Bothe, Niethammer, Pilz and Brenn2022), who identify a similar behaviour in raising droplets and identified a critical droplet volume for the kinematic stretching. They showed that the time scales due to (i) the bubble propagation in a viscoelastic medium and (ii) the relaxation time of the polymer give raise to a qualitatively different behaviour upon overcoming the critical value for the Weissenberg number
$\textit{Wi}_c$
. In our case, when
$\textit{Wi} \lt We_c$
, the polymer stretches near the shoulder of the trailing meniscus, hence the viscoelastic stresses are released in the liquid plug (see figure 7
a). If
$\textit{Wi}\gt We_c$
, the stretched polymeric chains will release their stresses in the trailing film, giving rise to (i) a hoop stress that increases the film thickness and (ii) an axial stress that leads to a speed up of the bubble (in Bothe et al. Reference Bothe, Niethammer, Pilz and Brenn2022) or plug (see figure 7
b). This critical kinematic stretching condition occurs in the trailing liquid film, which is stretched; conversely, the front liquid film experiences primarily compression, hence it does not undergo the same critical stress-intensification dynamics. Additionally, the hoop stress produced in the supercritical case has an indirect impact on the wall pressure. In fact, the hoop stress produces a deformation of the trailing film interface, which in turn affects the capillary stresses owing to a change of curvature that ultimately feeds back on the wall pressure.
Figure 6. Time evolution of polymer extra stress along the interface in plug propagation and rupture for low-
$\textit{Wi}$
(a) (
$\textit{Wi}=10$
) and high-
$\textit{Wi}$
(b) (
$\textit{Wi}=100$
). The rest of non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\mu_{S}=0.5$
,
${La}=300$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
Figure 7. Schematic of subcritical (a) and supercritical (b) scenarios for viscoelastic solutions with low (a) and high (b) relaxation time, respectively. The red lines and arrows illustrate the potential polymer chains and the direction of polymer stress release, respectively.
In addition, figure 6(b) depicts localized stress oscillations caused by a dynamics due to the high
$\textit{Wi}$
that will be discussed in the following section. When
$t=150$
, there are multiple
$S_{i,zz}$
peaks and a precursor peak emerges at a position far away from the meniscus shoulder where curvature is smooth. This is the location where the stresses get released. As time progresses (
$t=200$
), the continued stretching of the polymer increases
$S_{i,zz}$
and the fluctuation region of
$S_{i,zz}$
expands, but remains localized to the thin film. Following plug rupture (
$t=280$
),
$S_{i,zz}$
amplitude diminishes progressively.
Based on the above observation, we identify the critical conditions for which such a qualitative behavioural change occurs as shown in figure 8. The left-hand and bottom axes represent the Weissenberg number (
$\textit{Wi}$
) and the polymer concentration (
$\mu_{S}$
) for each of our simulations, respectively. The colour of markers represents the supercritical (red) and subcritical (blue) conditions for
${La}=50$
(star),
${La}=100$
(circle),
${La}=200$
(square),
${La}=300$
(triangle) and
${La}=500$
(pentagon). Figure 8 shows that
$\textit{Wi}$
is the only important parameter affecting the qualitative behaviour of the viscoelastic dynamics, which is consistent with our interpretation based on the time scales of stretching and stress release. The reference critical conditions are solely produced by the polymeric relaxation behaviour and are not affected by Newtonian inertia (
${La}$
) and concentration (
$\mu_{S}$
). In addition, we anticipate that such a critical threshold (
$\textit{Wi}_c$
) persists even when different constitutive models are utilized for the extra-stress tensor (see § 4.3). To the best authors’ knowledge, this is the first time that critical viscoelastic behaviour for polymeric stress intensification has been reported for airway reopening.
Figure 8. Critical kinematic stretching diagram for the plug propagation elastoinertial relaxation. The left and bottom axes represent the Weissenberg number (
$\textit{Wi}$
) and the solvent-to-total dynamic viscosity ratio (
$\mu_{S}$
) of each simulation, respectively. The colour of markers represents the supercritical (red) and subcritical (blue) conditions for
${La}=50$
(star),
${La}=100$
(circles),
${La}=200$
(squares),
${La}=300$
(triangles) and
${La}=500$
(pentagon). The other non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
4.1.4. Resonance diagram
Based on the previous analysis, the supercritical viscoelastic behaviour is solely determined by threshold in Weissenberg number, i.e.
$\textit{Wi}_c$
. However, a comparative examination of figure 3(a) reveals that the amplitude of the wall pressure oscillation is also affected by the Laplace number (cf. solid lines for
$\textit{Wi}=100$
,
$\mu_{S}=0.5$
and different
${La}$
). For comparison, let us consider
$\textit{Wi}=100$
,
$\mu_{S}=0.5$
: the
$\Delta p_w$
for
${La}=\lbrace 200,\ 300\rbrace$
shows a flat evolution during plug propagation, while for
${La}=100$
we observed a pronounced peak in
$\Delta p_w$
even though they all reside in supercritical region of figure 8. The pressure oscillations observed in figure 3(a) exhibit remarkable qualitative and quantitative characteristics compared with other supercritical flow cases, hence they deserve a dedicated analysis.
To understand such a new dynamic behaviour, the time evolution of the wall polymer shear stress
$S_{w}$
distribution is compared for two Laplace numbers in figure 9. Figures 9(a) and 9(b) represent
${La}=300$
and
${La}=100$
, respectively, while the other simulation parameters are kept constant to
$\textit{Wi}=100$
,
$\mu_{S}=0.5$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
. The black line indicates the gas–liquid interface, while the blue line illustrates the wall polymer shear stress at five different times. Despite both cases being in the supercritical regime, their
$S_w$
distributions differ significantly. For case
${La}=300$
,
$S_w$
remains small without significant fluctuations. However, for the case
${La}=100$
, beginning at
$t=150$
,
$S_w$
exhibits pronounced oscillations in the trailing film. The oscillation amplitude initially increases, then decreases during the time interval from
$t=150$
to
$t=180$
. Here
$\textit{Wi}$
reflects how rapidly the polymer responds to stretching, indicating the elastic response of the fluid to an external force. Besides, the Laplace number pertains to the Newtonian inertia. Relating such dynamic to a classic spring–mass–damper forced system (see figure 10), it results the following dynamical analogy:
\begin{equation} \begin{matrix} m \ddot {x}& + & k_e x & + &\eta \dot {x} & = & F\\ \updownarrow & & \updownarrow & &\updownarrow &&\updownarrow \\ \rho \frac {\partial \boldsymbol{u}}{\partial t}& & \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{S} & &\mu \nabla^{2} \boldsymbol{u} && \frac {\sigma }{a^2}\sin (\frac {t}{\varLambda })\\ \updownarrow & & \updownarrow & &\updownarrow &&\updownarrow \\ \sim \rho \frac {\sigma ^2}{\mu ^2a}& & \sim \frac {\mu }{a\varLambda }& & \sim \frac {\sigma }{a^2} &&\sim \frac {\sigma }{a^2} \end{matrix} \end{equation}
Figure 9. Comparison of the time evolution of polymer wall stress in plug propagation and rupture for non-resonance case (a) (
$\textit{Wi}=100$
,
${La}=300$
) and resonance case (b) (
$\textit{Wi}=100$
,
${La}=100$
). The rest of non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\mu_{S}=0.5$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
Figure 10. Schematic of the mass–spring damper model.
where
$m$
,
$k_e$
,
$\eta$
and
$F$
represent the mass, spring constant, damping coefficient and external force of the fluid volume whose motion is being analysed,
$\varLambda$
is the time scale of the polymeric response that controls the supercritical behaviour, hence it enters the forcing term spring–mass–damper system. The last row of (4.1) denotes the scaling for each of the four terms, while
$x$
,
$\dot {x}$
and
$\ddot {x}$
are the displacement, velocity and acceleration of the mass, respectively. Carrying out a dimensional analysis, it yields
\begin{equation} m \sim \frac {\rho \sigma ^2}{\mu ^2a} \ddot {x}^{-1}\sim \frac {\rho \sigma ^2}{\mu ^2a} \left (\frac {\sigma ^2a}{\mu ^2a^2}\right )^{-1}=\rho , \qquad k_e \sim \frac {\mu }{a\varLambda } x^{-1} \sim \frac {\mu }{a\varLambda } \frac {1}{a} =\frac {\mu }{a^2\varLambda } . \end{equation}
Hence the natural frequency of the analogous spring–mass–damper system is
$\omega _n=\sqrt {k/m}\sim \sqrt {\mu /\rho \varLambda a^2}$
. Combining it with the forcing frequency coming from the polymeric stress, i.e.
$\omega _f=\varLambda ^{-1}$
, it yields that a resonance is expected for
\begin{equation} \frac {\omega _n}{\omega _f}\sim \sqrt {\frac {\mu }{\rho \varLambda a^2}}\frac {1}{\varLambda ^{-1}}=\sqrt {\frac {\mu \varLambda }{\rho a^2}}=\sqrt {\frac {\textit{Wi}}{{La}}} = 1. \end{equation}
To verify this prediction of our analogy, we compared several simulation results based on the average time derivative of wall extra stress,
$\overline {d_t \Delta S_w}=\delta \Delta S_w/ \delta T$
, where
$\delta \Delta S_w = [\max (\Delta S_w) - \Delta S_{w,rg}]$
is the wall extra stress excursion between the peak for
$\Delta S_w$
, i.e.
$\max (\Delta S_w)$
, and the beginning of the rapid growth for
$\Delta S_w$
, i.e.
$\Delta S_{w,rg}$
, while
$\Delta T$
denotes the corresponding time interval.
Figure 11 demonstrates the actual occurrence of a resonance in the trailing film. The left-hand and bottom axes represent
$\overline {d_t \Delta S_w}$
and
$\sqrt {\textit{Wi}/{La}}$
, respectively, for three polymer concentrations. The inset in figure 11 demonstrates how
$\overline {d_t \Delta S_w}$
is computed and a corresponding mesh-independence study is reported in Appendix B. The dotted lines connecting the markers are guide to the eyes used to show that the resonance is expected for
$\sqrt {\textit{Wi}/{La}}\sim 1$
, as predicted by our dimensional analysis. In fact, we notice that
$\overline {d_t \Delta S_w}$
decays rapidly for
$\sqrt {\textit{Wi}/{La}}\lesssim 1$
and
$\sqrt {\textit{Wi}/{La}}\gtrsim 1$
. We can then understand the oscillations observed on the wall extra stresses depicted in figure 9. Under resonant conditions, the forcing driven by viscoelastic supercritical conditions gets amplified, leading to significant stress intensification. This feeds back on the wall pressure, as noticed before, due to the indirect link between polymeric stresses, interfacial deformation and wall pressure. We further notice that the effect of the polymeric concentration on resonance intensity is monotonic as the intensity of
$\overline {d_t \Delta S_w}$
declines by increasing
$\mu_{S}$
. This is understood by considering that higher polymer concentrations feed back stronger elastic stresses, and hence the force in our analogy is stronger for lower
$\mu_{S}$
. From the perspective of wall stress evolution, the occurrence of such resonance phenomena increases the risk of epithelial cell damage, even without necessitating plug rupture.
Figure 11. Resonance diagram for the plug propagation under supercritical conditions. The black dashed line indicates the expected resonance conditions. The colour of markers represents the polymer concentration,
$\mu_{S}=0.25$
(red),
$\mu_{S}=0.5$
(blue) and
$\mu_{S}=0.9$
(green) for
${La}=50$
(pentagram),
${La}=100$
(circles),
${La}=200$
(squares),
${La}=300$
(triangles) and
${La}=500$
(pentagon). An example for calculating the resonance intensity (
$\overline {d_t \Delta S_w}$
) is reported in the inset. The other non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
4.1.5. Effect of the viscoelastic parameters:
$\textit{Wi}$
and
$\mu_{S}$
In this section we explore the effects of two viscoelastic parameters across three polymer concentrations (
$\mu_{S}$
= 0.25, 0.5 and 0.9) and six Weissenberg numbers (
$\textit{Wi}$
= 5, 10, 50, 100, 500 and 1000). We first vary the Weissenberg number while keeping solvent-to-total viscosity ratio
$\mu_{S}=0.5$
. Newtonian parameters are kept constant to isolate the influence of viscoelastic properties, specifically
${La}=100$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
. Figure 12 shows the time evolution of the plug length
$L_p$
excursion (figure 12
a), plug propagation velocity
$u_p$
(solid lines) and front interface velocity
$u_{f}$
(dashed–dotted lines) (figure 12
b), the wall normal stress
$\Delta p_w = \max [p(r=1)] - \min [p(r=1)]$
excursion (figure 12
c) and the local trailing film thickness
$\epsilon _r$
measured at
$z = z_p - 2$
(figure 12
d). By comparing figure 12(a) and figure 12(c), it is found that once the plug ruptures,
$\Delta p_w$
shows a sharp rise followed by a rapid decline, which is qualitatively in agreement with studies reported for Newtonian fluids (Hassan et al. Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011; Muradoglu et al. Reference Muradoglu, Romanò, Fujioka and Grotberg2019). This is also the case for
$\Delta \tau _w$
(not shown). The semitransparent orange band in figure 12(b) corresponds to velocity measured by Viola et al. (Reference Viola2024) in their Newtonian experiment. They used a microfluidic channel with semicircular cross-section of hydraulic radius
$a_{{exp}} = 0.23$
mm with driving pressure
$\Delta p_{{exp}} \approx 0.15$
psi (see plateau of figure 3 in Viola et al. Reference Viola2024). The close agreement between simulated
$\langle {u}_p\rangle$
and experimentally measured
$\langle {u}_{p,{exp}}\rangle$
average propagation velocities reported in figure 4(b) of Viola et al. (Reference Viola2024) confirms the model’s fidelity. Moreover, the comparison with the current results highlights the comparable order of magnitude of the non-dimensional plug propagation velocity, despite the differences in cross-section.
Figure 12. The effects of Weissenberg number (
$Wi$
). Time evolution of (a) plug length (
$L_p$
), (b) plug propagation velocity
$u_p$
(solid lines) and front interface velocity
$u_{f}$
(dashed–dotted lines), (c) wall pressure (
$\Delta p_w$
) and (d) local trailing film thickness (
$\epsilon _r$
). The orange semitransparent area is the experimental result reported in Viola et al. (Reference Viola2024) obtained with microfluidic equipment. The non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\mu_{S}=0.5$
,
$\textit{Wi}= \{ 5,10,50,100,500,1000 \}$
,
${La}=100$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
Figure 12(a) shows that, at small
$Wi$
(
$Wi \leqslant 10$
), viscoelastic plug ruptures later than Newtonian plug, but at large
$Wi$
(
$Wi \gt 50$
), viscoelasticity early than Newtonian, which is consistent with the analysis of subcritical and supercritical conditions in figure 7. Increasing
$\textit{Wi}$
accelerates the viscoelastic plug rupture. Two potential reasons can contribute to it: (i) an increase in plug propagation velocity and (ii) an increase in trailing liquid film thickness. Let us first consider the liquid plug velocity
$u_p$
(see figure 12
b). For
$u_p$
, the initial (
$t \lt 10$
) acceleration is almost the same regardless of
$\textit{Wi}$
. As the plug propagation progresses, the plug velocity reaches a plateau for
$\textit{Wi} \leqslant 10$
. This is typical also of Newtonian plug propagation as the plateau is controlled by the equilibrium between the pressure-difference driving and the net capillary resistance resulting from the two menisci (Bahrani et al. Reference Bahrani, Hamidouche, Moazzen, Seck, Duc, Muradoglu, Grotberg and Romanò2022; Viola et al. Reference Viola2024). Upon an increase of
$\textit{Wi}$
, we run into the supercritical regime, which severely affects the qualitative behaviour of
$u_p$
. The impact of the kinematic polymeric stretching for supercritical conditions tends to accelerate the plug because polymers stretched from the liquid plug back into the trailing film will feed back stress that pulls the plug towards the advancing propagation direction. Hence, upon an increase of
$\textit{Wi}$
, the wall-tangential polymeric stresses in the trailing film become stronger and stronger (see horizontal red arrows in figure 7
b), which explains why
$u_p$
grows with
$\textit{Wi}$
. This is in agreement with the results of Bothe et al. (Reference Bothe, Niethammer, Pilz and Brenn2022) and França et al. (Reference França, Jalaal and Oishi2024). The difference between the plug propagation velocity
$u_p$
(solid lines) and front interface velocity
$u_{f}$
(dashed–dotted lines) in figure 12 is connected to the plug draining rate, which affects the rupture time. In particular, one can reason in terms of the rear meniscus velocity
$u_r = 2u_p-u_{f}$
that increases significantly with the increase of
$Wi$
. This means that the drainage rate increases with
$\textit{Wi}$
, which effect is further enhanced by the increase of the film thickness at the rear meniscus
$\epsilon _r$
(Bahrani et al. Reference Bahrani, Hamidouche, Moazzen, Seck, Duc, Muradoglu, Grotberg and Romanò2022) due to the supercritical dynamics. To have a corresponding analytical description of the plug drainage rate
$-{\textrm d} L_p/{\textrm d}t$
, one can subtract the liquid fed to plug by advancing the front meniscus
$\epsilon _f(2-\epsilon _f)u_f$
to the one deposited by the rear meniscus
$\epsilon _r(2-\epsilon _r)u_r$
per unit cross section. It results that the drainage rate
$-{\textrm d} L_p/{\textrm d}t=\epsilon _r(2-\epsilon _r)u_r - \epsilon _f(2-\epsilon _f)u_f$
is non-trivially controlled by the variations of
$u_r$
and
$\epsilon _r$
, with an increase of
$u_r$
and
$\epsilon _r$
nonlinearly determining the plug drainage rate. A further confirmation of such an interpretation is well illustrated in figure 6, where
$S_{i,zz}$
is pulling the plug along the interface in the direction of propagation. For low
$\textit{Wi}$
, the polymer undergoes minimal stretching before the elastic stress is released in the liquid plug and the plug dynamics is quasi-Newtonian (see figure 6 case
$\textit{Wi}=10$
).
Considering now the local trailing film thickness
$\epsilon _r$
, we observe that when the hoop stress reaches a critical threshold, it induces localized wrinkling at the interface, as highlighted in figure 7(b) and as observed in figure 6(b). This effect is far more pronounced when the resonance occurring in the trailing film tends to intensify the hoop stress for
$\sqrt {\textit{Wi}/{La}}\sim 1$
, hence the film deforms (see purple line in figure 12
d for
$\textit{Wi} = 100$
and
${La} = 100$
). The combined effect of the polymeric axial pulling and the hoop-stress enhancing of the trailing film thickness explain the acceleration of the plug rupture depicted in figure 12(a).
We further stress that capillary effects on the film thickening are negligible in our case. According to the thin film theory (Aussillous & Quéré Reference Aussillous and Quéré2000), the liquid film thickness scale as
$\epsilon _r \approx {1.34 u_p^{2 / 3}\times (1+1.34 \times 2.5 u_p^{2 / 3})^{-1}}$
, where
$\epsilon _r$
is the trailing edge thickness, non-dimensionalized by the airway radius,
$u_p$
is employed to approximate the trailing meniscus dimensionless velocity, that in our scaling corresponds to the trailing meniscus capillary number
$Ca_r$
. Considering the
$u_p$
depicted in figure 12(b),
$u_p \in \sim [0.05,0.1]$
, hence
$\epsilon _r \in \sim [0.08,0.1]$
in a Newtonian plug. The impact of capillary effects on
$\epsilon _r$
is therefore expected to be negligible.
Finally, considering the wall pressure excursion
$\Delta p_w = \max [p(r=1)] - \min [p(r=1)]$
(see figure 12
c), its peak resulting from the plug rupture shows a weak dependency on
$\textit{Wi}$
, and uniformly decays to
$\Delta p_w \approx 2$
. This is qualitatively and quantitatively consistent with the results of Newtonian fluids, therefore, we can conclude that the dynamics after and immediately before rupture are dominated by Newtonian-like behaviour.
In figure 13, the wall-tangential stress excursion
$\Delta \tau _w$
(solid line) is depicted, together with its two components: (i) the Newtonian stress contribution due to the solvent
$\Delta \tau _w^N$
(dashed line), and (ii) the extra stress excursion
$\Delta S_w$
(dashed–dotted line) caused by the presence of polymers. Most of the non-dimensional groups are fixed, i.e.
${La}=100$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
, while we vary the two viscoelastic parameters
$\mu_{S}$
and
$\textit{Wi}$
.
Figure 13. Time evolution of tangential stress in plug propagation and rupture. The tangential stress (
$\Delta \tau _w$
, solid line) consists of Newtonian stress (
$\Delta \tau _w^N$
, dashed line) and polymer extra stress (
$\Delta S_w$
, dashed–dotted line). The non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\mu_{S}= \{ 0.25, 0.5, 0.9 \}$
,
$\textit{Wi}= \{ 5,10,50,100,500,1000 \}$
,
${La}=100$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
When
$\mu_{S}=0.9$
, the increase of Weissenberg number has negligible effects on rupture time and the characteristics of tangential stress. Besides, even as
$\textit{Wi}$
increases by three orders of magnitude,
$\Delta S_w$
(dashed–dotted point line) always accounts for less than
$1\,\%$
of total tangential stress (
$\Delta \tau _w$
), irrespective of plug propagation or rupture. This is well understood considering that the polymer concentration is very low, hence the plug behaviour does not get significantly impacted by viscoelasticity.
For
$\mu_{S}=0.5$
and
$\textit{Wi} \leqslant 10$
in figure 13, the viscoelastic stress is significantly increased during liquid plug propagation compared with
$\mu_{S}=0.9$
. This observation aligns with the understanding that, at sufficiently low Weissenberg numbers, the extra stresses in the Oldroyd-B model exhibit behaviour akin to that of a Newtonian fluid. Taking the limit for
$Wi \to 0$
in (2.3), the viscoelastic stress tensor converges to the Newtonian stress tensor, i.e.
$\boldsymbol{S}=\mu_{P}(\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla} ^T \boldsymbol{u})$
. Consequently, the ratio of
$\tau ^{N}_{w}$
and
$\Delta S_w$
becomes the viscosity ratio of the solvent to the polymer. This also means that the limiting proportion of
$\Delta S_w$
in the
$\Delta \tau _N$
is
$1-\mu_{S}$
(when
$\textit{Wi}\to 0$
). As
$\textit{Wi}$
increases, the polymer extra stress
$\Delta S_w$
decreases and the proportion of Newtonian stress
$\Delta \tau _w^N$
of solvent increases. To understand it, we need to consider the physical interpretation of the Weissenberg number. Since
$\textit{Wi}$
represents a non-dimensional measure of the polymer relaxation time, a high relaxation time indicates a delayed polymer response to flow stretching, allowing to build up elastic stress. As the Weissenberg number increases significantly, the left-hand side of (2.3) becomes dominant, while the right-hand side of (2.3) remains constant since the velocity gradient retains its order of magnitude. Consequently, the only way to balance the momentum in the thin-film equation is by reducing
$S_{w}$
for
$\textit{Wi}\gt 100$
. For
$\mu_{S}=0.25$
in figure 13, the
$\Delta S_w$
distribution first exceeds
$\tau _w^{N}$
for case
$\textit{Wi}=5$
and
$\textit{Wi}=10$
. Moreover, the accelerated rupture due to
$\textit{Wi}$
is amplified as the polymer concentration increases. When
$\textit{Wi}$
is held constant, higher polymer concentrations amplify the viscoelastic effect, causing an increase in
$\Delta S_w$
and accelerated rupture. Once the plug ruptures, the peak of wall pressure and tangential stress due to the rapid radial retraction of the plug fluid is mainly driven by Newtonian dynamics, so the peak of
$\Delta p_w$
and
$\Delta \tau _w$
are weak functions with respect to
$\textit{Wi}$
. This is consistent with the previous observations done for
$\mu_{S}=0.5$
.
By comparing the evolution of
$\tau _w^{N}$
and
$S_w$
, it is evident that the oscillation in total wall shear stress
$\tau _w$
is primarily driven by the elastic stress component. For instance, in the case where
$\mu_{S} = 0.25$
and
$\textit{Wi} = 100$
, this influence is particularly noticeable (see purple curve in figure 13
a). As discussed in § 4.1.4, such elastic stress oscillations are due to resonance and
$\Delta S_w$
can reach amplitudes comparable to the peak of
$\Delta \tau$
caused by the plug rupture.
4.1.6. Effect of the Laplace number
The effect of
${La}$
on plug propagation and airway reopening has been studied for the Newtonian case by Hassan et al. (Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011) and Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019), and they found that the increase in
${La}$
delays plug rupture and increases the pressure and shear stress levels on the wall. In viscoelastic scenarios, resonance emerges form the coupling effect of
${La}$
and
$\textit{Wi}$
. The effect of the Laplace number is investigated by varying
${La}$
form 50–500, for the Weissenberg number ranging from 10 to 1000 and for three different solvent-to-total liquid dynamic viscosity ratios, i.e.
$\mu_{S} =$
0.25, 0.5 and 0.9.
Figure 14. Time evolution of polymer extra stress in plug propagation and rupture. The non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\mu_{S}= \{ 0.25, 0.5, 0.9 \}$
,
$\textit{Wi}= \{ 10,100,1000 \}$
,
${La}= \{ 100, 200, 300 \}$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
The time evolution of the polymer extra stress
$\Delta S_w$
on the wall is shown in figure 14, in order to highlight the role of
${La}$
in viscoelastic fluids. For both low (
$\textit{Wi} = 10$
) and high (
$\textit{Wi} = 1000$
) Weissenberg numbers, the influence of
${La}$
on the magnitude and behaviour of polymeric stress is negligible across all polymer concentration levels. Under the conditions of medium Weissenberg number
$\textit{Wi}=100$
and high polymer concentration
$\mu_{S} \leqslant 0.5$
, Laplace number shows significant influence on
$\Delta S_w$
. For the case of
$\mu_{S}=0.25$
and
$\textit{Wi}=100$
(purple lines),
$\Delta S_w$
expresses a strong resonance-induced extra-stress intensification leading to amplitudes comparable to the total wall tangential stress in figure 13. With the increase of
${La}$
, the resonance phenomenon is first enhanced and then suppressed, and the peak of
$\Delta S_w$
also first increases and then decreases, until
${La}=500$
, when the resonance is no longer significant. This inhibitory effect is very clear in the middle panel of figure 14 for the purple lines that depict the
$\mu_{S}=0.5$
and
$\textit{Wi}=100$
case.
4.2. Viscoplastic airway reopening
The effect of viscoplasticity is investigated by varying Bingham number (
$\textit{Bi}$
) and power-law index (
$n$
) while keeping
$\textit{Wi}=\mathcal{O}(1)$
and
$\mu_{S}=0.5$
in the Saramito–HB model. As discussed in § 4.1, when
$\textit{Wi}$
is small enough (
$\textit{Wi} \lt 10$
), the viscoelastic fluid approximates the behaviour of a Newtonian fluid, as shown by the evolution of
$L_p$
,
$\Delta p_w$
and
$\Delta \tau _w$
across different
$\textit{Wi}$
in figures 12 and 13 when polymer concentration is low (
$\mu_{S} \geqslant 0.5$
) and
$\textit{Wi} \leqslant 10$
. Therefore, setting
$\textit{Wi}=5$
is sufficient to rule out significant qualitative and quantitative effects of viscoelasticity in order to investigate viscoplasticity while maintaining numerical stability in solving the Saramito–HB model equations. The other parameters are fixed as:
${La}=100$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
,
$\Delta p=1$
and
$\epsilon =0.05$
.
Figure 15. Comparison of the time evolution of the wall pressure for two constitutive models, Newtonian (dashed line,
$\textit{Wi}=0, \mu_{S}=0$
) and viscoplasticity (solid line,
$\textit{Wi}=5, \textit{Bi}=0.01, n=0.5, \mu_{S}=0.5$
), where
${La}= \{100, 200, 300\}$
.
4.2.1. Comparison of the Newtonian and viscoplastic airway reopening models
The comparison between the viscoplastic (solid lines) and the Newtonian (dashed lines) airway reopenings is presented in figure 15 in terms of the excursion of the wall pressure (
$\Delta p_w$
) as a function of time. Three different Laplace numbers are considered, i.e.
${La} = 100$
(blue), 200 (green) and 300 (red). For the Saramito–HB model, we use
$\textit{Bi}=0.01$
and
$n=0.5$
in figure 15. For both fluids, increasing the Laplace number results in an increase in peak wall pressure, rising from max(
$\Delta p_w$
)
$\approx 4.1$
for
${La}=100$
to max(
$\Delta p_w$
)
$\approx 5.5$
for
${La}=300$
. As discussed in § 4.1.1, this peak wall pressure arises from the abrupt topological change following the rupture of the plug and is primarily driven by surface tension, which becomes more dominant upon an increase of
${La}$
. Consequently, the max(
$\Delta p_w$
) is not very sensitive to the constitutive model. The primary distinction between viscoplastic and Newtonian airway reopening lies in the substantial delay caused by viscoplasticity, leading to a significant increase in rupture time
$t_r$
.
4.2.2. Analysis of the yield zone
Figure 16 shows how yielded (yellow) and unyielded (dark blue) zones evolve during the airway reopening. In each snapshot, the upper halves depict
$\textit{Bi}=0.1$
, and the lower halves
$\textit{Bi}=0.001$
, both of which are simulated for
$n=0.5$
. The yielded (unyielded) regions are separated based on the von Mises criterion given by
$\left | \boldsymbol{S}^d \right | \gt \textit{Bi}$
(yield region) and
$\left | \boldsymbol{S}^d \right | \lt \textit{Bi}$
(unyield region), where
$\left | \boldsymbol{S}^d \right |$
denotes the second norm of the extra stress. The gas–liquid interface is indicated by the magenta line. For
$\textit{Bi}=0.1$
, the yielded region becomes more prominent as the plug propagates. The corresponding location remains near the meniscus and the trailing film all across the plug propagation. This is consistent with the distribution of stress concentration zones shown in figure 5.
Figure 16. Time evolution of yielded/unyielded (yellow/dark blue) region for
$\textit{Bi}=0.1$
(top halves) and
$\textit{Bi}=0.001$
(bottom halves). The other non-dimensional groups are
$n=0.5$
,
$\textit{Wi}=5$
,
${La}=100$
,
$\mu_{S}=0.5$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
For
$\textit{Bi}=0.001$
, the plug region is always yielding except for a small patch in the centre of the plug for
$t\lesssim 10$
. This is due to the fact that
$\textit{Bi}=0.001$
implies a low yield stress, hence a low level of stress is enough to yield the liquid plug. On the contrary, for
$\textit{Bi}=0.1$
unyielded regions are persistent in the plug owing to the large yield stress. Besides, an increase in Bingham number always delays the propagation, such that the
$\textit{Bi}=0.1$
case did not rupture until the end of the airway model. This can be well understood due to the increased resistance resulting from unyielded region.
It is noteworthy that the rupture in the
$\textit{Bi}=0.001$
,
$n \leqslant 0.5$
case (
$t_r \geqslant 430$
) is significantly delayed compared with the Newtonian case (
$t_r \approx 200$
) shown in figure 15. This delay suggests that the power law index
$n$
, i.e. the shear-thinning behaviour, rather than the yield stress, that significantly influences the reopening of viscoplastic mucus, as the effect of
$\textit{Bi}$
is negligible because the mucus is yielded almost everywhere for
$\textit{Bi}=0.001$
. Thus, the next section will investigate the effects of
$\textit{Bi}$
and
$n$
by setting
$\textit{Bi}=0.001,\ 0.01$
and 0.1,
$n=$
0.3, 0.5 and 1.
4.2.3. Effect of the viscoplastic parameters:
$\textit{Bi}$
and
$n$
Figure 17 shows the time evolution of plug length
$L_p$
(figure 17
a), wall pressure
$\Delta p_w$
(figure 17
b) and wall extra stress
$\Delta S_w$
(figure 17
c) during the plug propagation up until rupture and beyond for three power-law indices
$n=0.3$
(green),
$n=0.5$
(blue) and
$n=1$
(red). Three
$\textit{Bi}$
are considered:
$\textit{Bi}=0.1$
(solid line),
$\textit{Bi}=0.01$
(dashed line) and
$\textit{Bi}=0.001$
(dashed–dotted line) for a total of nine viscoplastic fluid simulations compared with their Oldroyd-B counterpart (black dotted line,
$\textit{Bi}=0$
,
$n=1$
,
$\textit{Wi} = 5$
and
$\mu_{S}=0.5$
).
Figure 17. Time evolution of plug length
$L_p$
(a), wall pressure
$\Delta p_w$
(b) and wall tangential stress
$\Delta \tau _w$
(c) in plug propagation and rupture for three power-law index
$n=0.3$
(green),
$n=0.5$
(blue) and
$n=1$
(red). The yield stress is varying as
$\textit{Bi}=0.1$
(solid line),
$\textit{Bi}=0.01$
(dashed line) and
$\textit{Bi}=0.001$
(dashed–dotted line) for viscoplastic fluid. The dashed line represents the viscoelastic fluid using Oldroyd-B model where
$\textit{Bi}=0$
and
$n=1$
. The rest non-dimensional groups are
$\textit{Wi}=5$
,
${La}=100$
,
$\mu_{S}=0.5$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
In the evolution of
$L_p$
, the results clearly show that the rupture time increases as
$\textit{Bi}$
increases and
$n$
decreases, i.e. the stronger the viscoplastic characteristic, the more delayed the rupture. Figure 17 also includes four cases where rupture does not occur. They carry the strongest viscoplastic character of all, which include the three cases with
$n=0.3$
(green lines) and one case with
$\textit{Bi}=0.1, n=0.5$
(solid blue line). As expected, when
$n=1, \textit{Bi}=0.001$
(red dashed–dotted line),
$L_p$
converges to the result of Oldroyd-B as the Saramito–HB model converges to the Oldroyd-B for weakly viscoplastic parameters. Remarkably, upon a one-order-of-magnitude increase of
$\textit{Bi}$
, passing from
$\textit{Bi}=0.001$
to 0.01, the rupture time increases only slightly. This indicates that the rupture time in viscoplastic fluids is only weakly influenced by
$\textit{Bi}$
within physiological regimes.
It seems counterintuitive that a decrease of
$n$
leads to a longer rupture time, as decreasing
$n$
corresponds to a shear-thinning fluid. The interpretation of the results upon a variation of
$n$
deserves therefore an additional explanation. To single out the role of
$n$
, we employed a constant reference polymeric viscosity
$1-\mu_{S}$
that is here set to 0.5 for all simulations. In terms of the one-dimensional effective polymeric viscosity
$\mu_{\textit{eff}}$
versus shear rate
$\dot \gamma$
curve for an HB fluid (
$\mu_{\textit{eff}}=\textit{Bi}\ |\dot \gamma |^{-1}+K_c|\dot \gamma |^{n-1}$
), it means that all considered polymeric viscosities must intersect at
$\mu_{\textit{eff}}=1-\mu_{S} = 0.5$
. Figure 18(a) demonstrates the three flow indices
$n$
and three Bingham numbers
$\textit{Bi}$
. For the
$n=1$
case, if mucus is not yielded, the viscosity approaches infinity with an asymptote at
$\dot \gamma = \textit{Bi}$
making mucus more like a solid. On the other hand, when mucus is yielded for
$\dot \gamma \gt \textit{Bi}$
, the viscosity nearly stays constant at 0.5 which is our reference effective polymer viscosity. For
$n=0.5$
and even more for
$n=0.3$
,
$\mu_{\textit{eff}}$
remains at values significantly larger than 0.5 for
$\dot \gamma \lt 1$
(reference crossing point), even after yielding. Since the shear rate at the wall is always
$\dot \gamma \lt 1$
in our simulations (see dotted lines in figure 18
b), the effective viscosity
$\mu_{\textit{eff}}(n=0.3) \gt \mu_{\textit{eff}}(n=0.5) \gt \mu_{\textit{eff}}(n=1)$
leads the mucus to a thicker state, hence less prone to flow for lower
$n$
. Of course, the shear-thinning behaviour of the mucus is still consistently included in our simulations, as for
$\dot \gamma \gt 1$
, lower
$n$
would rather oppose a lower resistance (see
$\dot \gamma \gt 1$
in figure 18
a). Figures 18(b) and 18(c) show the distribution of shear rate
$\dot {\gamma }$
and tangential wall stress
$\tau _w$
when
$L_p=0.7$
. First, the shear rate near the wall reaches a maximum near the front plug meniscus and decays on both sides, which is consistent with what has been observed for a Newtonian plug (see Hassan et al. Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011; Muradoglu et al. Reference Muradoglu, Romanò, Fujioka and Grotberg2019). As shown in figure 18(b), low
$\textit{Bi}$
(
$\textit{Bi}=0.001$
and
$\textit{Bi}=0.01$
) has minimal influence on the distribution of
$\dot \gamma$
. Notably, only at
$\textit{Bi}=0.1$
the absolute value of
$\dot {\gamma }$
shows a minor decrease. At this value,
$\tau _w$
reaches its maximum, indicating increased resistance that delays rupture, aligning with the trends observed in figure 17. On the other hand, the effects of
$n$
are significant, as shown in figure 18(c). As
$n$
decreases,
$\dot {\gamma }$
also decreases while
$\tau _w$
increases. Consequently, the trailing film resistance rises upon an decrease of
$n$
, leading to a slower plug propagation. The reason for the delayed rupture time for
$n\lt 1$
is thus understood: shear thinning is always present, but the weak shear rate keeps the mucus in the high-effective-polymeric-viscosity region
$\dot {\gamma }\lt 1$
(see figure 18
a). This explains why a decrease in
$n$
leads to increased resistance.
Figure 18. One-dimensional analysis of the HB model: (a) effective viscosity; (b,c) distribution of
$\tau _w$
and
$\dot {\gamma }$
for different combinations of
$\textit{Bi}$
and
$n$
.
In the excursion of
$\Delta p_w$
and
$\Delta \tau _w$
depicted in figure 17, no significant differences are observed with respect to the viscoelastic case for the case
$\textit{Bi}=0.001, n=1$
(red dashed–dotted line). Besides,
$\max (\Delta p_w)$
and
$\max (\Delta \tau _w)$
caused by plug rupture exhibit weak sensitivity to the constitutive model, regardless of
$\textit{Bi}$
and
$n$
. This is consistent with the analysis in § 4.1, which demonstrates that the rupture behaviour is predominantly influenced by Newtonian properties, irrespective of rheology characteristics. It is worth noting that, the wall excursions
$\Delta \tau _w$
and
$\Delta p_w$
during plug propagation are a weak function of
$\textit{Bi}$
and
$n$
, once the time axis is rescaled by the rupture time.
4.3. Elastoviscoplastic airway reopening
Building upon the understanding developed in the previous sections, in this section we consider all the rheological characteristics of mucus, i.e. viscoelasticity, yield stress and shear thinning. This allows us to study airway reopening taking into account the complexity of the macroscopic rheological properties of mucus, including their interplay. The effect of elastoviscoplasticity is investigated by varying
$0.001\leqslant \textit{Bi} \leqslant 0.1$
,
$0.3\leqslant n \leqslant 1$
and
$50 \leqslant \textit{Wi} \leqslant 1000$
. Unless stated otherwise, the rest of the parameters are fixed as
${La}=100$
,
$\mu_{S}=0.5$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
,
$\Delta p=1$
and
$\epsilon =0.05$
.
4.3.1. Comparison of the viscoelastic and EVP models
The comparison between the Saramito–HB (EVP with
$\textit{Bi}=0.01$
,
$n=0.5$
, solid lines) and the Oldroyd-B (EVP with
$n=1$
and
$\textit{Bi} = 0$
, dashed lines) is presented in figure 19 focusing on the excursion of the wall pressure
$\Delta p_w$
as a function of time. All the results in figure 19 are obtained for three different Laplace numbers, i.e.
${La} = 100$
(blue), 200 (green) and 300 (red). For both constitutive models, increasing the Laplace number results in an increase in peak wall pressure, which follows the results of other constitutive models as it is due to a capillary effect. For
${La}=100$
, the EVP model exhibits similar fluctuations during propagation as the ones observed in the Oldroyd-B model. In addition, it can be found that the rupture times of the EVP model is consistently larger than that of Oldroyd-B at all
${La}$
, since it incorporates viscoplastic delaying.
Figure 19. Comparison of the time evolution of the wall pressure for two constitutive models, VE (dashed line,
$\textit{Wi}=100, \mu_{S}=0.5$
) and EVP (solid line,
$\textit{Wi}=100, \textit{Bi}=0.01, n=0.5, \mu_{S}=0.5$
), where
${La}= \{100, 200, 300\}$
.
4.3.2. Robustness of the elastic supercritical behaviour for EVP fluids
In order to investigate the robustness of the elastically driven supercritical behaviour found for the Olyroyd-B model, we consider it within the context of EVP. Figure 20 shows the critical kinematic stretching diagram for the EVP plug propagation elastoinertial relaxation. Critical conditions are determined using the same method as figure 8. Figure 20 demonstrates that the behaviour of all extra stresses aligns qualitatively with the viscoelastic case (see figure 6) at the same
$\textit{Wi}$
. This consistency underscores the similar stress response under comparable
$\textit{Wi}$
conditions, regardless of the yield stress (
$\textit{Bi}$
) and the strength of the shear thinning (
$n$
). We therefore conclude that elastic critical condition, due to
$\textit{Wi}$
, is robust upon the introduction of EVP effect. As a result, the boundaries of the stability diagram reported in figure 8 do not change significantly due to the inclusion of viscoplasticity, hence the conclusions derived for the Oldroyd-B model also apply to the EVP cases under the supercritical conditions.
Figure 20. Critical kinematic stretching diagram for the EVP plug propagation elastoinertial relaxation. The colour of markers represents the supercritical (red) and subcritical (blue) conditions for
$\textit{Bi}=0.001$
(circles),
$\textit{Bi}=0.01$
(squares) and
$\textit{Bi}=0.1$
(triangles). The other parameters are fixed as
${La}=100$
,
$\mu_{S}=0.5$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
,
$\Delta p=1$
and
$\epsilon =0.05$
.
4.3.3. Robustness of the resonance behaviour in EVP fluids
As mentioned in § 4.1.4, the resonance is triggered by elastically driven supercritical behaviours. In the previous section, we showed that critical elastic behaviours are robust upon a change of
$\textit{Bi}$
and
$n$
. In the following, we inquire if the resonance phenomenon is similarly insensitive to significant viscoplastic effects.
As already shown for the Oldroyd-B case in figure 11,
$\Delta S_w$
also oscillates at the trailing liquid film for
$\sqrt {\textit{Wi}/{La}} \approx 1$
, indicating a resonance. The corresponding figure that demonstrates the robustness of the resonance phenomenon in EVP fluids is figure 21, where the measurement of the resonance intensity (
$\overline {d_t \Delta S_w}$
) is using the same method as in figure 11. We see an amplification of the resonance intensity for
$\textit{Wi}/{La} \approx 1$
, while in both orientations away from
$\sqrt {\textit{Wi}/{La}} \approx 1$
, the intensity of the resonance decays, where the damping of resonance as also observed for the Oldroyd-B model (see figure 11). It is important to note that the yield stress and shear-thinning do not have a major impact on the onset of the resonance, despite quantitative differences in the amplification of the resonance intensity.
Figure 21. Resonance diagram for the EVP plug propagation under supercritical conditions. The black dashed line indicates the expected resonance conditions. The colour of markers represents the
$n=0.3$
(green),
$n=0.5$
(blue) and
$n=1.0$
(red) for
$\textit{Bi}=0.001$
(circles),
$\textit{Bi}=0.01$
(squares) and
$\textit{Bi}=0.1$
(triangles). The other parameters are fixed as
${La}=100$
,
$\mu_{S}=0.5$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
,
$\Delta p=1$
and
$\epsilon =0.05$
.
4.3.4. Effect of the EVP parameters:
$\textit{Wi}$
,
$\textit{Bi}$
and
$n$
The last parametric analysis carried out in this study aims to investigate the impact of viscoplastic effects using a significantly elastic parameter set, i.e. for supercritical conditions. Figure 22(a,b) show the time evolution of the plug length
$L_p$
and wall pressure
$\Delta p_w$
, respectively. In each panel, the values of
$\textit{Wi}$
are indicated, i.e.
$\textit{Wi}=50, 100, 500, 1000$
. For each
$\textit{Wi}$
, we carry out nine simulations considering all the combinations of the two viscoplastic parameters within the tuples
$(n,\ \textit{Bi}) = \lbrace 0.3,\ 0.5,\ 1\rbrace \times \lbrace 0.001,\ 0.01,\ 0.1\rbrace$
. The colour coding refers to a change of
$n$
, while the line style depicts a change of
$\textit{Bi}$
. The black dotted line represents the results obtained using the Oldroyd-B model.
Figure 22. Time evolution of plug length
$L_p$
(a) and wall pressure
$\Delta p_w$
(b). The colour of solid lines represents
$n=0.3$
(green),
$n=0.5$
(blue) and
$n=1$
(red), each of which corresponds to Bingham numbers,
$\textit{Bi}=0.1$
(solid line),
$\textit{Bi}=0.01$
(dashed line) and
$\textit{Bi}=0.001$
(dashed–dotted line). In each panel, the values of
$\textit{Wi}$
are indicated:
$\textit{Wi}=50,100,500, 1000$
. The black dotted line represents the viscoelastic fluid using Oldroyd-B model , while the other parameters are fixed at
${La}=100$
,
$\mu_{S}=0.5$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
As expected, the results of Saramito–HB model tend to converge to the Oldroyd-B curves for
$L_p$
and
$\Delta p_w$
as
$\textit{Bi}$
decreases at
$n = 1$
, which is consistent with the results in § 4.2. At moderate
$\textit{Wi}$
, i.e. for
$\textit{Wi} \leqslant 100$
, both
$L_p$
and
$\Delta p_w$
are affected by
$\textit{Bi}$
and
$n$
, even if the qualitative behaviours remain consistent. For large
$\textit{Wi}$
, the strong elasticity overshadows the viscoplastic effects, leading all results to overlap with the purely viscoelastic behaviour (see
$\textit{Wi} \geqslant 500$
, figure 22
a iii,iv, b iii,iv). By combining the rupture time results of § 4.2, we notice that
$\textit{Wi}$
has a more pronounced impact on rupture time compared with
$\textit{Bi}$
and
$n$
. Consistent with the viscoelastic results, increasing
$\textit{Wi}$
results in a faster rupture, while a higher yield stress and shear thinning lead to a modest delay in the EVP rupture. We further stress that this is in agreement with the interpretation given on the role of the elastic extra stresses under the supercritical conditions for large
$\textit{Wi}$
. We recall that the relaxation time associated with large
$\textit{Wi}$
is so high that the elastic stresses get released in the thin film far from the liquid plug, giving a net elastic pulling to the liquid plug. It results in an increase in the plug velocity that leads to an acceleration of the liquid plug rupture. As this effect relies neither on yield stress nor on the power-law index, the corresponding dynamics of
$L_p$
and
$\Delta p_w$
is a weak function of
$\textit{Bi}$
and
$n$
.
Let us focus now on the effects of
$\textit{Bi}$
and
$n$
on
$L_p$
and
$\Delta p_w$
for
$\textit{Wi} = 50$
. For
$L_p$
, the effect of
$\textit{Bi}$
on rupture time is more significant when
$n=1$
. On the other hand,
$n$
has always a significant impact on the plug dynamics. The major impact of
$n$
is consistent with the results presented in figure 18. The non-monotonic behaviour of
$L_p$
upon a change of
$n$
results therefore from the two opposite effects: (i) on the one hand, the effective viscosity opposes a higher trailing film resistance, while (ii) on the other hand, the plug deposits more liquid upon a decrease of
$n$
. Additionally, it is worth noting that the effect of the Bingham number seems to be significantly tamed down by elastic effects for
$\textit{Bi}\leqslant 0.01$
. This is due to the enhanced extra stresses due to the resonance amplification that unyields almost the whole liquid plug. Conversely, for
$\textit{Bi}\leqslant 0.1$
, a significant region of the liquid film and plug remains unyielded.
Figure 23. Time evolution of tangential stress (
$\Delta \tau _w$
, solid line), which consists of Newtonian stress (
$\Delta \tau _w^N$
, dashed line) and polymer extra stress (
$\Delta S_w$
, dashed–dotted line). The colour of lines represents
$\textit{Bi}=0.1$
(green),
$\textit{Bi}=0.01$
(blue) and
$\textit{Bi}=0.001$
(red). In each panel, the values of
$n$
are indicated:
$n=0.3, 0.5 ,1.0$
. The black line represents the viscoelastic fluid using Oldroyd-B model, while the other parameters are fixed at
$\textit{Wi}=100$
,
${La}=100$
,
$\mu_{S}=0.5$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
Figure 23 shows the time evolution of tangential stress (
$\Delta \tau _w$
, solid line), which consists of Newtonian stress (
$\Delta \tau _w^N$
, dashed line) and polymer extra stress (
$\Delta S_w$
, dashed–dotted line). The line colours represent
$\textit{Bi}=0.1$
(green),
$\textit{Bi}=0.01$
(blue) and
$\textit{Bi}=0.001$
(red). In each panel, the values of
$n$
are indicated:
$n=0.3, \ 0.5, \ 1.0$
, while the other parameters are fixed at
$\textit{Wi}=100$
and
${La}=100$
. The black dotted-line represents the viscoelastic fluid simulated using the Oldroyd-B model.
Once again, we confirm that the peak wall shear stress due to the rupture event is insensitive to
$\textit{Bi}$
and
$n$
, which is consistent with such a phenomenon being due to capillary stresses. For the cases
$n=0.3$
and
$n=0.5$
, the amplitude of
$\Delta S_w$
decreases compared with the purely viscoelastic case due to the presence of shear thinning. This indicates that shear thinning partially suppresses the resonance strength.
4.4. Considerations on the rupture time
To characterize the EVP effects on the rupture time
$t_r$
, we first consider the viscoelastic liquid plugs. Figure 24 shows the plug rupture time diagram that combines the effects of viscoelastic parameters (
$\textit{Wi}$
and
$\mu_{S}$
) and Laplace number (
${La}$
) in this study, where the dashed line represents the trend of rupture time for each
$\mu_{S}$
. As consistent with the results analysed previously, augmenting the Laplace number postpones rupture, which is consistent with the results and analysis of figure 3(a), as well as with the literature (Hassan et al. Reference Hassan, Uzgoren, Fujioka, Grotberg and Shyy2011; Muradoglu et al. Reference Muradoglu, Romanò, Fujioka and Grotberg2019). An increase of
$\mu_{S}$
expedites rupture across the whole range of Weissenberg numbers.
Figure 24. Rupture time diagram for all parameters in this study. The colour of markers represents the
$\mu_{S}=0.25$
(red),
$\mu_{S}=0.5$
(green) and
$\mu_{S}=0.9$
(blue) conditions for
${La}=100$
(circles),
${La}=200$
(squares) and
${La}=300$
(triangles). The other non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
Generally, increasing
$\textit{Wi}$
accelerates rupture. Except for all the resonance cases with
$\mu_{S}=0.25$
, their rupture time decreases significantly, leading to a minimum in the rupture time (red dash lines). This is attributed to the intensification of the extra stresses as shown in figure 13. Finally, at low-enough Weissenberg number (
$\textit{Wi} \leqslant 10$
), the effects of viscoelastic parameters (
$\textit{Wi}$
and
$\mu_{S}$
) on the plug rupture time are not significant. We further remark that under certain conditions (
$\textit{Wi}=1000$
,
$\mu_{S}=0.25$
,
${La} \geqslant 200$
), the plug does not rupture until reaching the domain outlet, and the corresponding rupture time was estimated by extrapolating the time evolution curve of the plug length.
Figure 25. Rupture time diagram for EVP airway reopening. The colour of markers represents the
$n=0.3$
(green),
$n=0.5$
(blue) and
$n=1.0$
(red) conditions for
$\textit{Bi}=0.001$
(circles),
$\textit{Bi}=0.01$
(squares) and
$\textit{Bi}=0.1$
(triangles). The other non-dimensional groups are
${La}=100$
,
$\mu_{S}=0.5$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
Including now the viscoplastic properties of the liquid, the rupture time is computed upon a variation of
$\textit{Bi}$
(different markers),
$n$
(different colours) and
$\textit{Wi}$
(abscissa) as shown in figure 25. First,
$t_r$
decreases significantly as
$\textit{Wi}$
increases, and eventually converges to
$t_r \approx 150$
. At a low Weissenberg number (
$\textit{Wi} \leqslant 100$
),
$\textit{Bi}$
and
$n$
have a non-negligible impact on
$t_r$
, mainly manifested in that increasing
$\textit{Bi}$
and decreasing
$n$
delay the rupture. We remark that under certain conditions the plug does not rupture until reaching the domain outlet. This is observed for
$\textit{Wi}=5$
and
$n\leqslant 0.3$
, as well as for
$\textit{Wi}=5$
,
$n=0.5$
and
$\textit{Bi}=0.1$
. The reason has been analysed in 4.2 and is consistent with the trend in figure 25, so they are not shown here. At large Weissenberg number (
$\textit{Wi} \geqslant 500$
), the influence of
$\textit{Bi}$
and
$n$
is negligible due to the strong viscoelasticity.
4.5. Effects of surfactant
The role of surfactant in airway reopening is a subject of critical interest due to its significant implications for pulmonary mechanics, particularly under pathological conditions. Surfactant’s effect on the dynamics of airway reopening not only modifies surface tension, but also introduces interfacial solutocapillary Marangoni stresses, and it can hinder or even preclude the reopening process. The effect of surfactant on Newtonian airway reopening has been thoroughly studied by Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019), who demonstrated that the presence of surfactant leads to substantial deviations from clean cases. They identified the surfactant strength and the interfacial surfactant concentration ratio as the two non-dimensional groups that most severely affect the plug propagation and rupture dynamics. They also pointed out that a change of orders of magnitude in adsorption and desorption coefficients, penetration depth as well as bulk and interfacial Schmidt numbers, has a minor impact on airway reopening. This leads to the conclusion that considering the dynamics of insoluble surfactant and limiting our study to the evolution of the surfactant distribution for different initial surfactant concentrations and different surfactant strengths would be enough to make a conclusive statement on the interaction between the interfacial solutocapillary stresses and the bulk rheological mechanisms studied so far in this work. This subsection underscores, therefore, the need to account for the complex interaction between surfactant dynamics (modelled as insoluble) and the complex elastic supercritical behaviour identified for the first time in this study.
Following Farsoiya et al. (Reference Farsoiya, Popinet, Stone and Deike2024), we employ the module already implemented and tested in Basilisk that solves the governing equation derived by Stone (Reference Stone1990) using a level-set method. Insoluble interfacial surfactant concentration evolves by
\begin{equation} \cfrac {\partial \varGamma }{\partial t} + \boldsymbol{\nabla}_{s} \boldsymbol{\cdot }\left (\varGamma \boldsymbol{u}_{s} \right ) + \varGamma \left (\boldsymbol{\nabla}_{s}\boldsymbol{\cdot }\boldsymbol{n}\right )\left (\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\right ) = \cfrac {1}{{Pe}}{\nabla} ^{2}_{s}\varGamma , \end{equation}
where
$\varGamma$
denotes the interfacial surfactant ratio (normalized with the maximum equilibrium concentration) and the subscript ‘
$s$
’ denotes that the nabla operator is computed on the interface. The Péclet number (
${Pe}$
) is obtained by the product of the Laplace number and the surface Schmidt number, and is set to
$10^4$
in the following (Muradoglu et al. Reference Muradoglu, Romanò, Fujioka and Grotberg2019). As
${Pe}\gg 1$
, the impact of surfactant diffusion is assumed to be negligible with respect to the surfactant convective transport along the interface.
The surface tension variation due to the local surfactant concentration is modelled by the a linear equation of state in the form
$\sigma = \sigma _0(1 - \beta \varGamma )$
, where
$\sigma _0$
is the reference surface tension and is also used to evaluate the non-dimensional parameters
${La}$
,
$\textit{Bi}$
and
$\Delta p$
, while
$\beta$
denotes the surfacant strength. The coupling between the surfactant dynamics and the momentum equation is done at the interface via the variation of the surface tension term
$\chi \boldsymbol{n} \delta_{s}$
, which is now premultiplied by the local-to-reference surface tension ratio
$\sigma /\sigma _0$
. Moreover, the Marangoni stresses are included as explained in Farsoiya et al. (Reference Farsoiya, Popinet, Stone and Deike2024).
A parametric study is carried out to investigate the effect of initial interfacial surfactant concentrations and surfactant strengths on the boundary of the critical elastic behaviour identified for the liquid plug. We select
${La} = 100$
and test the interaction between the dynamics of the Oldroyd-B model and the surfactant for subcritical conditions (
$\textit{Wi} = 10$
) and supercritical (
$\textit{Wi} = 100$
) conditions. Four initial concentrations are considered, setting a uniform interfacial surfactant distribution equal to
$\varGamma (t=0) = 0$
(clean case), 0.005, 0.01 and 0.015 and keeping the surfactant strength constant at
$\beta = 0.7$
(in agreement with the baseline case selected by Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019)). In addition, we tested the strength of the surfactant for three different values of
$\beta$
, that is,
$\beta = 0$
(clean case), 0.2 and 0.7, while keeping the initial surfactant distribution constant at
$\varGamma (t=0) = 0.01$
.
Figure 26 depicts the location of interface, the extra-stress component
$S_{i,zz}$
evaluated at the interface, and the interfacial surfactant concentration ratio
$\varGamma$
for five snapshots in time and the several parameters tested in this subsection. In order to facilitate the comparison, we compare different surfactant parameters for constant values of
$L_p(t)$
. Figures 26(a i) and 26(b i) depict the results for
$\varGamma$
, while figures 26(a i) and 26(b i) refer to
$S_{i,zz}$
. The surfactant distribution at the interface is almost always uniform, except for the shoulder of the front meniscus where Marangoni stress contributes to a rigidification of the interface that affects the local curvature. This is in agreement with the results of Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019) and further confirms that it is not required to include bulk surfactant dynamics to obtain a leading-order approximation of
$\varGamma$
. We further observe that the presence of surfactant has a visible effect on the trailing film thickness, which becomes smoother thanks to the interface rigidification effect. Most importantly, we stress that the key observations reported in our work for the first time are robust: the same critical viscoelastic behaviour is observed in the thin film between
$\textit{Wi} = 10$
and
$\textit{Wi} = 100$
regardless of the initial surfactant concentration
$\varGamma (t=0)$
and of the surfactant strength
$\beta$
. Therefore, the presence of surfactant does not qualitatively affect the plug dynamics, even if it impacts the distribution of extra stress, yet it does not stabilizes the thin-film elastically driven supercritical behaviour found for clean cases.
Figure 26. Time evolution of surfactant concentration (i) and polymer extra stress (ii) along the interface for
$\textit{Wi}=10$
(a) and
$\textit{Wi}=100$
(b). The rest of non-dimensional groups are
$\textit{Bi}=0$
,
$n=1$
,
$\mu_{S}=0.5$
,
${La}=100$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
Figure 27(a,b) show the time evolution of the plug length
$L_p$
and wall pressure
$\Delta p_w$
, respectively. In particular, figure 27(a) shows that a significant plug rupture delay is due to interfacial contamination, with
$t_r$
increasing upon an increase of
$\varGamma (t=0)$
and
$\beta$
, in agreement with the study of Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019). The impact of the surfactant is far more significant when subcritical conditions are considered. This is consistent with what was reported for the EVP parametric study, where we showed that the elastic supercritical behaviour at the core of this study determines the dominant dynamics. Moreover, the resonance phenomenon identified for clean interfaces seems to be persistent, even if occurring at earlier times. This can be appreciated by considering the extra-stress excursions for supercritical conditions at
$\sqrt {\textit{Wi}/{La}} \approx 1$
. We stress, however, that the presence of surfactant renders the analysis more complex as the surface tension used to determine the reference Weissenberg and Laplace numbers is an overestimate of the average
$\sigma$
. Finally, the major impact of the surfactant can be observed on the pressure excursion at the wall
$\Delta p_w = \max [p(r=1)] - \min [p(r=1)]$
. When contaminated surfaces are considered, the impact of the supercritical elastic phenomenon on the wall pressure becomes negligible. This is well understood considering figure 26. The rigidification of the rear-film interface led to negligible
$\partial _z \epsilon _r$
, hence to negligible pressure gradients in the thin film. As a consequence, the presence of surfactant breaks the link between the resonance and the wall pressure, resulting in an effective annihilation of the wall pressure oscillations.
Figure 27. Time evolution of plug length
$L_p$
, wall pressure
$\Delta p_w$
and wall elastic stress
$\Delta S_w$
for
$\textit{Wi}=10$
(a) and
$\textit{Wi}=100$
(b). The coloured solid lines represent different surfactant concentrations and strengths, and the black dashed line represents the cleaning case, while the other parameters are fixed at
$\textit{Wi}=100$
,
${La}=100$
,
$\mu_{S}=0.5$
,
$\textit{Bi}=0$
,
$n=1$
,
$\Delta p=1$
,
$\mu =1.5 \times 10^{-3}$
,
$\rho =10^{-3}$
and
$\lambda =16$
.
5. Conclusions
The non-Newtonian effects on the liquid plug rupture have been investigated by computationally studying an airway reopening model, where a pressure-driven liquid plug ruptures inside a rigid pipe. The pressure difference induced by respiration pushes the liquid plug that tends to deposit a trailing film thicker than the leading film. This leads to effective drainage of the liquid plug, which eventually leads to the plug rupture. The simulation parameters are chosen to mimic physiologically relevant regimes for the ninth or tenth generation of a typical adult lung, where airway diameter is at millimetre scale and gravitational effects are negligible. Non-Newtonian effects are modelled using the Saramito–HB constitutive law. Elastoviscoplasticity has been thoroughly investigated exploring a parametric space that includes shear-thinning (
$n\in \lbrace 0.3,\ 0.5,\ 1\rbrace$
) viscoelastic (
$\textit{Wi}\in [0,\ 1000]$
), polymeric viscosity (
$1-\mu_{S}\in [0,\ 0.75]$
) and yield-stress (
$\textit{Bi}\in [0,\ 0.1]$
) effects. All parameter ranges have been selected to represent the healthy and the pathological conditions.
This work has been designed with the main goal of identifying the impact of the viscoelasticity and viscoplasticity, at first, and then their non-trivial interplays. Hence, we structured this paper into three main parts: (i) the first part considered a comparison between Newtonian and viscoelastic airway reopening; (ii) the second part compared Newtonian and viscoplastic plug rupture; (iii) the last part combined viscoelastic and viscoplastic effects into the EVP Saramito–HB model and compared the plug propagation and rupture with (i) and (ii) to identify which standalone features are robust and what significant leading-order effects are due to the interplay between viscoelasticity and viscoplasticity. Finally, the robustness of the analysis with respect to surfactant effects at the air–mucus interface is verified.
For viscoelastic airway reopening, Oldroyd-B model is used by setting
$\textit{Bi}=0$
and
$n=1$
inside of Saramito–HB model. We identified an elastically driven supercritical behaviour that is reported here for the first time. The critical kinematic stretching diagram reported in figure 8 showed that
$\textit{Wi}$
is the only parameter affecting the onset of critical conditions, whereas
${La}$
and
$\mu_{S}$
are not found to be significant for the identification of a critical threshold. This led us to conclude that the mechanism observed in our system is similar to the one reported for droplets raising in viscoelastic media (see Bothe et al. Reference Bothe, Niethammer, Pilz and Brenn2022). Furthermore, at supercritical conditions, a resonance can be triggered, leading to an amplification of the extra stress oscillations in the trailing film, where the polymer is stretched along the streamlines. We proposed a mechanical analogy with a forced mass–spring–damper system and deduced that the characteristic frequency ratio is represented by
$\sqrt {\textit{Wi}/{La}}$
. A verification of such an analogy has been reported in the resonance diagram of figure 11, which demonstrates that the resonance occurs when
$\sqrt {\textit{Wi}/{La}}\approx 1$
and is robust upon a variation of
$\textit{Wi}$
,
${La}$
and
$\mu_{S}$
. For a resonant airway reopening, the corresponding
$\Delta p_w$
and
$\Delta \tau _w$
values are found to be of comparable amplitude with the one reported at liquid plug closure (Romanò et al. Reference Romanò, Muradoglu, Fujioka and Grotberg2021). This stresses the significance of our study in physiological terms: the new elastocapillary mechanism at the core of the trailing film resonance is potentially as damaging for the epithelium as the fast capillary recoil of the interface upon plug rupture. Moreover, as the resonance occurs during plug propagation, our study further points out that such a novel elastocapillary mechanism increases the risk of epithelial cell damage regardless of the occurrence of the plug rupture.
We also have observed that viscoelasticity accelerates plug rupture by approximately
$20\,\%$
in comparison with Newtonian plug, i.e.
$t_r\downarrow$
for
$\textit{Wi} \uparrow$
and
$\mu_{S} \downarrow$
. A corresponding interpretation of the results has been given in terms of the release of elastic stresses in the trailing film, that tends to accelerate plug propagation upon an increase of
$\textit{Wi}$
and
$\mu_{S}$
. Combining the increase of
$t_r$
due to viscoelastic effects with a consistent result reported by Romanò et al. (Reference Romanò, Muradoglu, Fujioka and Grotberg2021) for viscoelastic plug formation, a physiologically relevant conclusion follows: viscoelasticity promotes the formation and the rupture of liquid plugs. We therefore speculate that the increased frequency resulting from shorter closure and rupture times of viscoelastic fluids could increase the risk of damaging the epithelium.
The effect of viscoplasticity is investigated by varying Bingham number (
$\textit{Bi}$
) and power-law index (
$n$
) while keeping
$\textit{Wi}={O}(1)$
and
$\mu_{S}=0.5$
within the Saramito–HB model. The primary distinction between viscoplastic and weakly viscoelastic (for
$\textit{Wi}=5$
and
$\mu_{S}=0.5$
) airway reopening lies in the significant increase in rupture time resulting form an increase of yield stress
$\textit{Bi}$
and a decrease of power-law index
$n$
. Increasing
$\textit{Bi}$
leads to increasing rupture time by increasing the resistance resulting from the unyielded region. The increase of
$t_r$
due to a decrease of
$n$
is rather due to the low level of non-dimensional shear rate in the trailing film (
$\dot \gamma \lt 1$
). As our comparison aims to single out the shear-thinning effect by keeping the equivalent reference polymeric viscosity fixed, i.e.
$\mu_{\textit{eff}} = 0.5$
at
$\dot \gamma =1$
, the flow regions at shear stress lower than the characteristic capillary stress (
$\sigma /a$
) have higher effective viscosity for lower
$n$
(see figure 18), hence they oppose a higher resistance that does not even out with the regions at higher stress. Apart from the effect on the reopening time, no other major effects are found for viscoplastic airway reopening, not even on the wall stresses.
The third part of our study includes all relevant rheological properties of the Saramito–HB model in the EVP liquid to corroborate the previous findings from viscoelastic and viscoplastic airway reopening. Our EVP simulations confirm that the elastic supercritical behaviour, as well as the resonance mechanisms identified for the Oldroyd-B model persist within the EVP, leading to an essentially unchanged critical threshold (
$\textit{Wi}_c\in [10,\ 50]$
) and resonance diagram trends. As previously discussed, viscoelasticity facilitates rupture while viscoplasticity impedes it. Their interplay leads to intermediate rupture times for the EVP up until viscoelastic effects do not become dominant. If
$\textit{Wi}\geqslant 500$
, the EVP plug rupture can be approximated to a purely viscoelastic rupture in terms of
$L_p$
,
$\Delta \tau _w$
and
$\Delta p_w$
, regardless of the persistence of unyielded regions. This could have important implications in future investigations, of either numerical or experimental kind, as it represents, de facto, a significant reduction of the non-dimensional parametric space. Hence, for accurately predicting physiological indicators such as
$L_p(t)$
,
$\Delta \tau _w$
and
$\Delta p_w$
, viscoplastic effects can be neglected if
$\textit{Wi}$
is large enough.
The last part of our investigation included the effect of surfactant. A previous study on their impact on a Newtonian airway reopening has been carried out by Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019), hence we focused on the interplay between the surfactant interfacial dynamics and the novel supercritical phenomena identified in this work. We showed that the supercritical elastic behaviour unravelled in our study is robust upon the increasing values of initial surfactant concentration ratios
$\varGamma (t=0)\in \lbrace 0,\ 0.005,\ 0.01,\ 0.015\rbrace$
and for different surfactant strengths
$\beta \in \lbrace 0,\ 0.2,\ 0.7\rbrace$
. We further confirmed that the presence of surfactant significantly delays rupture, especially for subcritical conditions, in agreement with the study of Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019).
These conclusive considerations indicate that the EVP airway reopening model presented in this study is capable of capturing additional physical phenomena that have not been previously documented. Future research will focus on refining the airway reopening model by incorporating two-layer coating and deformable wall. Moreover, as we unravelled a purely elastic supercritical behaviour, future investigations are in plan for extending the present viscoelastic model and include finite polymer extensibility for studying the entire range of
$\mathcal{L}\lt 150$
(e.g. by the FENE-P (finitely extensible nonlinear elastic – Peterlin) model), as well as high-concentration shear-thinning viscoelastic effects (e.g. by a multimode Giesekus model).
Funding
This work is supported by a scholarship from the China Scholarship Council (CSC) under grant no. 202106240007. We acknowledge financial support from the Scientific and Technical Research Council of Türkiye (TUBITAK) (grant no. 124M335) and the Research Council of Finland (grant no. 354620).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation
The implementation of the Newtonian solver has already been validated successfully for the airway closure problem (Romanò et al. Reference Romanò, Fujioka, Muradoglu and Grotberg2019, Reference Romanò, Muradoglu, Fujioka and Grotberg2021), besides Romanò et al. (Reference Romanò, Muradoglu, Fujioka and Grotberg2021) also verified the viscoelastic model with solver of Izbassarov & Muradoglu (Reference Izbassarov and Muradoglu2015) using the Oldroy-B constitutive law. On the other hand, the viscoplastic model has been verified by Deka, Pierson & Soares (Reference Deka, Pierson and Soares2019), Deka, Pierson & Soares (Reference Deka, Pierson and Soares2020) and Deoclecio, Soares & Popinet (Reference Deoclecio, Soares and Popinet2023), who used Bingham materials. Following the approach of Deoclecio et al. (Reference Deoclecio, Soares and Popinet2023), we use the HB model to evolve the polymeric viscosity based on the Oldroyd-B model, i.e. we end up with Saramito–HB material. In order to verify the applicability of the numerical framework to the dynamics of airway reopening, the process of plug propagation and rupture was first verified in the absence of polymer, i.e.
$\mu_{S}=1$
. Figure 28 shows the comparison of the results of this study with Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019). Figures 28(a) and 28(b), respectively, represent the evolution of the plug length and wall pressure
$\Delta p_w$
excursion over time. The flow parameters are set to
$\Delta p= \{1.5,2 \}$
as driving pressure difference,
$\lambda = 16$
,
$L_p(t=0)=1.0$
as initial plug length, Laplace number
${La}=100$
and dimensionless initial film thickness
$\epsilon =0.05$
. The present
$\Delta p_w$
results are shifted to align with the reference pressure in Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019). Despite the different numerical approaches (VOF/finite-volumes for basilisk and front tracking/finite-difference for Muradoglu et al. (Reference Muradoglu, Romanò, Fujioka and Grotberg2019)), our results compare well with the literature.
Figure 28. Validation for the Newtonian case. Evolution of plug length
$L_p$
(a) and wall pressure
$\Delta p_{w}$
(b) under the conditions of driving pressure difference
$\Delta p=1.5, 2$
,
$L_p(t=0)=1$
,
${La}=100$
,
$\lambda = 16$
and
$\epsilon =0.05$
.
Appendix B. Mesh convergence
Besides the validation with the literature results, we carry out a mesh independence study by varying the level of adaptive mesh refinement (Popinet Reference Popinet2009, Reference Popinet2015). The grid refinement level
$l$
leads to the radial number of cells
$N_r=2^{l}/\lambda$
and an axial number of cells
$N_z=2^{l}$
. The maximum mesh refinement level
$l_{max}$
is used to ensure that the grid does not become overly fine, while the minimum refinement level
$l_{min}$
is ensuring that the grid remains sufficiently resolved even in coarse regions. We tested
$l_{max}=\left \{9,10,11,12 \right \}$
respectively, while maintaining
$l_{min}=9$
across all cases for grid independence investigations. Figure 29 shows the results of the grid independence study in terms of the evolution of the plug length
$L_p(t)$
and the maximum wall pressure excursion over time
$\Delta p_w = \max p(r=1) - \min p(r=1)$
. Considering the negligible difference between the results for
$l_{max}=11$
and
$l_{max}=12$
, both in plug length, rupture time and maximum wall pressure,
$l_{max}=11$
is chosen for this study to balance accuracy and computational cost-effectively.
Figure 29. Grid convergence. Evolution of plug length
$L_p$
(a) and wall pressure excursion
$\Delta p_w = \max [p(r=1)] - \min [p(r=1)]$
(b) under the conditions of
$\Delta p=1.5$
,
$L_p(t=0)=1$
,
${La}=100$
,
$\textit{Wi}=10$
,
$\mu_{S}=0.5$
,
$\lambda = 16$
and
$\epsilon =0.05$
.
Figure 30. Grid convergence with resonance. Evolution of plug length
$lp$
(a) and wall pressure excursion
$\Delta p_w = \max [p(r=1)] - \min [p(r=1)]$
(b) under the conditions of
$\Delta p=1$
,
$L_p(t=0)=1$
,
${La}=100$
,
$\textit{Wi}=100$
,
$\mu_{S}=0.25$
,
$\lambda = 16$
and
$\epsilon =0.05$
.
Figure 30 depicts a mesh convergence study for the resonance phenomenon discussed in § 4.1.4. Our convergence study shows that resolved-enough simulations allow for the fine structure of the stresses to develop; hence they allow to capture the resonance phenomenon. We further stress that the qualitative resonant behaviour is grid independent starting from
$l_{{max}} = 11$
, hence we consider it as a robust feature of our simulations. Owing to the sensitivity of such a phenomenon to the initial stress growing conditions, we cannot expect a very good quantitative grid convergence. However, our results are used to provide an order of magnitude of the stresses, and the accuracy granted by the selected grid (
$l_{{max}} = 11$
) is surely enough to draw all the conclusions reported in our study.
Open access
