Mean streaming and the onset of turbulence in the reciprocating flow in a double bifurcation airway model: Insights for high-frequency ventilation

This paper reports on a numerical study of the reciprocating flow through a double bifurcation. Its aim is to quantify the potential gas transport in flows similar to those in the human airway during high-frequency ventilation, a medical ventilation technique that uses fast, yet shallow, inflations to minimize pressure and volume variation in the lungs and therefore protect lungs from ventilator-induced lung injury. The shallow inflations mean that gas transport is not achieved by bulk advection and other mechanisms must be used. Here we focus on nonlinear mean streaming and turbulent diffusion. We show that both of these mechanisms are driven by the formation of Dean vortices due to centrifugal instability in the curved sections of the bifurcation. We report the impact of the upstream and downstream conditions on both of these during both inhalation and exhalation portions of the reciprocating flow; in particular, we show that the development of turbulence in a given airway vessel is influenced by downstream conditions and not simply by the local flow conditions. This result highlights the importance of using a geometry of appropriate complexity for the modeling of physiological flows, particularly if the results are to be parameterised for use in larger-scale but lower-order models.

Figure 1

Attributes of the bifurcation model used in the paper. (a) Dimensions of the planar 1:2:4 bifurcation. (b) Three-dimensional CAD model of the 1:2:4 bifurcation. Note that the CAD model shown here is not to scale and the actual computational domain has much longer inlet and outlet sections to minimize the effect of boundary conditions on the flow generated by the bifurcation. $D_0$ corresponds to the diameter of the endotracheal tube of $3\text{mm}$.

Figure 2

Visualization of dimensionless mean streaming velocity ${\overline{\mathbf{U}}}_{\mathrm{ax}}/{\widehat{U}}_{\text{max}}$ fields for (a) ${}^{0}G_2$, (${\text{Re}}_{\text{max}}=9600, \alpha =3.36, L/a=849$); (b) ${}^{1}G_3$, (${\text{Re}}_{\text{max}}=6076, \alpha =2.66, L/a=861$); (c) ${}^{2}G_4$, (${\text{Re}}_{\text{max}}=3846, \alpha =2.1, L/a=873$); and (d) ${}^{8}G_{10}$, (${\text{Re}}_{\text{max}}=247, \alpha =0.51, L/a=950$). The first panel of each case shows the longitudinal symmetry plane along with several transverse planes. The second panel shows the same field without the longitudinal symmetry plane.

Figure 3

Visualization of dimensionless axial vorticity ${\overline{\omega }}_{\mathrm{ax}}/{\overline{\omega }}_{\max }$ (where ${\overline{\omega }}_{\max }={\widehat{U}}_{\text{max}}/D_0$), computed based on the mean streaming velocity ${\overline{\mathbf{U}}}_{\mathrm{ax}}$ fields for (a) ${}^{0}G_2$, (${\text{Re}}_{\text{max}}=9600, \alpha =3.36, L/a=849$); (b) ${}^{1}G_3$, (${\text{Re}}_{\text{max}}=6076, \alpha =2.66, L/a=861$); (c) ${}^{2}G_4$, (${\text{Re}}_{\text{max}}=3846, \alpha =2.1, L/a=873$); and (d) ${}^{8}G_{10}$, (${\text{Re}}_{\text{max}}=247, \alpha =0.51, L/a=950$).

Figure 4

Contours of dimensionless streaming velocity ${\overline{\mathbf{U}}}_{\mathrm{ax}}/{\widehat{U}}_{\text{max}}$ fields on cross sections A-${\mathrm{A}}^{\prime }$, B-${\mathrm{B}}^{\prime }$, C-${\mathrm{C}}^{\prime }$, and D-${\mathrm{D}}^{\prime }$ superposed with secondary flow vectors ${\mathbf{U}}_{\mathrm{tr}}$. Note the vectors are uniform length, not scaled by magnitude, and so indicate the flow direction only. Cross sections A-${\mathrm{A}}^{\prime }$, B-${\mathrm{B}}^{\prime }$, C-${\mathrm{C}}^{\prime }$, and D-${\mathrm{D}}^{\prime }$ are located at the middle of the straight section at each generation level.

Figure 5

Dimensionless recirculating flux computed at the borderline between the bifurcating region 1 and intermediate (daughter) branch as a function of ${\text{Re}}_{\text{max}}$ for the first nine intermediate generations modeling the airway. (a) The red curve shows a power law fit for the first five intermediate generations: $Q_r=Q_{\text{max}}[0.0034{\left(R{e}_{\max }\right)}^{1/4}-0.02]$. (b) The geometric location of the measurements (referred to as location 1 hereafter).

Figure 6

Dimensionless recirculating flux computed at the borderline between the intermediate (daughter) branch and bifurcating region 2 as a function of ${\text{Re}}_{\text{max}}$ for the first nine intermediate generations modeling the airway. (a) The red curve shows a power law fit for the first five intermediate generations: $Q_r=Q_{\text{max}}[0.0012{\left({\text{Re}}_{\text{max}}\right)}^{1/4}-0.0072]$. (b) The geometric location of the measurements (referred to as location 2 hereafter).

