Effects of varying inhalation duration and respiratory rate on human airway flow

Studies of flow through the human airway have shown that inhalation time (IT) and secondary flow structures can play important roles in particle deposition. However, the effects of varying IT in conjunction with respiratory rate (RR) on airway flow remain unknown. Using 3D numerical simulations of oscillatory flow through an idealized airway model consisting of a mouth inlet, glottis, trachea and symmetric double bifurcation at trachea Reynolds number ($Re$) of 4,200, we investigated how varying the ratio of IT to breathing time (BT) from 25% to 50% and RR from 10 breaths per minute (bpm) corresponding to Womersley number ($Wo$) of 2.37 to 1,000 bpm ($Wo$=23.7) impacts airway flow characteristics. Irrespective of IT/BT, axial flow during inhalation at tracheal cross-sections was non-uniform for $Wo$=2.37 as compared to centrally concentrated distribution for $Wo$=23.7. For a given $Wo$ and IT/BT, both axial and secondary (lateral) flow components unevenly split between left and right branches of a bifurcation. Irrespective of $Wo$, IT/BT and airway generation, lateral dispersion was stronger than axial flow streaming. Despite left-right symmetry of the lower airway in our model, the right-sided mouth-to-glottis portion generated turbulence in the upper airway. Varying IT/BT for a given $Wo$ did not noticeably change flow characteristics. Discrepancy in the oscillatory flow relation $Re$/$Wo^2$=2$L$/$D$ ($L$=stroke length; $D$=trachea diameter) was observed for IT/BT$\neq$50%, as $L$ changed with IT/BT. We developed a modified dimensionless stroke length term including IT/BT. While viscous forces and convective acceleration were dominant for lower $Wo$, unsteady acceleration was dominant for higher $Wo$.


