The reopening of an occluded airway can lead to the formation of droplets and aerosols, which can be released during exhalation, providing a possible mechanism of disease transmission. In this study, the flow behavior of airway occlusions (“plugs”) close to their point of rupture is examined using a free-surface model (volume of fluid), such that factors influencing the formation of droplets during their reopening can be identified. The propagation of airway occlusions is highly influenced by recirculating flow at the edge of the front interface, where significant fluctuations of wall shear stresses occur. The resulting drag force causes the rear interface to advance at a greater rate, destabilizing the plug. As the plug thickness decreases, a thin film with uniform thickness forms, resulting in a disk-like structure around the centerline. Rupture occurs around the disk formation largely due to surface tension instability. At lower pressures, smaller disks form causing the rupture to occur through a puncture point (forming no droplets); at higher pressures, a larger disk forms, with rupture occurring along the disk edge and at the center (forming multiple droplets). Upon reopening, a jet of air is produced, causing a temporary increase in shear stress along the wall. However, the magnitude and duration of this increase do not scale directly to the applied pressure, as the formation of droplets and irregularities in airway lining were found to disrupt the flow field and the shear stresses at the wall.

Due to the recent outbreak of SARS-CoV-2, the mechanisms of transport of viruses and diseases through air have gained considerable attention. Substantial effort has been put toward the minimization of droplet transport and disease transmission through social distancing and mask wearing. Various flow visualization and numerical studies have been conducted to determine the transfer of droplets during coughing and sneezing.1,2 Such processes typically involve highly turbulent flows, with the formation of a large distribution of droplet sizes.3–5 Bioaerosols, however, can be formed even in the absence of turbulence, e.g., during a normal breathing process.6 This is caused by the instabilities of the respiratory tract lining fluid (RTLF), which close the small airways (e.g., bronchioles) during contraction, forming mucus/liquid plugs. Upon re-expansion, the closure/plug decreases in thickness and bursts open, forming aerosol droplets that are then released alongside the exhaled breath.7–9 This phenomenon has been termed as bronchiole fluid film burst (BFFB), where an increase in the concentration of aerosol following a rapid inhalation and deep exhalation cycle has been observed.10 In contrast, no significant increase in aerosol production was observed with an increase in exhalation rate, which indicates an independence of the phenomenon from the level of turbulence within the airways. Furthermore, various studies have been conducted on the effects of various physical activities on aerosol concentration of exhaled breath,11–13 finding increased aerosol production rates with increased exercise intensity.

The study of the formation of aerosols and droplets due to the closure and reopening of small airways is interesting due to several factors. First, due to the small sizes of the droplets, they typically stay for extended periods of time in the atmosphere, increasing the transmission probability of viruses and diseases.14,15 Second, there has also been considerable interest in the sampling and measurement of exhaled breath condensates, as it can potentially provide a convenient and noninvasive approach for sampling RTLF, e.g., for the detection and measurement of inflammatory indicators in asthma and lung diseases. Unfortunately, while this is a promising method, various problems still exist, one of which is the uncertainties associated with the source and mechanisms of the droplet formation.16 A greater understanding of factors affecting the formation of droplets upon the reopening of closed airways can, therefore, contribute toward the development of this method.

The propagation of a liquid plug in a flexible lined tube in the Stokes flow regime has been studied analytically by Howell et al.17 It was found that the liquid plug behavior can be generally described by the Capillary number, C a = 3 μ U / σ. Above a certain threshold (CaC), the thickness of the plug decreases and vice versa. CaC was found to scale with the thickness of the mucus film ahead of the plug. The elasticity of the wall was found to be a contributing factor, wherein a more compliant tube was found to be harder to re-open. Fujioka and Grotberg18 analyzed the effect of inertia ( R e = ρ U R / μ 40 80) on the dynamics of the mucus plug in a two-dimensional geometry using a single-fluid Lagrangian approach. Recirculation flow was found to occur within the plug, with its center being skewed further toward the rear meniscus, resulting in a fore-aft asymmetry in the plug (particularly at high Re values). Fujioka et al.19 found a dependency of the resistance experienced by the plug on the thickness of the mucus lining, whereby with decreasing lining thickness, an increase in resistance/drag is obtained at the interface front. This resulted in a condition where the plug decreases in thickness and ruptures. Magniez et al.20 proposed an analytical model for the critical pressure required for the destabilization of a mucus plug at low values of Reynolds number (Re < 10). Various studies have also been conducted on other factors affecting the propagation, growth, and thinning of the plug, e.g., surfactants and surface tension gradients,21,22 non-Newtonian mucus behavior,23,24 non-constant pressure-driving force,25 and airway geometry and bifurcations.26 

The rupture process of a mucus plug involves significant topological changes that can be difficult to capture accurately in a numerical and computational fluid dynamics (CFD) framework. The works of Fujioka et al.18,19 (and other associated works) were mostly done using a Lagrangian approach, whereby only the dynamics of the liquid (mucus) phase is considered, and the numerical domain and grid are moved dynamically according to the fluid flow (as described by the Navier–Stokes equations). While this approach is able to capture the dynamics of the plug during propagation, it cannot capture its rupture dynamics. Hassan et al.27 bridged this gap using an adaptive Lagrangian–Eulerian approach and found that the front meniscus carries the maximum shear and pressure stresses at the wall and that peak values are obtained during the rupture, possibly causing an injury on the airway wall during the reopening process. This is confirmed by the study of Muradoglu et al.,22 using the Front Tracking method to capture the interface during rupture. Malashenko et al.28 studied the rupture of the mucus plug using an axisymmetric volume-of-fluid (VOF) model and observed three possible modes of rupture: (i) a simple collapse yielding no droplets; (ii) the formation of a large droplet containing most of the mass of the disintegrated plug; and (iii) the formation of several small droplets. He et al.29 employed a two-dimensional Lattice Boltzmann method and observed the first and second modes of rupture. It was observed that the size of the aerosol formed upon the rupture of the mucus plug depends greatly on the thickness of the liquid film; increased thickness in mucus lining results in a decrease in aerosol size (and vice versa).

In the current study, a 3D analysis of the propagation of a mucus plug in a pulmonary airway is presented. This is in contrast to most numerical studies that have been done to date on the dynamics of airway closure and reopening, which generally involve 2D geometries with inherent symmetrical/axi-symmetrical assumptions (as discussed above). The effects of fluid properties (mucus viscosity and surface tension) and flow parameters on the stability of the mucus plugs will be assessed, so as to obtain a greater understanding of determinant parameters for the reopening of blocked airways. Furthermore, in the case where the mucus plug undergoes a reduction in thickness, its behavior upon rupture will be analyzed to gain a better understanding of the generation of droplets and aerosols upon the reopening of blocked airways.