Figure 7

Strength of the recirculating flux measured at the start and end of the intermediate generation level as a function of the generation number. Blue curves correspond to the strength of recirculating flux directly measured at a location of interest in the computational domain. Orange curves correspond to the total recirculating flux computed at each generation level. Dashed lines correspond to the measurements obtained at the borderline between the bifurcating region 1 and intermediate branch (location 1) and solid lines correspond to the measurements obtained at the borderline between the intermediate branch and bifurcating region 2 (location 2). The horizontal gray line corresponds to the flow rate required to meet the oxygen demand of a neonate in the clinical conditions based on the data reported by Hill and Robinson [42].

Figure 8

Isosurfaces of ${\lambda }_2$ vortical structures based on the threshold value of $-0.0025$. (a1), (a2), and (a3) correspond to isosurfaces generated based on the streaming velocity field of ${}^{0}G_2, {}^{1}G_3$, and ${}^{2}G_4$ cases, respectively. Similarly, (b1), (b2), and (b3) correspond to isosurfaces generated based on the peak inhalation flow field of ${}^{0}G_2, {}^{1}G_3$, and ${}^{2}G_4$ cases, respectively. Further, (c1), (c2), and (c3) correspond to isosurfaces generated based on the peak exhalation flow field of ${}^{0}G_2, {}^{1}G_3$, and ${}^{2}G_4$ cases, respectively. The animation |lam_0G2.avi| provided in the Supplemental Material shows the evolution of isosurfaces of ${\lambda }_2$ vortical structures pertinent to ${}^{0}G_2$ case.

Figure 9

Evolution of the instantaneous velocity magnitude, $\left|\mathbf{U}\right|$ for ${}^{0}G_2$ (${\text{Re}}_{\text{max}}=9600, \alpha =3.36$). The dotted lines highlight areas where turbulent bursts have developed in the flow. The animation |vel_0G2.avi| provided in the Supplemental Material shows the evolution of $\left|\mathbf{U}\right|$.

Figure 10

Evolution of the instantaneous velocity magnitude, $\left|\mathbf{U}\right|$ for ${}^{1}G_3$ (${\text{Re}}_{\text{max}}=6076, \alpha =2.66$).

Figure 11

A decomposition of the time history of axial velocity ${\mathbf{U}}_{\text{ax}}$ of the ${}^{0}G_2$ case. (a) Decomposition of temporal axial velocity (${\mathbf{U}}_{\text{ax}}$) at the center of the cross section A-${\mathrm{A}}^{\prime }$. (b) Decomposition of temporal axial velocity (${\mathbf{U}}_{\text{ax}}$) at the center of the cross section B-${\mathrm{B}}^{\prime }$. (c) Decomposition of temporal axial velocity (${\mathbf{U}}_{\text{ax}}$) at the center of the cross section C-${\mathrm{C}}^{\prime }$. (d) Decomposition of temporal axial velocity (${\mathbf{U}}_{\text{ax}}$) at the center of the cross section D-${\mathrm{D}}^{\prime }$. These cross sections are shown in Fig. 4. Solid black line shows the measured temporal data from the instantaneous velocity field; solid green line shows the velocity scale associated with the filtered low-frequency components while the solid red line shows the velocity scale associated with the filtered high-frequency components. The cutoff frequency for the low- and high-pass filters was $f_{\text{thr}}{D}_0/{\widehat{U}}_{\text{max}}=0.05$.

Figure 12

A decomposition of the time history of axial velocity ${\mathbf{U}}_{\text{ax}}$ of the ${}^{1}G_3$ case. (a) Decomposition of temporal axial velocity (${\mathbf{U}}_{\text{ax}}$) at the center of the cross section A-${\mathrm{A}}^{\prime }$. (b) Decomposition of temporal axial velocity (${\mathbf{U}}_{\text{ax}}$) at the center of the cross section B-${\mathrm{B}}^{\prime }$. (c) Decomposition of temporal axial velocity (${\mathbf{U}}_{\text{ax}}$) at the center of the cross section C-${\mathrm{C}}^{\prime }$. (d) Decomposition of temporal axial velocity (${\mathbf{U}}_{\text{ax}}$) at the center of the cross section D-${\mathrm{D}}^{\prime }$. Solid black line shows the measured temporal data from the instantaneous velocity field; solid green line shows the velocity scale associated with the filtered low-frequency components while the solid red line shows the velocity scale associated with the filtered high-frequency components. The cutoff frequency for the low- and high-pass filters was $f_{\text{thr}}{D}_0/{\widehat{U}}_{\text{max}}=0.05$.

Figure 13

