Liquid plug formation in an airway closure model

The closure of a human lung airway is modeled as an instability of a two-phase flow in a pipe coated internally with a Newtonian liquid. For a thick enough coating, the Plateau-Rayleigh instability creates a liquid plug which blocks the airway, halting distal gas exchange. Owing to a bifrontal plug growth, this airway closure flow induces high stress levels on the wall, which is the location of airway epithelial cells. A parametric numerical study is carried out simulating relevant conditions for human lungs, in either ordinary or pathological situations. Our simulations can represent the physical process from pre- to postcoalescence phases. Previous studies have been limited to precoalescence only. The topological change during coalescence induces a high level of stress and stress gradients on the epithelial cells, which are large enough to damage them, causing sublethal or lethal responses. We find that postcoalescence wall stresses can be in the range of 300% to 600% greater than precoalescence values and so introduce an important source of mechanical perturbation to the cells.

Figure 1

Schematic of the airway closure process. The air-liquid interface is depicted at four instants of time: for $t=t_1$ the interface is slightly deformed, for $t=t_2$ the deformation has increased because of the Plateau-Rayleigh instability, for $t=t_c$ it coalesces while undergoing a topological change, and for $t=t_4$ the newly formed liquid plug undergoes a bifrontal growth.

Figure 2

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

Figure 3

Comparison of the two-phase distribution at different times: first (a)–(c) and third (g)–(i) rows from the top are the experimental visualizations reported in Ref. [37] using fluorescent particles, whereas second (d)–(f) and fourth (j)–(l) rows from the top are our numerical results for the corresponding time relative to the coalescence event. The annular film is in red, while the core fluid is in black. The nondimensional parameters of interest are $\varepsilon =0.23, \text{La}=0.072, \mu =0.01, \varrho =0.95$, and $\lambda =9.2$ and the closure event is depicted in the axial range $-1.6\le z\le 1.6$. Time is scaled as in Ref. [37] by ${\mu }_{\text{L}}a/{\varepsilon }^3\sigma$. (a)–(f) The precoalescence phase due to the interfacial instability; (g)–(l) the postcoalescence phase, where a bifrontal plug growth is observed.

Figure 4

Velocity vectors (arrows) and streamlines for $\varepsilon =0.23, \text{La}=0.072, \mu =0.01, \varrho =0.95$, and $\lambda =9.2$: first (a)–(b) and third (e)–(f) rows from the top are the micro-PIV measurements of Ref. [37], whereas second (c)–(d) and fourth (g)–(h) rows from the top are numerical results. Time is scaled by ${\mu }_{\text{L}}a/{\varepsilon }^3\sigma$.

Figure 5

Temporal evolution of the wall shear stress between $-0.8\le z\le 0.8$: the solid black lines denote the numerical results for $z=\pm 0.8$, the dashed black lines for $z=\pm 0.4$, while the dashed-dotted line refers to $z=0$; markers refer to the experimental results of Ref. [37] as indicated in the legend. The two colored dashed lines denote the closure time ($t_c$) and the numerical limit ($t_s$) of Ref. [38]. In these results the numerical closing time $t_c$ is used as the reference time, and the experimental results are thus shifted in time to match with it, considering the scaling employed in (2).

Figure 6

Pressure field (colors), streamlines (yellow lines) and axial velocity in the liquid phase (arrows) at two precoalescence [(a) $t=600$ and (b) $t=736$] and two postcoalescence instants [(c) $t=737$ and (d) $t=800$]. The turquoise line denotes the interface location. The phenomenon is symmetric in $z$.

Figure 7

Bottom panel: Location of the liquid-gas interface at significant instants in time. Top-right panel: wall shear stress ${\tau }_{\text{w}}$ and wall pressure gradient ${\partial }_z{p}_{\text{w}}$ depicted with the same color coding of the time legend. Top-left panel: minimum of the interface radial location $R_{\text{min}}=min\left(R_{\text{I}}\right)$ (solid black line), asymptotic power laws predicted in Ref. [38] (dashed black line), excursions of wall pressure $\mathrm{\Delta }{p}_{\text{w}}=max\left(p_{\text{w}}\right)-min\left(p_{\text{w}}\right)$ and shear stress $\mathrm{\Delta }{\tau }_{\text{w}}=max\left({\tau }_{\text{w}}\right)-min\left({\tau }_{\text{w}}\right)$ depicted as solid and dashed lines without markers, maximum absolute values of wall pressure and wall shear stress gradients, $|{\partial }_z{p}_{\text{w}}{|}_{\text{max}}$ and $|{\partial }_z{\tau }_{\text{w}}{|}_{\text{max}}$, respectively, depicted by solid and dashed red lines marked with bullets. The nondimensional groups are $\text{La}=200, \mu =1.5\times {10}^{-3}, \varrho ={10}^{-3}, \lambda =6$, and $\varepsilon =0.25$.