We consider a simple airway geometry, represented by a rigid cylinder/tube lined with a thin layer of liquid, which is a mixture of mucus and serous layers (see Fig. 1). Mucus is known to have a highly variable viscosity ( 0.01 100 Pa s), varying greatly from one individual to the next, as well as between different physiological functions within the body. Mucus viscosity also depends on the composition of mucins and their glycosylation, thus varying with age, diet, the presence of specific antigens, commensals, and pathogens.30 It is also known to be highly non-Newtonian, having shear-rate dependent viscosity as well as some elastic responses in deformation. The serous layer, on the other hand, has a viscosity that is similar to water (∼0.001 Pa s). In the current study, the contributions of mucus and serous layer toward the dynamics of the mucus plug and lining layer are simplified, such that the liquid phase is assumed to have a constant viscosity (Newtonian). A base case scenario is selected, such that the viscosity of the liquid phase is μ L = 0.01 Pa s and the surface tension between the two phases is σ = 0.167 N/m. Furthermore, the viscosity and density of the gas phase are assumed constant ( ρ G = 1.3 kg / m 3 and μ G = 4 × 10 6 Pa s). Based on these parameters, the characteristic Laplace number of the problem (assuming constant liquid density of ρ L = 1000 kg / m 3) is L a = ρ l σ R / μ L 2 = 1000.

FIG. 1.

Schematic (cross-sectional) representation of the occluded airway.

FIG. 1.

Schematic (cross-sectional) representation of the occluded airway.

Close modal

1. Initial and boundary conditions

The wall of the airway/tube is modeled with a no-slip boundary condition that is fully wetted with the liquid phase. A typical adult lung airway at tenth generation is 0.6 mm in radius; this figure is, therefore, selected as the radius for the tube, a. The length of the tube, L, is taken to be 24a to allow for the development of flow field around the highly idealized (and initially stationary) mucus plug and for it to increase/reduce in thickness. The annular space inside the tube is occupied by a layer of liquid with a thickness of h = 0.03 0.06 mm. The airway occlusion is modeled as a liquid plug that completely blocks airflow within the tube; its initial thickness is equivalent to the radius of the tube, i.e., Lb = a, and the curvature of the front and back interface of the plug is a. The liquid plug travels through the tube under a pressure-driving force. Considering that there is minimal variation in the pressure field within the gas phase, a differential pressure can then be assigned at the back and the front of the plug ( Δ P = P 1 P 2). Zero-gradient boundary conditions are assigned for the phase volumetric fraction (c) at each end of the tube.

The flow of the liquid and gas phases within the airway is modeled through the incompressible Navier–Stokes and continuity equations with variable density and surface tension force
ρ ( u t + u · u ) = P + μ 2 u + σ κ δ s n ,
(1)
ρ t + · ( ρ u ) = 0 ,
(2)
· u = 0.
(3)
Equation (1) describes the conservation of momentum, (2) shows the advection of the two phases (as described by their densities), and (3) completes the mass conservation condition of the system and applies the incompressibility assumption for the two phases. u = ( u , v , w ) is the velocity vector, P is the pressure, and the last term on the RHS of Eq. (1) represents the surface tension force between the liquid and gas phases. σ is the surface tension coefficient between the two phases, κ is the interface curvature, and n is the normal vector to the interface. δs is the Dirac function, here implemented to confine the surface tension force within regions surrounding the interface. The physical properties used in Eqs. (1) and (2), i.e., the density and viscosity (ρ and μ, respectively), are calculated based on a volume fraction scalar field (c), such that
ρ = c ρ L + ( 1 c ) ρ G ,
(4)
μ = c μ L + ( 1 c ) μ G .
(5)
ρL, ρG, muL, and muG represent the density and viscosity of the liquid and gas phases, respectively. c represents the fraction of volume occupied by the liquid phase in each computational grid/mesh, such that a value of 0 represents a computational grid fully occupied by the gas phase (and vice versa for a cell c value of 1). The advection of the fluids [Eq. (2)] is evaluated based c, such that
c t + · ( c u ) = 0.
(6)

The Navier–Stokes and continuity equations [Eqs. (1)–(3)] are non-dimensionalized against the characteristic velocity of the flow due to the capillary force between the two phases, i.e., u ref = σ / ( ρ L a ). Furthermore, considering a as the characteristic length of the problem, the temporal scaling then becomes t ref = ( a 1.5 ρ L 0.5 / σ 0.5 ). The reference pressure due to capillary force is P ref = σ / a. Furthermore, the density and viscosity of the two phases are normalized against the respective values of the liquid phase, i.e., ρ * = ρ / ρ L and μ * = μ / μ L. This results in the scaling of the viscous term of Eq. (1) with respect to the Reynolds number, R e = ρ L U ref a / μ L. The Reynolds number of the base case scenario as detailed in Sec. II A is 31.6.

In the following, dimensionless terms within the manuscript are indicated by a superscript *, i.e., P * = P / P ref, etc. Further details of the non-dimensional form of the equations have been provided in  Appendix A.