I. INTRODUCTION
Flow through the human airway is characterized by regions of flow separation and the formation of secondary flow structures 1,2 . The oscillatory nature of airflow can facilitate gas exchange in higher generations of the human airway 1,3 . Normal respiratory rate (RR) in humans ranges between 10-15 breaths per minute (bpm), corresponding to a Womersley number (Wo) range of 2.36-2.89 as defined using the following relation: where ω = 2π/BT is the angular frequency based on breathing time (BT), D is the tracheal diameter, and ν is the kinematic viscosity of air. RR increases during exercise and also in mechanical ventilation strategies such as high frequency oscillatory ventilation (HFOV) [4][5][6] . Further, inhalation time (IT) is about 45% of BT under normal breathing conditions 4 . Characterizing the airflow inside the human airway during the entire breathing cycle is important for understanding particle deposition characteristics in lung diseases (e.g., chronic obstructive pulmonary disease or COPD), during the use of e-cigarettes, and in inhaled drug therapy 7,8 . Previous studies of inhaled drug delivery in human subjects with COPD and asthma have shown improved respiratory outcomes 9 and increased aerosol deposition with increased inhalation time 10 . However, computational fluid dynamics (CFD) studies examining particle deposition relevant to aerosols 11,12 and e-cigarettes 11,13,14 typically consider either steady flow rates or only the inhalation phase for simulations. Studies examining unsteady breathing patterns over an entire breathing cycle are limited.
Previous studies examining the fluid dynamics of deep inspiration were conducted by prescribing steady inhalation flow through idealized 15 and anatomically accurate physical models 16 . These studies identified both axial (streamwise) and secondary (transverse) flow dispersion to be effective transport mechanisms. Unsteady physiological flow through subject-specific airways has also been experimentally investigated [17][18][19][20] . Jalal et al. 20 found that with respect to idealized airway models 15 , a realistic airway geometry produced stronger secondary flows as well as increased axial and secondary flow dispersion. Secondary flows in realistic airways propagated deeper in the bronchial tree, and were stronger during exhalation as compared to inhalation 20 . Jalal et al. 21 conducted magnetic resonance velocimetry (MRV) experiments on an idealized double bifurcation airway model across a range of Re and Wo in the convective region of the flow regime diagram of Jan et al. 3 . Strong secondary flows were observed for Wo≥6 during crossover from inhalation to exhalation phase. Große et al. 18 conducted particle image velocimetry (PIV) measurements of steady and oscillating flows at the first bifurcation of an anatomical silicone model of the human airway. At steady inhalation, they found that the asymmetrical bifurcation promoted the formation of flow structures that were responsible for continuous streamwise transport to the lung 18 .
Oscillating flow studies showed that the size of secondary flow structures strongly depended on the instantaneous values of Re and Wo 18 . 2D PIV measurements of oscillating flow through an asymmetric idealized model (based on Weibel et al. 22 and Horsfield et al. 23 ) at normal and HFOV conditions showed mass exchange at higher frequencies 17 . Experimental studies of HFOV flow through a subject-specific airway model 24,25 reported homogeneous ventilation at higher generations with increasing RR (and hence Wo). However, the above studies of oscillatory flow through the human airway only considered IT/BT= 50%. While IT/BT varies in normal conditions and is not strictly equal to 50% (e.g., 46% 4  HFOV was proposed by Lunkenheimer et al. 26 to generate both active inspiration and expiration for eliminating entrained gas and gas decompression in airway models. HFOV technique has since been used to overcome the injuries caused due to the use of conventional ventilation for patients with acute lung injury and acute respiratory distress syndrome (ARDS) 1 . Fredburg et al. 27 observed changes in local air distribution and pressure variations at resonant frequency with the use of HFOV, and speculated that lung ventilation can be controlled by varying frequencies. Clinical oxygen delivery 28

II. METHODS
An idealized airway model 22 with the addition of idealized mouth-glottis (previously used in a particle deposition study 29 ) was designed in Solidworks (Dassault Systèmes SolidWorks Corporation, Waltham, MA) ( Figure 1A). The airway includes the following: (i) a circular mouth-inlet of 20 mm diameter; (ii) a circular glottis of 8 mm diameter; (iii) a smooth trachea with a uniform, circular cross-section of 120 mm length and diameter (D)=18 mm; and (iv) second generation (G) with G1 and G2 diameters of 12.2 and 8.3 mm, respectively. The lengths of G1 and G2 were equal to 47 mm and 19 mm, respectively. All the bifurcation angles of the airway geometry were identically equal to 70 • . The airway geometry was symmetric with respect to the coronal plane (defined along x-y plane at z=0 m, Figure 1A). Meshing was performed in ANSYS 2019 R3 (AN-SYS, Inc., PA, USA). The geometry file was imported to Fluent solver and three meshes were generated consisting of tetrahedral cells with prism layers on walls. y + for all the meshes were <5, corresponding to peak velocity condition in the trachea. k-ω SST turbulence model was used for all simulations performed in this study.
All simulations in ANSYS Fluent were set to transient conditions. Nine oscillatory profiles of time-varying inflow velocity (u), prescribed perpendicular to the inlet ( Figure 1A), were developed and prescribed as initial conditions. Sinusoidal profiles were used as a simplified representation of realistic breathing waveforms 5 where U T (=3.39 m s −1 ) denotes mean flow speed in the trachea ( Figure 1B). Inlet velocity magnitude during peak inhalation and peak exhalation phases was maintained constant at 3.39 m s −1 across all test conditions (i.e., Wo and IT/BT). After prescribing the velocity profile, each simulation was standard initialized with initial x-velocity, y-velocity, z-velocity to identically be equal to 0 m s −1 . Simulations were conducted at the High Performance Computing Center at Oklahoma State University using 16 processor units. The solution method included a coupled scheme, second-order for spatial discretization, and first-order upwind scheme for turbulence kinetic energy. Uniform time step of 0.001 s was used for running simulations for one breathing cycle. All the processing results were auto-saved at 12.5% inhalation time to maintain uniformity across all breathing frequencies and IT/BT ratios.
Mesh independence tests were performed on the airway model for different mesh sizes (Table I) at a time step of 0.001 s for HFOV velocity profiles (Re=4200, Wo=23.7, IT/BT=50%) for one breathing cycle. This particular case (1,000 bpm, Wo=23.7) with IT/BT = 50% was used in mesh independence tests to decrease the solution time. 3D velocity magnitude (|u|) was extracted from each converged simulation along a line 1-1 at the glottis in the coronal plane (z=0 m; see inset of Figure 1C). Figure 1C shows the extracted velocity magnitude along the line 1-1' for mesh 1 (blue), mesh 2 (red), and mesh 3 (black) as a function of non-dimensional diameter y/d, where d=8 mm is the glottis diameter. The velocity of mesh 2 was different compared to mesh 1 as the smaller element size increased the spatial resolution. However, the spatial resolution of mesh 2 was nearly the same as that of mesh 3, but the simulation time of mesh 3 was four times longer than mesh 2 due to the larger number of cells. Since the velocity profiles at 1-1 for mesh 2 and 3 were nearly the same ( Figure 1A), mesh 2 was selected for use in this study.    Integral parameters (I 1 , I 2 ) have been used in previous studies of flow through the human airway to examine the relative importance of axial flow streaming and lateral dispersion mechanisms 15,19 .
Axial (i.e., streamwise direction) flow uniformity at a cross-sectional plane was quantified using integral parameter I 1 15,19 : where u is the 3D velocity field for a specific phase-point in the cycle, n is the unit normal vector at a given cross-section of the airway, and A is the cross-sectional area. Lateral dispersion arising on account of secondary (i.e., transverse direction) flow at a cross-sectional plane was quantified using integral parameter I 2 15,19 : I 1 and I 2 were calculated to examine the effects of varying Wo and IT/BT at G0, G1 and G2 ( Figure 1A, Table II).
Turbulence statistics, including turbulence kinetic energy (k) and turbulence intensity (I), were characterized along the sagittal plane for the different test conditions at multiple phases per condition. Starting from the instantaneous 3D velocity vector field u(x,y,z,t), defined as: whereî,ĵ andk are the unit normal vectors along x, y, z coordinates, respectively, Reynolds decomposition was performed to separate time-averaged and fluctuating velocities as follows: w(x, y, z,t) = w(x, y, z) + w (x, y, z,t) where u, v and w represent time-averaged velocity components along x, y and z coordinates, respectively, and u , v , and w are the fluctuating velocity components along x, y and z coordinates, respectively. At each phase, k was calculated using the following equation: where u , v and w represent x, y and z components, respectively, of the fluctuating velocity. In addition, I for a specific phase of a breathing cycle was defined as k and I were both calculated using the above equations in ANSYS Fluent. Time-averaging needed for k and I calculations was performed across multiple iterations of a specific phase of the breathing cycle, after the convergence of the solution for that phase. Visualization of k and I contours were performed using Tecplot 360 software after averaging over one complete cycle.

III. RESULTS
A. Inhalation flow fields with varying breathing frequency Figure 2 shows the detailed view of the flow field at different planes along the trachea (G0) and G1 specified in Table II (plane 1). Secondary flow was more uniform at plane 1 for Wo=23.7 as compared to the pat- shown in Figure 5 for Wo=2.37 (top), and Wo=23.7 (bottom). For Wo=2.37 and IT/BT=25%, two counter-rotating vortex pairs were observed at plane 3 during phase D (θ =225 • , 25% ET) along with strong axial flow near the wall. Axial velocity increased in magnitude at peak exhalation (phase E) and later decreased at phase F for Wo=2.37 (IT/BT=25%). Secondary flow pattern at plane 3 remained the same across phases D-F for Wo=2.37 (IT/BT=25%). Similar flow fields were observed at plane 4 during exhalation at Wo=2.37 (IT/BT=25%).
For Wo=23.7 (IT/BT=25%), axial flow increased at plane 3 from phases D to E and then decreased at phase F. Secondary flow was observed at plane 3 for Wo=23.7 (IT/BT=25%) from center to outward direction at phase D, followed by formation of two counter-rotating vortex pairs at phases E and F. At peak exhalation (phase E), axial velocity slightly increased in magnitude at plane 4 as compared to that of plane 3. Similar to plane 3, two counter-rotating vortex pairs were also observed near the wall at plane 4 at peak exhalation (phases E and F) for Wo=23.7 (IT/BT=25%).
During exhalation at Wo=2.37 (IT/BT=25%) for a given phase point, the zone with higher axial flow streaming increased with increase in generation number (e.g., compare across planes 1-3).    Figure 5) to 50% (bottom half of Figure S4 in supplementary material).