A decomposition of the time history of axial velocity ${\mathbf{U}}_{\text{ax}}$ of the single bifurcation ${}^{0}G_1$ and ${}^{1}G_2$ cases. ${}^{0}G_1$ generation 0: Decomposition of ${\mathbf{U}}_{\text{ax}}$ at the geometric center of the mother branch of ${}^{0}G_1$ case. ${}^{0}G_1$ generation 1: Decomposition of ${\mathbf{U}}_{\text{ax}}$ at the geometric center of a daughter branch of ${}^{0}G_1$ case. ${}^{1}G_2$ generation 1: Decomposition of ${\mathbf{U}}_{\text{ax}}$ at the geometric center of the mother branch of ${}^{1}G_2$ case. ${}^{1}G_2$ generation 2: Decomposition of ${\mathbf{U}}_{\text{ax}}$ at the geometric center of a daughter branch of ${}^{1}G_2$ case. Solid black line shows the measured temporal data from the instantaneous velocity field; solid green line shows the velocity scale associated with the filtered low-frequency components, while the solid red line shows the velocity scale associated with the filtered high-frequency components. The cutoff frequency for the low- and high-pass filters was $f_{\text{thr}}{D}_0/{\widehat{U}}_{\text{max}}=0.05$.

Figure 14

Geometric locations compared in various simulation configurations shown in Fig. 15. (a1), (a2), (b1), and (b2) are dynamically equivalent location pairs located at the geometric center of the mother branch of ${}^{1}G_2$ case, the left daughter branch (section B-${\mathrm{B}}^{\prime }$ of Fig. 4) of ${}^{0}G_2$ case, the left daughter branch of ${}^{1}G_2$ case, and the upper left granddaughter branch (section D-${\mathrm{D}}^{\prime }$ of Fig. 4) of ${}^{0}G_1$ case, respectively. The scaling factors used to scale up each cross section to make them identical are shown.

Figure 15

A comparison of the evolution of the instantaneous velocity magnitude, $\left|\mathbf{U}\right|$, among dynamically equivalent flow configurations shown in Fig. 14. Rows (a1), (a2), (b1), and (b2) display 12 snapshots of $\left|\mathbf{U}\right|$ at cross sections located at the geometric center of the mother branch of ${}^{1}G_2$ case, the left daughter branch (section B-${\mathrm{B}}^{\prime }$ of Fig. 4) of ${}^{0}G_2$ case, the left daughter branch of ${}^{1}G_2$ case, and the upper left granddaughter branch (section D-${\mathrm{D}}^{\prime }$ of Fig. 4) of ${}^{0}G_1$ case, respectively. Note that the diameter of daughter and granddaughter branches are scaled up by a factor of 0.79 and $0.{79}^2$ to make them identical for comparison purposes. The top panel shows the variation of flow rate as a function of time.

Figure 16

The mesh topology of a bifurcating junction used in the study. (a) Macro mesh topology of a bifurcating junction. (b) Resulting mesh topology of hexahedral elements associated with Gauss-Legendre-Lobatto quadrature points of eighth-order tensor-product Lagrange polynomials.

Figure 17

The variation of energy norm ($E_{\text{norm}}$) as the global indicator of convergence. (a) The temporal evolution of $E_{\text{norm}}$ over a portion of the reciprocating flow cycle, i.e., $0\le t\le \frac{2T}{3}$, where $T$ is the reciprocation period; for varying interpolation polynomial order, i.e., $P=8,9,10,12$. (b) The variation of maximum energy norm $\text{max}|E_{\text{norm}}|$ as a function of interpolation polynomial order $P$. The red curve shows an exponential convergence of $\text{max}|E_{\text{norm}}|$ as a function of $P$: $\text{max}|E_{\text{norm}}|=0.464+\frac{1.216}{e^{0.441P}}$.

Figure 18

Frequency spectrum of the instantaneous axial (streamwise) velocity ${\mathbf{U}}_{\text{ax}}$ time history. (a)–(c) Frequency spectrum of the ${\mathbf{U}}_{\text{ax}}$ time history measured at the geometric center of generation 1 (geometric center of section B-${\mathrm{B}}^{\prime }$) for $P=8,10,12$. The cutoff frequency for the high-pass filter was $f_{\text{thr}}{D}_0/{\widehat{U}}_{\text{max}}=0.05$. The peak dimensionless oscillation frequency $f_{\text{max}}{D}_0/{\widehat{U}}_{\text{max}}$ is marked with a red dashed line.

Figure 19

A comparison of instantaneous velocity magnitude $\left|\mathbf{U}\right|$ and instantaneous axial vorticity ${\omega }_{\mathrm{ax}}$ at the peak flow rate indicated by the vertical grey line of Fig. 17 for varying interpolation polynomial order ($P$). (a1)–(a3) Snapshots of $\left|\mathbf{U}\right|$ at the peak flow rate for $P=8,9,10$, respectively. (b1)–(b3) Snapshots of ${\omega }_{\mathrm{ax}}$ at the peak flow rate for $P=8,9,10$, respectively. All the snapshots are taken from the vertical plane passing through the geometric center of generation 1 (section B-${\mathrm{B}}^{\prime }$).