The set of governing equations is solved through the VOF implementation currently available in the open-source BASILISK flow solver (http://basilisk.fr/), which features an efficient adaptive mesh capability, allowing for a reasonable level of grid resolution to be implemented on thin films and interfaces within the three-dimensional airway geometry. The solver is based on the implementation of the projection method of Chorin31 on a quad/octree computational grid/mesh,32 along with a multilevel Poisson solver to allow for pressure and continuity correction to be applied across a dynamically adapted mesh.32 Furthermore, the geometric VOF method (piecewise-linear) is implemented, with the interface segments being reconstructed using the mixed-Young's-centered (MYC) approximation,33 allowing for an explicit tracking of the interface. The surface tension between the two phases is modeled as a bulk volumetric force based on the continuum surface force (CSF) approach.34 To minimize the well-known issue of parasitic currents associated with the VOF approach, a balanced-force description of the surface tension and pressure gradient is employed. This is aided with the application of a height function scheme based on the VOF field to obtain a second-order representation of the interface curvature. As a result of the adaptive mesh refinement, the BASILISK solver has been used to simulate problems with large distribution of length scales and density ratios, e.g., the generation of up to 1000 droplets (100–400 μm) in a cough/sneeze.35 Gumulya et al.36 have demonstrated the validity of the code for the case of a bubble rising in a quiescent liquid against experimental data, whereas Wu et al.37 demonstrated a good correspondence on the droplet impact and splashing dynamics onto a deep pool to experimental data. Furthermore, Romano et al.38 employed the BASILISK solver to present an axisymmetric analysis of the closure of an airway and demonstrated a good correspondence with experimental observations.

To initialize the simulation, a cubic domain of size ( L 0 * = 24 3) is constructed, with an initial mesh of 1283 cells. The cylindrical geometry of the airway is then applied using embedded boundaries, with cells intersecting the tube being refined to a size of L 0 * / 2 10. These boundaries are defined through a separate volume-fraction field, denoting the ratio of volume occupied by the fluid relative to the total volume of the cell.39 Following the initialization of the pressure, velocity, and phase fraction fields, the computational mesh is adapted along the liquid/gas interface, such that cells along the interface are refined to a size of L 0 * / 2 10, and coarser mesh can be obtained in the bulk of the phases. The grid adaption step is performed at each time step, based on estimated errors associated with the volume fraction and velocity fields, with maximum tolerated errors of 0.001 and 0.005, respectively. The number of computational cells resolving the computational domain, therefore, varies between 3–6 × 106 cells, depending on the shape of the interface and flow field. To assure numerical stability, the simulation is progressed with (non-dimensional) time step of 10 4, and the tolerance on the solver is set to 10 6 in order to minimize mass conservation errors throughout the simulation. As a result of the grid refinement, the lining layer inside the airway is resolved by 2.1 cells for the thinnest lining case considered in this study ( h * = 0.05). In cases where the plug ruptures, the size and number of droplets generated at the reopening are determined by tagging connected regions in the domain, identified through neighboring cells with c value > 10 3. Further discussion on the effect of grid resolution and convergence of this framework is presented in  Appendix B.

To check the validity of the current numerical setup, the propagation behaviors of the mucus plugs are compared against the published studies of Fujioka et al.19, Figure 2 presents the evolution of the thickness of the mucus plugs at the centerline ( L b *) and its propagation velocity ( U *) under different values of pressure-driving force and lining thicknesses. At h * = 0.1 and Δ P * = 1, the mucus plug is stable and grows in thickness over time. The propagation rate of the plug increases as it accelerates from stationary. A maximum value of propagation velocity is obtained at t * > 5.7 ( U max * = 0.468), after which the mucus plug gradually decelerates (by t * = 11, the propagation velocity has decreased by 3 %). The growth in the thickness of the plug and its propagation velocity at Δ P * = 1 are found to be in reasonable agreement to the axisymmetric results of Fujioka et al.,19 supporting the validity of the current numerical method. With higher values of pressure differential, a departure from stability is obtained. At Δ P * = 4 (and h * = 0.1), the mucus plug decreases in thickness over time although rupture is not obtained over the length of the tube. As can be expected, the propagation rate of the plug is higher than in the case of Δ P * = 1, with U max * = 1.83. After reaching Umax, the mucus plug progresses to decelerate rapidly, indicating an increased drag force as the plug thickness decreases.

FIG. 2.

The dynamics of the mucus plug at various values of Δ P *: (a) the thickness of the plug at the centerline ( h * = 0.1); (b) propagation velocity of the plug ( h * = 0.1); (c) the thickness of the plug at the centerline ( h * = 0.05); (d) propagation velocity of the plug ( h * = 0.05).

FIG. 2.

The dynamics of the mucus plug at various values of Δ P *: (a) the thickness of the plug at the centerline ( h * = 0.1); (b) propagation velocity of the plug ( h * = 0.1); (c) the thickness of the plug at the centerline ( h * = 0.05); (d) propagation velocity of the plug ( h * = 0.05).

Close modal

At h * = 0.05, the mucus plug is unstable even at Δ P * = 1. This is in agreement with the results of Fujioka et al.19 although the current data indicates a lower rate of thickness reduction and acceleration in plug propagation rate than in the case of the axisymmetric study. In comparison to the case of h * = 0.1, higher rates of plug thickness reduction are obtained, indicating that a lower value of pressure-driving force is required to destabilize mucus plugs in airways with thinner lining layers. Furthermore, rupture is more likely to be obtained over a lower propagation distance. This is in agreement with the steady-state analytical study of Howell et al.,17 who postulated a critical parameter of differential pressure, above which the mucus plug would eventually rupture (and below which the plug is stable). It was determined that the critical value of pressure differential is largely dependent on the film thickness ahead of the mucus plug. Furthermore, the plug propagation rate in the case of h * = 0.05 appears to follow a different pattern than in the case of h * = 0.1—instead of reaching a maximum value and decelerating, the mucus plugs accelerate continuously toward the point of rupture. Both of these propagation behaviors have previously been observed19,27 and are in relatively good agreement with their results, supporting the validity of the current numerical model.

The times and locations (travel distance) of the mucus plugs at rupture have been presented in Table I. It is evident that mucus plugs with thicker lining fluids require greater pressure drops for their destabilization. Furthermore, the destabilization and thickness reduction process take place over longer distances, and the plug has a greater chance of traveling to smaller airways before rupturing.

TABLE I.

Rupture times ( t R *), rupture locations ( x R *), and droplet formation of mucus plugs with lining thicknesses of h * = 0.05 0.1 under different pressure-driving forces.

h * Δ P * t R * (ms) x R * (mm) No. of droplets D max * (μm) D mean * (μm)
0.1  4.90 [5.58]  20.6 [12.4]  0.231 [138.8]  0.231 [138.8] 
0.1  2.40 [2.73]  15.4 [9.2]  0.261 [156.8]  0.212 [127.0] 
0.05  14.7 [16.7]  17.5 [10.5]  ⋯  ⋯  ⋯ 
0.05  4.40 [5.01]  9.7 [5.8]  ⋯  ⋯  ⋯ 
0.05  3.41 [3.88]  9.3 [5.6]  0.223 [134.1]  0.223 [134.1] 
0.05  2.50 [2.85]  8.8 [5.3]  0.227 [136.1]  0.160 [96.3] 
0.05  1.81 [2.06]  8.5 [5.1]  0.281 [168.8]  0.165 [99.2] 
h * Δ P * t R * (ms) x R * (mm) No. of droplets D max * (μm) D mean * (μm)
0.1  4.90 [5.58]  20.6 [12.4]  0.231 [138.8]  0.231 [138.8] 
0.1  2.40 [2.73]  15.4 [9.2]  0.261 [156.8]  0.212 [127.0] 
0.05  14.7 [16.7]  17.5 [10.5]  ⋯  ⋯  ⋯ 
0.05  4.40 [5.01]  9.7 [5.8]  ⋯  ⋯  ⋯ 
0.05  3.41 [3.88]  9.3 [5.6]  0.223 [134.1]  0.223 [134.1] 
0.05  2.50 [2.85]  8.8 [5.3]  0.227 [136.1]  0.160 [96.3] 
0.05  1.81 [2.06]  8.5 [5.1]  0.281 [168.8]  0.165 [99.2] 

Figure 3 shows the shape of the mucus plug and lining layer during propagation. At Δ P * = 1, it can be seen that the mucus plug generally retains strong fore-aft symmetry, particularly in the case where the mucus plug is stable at h * = 0.1. With higher pressure-driving force, the front interface adopts a flatter curvature, and the back of the plug features a thicker trailing liquid layer (thus a higher curvature at the back interface). The increased thickness of the lining layer at the back of an unstable mucus plug has previously been observed in various studies,19,27 attributed toward the mass balance of the problem—an unsteady plug that is decreasing in size would “leave behind” a trail of liquid during propagation, resulting in a thicker lining layer at the back of the plug. The rapid decrease in plug thickness in cases of high pressure drops, therefore, results in the formation of thicker lining liquids behind the mucus plugs. At even higher pressure drops, the curvature of the front interface is seen to gradually reverse in direction, toward the direction of the propagation [see Figs. 3(d) and 3(g)]. Furthermore, it can also be seen that different outcomes are obtained upon the rupture of the mucus plug. In the cases of h * = 0.05 and Δ P * = 1, no droplet formation occurs, whereas at Δ P * 3, droplets are formed upon rupture. Further details have been presented in Table I, where the number of droplets formed, along with their mean diameters, are shown. With increasing pressure-driving force, smaller (and more numerous) droplets are generally formed upon plug rupture. Furthermore, cases with thicker lining layers ( h * = 0.1) generally produce less number of droplets than in the case of thinner lining layers ( h * = 0.05). In Sec. III B, the flow field within the propagating plug is analyzed. The rupture behavior of the plugs is then discussed in Sec. III C, along with factors influencing the formation of droplets upon the reopening of the airway.

FIG. 3.

Cross-sectional view of the volume fraction field: (a) h * = 0.05 , Δ P * = 1; (b) h * = 0.05 , Δ P * = 3; (c) h * = 0.05 , Δ P * = 4; (d) h * = 0.05 , Δ P * = 5; (e) h * = 0.1 , Δ P * = 1; (f) h * = 0.1 , Δ P * = 4; (g) h * = 0.1 , Δ P * = 6.

FIG. 3.

Cross-sectional view of the volume fraction field: (a) h * = 0.05 , Δ P * = 1; (b) h * = 0.05 , Δ P * = 3; (c) h * = 0.05 , Δ P * = 4; (d) h * = 0.05 , Δ P * = 5; (e) h * = 0.1 , Δ P * = 1; (f) h * = 0.1 , Δ P * = 4; (g) h * = 0.1 , Δ P * = 6.

Close modal

The propagation and stability of a mucus plug in a lined airway tube has previously been found to be largely influenced by the presence of a recirculation region that occurs ahead of the front interface of the plug, around the region of the lining layer.19,27 This can be seen in Fig. 4(a), where the distribution of axial velocity at both interfaces has been presented along the radial direction of the tube for the case of h * = 0.05 and Δ P * = 1, with the recirculation at the front edge of the mucus plug causing the axial velocity of the front interface close to the tube wall to be negative. On the other hand, the axial velocity along the back interface at | y * | 0.695 is greater than that along the front interface (resulting in the reduction of the plug thickness). Furthermore, the maximum axial velocity at the back interface occurs at | y * | 0.56, with a local minimum obtained at the centerline of the tube. A similar velocity distribution is observed with the front interface, although the maxima in u x * in this case occur further away from the centerline (at | y * | 0.71), due to the thinner lining layer at the front of the mucus plug. The distribution of interfacial axial velocity for the case of h * = 0.1 and Δ P * = 1 [see Fig. 4(b)] can be seen to follow a similar pattern although, in this case, the axial velocity of the back interface around the centerline of the tube | y * | 0.56 is generally lower than that at the front interface (thus causing the growth of the plug).

FIG. 4.

Flow field behavior for the case of Δ P * = 1 at t * = 7: (a) axial velocity at the interface ( h / a = 0.05); (b) axial velocity at the interface ( h / a = 0.1); (c) pressure contour (top half), axial velocity contour, and flow streamlines (bottom half) ( h / a = 0.05); (d) distribution of pressure along the wall; (e) distribution of shear stress along the wall.

FIG. 4.

Flow field behavior for the case of Δ P * = 1 at t * = 7: (a) axial velocity at the interface ( h / a = 0.05); (b) axial velocity at the interface ( h / a = 0.1); (c) pressure contour (top half), axial velocity contour, and flow streamlines (bottom half) ( h / a = 0.05); (d) distribution of pressure along the wall; (e) distribution of shear stress along the wall.

Close modal

The local minima in axial velocity seen at the center of the tube in the case of Figs. 4(a) and 4(b) can be inspected further in Fig. 4(c) where the streamlines, pressure, and velocity contours in and around the plug at h * = 0.05 and Δ P * = 1 have been presented. The pressure field at either side of the plug can be seen to be relatively uniform, in reflection of the adopted assumption that the pressure drop in the gaseous phase is minimal relative to the pressure variation within the liquid phase. Behind the plug, a region of local minimum in axial velocity occurs (still in the positive direction) as the flow diverts away from the centerline due to the presence and the curved shape of the interface. A similar observation is made with the flow field ahead of the front interface, as the flow diverts toward the centerline after following the curve of the interface and the recirculation region. The flow within the plug itself largely follows a linear pathline, apart from the front edge of the mucus plug, where a recirculation region occurs. As a result, a region of velocity minimum occurs, which coincides with a sharp fluctuation of the pressure field at the wall.

Figures 4(d) and 4(e) present the distribution of pressure and shear stress [ τ x y = μ L ( u x / y )] fields along the wall for the cases of h * = 0.05 and 0.1 (at Δ P * = 1). The x * values indicated in these graphs represent the distance from the center of mass of the mucus plug along the length of the airway tube. Figure 4(d) shows that along the mucus plug, regions of low pressure are obtained, whereas ahead of the plug, a region of local pressure maximum occurs due to the recirculation of the flow field in the area. This is accompanied by a region of high shear stress values along the wall, in reflection of the viscous forces affected upon the propagation of the plug, which is followed by a sharp reversal in the direction of the shear stress. These distributions of pressure and shear stress along the wall are in good agreement with the findings of Fujioka et al.19 In comparing the two cases presented here, it is seen that the fluctuations in both the pressure and shear stress fields in the case of h * = 0.1 are not as significant to those at h * = 0.05. This is caused by the position of the center of flow recirculation, which occurs closer to the tube wall (at r * = 0.85 and 0.83 for h * = 0.05 and 0.1, respectively). The smaller recirculation region in the case of h * = 0.05 intensifies the flow reversal at the edge of the plug, thus causing the sharp fluctuation in pressure and shear stress in this region.

1. Base case

The dynamics of the mucus plug ( h * = 0.05 and Δ P * = 1) around the point of rupture is examined in Fig. 5. Figure 5(a) shows the flow streamlines in and around the mucus plug (along with three-dimensional representation of the interface), whereas the distribution of u x * along the radial direction at the front interface can be seen in Fig. 5(b). It is seen that the distribution of interfacial velocity as shown earlier in Fig. 4(a) is maintained, even as the central thickness of the plug reaches L b * = 0.07 (at t * = 14). At t * = 15, a region of pressure minimum develops at the center of the mucus plug, translating to a sharp increase in axial velocity at the front interface. This is followed by the development of a local pressure maximum at the same location at t * = 15.1. This fluctuation in pressure results in a break from axi-symmetry in the flow field at the centerline of the tube, and a dimple is seen to develop in the interface at t * = 15.2. The dimple formation at the front interface causes the development of vortices ahead of the mucus plug [see the flow field in Fig. 5(a) at t * = 15.3], in turn, causing the formation of a disk-like region around the tube centerline with a (dimensionless) radius of r d * 0.2 with uneven surface, but with a relatively uniform film thickness of 0.0513 at t * = 15.4. It is noted that within the framework of the current numerical simulation, this thin film is resolved with approximately two numerical cells. Future work will seek to further establish the effects of this numerical limitation on the dynamics of the rupture behavior. It is also noted that the rupture behavior of the mucus plug, in this case, is quite disimilar to the case of film rupture due to bubble film drainage, in which the thinnest part of the film is at the axis of symmetry. In the current case, rupture tends to occur at the edge of the disk-like thin-film structure due to pressure and surface tension instabilities in the thin film. For the case of h * = 0.05 and Δ P * = 1, rupture occurs following a puncture in the thin film at r * = 0.22 away from the centerline.

FIG. 5.

Rupture behavior for the case of h * = 0.05 and Δ P * = 1: (a) pressure contour and flow streamlines, with 3D representation of the interface around the time of rupture; (b) axial velocity at the front interface; (c) distribution of shear stress along the wall before and after rupture.

FIG. 5.

Rupture behavior for the case of h * = 0.05 and Δ P * = 1: (a) pressure contour and flow streamlines, with 3D representation of the interface around the time of rupture; (b) axial velocity at the front interface; (c) distribution of shear stress along the wall before and after rupture.

Close modal

The distribution of shear stress along the bottom wall can be seen in Fig. 5(c). The x * values indicated in these graphs represent the distance from the center of mass of the mucus plug at t * = 15; subsequent shear stress curves were adjusted around this reference position for clarity of their presentation. Prior to the rupture of the mucus plug, it is seen that the shear stress along the wall maintains a similar distribution trend to earlier timelines [c.f. Fig. 4(e) and the shear stress distribution at t * 15.6], indicating minimal disruptions in the viscous flow field along the wall despite the three-dimensional flow structure at the centerline. Upon rupture, a jet of air is formed at the rupture opening, resulting in a sharp increase in the maximum shear stress along the wall (at t * = 15.7). Over time, the intensity of the jet progresses to dissipate, and the flow field returns toward regular streamline flow, with more uniform distribution of lining thickness within the tube.

2. Effects of pressure

With increasing pressure-driving forces, an increase in interfacial velocity is observed, particularly at the centerline of the tube. Figures 6 and 7 show the flow fields in the liquid and gas phases around the point of rupture at Δ P * = 3 and 5, respectively, with h * = 0.05. Figure 6(b) shows that a local maximum is developed at the centerline of the back interface at Δ P * = 3. This continues to develop with increasing pressure-driving force, and at Δ P * = 5 [see Fig. 7(b)], global maximum axial velocity is obtained at the centerline. This is due to the increased thickness of the film at the back of the mucus plug with increasing pressure drop, which allows for the flow field in the gas phase behind the back interface to be more streamlined, requiring less diversion away from the centerline to follow the curvature of the interface [Fig. 7(a)]. Meanwhile, the velocity field at the front interface can be seen to attain a greater degree of uniformity with increasing pressure-driving force. The local velocity minimum earlier seen at the centerline of the front interface in the case of Δ P * = 1 becomes less apparent in the case of 3 Δ P * 4. Furthermore, the locations of the maxima in axial velocity at the front interface shift further out in the radial direction (at | y * | 0.735 and 0.75 for cases of Δ P * = 3 and 4, respectively), as a result of the decreased lining thickness ahead of the plug with increasing pressure-driving forces. The velocity maximum observed at the centerline of the back interface in the case of Δ P * = 5 is also seen to occur at the front interface. This is caused by the flat curvature of the front interface at the center of the tube and the sharp curvature of the interface close to the tube wall, which pushes the flow field toward the centerline [see Fig. 7(a)].

FIG. 6.

Rupture behavior at h * = 0.05 and Δ P * = 3: (a) pressure contour and flow streamlines, with 3D representation of the interface around the time of rupture; (b) axial velocity of the front and back interfaces at t * = 3.2; (c) distribution of shear stress along the wall before and after rupture.

FIG. 6.

Rupture behavior at h * = 0.05 and Δ P * = 3: (a) pressure contour and flow streamlines, with 3D representation of the interface around the time of rupture; (b) axial velocity of the front and back interfaces at t * = 3.2; (c) distribution of shear stress along the wall before and after rupture.

Close modal
FIG. 7.

Rupture behavior at h * = 0.05 and Δ P * = 5: (a) pressure contour and flow streamlines, with 3D representation of the interface around the time of rupture; (b) axial velocity of the front and back interfaces at t * = 1.7; (c) distribution of shear stress along the wall before and after rupture.

FIG. 7.

Rupture behavior at h * = 0.05 and Δ P * = 5: (a) pressure contour and flow streamlines, with 3D representation of the interface around the time of rupture; (b) axial velocity of the front and back interfaces at t * = 1.7; (c) distribution of shear stress along the wall before and after rupture.

Close modal

The differences in the axial velocity of the interfaces cause different rupture behaviors within the mucus plugs. The disk-like formation of thin film around the centerline seen earlier in the case of Δ P * = 1 (and h * = 0.05) is again observed in Fig. 6(a) (at t * = 3.42–3.44), with a slightly larger radius of r d * 0.25. At t * = 3.46, rupture occurs along the edges of the disk formation, and a section of thin film is dislodged from the bulk of the mucus plug. This results in the formation of a jet of air around the disloged segment and a large droplet within the gap.

With increasing pressure drive, the development of thin-film sections around the axis of symmetry is observed to occur across larger radial distances ( r d * 0.274 and 0.358 for the cases of Δ P * = 4 and 5, respectively). Furthermore, due to the increased axial velocity at the center of the back interface, it is seen that the thickness of the film at the centerline of the tube decreases with increasing pressure drive, resulting in the fragmentation of the thin film at multiple locations, i.e., at the edges of the disk formation as well as at the central axis of the tube [see Fig. 7(a)]. This results in the formation of multiple droplets at higher values of pressure drops.

In comparing the distribution of shear stress along the wall [see Figs. 5–7(c)], it can be seen that a considerable increase in shear stress is obtained in the lining layer behind the mucus plug with increasing pressure-driving forces. This is as expected due to the increased velocity gradients presented in the lining layer in cases of high pressure-driving forces and the increased thickness of lining layer behind the plugs. Ahead of the plug, a minimum in shear stress is observed, as a result of the flow recirculation at the edge of the lining layer. The effect of the decreased lining thickness ahead of the plug is evident here, as the area affected by the shear stress minimum (where τ w 0) can be seen to decrease with increasing pressure drops.

For all the cases considered in the current study, it is seen that the distribution of shear stress along the wall remains comparatively steady up to the point of rupture. Upon rupture, the formation of a jet of air in the film opening results in a temporary increase in the maximum shear stress at the wall, caused by instabilities of flow fields. Over time, local regions of high velocities dissipate and the flow returns to regular flow, here indicated by subsequent decrease in wall shear stress. With increasing pressure-driving force, it is observed that the decrease in shear stress after the point of rupture occurs over shorter periods of time ( Δ t *  0.54, 0.38, and 0.14 for cases of Δ P * = 3, 4, and 5, respectively, with Δ t * in this case being the time lapse since the point of rupture to the instance where the maximum shear stress along the wall is less or equal to the maximum shear stress prior to rupture). This could be caused by the formation of droplets in cases of higher pressure drops, which carry some of the momentum borne by the mucus plug prior to the point of rupture and present disruptions in the air flow at the rupture opening.

3. Effects of lining layer thickness

In comparing the rupture behaviors of plugs with thinner/thicker lining layers, it is seen that the formation of disk-like structure of thin liquid at the center of the plug is generally smaller in cases of h * = 0.1 (c.f. h * = 0.05). For the case of Δ P * = 4, the radius of the disk formation is found to be r d *  0.265, whereas for Δ P * = 6 , r d *  0.3. For the case of Δ P * = 4, rupture occurs along the edge of the disk formation, thus producing a single droplet in the process. For the case of Δ P * = 6, on the other hand, rupture occurs along the edge of the disk formation, as well as the centerline of the tube, thus creating multiple droplets in the process.

The smaller thin-film disk formation in cases of h * = 0.1 can be attributed toward the decrease in plug propagation rate over time [see Fig. 2(b)]. This is caused by the formation of multiple recirculating regions in the lining layer ahead of the plug [see Fig. 8(a)], which enlarges the region of pressure minimum at the edge of the front interface and intensifies the shear stress minimum and flow reversal in this region [see Fig. 8(c)]. The resulting increase in drag force at the edge of the front interface is observed in the curvature of the interface, which is generally found to be flatter in the case of h * = 0.1 in comparison to h * = 0.05 [see Figs. 3(f) and 3(g), with the curvature of the front interface at Δ P * = 6 being reversed along the direction of the plug propagation]. Furthermore, it is seen that this increase in drag force affects the propagation rate of the interface. For the case of Δ P * = 4, the values of axial velocity at the front and back interfaces noticeably decrease between t * = 2.3 and 4.5 (closer to the point of rupture) [see Fig. 8(b)], with the values of u x * at the back interface at t * = 4.5 closely matching those at the front interface, particularly at | y * | 0.7. The formation of the thin-film disk section that occurs prior to the rupture of the mucus plug for the case of h * = 0.1, therefore, appears to be highly affected by the increased drag force due to the formation of multiple circulation regions in the lining layer ahead of the plug.

FIG. 8.

Rupture behavior at h * = 0.1 and Δ P * = 4: (a) pressure contour and flow streamlines, with 3D representation of the interface around the time of rupture; (b) axial velocity of the front and back interfaces at t * = 2.2 and 4.5, compared against the equivalent parameters at h * = 0.05 , Δ P * = 4, and t * = 2.3 (shown as smooth and dashed lines); (c) distribution of shear stress along the wall before and after rupture.

FIG. 8.

Rupture behavior at h * = 0.1 and Δ P * = 4: (a) pressure contour and flow streamlines, with 3D representation of the interface around the time of rupture; (b) axial velocity of the front and back interfaces at t * = 2.2 and 4.5, compared against the equivalent parameters at h * = 0.05 , Δ P * = 4, and t * = 2.3 (shown as smooth and dashed lines); (c) distribution of shear stress along the wall before and after rupture.

Close modal

Figure 8(c) also shows that upon the rupture of the mucus plug, minimal increase in the maximum value of shear stress is obtained in the current case (unlike in cases of h * = 0.05). This can also be attributed toward the formation of multiple recirculation regions along the lining layer ahead of the point of rupture, which confine the direction of the jet flow toward the centerline of the tube, thus minimizing the effect of the temporal variations of velocity gradients at the tube wall.

The results of the current study show that the rupture of a mucus plug, particularly at higher pressure drops, is a three-dimensional phenomenon, possibly resulting in the formation of multiple droplets (cf. the two-dimensional analysis of He et al.)29 The dependency of the rupture behavior on the level of pressure drop indicates the effects of flow rate during inhalation upon droplet formation, consistent with the observation of Johnson and Morawska.10 On the other hand, the range of average and maximum droplet sizes obtained in the current study is large in comparison to those reported in experimental studies of exhaled breath (Morawska et al.6 indicated bimodal droplet size distribution of 0.8 and 1.8 μm in a variety of breathing and speech patterns). Larger droplets of such size generated in lower airways are more likely to settle within the airway and, therefore, do not get released during exhalation. However, the results do indicate that the formation of multiple (smaller) droplets is more likely in cases of high pressure drops, e.g., in cases where a person is undertaking rigorous physical activities with higher levels of tidal volume. This is consistent with observations of various experimental studies considering the effects of physical activity on exhaled droplet concentration,11–13 whereby increased concentration of aerosol generation rate has been observed with increased exercise intensity.

The reopening of obstructed airways has also previously been associated with injuries in the epitheleal cells18,19,25,29,40,41 due to the significant variations of shear stress and pressure along the tube wall during this process. The results of the current study demonstrate that heightened peaks in shear stress do occur during the propagation of a mucus plug, particularly due to the presence of flow recirculation regions at the edge of the front interface of the occlusion. In the rupture process, the peaks of the shear stress are increased due to the formation of jet flow that disrupts the flow field at the lining layer shortly after rupture. Unsurprisingly, the peaks in shear stress along the wall appear to scale with the applied pressure drop. However, the current findings suggest that the level and probability of injury in the epitheleal cells cannot be scaled directly against the applied pressure drop, as the formation of droplets and irregularities in the lining layer ahead of the rupture may present disruptions in the formation of the jet. These disruptions have been found to decrease the duration of the peaks in wall shear stress, therefore decreasing the time of exposure. This has previously been observed experimentally in the study of Bilek et al.,40 who found that epitheleal cell damage increased with decreasing reopening velocity, and may provide some rationalization toward this observation.

The current study examines the rupture and reopening dynamics of airway occlusions propagated under a pressure-driving force (such as in the inhalation phase of a breathing cycle) through numerical analysis. It is shown that the destabilization of the occlusion (decrease in thickness) is highly dependent on the flow fields in both the liquid and gas phases at the front of the plug, where the presence of recirculating flow at the edge of the interface, close to the tube wall, causes significant drag force in the region and the rear interface to advance at a greater rate than the front. Within the parameters of the current study, it is observed that, when the plug is close to rupturing, a disk-like structure of thin liquid film is formed at the centerline of the plug. The size of the disk formation is highly dependent on the applied pressure drop. In cases of low values of applied pressure drop, the area of the disk formation is found to be relatively small; rupture occurs at a point along the edge of the disk, with no droplets forming in the process. At higher values of pressure drop, rupture occurs along the edge of the disk, thus causing the formation of a droplet. With even higher levels of applied pressure drop, rupture is found to occur along the edge of the disk, as well as at the center of the disk, resulting in the formation of multiple droplets of varying sizes. Different formations of droplets can, therefore, be attained, depending on the pressure drop available within the airway during inhalation (i.e., flow rate).

Upon rupture, a jet of air is formed at the opening of the plug, resulting in a temporary increase in shear stress at the tube wall. The magnitude and duration of this increase is found to be highly dependent on droplet formation, which causes disruptions in the flow field of the jet formation. Thus, while higher peaks in shear stress are obtained in the case of higher (applied) pressure drops, the duration of the period of increased shear stress decreases. This finding is consistent with previous clinical observation on the extent of epitheleal cell damage in the process of airway reopening (due to significant shear stress at the wall), whereby it was observed that a lower reopening velocity (lower pressure drop) causes a greater extent of damage in epitheleal cells at the airway.

The authors would like to acknowledge the remarkable contribution of the Basilisk and Gerris Flow Solver community (http://basilisk.fr/).

This work was supported by resources provided by The Pawsey Supercomputing Center with funding from the Australian Government and the Government of Western Australia.

The authors have no conflicts to disclose.

Monica M. Gumulya: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Ryan Mead-Hunter: Conceptualization (equal); Formal analysis (equal); Investigation (supporting); Writing – review & editing (equal). Benjamin J. Mullins: Conceptualization (equal); Formal analysis (equal); Resources (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Taking the reference parameters for velocity, pressure, and length, the following form of Navier–Stokes equations took place:
ρ * ( u * t * + u * · * u * ) = * P * + 1 R e * 2 u * + 1 W e κ * δ s n * ,
(A1)
ρ * t * + * · ( ρ * u * ) = 0 ,
(A2)
* · u * = 0.
(A3)

The Weber number, W e = ρ L u ref 2 x ref / σ, has been set to unity due to the form of the reference velocity and length (uref and xref) adopted in this study, thus leaving the Reynolds number and pressure-driving force as the variables for the simulation studies.

The grid dependence of the numerical framework is assessed by assigning two different levels of refinement at the interface, L 0 / 2 9 (“Ref. 9”) and L 0 / 2 10 (“Ref. 10”), for cases of Δ P * = 2 and 5 at h * = 0.05. The shape and structure of the grid at Ref. 10 can be seen in Fig. 9(a).

FIG. 9.

Comparison of plug behavior at Δ P * = 2 and 5 (at h * = 0.05), using two levels of grid refinement: (a) the grid refinement applied on the gas–liquid interface at Ref. 10, (b) plug thickness at the centerline, (c) plug propagation rate, (d) 3D reconstruction of the interface ( h * = 0.05 and Δ P * = 5) at Refs. 9 (top row) and 10 (bottom row).

FIG. 9.

Comparison of plug behavior at Δ P * = 2 and 5 (at h * = 0.05), using two levels of grid refinement: (a) the grid refinement applied on the gas–liquid interface at Ref. 10, (b) plug thickness at the centerline, (c) plug propagation rate, (d) 3D reconstruction of the interface ( h * = 0.05 and Δ P * = 5) at Refs. 9 (top row) and 10 (bottom row).

Close modal

Figures 9(b) and 9(c) show that both levels of refinement yield similar propagation rates and thickness reduction for the case considered here. In the case of Ref. 9, the interface is seen to be less smooth in comparison to Ref. 10. Based on this assessment, Ref. 10 refinement is adopted for all the cases in this study.

1.
A.
Agrawal
and
R.
Bhardwaj
, “
Probability of covid-19 infection by cough of a normal person and a super-spreader
,”
Phys. Fluids
33
,
031704
(
2021
).
2.
H.
Wang
,
Z.
Li
,
X.
Zhang
,
L.
Zhu
,
Y.
Liu
, and
S.
Wang
, “
The motion of respiratory droplets produced by coughing
,”
Phys. Fluids
32
,
125102
(
2020
).
3.
M.-R.
Pendar
and
J. C.
Páscoa
, “
Numerical modeling of the distribution of virus carrying saliva droplets during sneeze and cough
,”
Phys. Fluids
32
,
083305
(
2020
).
4.
T.
Dbouk
and
D.
Drikakis
, “
On coughing and airborne droplet transmission to humans
,”
Phys. Fluids
32
,
053310
(
2020
).
5.
S.
Ren
,
M.
Cai
,
Y.
Shi
,
Z.
Luo
, and
T.
Wang
, “
Influence of cough airflow characteristics on respiratory mucus clearance
,”
Phys. Fluids
34
,
041911
(
2022
).
6.
L.
Morawska
,
G.
Johnson
,
Z.
Ristovski
,
M.
Hargreaves
,
K.
Mengersen
,
S.
Corbett
,
C. Y. H.
Chao
,
Y.
Li
, and
D.
Katoshevski
, “
Size distribution and sites of origin of droplets expelled from the human respiratory tract during expiratory activities
,”
J. Aerosol Sci.
40
,
256
269
(
2009
).
7.
P. T.
Macklem
, “
Airway obstruction and collateral ventilation
,”
Physiol. Rev.
51
,
368
436
(
1971
).
8.
R.
Kamm
and
R.
Schroter
, “
Is airway closure caused by a liquid film instability?
,”
Respir. Physiol.
75
,
141
156
(
1989
).
9.
B.
Bake
,
P.
Larsson
,
G.
Ljungkvist
,
E.
Ljungström
, and
A.-C.
Olin
, “
Exhaled particles and small airways
,”
Respir. Res.
20
,
1
14
(
2019
).
10.
G. R.
Johnson
and
L.
Morawska
, “
The mechanism of breath aerosol formation
,”
J. Aerosol Med. Pulm. Drug Delivery
22
,
229
237
(
2009
).
11.
P.
Sajgalik
,
A.
Garzona-Navas
,
I.
Csécs
,
J. W.
Askew
,
F.
Lopez-Jimenez
,
A. S.
Niven
,
B. D.
Johnson
, and
T. G.
Allison
, “
Characterization of aerosol generation during various intensities of exercise
,”
Chest
160
,
1377
1387
(
2021
).
12.
C. M.
Orton
,
H. E.
Symons
,
B.
Moseley
,
J.
Archer
,
N. A.
Watson
,
K. E.
Philip
,
S.
Sheikh
,
B.
Saccente-Kennedy
,
D.
Costello
,
W. J.
Browne
et al, “
A comparison of respiratory particle emission rates at rest and while speaking or exercising
,”
Commun. Med.
2
,
44
(
2022
).
13.
B.
Mutsch
,
M.
Heiber
,
F.
Grätz
,
R.
Hain
,
M.
Schönfelder
,
S.
Kaps
,
D.
Schranner
,
C. J.
Kähler
, and
H.
Wackerhage
, “
Aerosol particle emission increases exponentially above moderate exercise intensity resulting in superemission during maximal exercise
,”
Proc. Natl. Acad. Sci.
119
,
e2202521119
(
2022
).
14.
W. F.
Wells
and
M. W.
Wells
, “
Air-borne infection
,”
J. Am. Med. Assoc.
107
,
1698
1703
(
1936
).
15.
R. R.
Netz
, “
Mechanisms of airborne infection via evaporating and sedimenting droplets produced by speaking
,”
J. Phys. Chem. B
124
,
7093
7101
(
2020
).
16.
R. M.
Effros
,
M. B.
Dunning
III
,
J.
Biller
, and
R.
Shaker
, “
The promise and perils of exhaled breath condensates
,”
Am. J. Physiol. Lung Cell. Mol. Physiol.
287
,
L1073
L1080
(
2004
).
17.
P.
Howell
,
S.
Waters
, and
J.
Grotberg
, “
The propagation of a liquid bolus along a liquid-lined flexible tube
,”
J. Fluid Mech.
406
,
309
335
(
2000
).
18.
H.
Fujioka
and
J. B.
Grotberg
, “
Steady propagation of a liquid plug in a two-dimensional channel
,”
J. Biomech. Eng.
126
,
567
577
(
2004
).
19.
H.
Fujioka
,
S.
Takayama
, and
J. B.
Grotberg
, “
Unsteady propagation of a liquid plug in a liquid-lined straight tube
,”
Phys. Fluids
20
,
062104
(
2008
).
20.
J.
Magniez
,
M.
Baudoin
,
C.
Liu
, and
F.
Zoueshtiagh
, “
Dynamics of liquid plugs in prewetted capillary tubes: from acceleration and rupture to deceleration and airway obstruction
,”
Soft Matter
12
,
8710
8717
(
2016
).
21.
S.
Waters
and
J.
Grotberg
, “
The propagation of a surfactant laden liquid plug in a capillary tube
,”
Phys. Fluids
14
,
471
480
(
2002
).
22.
M.
Muradoglu
,
F.
Romanò
,
H.
Fujioka
, and
J.
Grotberg
, “
Effects of surfactant on propagation and rupture of a liquid plug in a tube
,”
J. Fluid Mech.
872
,
407
437
(
2019
).
23.
D.
Halpern
,
H.
Fujioka
, and
J. B.
Grotberg
, “
The effect of viscoelasticity on the stability of a pulmonary airway liquid layer
,”
Phys. Fluids
22
,
011901
(
2010
).
24.
Y.
Hu
,
F.
Romanò
, and
J. B.
Grotberg
, “
Effects of surface tension and yield stress on mucus plug rupture: A numerical study
,”
J. Biomech. Eng.
142
,
061007
(
2020
).
25.
S. S.
Mamba
,
J.
Magniez
,
F.
Zoueshtiagh
, and
M.
Baudoin
, “
Dynamics of a liquid plug in a capillary tube under cyclic forcing: Memory effects and airway reopening
,”
J. Fluid Mech.
838
,
165
191
(
2018
).
26.
B.
Munir
and
Y.
Xu
, “
Effects of gravity and surface tension on steady microbubble propagation in asymmetric bifurcating airways
,”
Phys. Fluids
32
,
072105
(
2020
).
27.
E. A.
Hassan
,
E.
Uzgoren
,
H.
Fujioka
,
J. B.
Grotberg
, and
W.
Shyy
, “
Adaptive Lagrangian–Eulerian computation of propagation and rupture of a liquid plug in a tube
,”
Int. J. Numer. Methods Fluids
67
,
1373
1392
(
2011
).
28.
A.
Malashenko
,
A.
Tsuda
, and
S.
Haber
, “
Propagation and breakup of liquid menisci and aerosol generation in small airways
,”
J. Aerosol Med. Pulm. Drug Delivery
22
,
341
353
(
2009
).
29.
B.
He
,
C.
Qin
,
W.
Chen
, and
B.
Wen
, “
Numerical simulation of pulmonary airway reopening by the multiphase lattice Boltzmann method
,”
Comput. Math. Appl.
108
,
196
205
(
2022
).
30.
S. K.
Lai
,
Y.-Y.
Wang
,
D.
Wirtz
, and
J.
Hanes
, “
Micro-and macrorheology of mucus
,”
Adv. Drug Delivery Rev.
61
,
86
100
(
2009
).
31.
A. J.
Chorin
, “
Numerical solution of the Navier-Stokes equations
,”
Math. Comput.
22
,
745
762
(
1968
).
32.
S.
Popinet
, “
Gerris: A tree-based adaptive solver for the incompressible Euler equations in complex geometries
,”
J. Comput. Phys.
190
,
572
600
(
2003
).
33.
S.
Popinet
, “
An accurate adaptive solver for surface-tension-driven interfacial flows
,”
J. Comput. Phys.
228
,
5838
5866
(
2009
).
34.
J. U.
Brackbill
,
D. B.
Kothe
, and
C.
Zemach
, “
A continuum method for modeling surface tension
,”
J. Comput. Phys.
100
,
335
354
(
1992
).
35.
C.
Pairetti
,
R.
Villiers
, and
S.
Zaleski
, “
A numerical cough machine
,” arXiv:2101.05662 (
2021
).
36.
M.
Gumulya
,
J. B.
Joshi
,
R. P.
Utikar
,
G. M.
Evans
, and
V.
Pareek
, “
Bubbles in viscous liquids: Time dependent behaviour and wake characteristics
,”
Chem. Eng. Sci.
144
,
298
309
(
2016
).
37.
S.
Wu
,
J.
Zhang
,
Q.
Xiao
, and
M.-J.
Ni
, “
Comparison of two interfacial flow solvers: Specific case of a single droplet impacting onto a deep pool
,”
Comput. Math. Appl.
81
,
664
678
(
2021
).
38.
F.
Romanò
,
H.
Fujioka
,
M.
Muradoglu
, and
J.
Grotberg
, “
Liquid plug formation in an airway closure model
,”
Phys. Rev. Fluids
4
,
093103
(
2019
).
39.
P.
Schwartz
,
M.
Barad
,
P.
Colella
, and
T.
Ligocki
, “
A cartesian grid embedded boundary method for the heat equation and Poisson's equation in three dimensions
,”
J. Comput. Phys.
211
,
531
550
(
2006
).
40.
A. M.
Bilek
,
K. C.
Dee
, and
D. P.
Gaver
III
, “
Mechanisms of surface-tension-induced epithelial cell damage in a model of pulmonary airway reopening
,”
J. Appl. Physiol.
94
,
770
783
(
2003
).
41.
D.
Huh
,
H.
Fujioka
,
Y.-C.
Tung
,
N.
Futai
,
R.
Paine
III
,
J. B.
Grotberg
, and
S.
Takayama
, “
Acoustically detectable cellular-level lung injury induced by fluid mechanical stresses in microfluidic airway systems
,”
Proc. Natl. Acad. Sci.
104
,
18886
18891
(
2007
).