D. Integral parameters
To quantify axial flow streaming and strength of secondary flow, integral parameters 19 were calculated using equations (3) and (4) for all Wo and IT/BT conditions at the planes in Table II.
During inhalation or exhalation for a given Wo and IT/BT, both the integral parameters I 1 (Figure 8A) and I 2 ( Figure 8B) at plane 1 followed a periodic trend such that they increased in value until peak inhalation or exhalation and subsequently decreased. Across the entire breathing cycle (θ =0-360 • ), maxima of I 1 and I 2 both occurred for Wo=23.7 at IT/BT=25% during peak inhalation (θ =90 • ). This is in line with the flow field observations discussed previously, where axial and secondary flows were generally stronger during inhalation as compared to exhalation. Specific to exhalation, highest axial flow streaming (I 1 ) was observed for Wo=23.7 at IT/BT=50%. Across all Wo and IT/BT, secondary flow (characterized by I 2 ) was generally stronger compared to axial streaming (compare I 1 and I 2 values in Figures 8A,B).
Integral parameter I 1 at planes 2-4 are shown in Figure 9. Unlike plane 1 where axial flow streaming (I 1 ) during inhalation was stronger than I 1 during exhalation (Figure 8), I 1 variation  right asymmetry of I 1 that was noted at planes 3 and 4 ( Figure 9) is also sustained in the higher generation G2 (planes 7-10). IT/BT for Wo=2.37 and 7.51. Figure 12 shows phase-variation of integral parameter I 2 at planes 5-10. In general, I 2 variation along planes 5-10 mirrored those of I 1 (Figure 10). Similar to I 1 , I 2 increased from the start of inhalation (or start of exhalation) to reach peak value at mid-inhalation (or mid-exhalation) and decreased until the end of inhalation (or end of exhalation). Also similar to I 1 (Figure 10 and IT/BT=50% (bottom row). Refer to Figure 9 for marker definitions. throughout the cycle, and I 2 was lowest at plane 9. This exact same pattern was also observed with respect to I 1 variation across planes on the right side of the second bifurcation ( Figure 10). On the left side of the second bifurcation, I 2 followed the same trend as I 1 throughout the cycle and for all Wo and IT/BT conditions, i.e., I 2, plane 10 > I 2, plane 8 > I 2, plane 6 . For IT/BT =25%, I 2 at the end of inhalation increased with increasing Wo, and this same trend was also observed in I 1 (Figure 10).
Overall, I 2 was greater than I 1 across planes 2-10 for all IT/BT and Wo, suggesting that secondary flows are stronger transport mechanisms as compared to axial streaming.

E. Turbulence characteristics
Turbulence levels in the airway flow at varying Wo and IT/BT were examined along the coronal plane using time-averaged contours of turbulence kinetic energy (k) calculated using equation (9) and dimensionless turbulence intensity (I) calculated using equation (10). Figure 13   Similar to the observations of Banko et al. 16,19 where Wo=7 and IT/BT=50%, a single-sided streamwise vortex (see plane 1B in Figures 2, 6) was observed during peak inhalation in the trachea for Wo=2.37 across all IT/BT ratios. Banko et al. 19 noted that the vortex at peak inhalation was not previously observed in studies using idealized symmetric geometries. Compared to idealized airway models, our airway geometry includes the mouth-to-glottis section in the coronal plane instead of the sagittal plane. Despite the left-right symmetry in our model starting from the glottis, the presence of a mouth-to-glottis section on the anatomical right side promoted singlesided streamwise swirls at lower Wo. The effect of the glottis can also be seen during inhalation, with higher I 1 and I 2 values at plane 1 as compared to planes 2-10. Also, Banko et al. 19 found weak secondary flow at peak exhalation followed by uniformly distributed axial flow at a location near the glottis. Our study showed the presence of axial and secondary flows in both Wo=2.37 and Wo=23.7 during exhalation.
At a given cross-section, I 1 and I 2 for Wo=7.51 and IT/BT=50% showed no variation at peak inhalation and peak exhalation across all generations. This was in disagreement with the Wo=7 findings of Banko et al. 19 , where I 1 was larger at peak inhalation and I 2 was larger at peak ex-halation. Our idealized symmetric airway geometry is likely the reason for this disagreement as compared to the asymmetric subject-specific airway of Banko et al. 19 . However, our I 1 and I 2 trends agreed with those reported by Banko et al. 19 for higher generations, such as I 1 and I 2 increasing and decreasing in a manner identical to the prescribed inflow velocity profile throughout a breathing cycle. Interestingly, I 1 and I 2 values showed noticeable variation between planes at all the generations (e.g., I 1, plane 5 = I 1, plane 7 > I 1, plane 9 ; I 1, plane 10 > I 1, plane 8 > I 1, plane 6 ). This suggests that local flow features can be responsible for lagging/leading axial flow streaming and lateral dispersion relative to the inflow profile, depending on the bronchial path and the generation number. This observation was also previously noted by Banko et al. 19 .
At the upper tracheal section (plane 1), both I 1 and I 2 were larger during peak inhalation as compared to their values in peak exhalation. In contrast, I 1 and I 2 individually reached their maximum value at both peak inhalation and peak exhalation for other cross-sections (planes 2-10).
The geometric asymmetry of the upper airway thus appears to have contributed to the following at plane 1: differences observed in a particular integral parameter (I 1 or I 2 ) between inhalation and exhalation; as well as differences between I 1 and I 2 . With increase in Wo at a given IT/BT ratio, I 1 and I 2 remained essentially the same for each plane at higher generations, showing no effect of Wo for given idealized geometry. Turbulence intensity levels in this study were pronounced in the mouth-to-glottis section, showing the effect of including realistic geometry even at higher RR. Our results point to the importance of including the upper airway structure and its orientation relative to the sagittal plane, which was also noted previously by Lin et al. 31 . Obstructions such as tongue and upper mouth geometry that were not included in our model can also enhance turbulence in the airway flow, as has been noted by Lin et al. 31 . Further studies are needed to isolate the roles of each of these anatomical structure in generating turbulence and asymmetric axial flow in the airway.
Standard treatment strategies of HFOV therapy uses mechanisms of bulk convection, molecular diffusion, and transport 32,33 , which are potentially affected due to turbulence levels. Our results showed the change in turbulence levels with the change in Wo and IT/BT ratios, suggesting for need more studies to examine the effect of each individual mechanism in HFOV. Although HFOV was proposed to overcome the acute lung injury (ALI) and adult respiratory distress syndrome, similar to airway turbulence induced by laryngeal jet 31   where Re=4,200 based on U T and D. It can be seen that Re/Wo 2 does not equal 2L/D when To include the effect of IT/BT =50%, we examined the use of a modified dimensionless stroke length, where the parameter (BT/2 IT) is included as a multiple of 2L/D (i.e., (2L/D)(BT/2 IT)).
The modified dimensionless stroke length matched Re/Wo 2 for all Wo and IT/BT conditions considered in this study (Table III), showing the importance of including IT/BT to accurately classify the operating flow regime. To examine the operating flow regimes at the cross-sectional locations where we analyzed the 3D CFD data (planes 1-10), we calculated local values of Reynolds number (Re L ) and Womersley number (Wo L ) as follows: where V L and D L denote axial velocity and in-plane airway diameter, respectively, at a given cross-sectional plane (Table S1 in supplementary material). Re L and Wo L are interrelated via the modified dimensionless stroke length as: Re L /Wo 2 L = (2L L /D L )(BT /2 IT ). The values of Re L and Wo L for a given test condition (i.e., Wo, IT/BT) were averaged over the different planes included within a particular generation (see Table II for plane locations).
The regime diagram for classifying the flow regime at different generations is shown in Figure 14, where Wo 2 L is plotted along the x-axis and the modified dimensionless stroke length ((2L L /D L )(BT/2 IT)) is plotted along the y-axis. For a given Wo and generation number, varying IT/BT did not noticeably alter flow regime location when using the modified dimensionless stroke length. For Wo=2.37, the trachea (G0) was in the turbulent zone of the regime diagram. This is in agreement with turbulence characteristics discussed previously (see Figure 13) where k and I values were found to be larger in the mouth-to-glottis section.
With increasing Wo beyond 2.37, flow through the trachea was in the unsteady-convective zone, meaning that both unsteady and convective acceleration terms are of importance in the Navier-Stokes' equations. With increase in generation number (i.e., moving further down the airway) for either Wo=2.37 or Wo=7.51, the flow regime tends towards the viscous-convective zone so that unsteady effects are not dominant (quasisteady). At Wo=23.7, all three generations are in the unsteady-convective zone, such that time-dependent effects on the flow field are the most dominant. Additionally, Wo=23.7 shows the farthest departure from the turbulent regime, which is in agreement with our k and I observations with increasing Wo (see Figure 13). Our flow regime results are in agreement with MRV studies on a subject-specific anatomical model at Wo=7 19 and on a planar double bifurcation model at Wo=6 and Wo=12 21 . At Wo relevant to HFOV, gas exchange and particle deposition at higher airway generations are expected to be most impacted by unsteady acceleration of the flow. By contrast, viscous forces and convective acceleration are expected to be more influential at lower Wo.

V. CONCLUSIONS
Using 3D CFD simulations, we examined the roles of varying respiratory rate and inhalation duration (IT/BT) on flow through an idealized human airway model consisting of a mouth-to- We developed a modified dimensionless stroke length including IT/BT to correct this discrepancy.
While lower Wo regime was dominated by viscous forces and convective acceleration, unsteady acceleration was dominant for higher Wo.
A central limitation of our study is the simplification of our model geometry. A range of morphological complexities observed from G0-G2, including (but not limited to) asymmetric branching and non-circular lumen of the airway branches, are expected to alter the respiratory flows.
The general trends reported here relative to Wo and IT/BT, specifically the importance of secondary flows with increasing Wo as well as the modified dimensionless stroke length accounting for IT/BT, are expected to be applicable in anatomically realistic airways. Additional model simplifications, such as the use of rigid walled vessels and a sinusoidal flow profile, should also be noted as limitations of this study. While changes to waveform shapes have shown to modestly impact the flow physics 5 , including airway wall motion has been reported to influence axial and secondary flows compared to rigid-walled airways 34 . Finally, RANS models (such as in this study) have well-known limitations in modeling unsteady turbulent flows, and comparisons of our study findings with higher-fidelity turbulence models (e.g., large eddy simulations) are needed to examine how the choice of turbulence model impacts time-varying flow physics.

SUPPLEMENTARY MATERIAL
See the supplementary material for: Figures S1-S4 of plane-normal velocity magnitude contours with superimposed in-plane velocity vector fields for inhalation at IT/BT=33% and exhalation at IT/BT=50% (Wo=2.37 and 23.7); Figure S5 showing time-averaged k and I contours for IT/BT=33%; Table S1 showing local Reynolds number (Re L ) and local Womersley number (Wo L ) for different planes across all the test conditions of Wo and IT/BT.