Figure 8

Minimum of the interface radial location $R_{\text{min}}=min\left(R_{\text{I}}\right)$ (solid black line), excursions of wall pressure $\mathrm{\Delta }{p}_{\text{w}}=max\left(p_{\text{w}}\right)-min\left(p_{\text{w}}\right)$ and shear stress $\mathrm{\Delta }{\tau }_{\text{w}}=max\left({\tau }_{\text{w}}\right)-min\left({\tau }_{\text{w}}\right)$ depicted as solid and dashed lines without markers, and maximum absolute values of wall pressure and wall shear stress gradients, $|{\partial }_z{p}_{\text{w}}{|}_{\text{max}}$ and $|{\partial }_z{\tau }_{\text{w}}{|}_{\text{max}}$, respectively, depicted by solid and dashed red lines marked with bullets. The nondimensional groups are $\text{La}=100, \mu =1.5\times {10}^{-3}, \varrho ={10}^{-3}, \lambda =6$, and $\varepsilon =0.25$.

Figure 9

The effects of the liquid viscosity. The time evolution of the minimum radius of the interface, $R_{\text{I}}, R_{\text{min}}$, is plotted for $\varrho ={10}^{-3}, \lambda =6, \varepsilon =0.25$, and $\text{La}=1$, and $\mu =1.5\times {10}^{-4}$ (solid line), $\text{La}=4$, and $\mu =3\times {10}^{-4}$ (dashed line), $\text{La}=100$, and $\mu =1.5\times {10}^{-3}$ (dashed-dotted line).

Figure 10

The effects of the liquid viscosity on the time evolution of the excursion of wall stresses $\mathrm{\Delta }{p}_{\text{w}}$ (higher-amplitude curves) and $\mathrm{\Delta }{\tau }_{\text{w}}$ (lower-amplitude curves) for $\varrho ={10}^{-3}, \lambda =6, \varepsilon =0.25$ and $\text{La}=1, \mu =1.5\times {10}^{-4}$ (solid line), $\text{La}=4$ and $\mu =3\times {10}^{-4}$ (dashed line), $\text{La}=100$ and $\mu =1.5\times {10}^{-3}$ (dashed-dotted line).

Figure 11

The effects of the initial liquid film thickness on the minimum radius of the interface $R_{\text{min}}$ (solid black line) and excursions of wall pressure $\mathrm{\Delta }p=max\left(p_{\text{w}}\right)-min\left(p_{\text{w}}\right)$ and shear stress $\mathrm{\Delta }\tau =max\left({\tau }_{\text{w}}\right)-min\left({\tau }_{\text{w}}\right)$ depicted as solid and dashed lines without markers, and the maximum absolute values of wall pressure and shear stress gradients, $|{\partial }_z{p}_{\text{w}}{|}_{\text{max}}$ and $|{\partial }_z{\tau }_{\text{w}}{|}_{\text{max}}$, respectively, depicted by solid and dashed red lines marked with bullets. The nondimensional groups in common are $\text{La}=1, \mu =1.5\times {10}^{-3}, \varrho ={10}^{-3}$, and $\lambda =6$; the three sets of curves refer to $\varepsilon =0.2$, 0.25, and 0.3.

Figure 12

The time evolutions of the minimum radius of the interface, $R_{\text{min}}$ (black line), the minimum axial distance between the fronts, $Z_{\text{min}}$ (black line with bullets), the minimum film thickness, $1-max\left(R_{\text{I}}\right)$ (turquoise solid line), its axial location for $z<3, Z_{\text{mt}}$ (turquoise solid line with bullets), and the maximum wall pressure gradient $|{\partial }_z{p}_{\text{w}}{|}_{\text{max}}$. The nondimensional groups are $\text{La}=1, \mu =1.5\times {10}^{-3}, \varrho ={10}^{-3}, \lambda =6$, and $\varepsilon =0.25$.

Figure 13

Normalized wall shear stress ${\tilde{\tau }}_{\text{w}}={\tau }_{\text{w}}/max\left({\tau }_{\text{w}}\right)$ at $t=t_{\text{s}}$ against the normalized coordinate $\eta =z/z_{{\tau }_{\text{w,max}}}$. The present results (solid line) are compared to the results in Ref. [38] (dashed line).

Figure 14

Minimum of the interface radial location $R_{\text{min}}=min\left(R_{\text{I}}\right)$, excursions of wall pressure $\mathrm{\Delta }{p}_{\text{w}}=max\left(p_{\text{w}}\right)-min\left(p_{\text{w}}\right)$ and shear stress $\mathrm{\Delta }{\tau }_{\text{w}}=max\left({\tau }_{\text{w}}\right)-min\left({\tau }_{\text{w}}\right)$ computed with basilisk (lines) and with the numerical solver of Ref. [55] (markers). The nondimensional groups are $\text{La}=100, \lambda =6$, and $\varepsilon =0.25$.