Effects of Varying Inhalation Duration and Respiratory Rate on Human Airway Flow

Abstract

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 the respiratory rate (RR) on airway flow remain unknown. Using three-dimensional numerical simulations of oscillatory flow through an idealized airway model (consisting of a mouth, glottis, trachea, and symmetric double bifurcation) at a trachea Reynolds number ($Re$) of 4200, 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 a Womersley number ($Wo$) of 2.41 to 1000 bpm ($Wo$ = 24.1) impacts airway flow characteristics. Irrespective of IT/BT, axial flow during inhalation at tracheal cross-sections was non-uniform for $Wo$ = 2.41, as compared to centrally concentrated distribution for $Wo$ = 24.1. 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 a stronger transport mechanism than axial flow streaming. Discrepancy in the oscillatory flow relation $Re$/$W{o}^2$ = 2 L/D (L = stroke length; D = trachea diameter) was observed for IT/BT ≠ 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$.

1. 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.41–2.95 based on the following relation:

$$Wo=\frac{D}{2}\sqrt{\frac{\omega }{\nu }},$$

(1)

where $\omega =2\pi /\mathrm{BT}$ is the angular frequency based on breathing time (BT, where BT = 60/RR seconds), D is the tracheal diameter, and $\nu$ is the kinematic viscosity of air. RR can be altered from the resting range of 10–15 bpm during exercise (such as walking RR = 20–60 bpm [4] and yoga RR = 3–7 bpm [5]), and also in mechanical ventilation strategies, such as high-frequency oscillatory ventilation (HFOV) [6,7,8]. Further, inhalation time (IT) is about 45% of BT under normal breathing conditions [6]. 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 [9,10]. Previous studies of inhaled drug delivery in human subjects with COPD and asthma have shown improved respiratory outcomes [11] and increased aerosol deposition with increased inhalation time [12]. However, computational fluid dynamics (CFD) studies examining particle deposition relevant to aerosols [13,14] and e-cigarettes [13,15,16] 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.

The respiratory system consists of airways with complex geometrical cross-sections, including the nasal cavity [17,18], mouth, pharynx, larynx, glottis [19], and the trachea that further bifurcates up to 16–23 Generations (G) [20] in addition to inter-subject geometry variability [18,21]. Due to the geometrical complexity, idealized geometries of the nasal cavity [22], mouth-to-glottis section [23], and airway generations [24,25] have been used in several studies. Weibel proposed an idealized airway geometry consisting of circular tubes at each generation [24]. This geometry bifurcated symmetrically along the sagittal plane and was able to capture key aspects of airflow and its features [26,27]. A review on using idealized airway models can be found elsewhere [21]. In addition to airway geometry, head rotation and posture can also influence airway fluid dynamics. Idealized airway model studies have included mouth-to-glottis geometry in the sagittal plane that incorporate head rotation [23], and clinical studies of HFOV have reported prone positions to improve gas exchange, as compared to the supine position [28]. Finally, voluntary exercises, such as freestyle swimming and some yoga breathing techniques, can be conducted with head rotation [29].

Previous studies examining the fluid dynamics of deep inspiration were conducted by prescribing steady inhalation flow through idealized [30] and anatomically accurate physical models [18,31]. Borojeni et al. [18] recently used computational fluid dynamics simulations to establish physiological ranges of nasal airflow variables, and documented inter-subject variability. Steady inhalation studies with varying airflow parameters (such as $Re$, $Wo$) identified both axial (streamwise) and secondary (transverse) flow dispersion to be effective transport mechanisms [30,31]. Unsteady physiological flow through subject-specific airways has also been experimentally investigated [32,33,34,35]. Jalal et al. [35] found that with respect to idealized airway models [30], 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 [35]. Jalal et al. [36] 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 the inhalation to exhalation phase [35]. Große et al. [33] 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 axial transport to the lung [33]. Oscillating flow studies showed that the size of secondary flow structures strongly depended on the instantaneous values of $Re$ and $Wo$ [33]. Two-dimensional PIV measurements of oscillating flow through an asymmetric idealized model (based on Weibel et al. [24] and Horsfield et al. [25]) at normal and HFOV conditions showed mass exchange at higher frequencies [32]. Experimental studies of HFOV flow through a subject-specific airway model [37,38] 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., IT/BT = 46% [6]), IT/BT as low as 20% has been used in mechanical ventilation of patients with chronic obstructive pulmonary disease (COPD) and asthma [4,39,40]. Further, breathing exercises that involve voluntary manipulation of the normal rhythm (e.g., hatha yoga) also result in changes to IT/BT. The effects of varying IT/BT on axial dispersion and secondary flows currently remains unknown.

In terms of the physical mechanisms driving flow and gas exchange in different portions of the human airway, Jan et al. [3] used order-of-magnitude analysis to examine the relative importance of various terms in axial flow and secondary momentum equations. Based on the regional domination of viscosity, unsteadiness, and convective acceleration, they developed a flow regime map by relating $W{o}^2$ with the product of the ratio of stroke length (L) to hydraulic diameter (D) of the trachea when using a sinusoidal breathing profile with IT/BT = 50%. This map was developed to: (a) identify what approximation can be used in modeling flow in different regions of the airway; and (b) investigate the fluid dynamic phenomena that affect gas exchange in HFOV with equal inhalation and exhalation durations. To extend the applicability of their map to realistic breathing patterns in normal and HFOV conditions, this regime chart needs to be modified to account for variations in IT/BT.

HFOV was proposed by Lunkenheimer et al. [41] to generate both active inspiration and expiration for eliminating entrained gas and gas decompression in airway models. The 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. [42] 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 studies on ARDS patients have shown improved gas exchange with HFOV without compromise in oxygen delivery [43]. Zhang et al. [6] performed computations using realistic breathing patterns under normal conditions ($Wo$ = 2.4) with IT/BT = 46% and HFOV ($Wo$ = 23.3) with IT/BT = 50%, and characterized fluid flow in the fourth-generation airway. However, the extent to which active inspiration and expiration at larger $Wo$ (characteristic to HFOV) affect the flow characteristics in upper airways remains fundamentally unclear.

In this study, we examine the effects of varying $Wo$ and IT/BT on airway flow characteristics. Three-dimensional (3D) computational fluid dynamics (CFD) simulations were performed by prescribing an oscillatory inflow profile in an idealized double bifurcation airway model (two generations) with a mouth-to-glottis section. The aims of this study were to: elucidate the flow physics under varying inhalation to breathing time, and modify the flow regime diagram of Jan et al. [3] to incorporate effects of varying IT/BT.

2. Methods

An idealized airway model [24] with the addition of an idealized mouth-to-glottis section (previously used in a particle deposition study [23]) 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) two generations (G) G1 and G2 of diameters 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 (ANSYS, Inc., Canonsburg, PA, USA) by importing an STL geometry file into ICEM-CFD. The surface was modified to create parts such as an inlet, walls, and outlets using the angle method, and a material region was created. Octree mesh was created using minimum and maximum element sizes, along with curvature refinement. Walls were selected to generate prism layers with growth rate = 1.1, and volume mesh was computed. A total of three different meshes, each consisting of tetrahedral cells with five prism layers on walls, were generated for mesh independence tests. $y_{+}$ for all the meshes were <5, corresponding to peak velocity condition in the trachea. Previous studies have reported that $y_{+}$ < 5 in the glottis region and $y_{+}$ < 1 in the trachea are sufficient for resolving flow in the near-wall region when using the k-$\omega$ SST turbulence model [44,45].

Three-dimensional CFD simulations of transient, incompressible flow through the above airway geometry were performed in ANSYS Fluent. k-$\omega$ SST turbulence model, used in previous studies of flow through the human airway [46,47], was chosen for all simulations in this study. Sinusoidal profiles were used as a simplified representation of realistic breathing waveforms [7]. The Reynolds number at the trachea during peak inhalation ($R{e}_{\mathrm{T}}$) was defined as

$$R{e}_{\mathrm{T}}=\left(\frac{U_{\mathrm{T}}D}{\nu }\right),$$

(2)

where $U_{\mathrm{T}}$ is the mean flow speed in the trachea, $\nu$ is the kinematic viscosity of air (1.46 × 10⁻⁵ m² s⁻¹), and D is the tracheal diameter = 18 mm. $R{e}_{\mathrm{T}}$ was maintained constant in all breathing waveforms at $R{e}_{\mathrm{T}}$ = 4200, which is in the physiological range [14,31]. The peak inhalation flow rate ($Q_{\mathrm{in},\mathrm{peak}}$) was determined using the relation:

$$Q_{\mathrm{in},\mathrm{peak}}=\frac{R{e}_{\mathrm{T}}\nu A_{\mathrm{T}}}{D},$$

(3)

where $\nu$ is the kinematic viscosity of air (1.46 × 10⁻⁵ m² s⁻¹) and $A_{\mathrm{T}}$ is the cross-sectional area along the transverse cross-sectional Plane 2 with diameter D = 18 mm. Inhalation and exhalation volumes were kept equal for a given RR and IT. The peak exhalation flow rate ($Q_{\mathrm{ex},\mathrm{peak}}$) was defined based on IT, BT, and $Q_{\mathrm{in},\mathrm{peak}}$ using the equation below:

$$Q_{\mathrm{ex},\mathrm{peak}}=Q_{\mathrm{in},\mathrm{peak}}\left(\frac{\mathrm{IT}}{\mathrm{BT}-\mathrm{IT}}\right),$$

(4)

such that $Q_{\mathrm{ex},\mathrm{peak}}$ = $Q_{\mathrm{in},\mathrm{peak}}$ when IT = 0.5 BT, and $Q_{\mathrm{ex},\mathrm{peak}}<Q_{\mathrm{in},\mathrm{peak}}$ when IT < 0.5 BT. Instantaneous profiles of inhalation velocity ($V_{\mathrm{in}}\left(t\right)$) and exhalation velocity ($V_{\mathrm{ex}}\left(t\right)$) were defined using the equations:

$$V_{\mathrm{in}}\left(t\right)=\frac{Q_{\mathrm{in},\mathrm{peak}}}{A_{\mathrm{in}}}sin\left(\frac{\pi t}{\mathrm{IT}}\right),0<t\le \mathrm{IT}$$

(5)

$$V_{\mathrm{ex}}\left(t\right)=\frac{Q_{\mathrm{ex},\mathrm{peak}}}{A_{\mathrm{in}}}sin\left(\frac{\pi (\mathrm{IT}-t)}{\mathrm{BT}-\mathrm{IT}}\right),\mathrm{IT}<t<\mathrm{BT},$$

(6)

where $A_{\mathrm{in}}$ is the area of the inlet. The above relations were used to develop nine oscillatory profiles of time-varying inflow velocity ($\mathbf{u}$), prescribed perpendicular to the inlet, and prescribed as initial conditions. The tests conducted here include three RRs (10 bpm, 100 bpm, and 1000 bpm (later being clinically impractical in adults)), each for IT/BT = 25%, 33%, and 50% [4,39,40,48] to examine fluid dynamics across several orders of magnitude variation in $Wo$. We expect that our findings can inform future studies of low Reynolds number ($Re$) airway flows (such as in children), where HFOV is clinically used at higher frequencies than in adults (e.g., RR = 960 bpm in neonates [49,50]). All the oscillatory velocity profiles consisted of positive velocities for inhalation phases ($\theta$ from 0° to 180°) and negative velocities for the exhalation phase ($\theta$ from 180° to 360°).

Boundary conditions included zero pressures at the outlets and no-slip walls. To characterize the oscillatory flow relative to viscous effects, $Wo$ was calculated for each condition using Equation (1), where $\nu$ = 1.46 × 10⁻⁵ m² s⁻¹ was used as the kinematic viscosity of air. $U_{\mathrm{T}}$ denotes mean flow speed in the trachea, and was calculated to be 3.4 m s⁻¹ by equating flow rates at the inlet and the trachea via mass conservation for an incompressible flow:

$$U_{\mathrm{T}}=\left(\frac{V_{\mathrm{in},\mathrm{peak}}}{A_{\mathrm{T}}}\right)A_{\mathrm{in}},$$

(7)

where $V_{\mathrm{in},\mathrm{peak}}$ is the peak inhalation velocity calculated from Equation (5).

After prescribing the velocity profile, each simulation was standard-initialized with initial x-velocity, y-velocity, and z-velocity to be equal to 0 m s⁻¹. Simulations were conducted at the High-Performance Computing Center at Oklahoma State University using 16 processor units. The solution method included a coupled scheme, and second-order spatial discretization. A uniform time-step of 10⁻³ 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 1. Mesh parameters used for mesh independence tests.

	Element Size [m]	Number of Cells	Simulation Time [min]
Mesh1	1 × 10−3	594,746	62
Mesh2	5 × 10−4	2,550,857	118
Mesh3	3 × 10−4	8,599,521	580

) at a time-step of 10${}^{-3}$ s for HFOV velocity profiles ($Re$ = 4200, $Wo$ = 24.1, IT/BT = 50%) for one breathing cycle. This particular case was used in mesh independence tests to decrease the solution time. Three-dimensional velocity ($\mathbf{u}$) was extracted from each converged simulation along a line 1-1’ at the upper trachea in the coronal plane (z = 0 m; see inset in Figure 1C). Figure 1C shows the extracted velocity magnitude along the line 1-1’ for Mesh 1 (dashed line), Mesh 2 (solid line), and Mesh 3 (dotted line with circle markers) as a function of non-dimensional diameter y/D, where D = 18 mm is the trachea 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 Meshes 2 and 3 were nearly identical (Figure 1C), Mesh 2 was selected for use in this study. In addition, the solution from each mesh study was analyzed further by measuring bulk velocity at a lower trachea (Plane 2). Average velocities of 3.4256 m s⁻¹, 3.4281 m s⁻¹ and 3.4298 m s⁻¹ was found, respectively, in mesh independence test simulations. Since the error in average velocity from Mesh 3 to Mesh 2 (0.05%) was less than the error between Mesh 3 and Mesh 1 (0.1%), Mesh 2 was selected for all simulations and time-step independence studies.

Time-step validation for higher $Wo$ was conducted on the finalized Mesh 2. For time-step size validation, three different time-step sizes (0.0005 s, 0.001 s, and 0.00375 s) were used with inlet velocity profile corresponding to $Re$ = 4200 and $Wo$ = 24.1 at IT/BT = 50% with the same computational setup as described earlier. Velocity profiles were extracted from the completed simulations along the coronal plane in the upper trachea (Plane 1 in

Table 2. Location of Planes used for data analysis and corresponding sectional airway diameters.

Plane	Description	$(\mathit{x},\mathit{y},\mathit{z})$ [mm]	$({\mathit{\theta }}_{\mathit{x}},{\mathit{\theta }}_{\mathit{y}},{\mathit{\theta }}_{\mathit{z}})$ [°]	Diameter [mm]
Plane 1	Upper trachea (G0)	(−78, 0, 0)	(90, 0, 0)	18
Plane 2	Lower trachea (G0)	(−170, 0, 0)	(90, 0, 0)	18
Plane 3	G1	(−200, 65, 0)	(55, 35, 0)	12.2
Plane 4	G1	(−200, 82, 0)	(−55, 35, 0)	12.2
Plane 5	G1	(−238, 35, 0)	(55, 35, 0)	12.2
Plane 6	G1	(−238, 108, 0)	(−55, 35, 0)	12.2
Plane 7	G2	(−244, 24, 0)	(20, 70, 0)	8.3
Plane 8	G2	(−244, 120, 0)	(−20, 70, 0)	8.3
Plane 9	G2	(−250, 32, 0)	(90, 0, 0)	8.3
Plane 10	G2	(−250, 110, 0)	(90, 0, 0)	8.3

). The three velocity profiles showed no major variation at the center of the airway, with minor differences in near-wall velocity (see Figure S1 in Supplementary Materials). As the error in velocity prediction was better in $\delta t$ = 0.001 s when comparing $\delta t$ = 0.00375 s and simulation time for $\delta t$ = 0.001 s was lower compared to $\delta t$ = 0.0005 s, the time-step of $\delta t$ = 0.001 s was used for all the simulations reported in this study.

Steady inhalation with similar boundary conditions, as described for transient simulations, were simulated at inlet flow rates of 10, 15, 30, 45, and 60 L min⁻¹. Pressure drop ($\mathsf{\Delta }p$) was estimated for both validation (steady flow simulations) and varying $Wo$ and IT/BT (transient flow simulations) in ANSYS Fluent.

Solution files from ANSYS Fluent were imported to ANSYS CFD-POST for further analysis and generating results. All $Wo$-tested, velocity magnitude contours were extracted along the planes in Table 2 and visualized in Tecplot 360 (Tecplot Inc., Bellevue, WA, USA). The Q-criterion previously used in a study of steady airflow through a similar geometry [46] was used to identify flow vortex structures at various breathing conditions. Q = 100,000 s⁻² was used for isosurfaces at peak inhalation, and Q = 15,000 s⁻² was used for isosurfaces at peak exhalation.

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 [30,34]. Axial (i.e., streamwise direction) flow uniformity at a cross-sectional plane was quantified using the integral parameter $I_1$ [30,34]:

$$I_1={\left(\frac{{\int \int }_A{(\mathbf{u}·\widehat{\mathbf{n}})}^{2}dA}{U_{\mathrm{T}}^{2}A}\right)}^{1/2},$$

(8)

where $\mathbf{u}$ is the 3D velocity field for a specific phase-point in the cycle, $\widehat{\mathbf{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$ [30,34]:

$$I_2={\left(\frac{{\int \int }_A{\parallel \mathbf{u}-(\mathbf{u}·\widehat{\mathbf{n}})\widehat{\mathbf{n}}\parallel }^{2}dA}{U_{\mathrm{T}}^{2}A}\right)}^{1/2}.$$

(9)

$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 2).

3. Results

3.1. Model Validation

We conducted steady flow simulations to compare with published results from studies of idealized, as well as subject-specific airway models [51,52,53]. A pressure drop in the airway as a function of the inhalation flow rate is shown in Figure 2A, comparing the results of our simulations with those of previous CFD studies that considered a whole airway model [51] and upper airway truncated models (only including the trachea and higher generations) [52,53]. From the steady flow simulations of our idealized airway model, a pressure drop was observed to increase with an increasing inhalation flow rate. The trend of pressure drop variation with inhalation flow rate was comparable to those reported in previous studies. When compared to previous studies, our results showed a higher pressure drop at larger inhalation flow rates. Variation in pressure drop values between the studies are most likely caused by the use of different airway geometries. Pressure drop ($\mathsf{\Delta }p$) in airways at peak inhalation ($\theta$ = 90°) and peak exhalation ($\theta$ = 270°) are shown in Figure 2B for different $Wo$ and IT/BT conditions. For a given $Wo$, the pressure drop remained almost the same in all IT/BT ratios at peak inhalation, while the pressure drop increased with increasing IT/BT at peak exhalation. In addition, with increase in $Wo$, a pressure drop at peak exhalation remained constant for the given IT/BT ratio. For a given IT/BT ratio, with an increase in $Wo$, the pressure drop at peak inhalation remained constant from $Wo$ = 2.61 to 7.61 and later increased at $Wo$ = 24.1.

3.2. Inhalation Flow Fields with Varying Breathing Frequency

Figure 3 shows the detailed view of the flow field at different planes along the trachea (G0) and G1 specified in Table 2. Flow at each planar location is visualized via contours of plane-normal (i.e., axial) velocity magnitude (non-dimensionalized by $U_{\mathrm{T}}$) overlaid with velocity vectors of secondary (i.e., transverse direction) flow at various instants of inhalation (A = early inhalation; B = peak inhalation; C = late inhalation). All the contours are shown from the top view, such that Planes 1 and 2 are parallel to Plane y–z, and Planes 3 and 4 are rotated locally to show the perpendicular views.

The presence of a jet in the trachea can be seen through the coronal, increasing in strength from the start of inhalation until peak inhalation ($\theta$ = 90°) and later decreasing until the end of inhalation ($\theta$ = 180°). A strong jet core was observed throughout inhalation (refer Movie 1). Additionally, the jet core was concentrated towards the ventral side of the trachea at Phase A, moved towards the dorsal side at Phase B, and returned towards the ventral side at Phase C. Compared to Plane 1, axial flow magnitude was lower at Planes 2 to 4. Similar to Plane 1, axial flow at Planes 2 to 4 increased in strength until peak inhalation and later decreased in strength. Irrespective of the specific inhalation phase, secondary flow was strongest at Plane 1. At Plane 4, secondary flow vectors opposite to that of Plane 3 were observed during all phases (A to C). Stronger axial velocity at Plane 3 compared to Plane 4 (both at peak inhalation) shows axial flow is asymmetrically distributed after the first bifurcation.

Non-dimensional axial velocity contours and secondary flow fields are shown in Figure 4 for $Wo$ = 24.1 at IT/BT = 25% at various instants of inhalation for Planes 1–4. Flow in Plane 1 showed uniform axial flow at the beginning of inhalation (Phase A), and centrally concentrated axial flow at peak inhalation (Phase B) and late inhalation (Phase C). The secondary flow at Plane 1 moved toward the left of the trachea from peak inhalation to late inhalation. Additionally, at Plane 1, axial flow in the coronal plane increased in strength from Phase A to Phase B, followed by a decrease in strength until Phase C. Axial flow at Plane 2 was uniform throughout inhalation with increased secondary flow velocity vectors from early inhalation until peak inhalation. Secondary flow at Plane 2 was directed from the dorsal and ventral sides and directed outward towards the left and right side of the trachea across inhalation. Planes 3 and 4 showed similar magnitudes of axial flow velocities during phases A to C. Flow separation zones at Planes 3 and 4 can be observed due to airway bifurcation at late inhalation (Phase C).

At IT/BT = 25%, for both $Wo$ = 2.41 and $Wo$ = 24.1 conditions, strong axial flow (refer Figure 3 and Figure 4) and strong secondary flow were observed in the trachea at Plane 1. Increase in $Wo$ weakened axial flow across Planes 1 and 2 during early inhalation (Phase A), while at peak inhalation, Plane 1 showed similar velocity magnitude even with an increase in $Wo$. During peak inhalation (Phase B) and late inhalation (Phase C) at Plane 1, axial flow was concentrated near the center of the airway at $Wo$ = 24.1, as opposed to asymmetric distribution at $Wo$ = 2.41. During peak inhalation (Phase B) at Planes 3 and 4, axial flow was stronger near the inner wall of the bifurcation for $Wo$ = 2.41, as compared to uniformly distributed axial flow for $Wo$ = 24.1 (discussed later under velocity profiles). Across the inhalation phases considered for analysis, non-uniform axial flow distribution was observed along Planes 1–4 for $Wo$ = 2.41. In contrast, axial flow appears to be uniformly distributed at Planes 1–4 during Phase A; centrally concentrated at Plane 1 during Phase B; and uniformly distributed for Planes 2–4 across phases B and C for $Wo$ = 24.1. For $Wo$ = 2.41 and IT/BT = 25%, axial flow at peak inhalation at Plane 3 was higher compared to Planes 2 and 4, while in $Wo$ = 24.1, flow through Planes 3 and 4 was equal, and was higher compared to Plane 2 (discussed later under integral parameters).

3.3. Exhalation Flow Fields with Varying $Wo$

Figure 5 shows non-dimensional axial velocity magnitude contours overlaid with secondary flow fields during exhalation at Planes 1–4 for $Wo$ = 2.41 at IT/BT = 25%. The coronal planes showed strong axial flow in mouth-to-glottis section with maximum axial velocity magnitude at Phase B. Uniform axial flow and non-uniform secondary flow can be observed across all planes during exhalation time (ET). Weak secondary flow was observed during exhalation at Plane 1 for $Wo$ = 2.41 at IT/BT = 25%. Multiple counter-rotating zones of secondary flow were observed near the wall at Plane 2 during all three phases of exhalation (D, E, F). Exhalation flow at Planes 3 and 4, located in a higher generation (G1) of the airway, show similar secondary flow fields from Phase D ($\theta$ = 225°, 25% ET) until Phase F for $Wo$ = 2.41 (IT/BT = 25%).

Figure 6 shows non-dimensional axial velocity magnitude contours overlaid with secondary flow fields during exhalation at Planes 1–4 for $Wo$ = 24.1 at IT/BT = 25%. Compared to exhalation flow in the mouth-to-glottis section for $Wo$ = 2.41 at IT/BT = 25%, weaker axial flow can be observed for $Wo$ = 24.1 at the same IT/BT (refer Movie 1,3). At Plane 1 of $Wo$ = 24.1 at IT/BT = 25%, the region of strong axial flow (near anatomical right) at Phase D decreased in size with increasing time until Phase F. In addition, two rotational secondary flow fields were observed at early exhalation (Phase D) until peak exhalation (Phase E), while at Plane 2, secondary flow was initially weak at Phase D, increased at Phase E, and later decreased until late exhalation (Phase F). As compared to Planes 1 and 2, secondary flow was weakened at Planes 3–4 throughout exhalation for $Wo$ = 24.1 at IT/BT = 25%.

3.4. Inhalation and Exhalation Flow Fields with Varying IT/BT

Figure 7 shows the detailed view of the flow field at Planes 1–2 (trachea, G0) and Planes 3–4 (G1) for $Wo$ = 2.41 and IT/BT = 50%. The coronal planes showed an increase in axial velocity from Phase A until Phase B, and later decreased at Phase C. Compared to Planes 2–4, Plane 1 showed stronger axial velocity magnitude and secondary flow at each inhalation phase. Further, the strong axial flow region was located near the right-dorsal side of Plane 1 during inhalation. At Plane 2, axial velocity magnitude increased until peak inhalation (Phase B) and later decreased at Phase C. The strong axial flow region at Plane 2 was located near the ventral side (bottom of contour) for phases A and C and uniformly distributed during peak inhalation (Phase B).

During peak inhalation (Phase B) at the coronal plane, axial flow was stronger in the anatomical right side of the airway compared to the left side of the airway. At Plane 3, with an increase in inhalation time, axial velocity magnitude increased until peak inhalation (Phase B) and later decreased at Phase C. Strong axial flow was observed at Plane 3 during Phase B as compared to Phase A, near the inner walls of the bifurcation (anatomical left side). At Plane 4, uniformly distributed flow with small flow separation near the outer airway walls at Phase A, a uniform axial flow at Phase B and axial flow similar to Phase A was observed in Phase C. At Planes 3 and 4, nearly similar secondary flow patterns can be seen throughout inhalation.

Figure S2 (Supplementary Materials) shows the detailed view of the flow field for $Wo$ = 2.41 at IT/BT = 33% during inhalation at Planes 1–2 (trachea, G0) and Planes 3–4 (G1). The magnitude of axial velocity at each plane increased until peak inhalation (Phase B) and later decreased at Phase C. At peak inhalation, secondary flow with a rotational flow pattern was observed near the wall at Plane 1 (phases A, B, and C), Plane 2 (Phase A), and Plane 3 (phases A, B, and C). Secondary flow was not observed during peak inhalation at Plane 4. coronal planes showed unequal flow distribution at first bifurcation during peak inhalation (Phase B), similar to those observed for $Wo$ = 2.41 at IT/BT = 25% and 50%.

Axial flow patterns were altered with increasing $Wo$ to 24.1 for IT/BT = 50% (Figure 8), as compared to those previously observed for $Wo$ = 2.41 for the same IT/BT. coronal planes showed increased axial flow from Phase A to Phase C during inhalation. Higher generations (G3) at peak inhalation had a stronger axial flow strength compared to G1 for $Wo$ = 24.1 and IT/BT = 50%. During inhalation, uniform axial flow can be observed at Planes 2–4 with increased secondary flow strength from G1 to G2. Across phases A to C, very weak secondary flow with uniformly distributed axial flow was observed during inhalation.

Figure S3 (Supplementary Materials) shows the detailed view of the flow field during inhalation at Planes 1–2 (trachea, G0) and Planes 3–4 (G1) for $Wo$ = 24.1 and IT/BT = 33%. At Plane 1, axial velocity magnitude increased until peak inhalation (Phase B) and later decreased at Phase C. Axial velocity was strongest at the center of Plane 1. Secondary flow patterns at Planes 2 to 4 for IT/BT = 33% were similar to those for IT/BT = 50% (at the same planes and $Wo$). When comparing IT/BT = 33% to IT/BT = 50%, only minor differences were observed between the cross-sectional distribution of axial flow at Planes 2 to 4. Additionally, a slight decrease in the jet can be noticed from the flow contours of the coronal plane.

With increasing IT/BT from 25% (Figure 3) to 50% (Figure 7) at the lower $Wo$ of 2.41, axial flow of similar magnitude was observed at Plane 1 at all phases, respectively. However, peak axial flow was concentrated towards the ventral side during early and late inhalation for 25% IT/BT compared to the dorsal side at the same phases in 33% (Figure S2 in the Supplementary Materials) and 50% IT/BT conditions. In addition, at lower $Wo$, all planes showed nearly similar velocity magnitudes, respectively at each phase, with a change in the IT/BT ratio. Secondary flow at lower $Wo$ from IT/BT = 25% to 50% showed increased velocity vectors during early inhalation and decreased secondary flow velocity vectors during peak until late inhalation. An increase in IT/BT slightly affected the axial flow in the upper trachea in the early inhalation (Phase A). Over all, the effect of IT/BT at lower $Wo$ showed a minimal effect on axial flow and secondary flow during the early inhalation phase.

With an increase in IT/BT at the higher $Wo$ of 24.1, axial flow magnitude in the trachea increased both axial flow (compare Figure 4, Figure S3 in the Supplementary Materials, Figure 8) and secondary at early inhalation. In addition, axial flow remained the same during peak and late inhalation with an increase in the IT/BT ratio, while the secondary flow decreased respectively at each plane. Contours from the coronal plane showed increased axial jet flow in the upper trachea with increasing IT/BT for both $Wo$, respectively.

Figure S4 shows non-dimensional axial velocity magnitude contours overlaid with secondary flow fields during exhalation at $Wo$ = 2.41 for IT/BT = 50%. Strong axial flow can be observed at Phase D until peak exhalation (Phase E) and later decreased at late exhalation (Phase F). Plane 1 showed nonuniform secondary flow throughout the exhalation. In addition, secondary flow increased in strength at peak exhalation (Phase E) compared to Phase D, and later decreased during 75% exhalation (Phase F). Multiple rotational flow patterns were observed near the wall of Plane 2 during all three phases of exhalation (D, E, F). Secondary flow patterns during the exhalation remained nearly same in Planes 3 and 4, respectively, with higher strength in Plane 4 compared to Plane 3.

Drastic differences in exhalation flow can be observed for $Wo$ = 23.7 with an increase in IT/BT = 25% to 50%. Axial flow increased with an increase in IT/BT ratio, as observed from all planes of Figure 6 and Figure S4 at each phase. Further, axial velocity increased in magnitude with increasing time from Phase D to Phase E and later decreased until Phase F in the mouth-to-glottis section. Secondary flow decreased at the early exhalation phase in the upper trachea with an increase in the IT/BT ratio. Additionally, the secondary flow strength increased in the lower trachea during Phase E (peak exhalation) and Phase F (late exhalation) with an increase in the IT/BT ratio.

To conclude, during exhalation at Planes 1 to 4 for both $Wo$, axial flow increased with increasing IT/BT from 25% (Figure 5 and Figure 6) to 50% (Figures S4 and S5 in the Supplementary Materials). This increase in axial flow was mainly due to increasing $Re$ with increasing IT/BT during exhalation. Further, increase in IT/BT increased the secondary flow strength in the lower portion of the trachea for each $Wo$, respectively.

3.5. Velocity Profiles

Three-dimensional velocity ($\mathbf{u}$) was extracted along the coronal plane during peak inhalation at Planes 1–10 for $Wo$ = 2.41 and 24.1 at IT/BT = 25% and 50%, and are presented in Figure 9, Figure 10 and Figure S6 (Supplementary Materials). Velocity profiles for $Wo$ = 2.41 at IT/BT = 25% and 50% showed narrow distribution, while a $Wo$ wider distribution of velocity was observed in the upper trachea (line 1) for $Wo$ = 24.1 at IT/BT = 25% and 50%. However, the flow in the trachea decreased at Plane 2 (line 2) compared to the upper trachea (line 1) for all conditions. Comparing across $Wo$, velocity magnitude at Plane 2 for $Wo$ = 2.41 was greater than that of $Wo$ = 24.1 for each IT/BT ratio. Velocity profiles showed variation in axial flow between Planes 3 and 4 at $Wo$ = 2.41 for both IT/BT. In comparison, axial flow from velocity profiles for $Wo$ = 24.1 at Planes 3 and 4 was uniform and of the same magnitude (Figure 10).

The right bifurcation (lines 5, 7 and 9 in Planes 5, 7 and 9, respectively) showed higher axial flow compared to the left bifurcation (lines 6, 8 and 10 in Planes 6, 8 and 10, respectively) at $Wo$ = 2.41 for IT/BT = 25% and 50% (refer to Figure S6 (Supplementary Materials)), while for $Wo$ = 24.1, nearly uniform axial flow was observed between left and right bifurcation for both IT/BT conditions.

3.6. Secondary Flow Vortices

Secondary flow vortices were identified using the Q criterion, and are shown in Figure 11 for $Wo$ = 2.41 and $Wo$ = 24.1 at IT/BT = 25% and 50%, during peak inhalation ($\theta$ = 90°) and peak exhalation ($\theta$ = 270°). Isosurfaces at Q = 100,000 s⁻² (peak inhalation) and Q = 15,000 s⁻² (peak exhalation) were used for visualization of secondary flows at Plane 1-1’ (at peak inhalation) and Plane 2-2’ (at peak exhalation). For $Wo$ = 2.41 at peak inhalation ($\theta$ = 90°), the secondary flow distribution in the trachea remained unaltered with an increase in IT/BT, while a small change in secondary flow pattern can be seen for $Wo$ = 24.1 with an increase in the IT/BT ratio. For $Wo$ = 2.41, isosurface contours at Plane 1 showed two regions with secondary flow vortices near the anatomical left-ventral side and a smaller region along the anatomical right-dorsal side. For $Wo$ = 24.1, secondary flow was observed along the walls at peak inhalation from Plane 1 for all IT/BT conditions.

At peak exhalation for $Wo$ = 2.41 at IT/BT = 25%, small secondary flow vortices were observed near the bifurcation of the parent generation, while for $Wo$ = 24.1 and IT/BT = 25%, secondary flow vortices were observed in G1 only. Increasing IT/BT for both $Wo$ increased the length of the secondary flow regions (axial dimension) vortex structures. Comparing across all conditions, $Wo$ = 2.41 at IT/BT = 50% showed a larger distribution of secondary flow regions (refer to Movie 1–3).

3.7. Integral Parameters

To quantify axial flow streaming and strength of secondary flow, integral parameters [34] were calculated using Equations (8) and (9) for all $Wo$ and IT/BT conditions at the planes in Table 2. Integral parameter $I_1$ at Planes 2–4 are shown in Figure 12 and Figure 13. Axial flow streaming ($I_1$) during inhalation was stronger than $I_1$ during exhalation (Figure 12) for IT/BT = 25%, while $I_1$ variation at Planes 2–4 was similar when comparing inhalation to exhalation for IT/BT = 50%. For $Wo$ = 2.41 and 7.61, axial flow streaming was asymmetric after the first bifurcation, and for $Wo$ = 24.1, axial flow streaming remained almost the same in both planes. At first bifurcation, within the G1, $I_1$ at the right bifurcation (Plane 3) was higher compared to the left bifurcation (Plane 4) in $Wo$ = 2.41 (both IT/BT) and 7.61 (IT/BT = 50%), while $I_1$ was higher in the left bifurcation for $Wo$ = 7.61 (IT/BT = 25%) and remained almost the same in $Wo$ = 24.1 in both IT/BT conditions. During inhalation or exhalation for a given $Wo$, $I_1$ increased from start to mid-phase (inhalation or exhalation) and later decreased until the end of the phase (inhalation or exhalation). With increase in IT/BT to 50%, a similar trend was observed for $I_1$. At a fixed cross-sectional location (Planes 2–4), changing either $Wo$ or IT/BT did not impact much of axial flow streaming ($I_1$) during inhalation. For each $Wo$ comparing across the IT/BT, the $I_1$ during exhalation was higher at IT/BT = 50% following the trend of an inlet velocity profile.

Figure 13 shows integral parameter $I_1$ at Planes 5–10. Across all conditions, strong axial flow streaming was observed throughout the breathing cycle. The phase-variation of $I_1$ during either inhalation or exhalation was identical to the corresponding phase-variation observed in Planes 2–4, that is, an increase from the start of inhalation (start of exhalation) to the peak value at mid-inhalation (or mid-exhalation) and decreasing until the end of inhalation (end of exhalation). For $Wo$ = 2.41 (both IT/BT) and $Wo$ = 7.61 (IT/BT = 50%), the right bifurcation (Planes 5, 7, 9) showed increased axial flow streaming compared to the left bifurcation (Planes 6, 8, 10), while for $Wo$ = 7.61, the right bifurcation (Planes 5, 7, 9) showed decreased axial flow streaming compared to the left bifurcation (Planes 6, 8, 10) at IT/BT= 25%, and $Wo$ = 24.1 (both IT/BT) showed equal axial flow streaming between the left (Planes 6, 8, 10) and right (Planes 5, 7, 9) bifurcation planes, respectively. For all the conditions, within the G2, $I_1$ at the right bifurcation showed to be higher or equal in Plane 9 compared to Plane 7, and at the left bifurcation, Plane 10 had higher or equal $I_1$ compared to Plane 8. The left–right asymmetry of $I_1$ that was noted at Planes 3 and 4 (Figure 12) was also sustained in the higher-generation G2 (Planes 7–10). Additionally, the $I_1$ during exhalation was higher at IT/BT = 50% in all conditions following the velocity profile.

Figure 14 shows phase-variation of the integral parameter $I_2$ at Planes 2–4, which is indicative of the lateral dispersion or the strength of the secondary flow. Across all conditions ($Wo$, IT/BT), $I_2$ variation throughout the cycle at Planes 2–4 followed the same order as that of $I_1$ (Figure 12). $I_2$ values were greater than $I_1$ values at any of these Planes, irrespective of $Wo$ and IT/BT, where $I_1$ and $I_2$ were compared. While axial flow streaming was present at Planes 2–4 for all $Wo$ and IT/BT examined here (Figure 12), secondary flows appear to be stronger. For a given IT/BT, increasing $Wo$ resulted in increasing $I_2$ at the end of inhalation ($\theta$ = 180°). Specific to $Wo$ = 23.7, the secondary flow was also observed at the end of exhalation ($\theta$ = 360°) for IT/BT = 50%. At a fixed cross-sectional location (Planes 2–4), individually changing either $Wo$ or IT/BT did not impact secondary flow streaming ($I_2$) during inhalation. For each $Wo$, comparing across the IT/BT, the $I_2$ during exhalation was higher at IT/BT = 50% following the velocity profile.

Figure 15 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 13). Similar to $I_1$, even $I_2$ increased from the start of inhalation (start of exhalation) to reach peak value at mid-inhalation (mid-exhalation) and decreased until the end of inhalation (end of exhalation). Additionally, similar to $I_1$ (Figure 13), for $Wo$ = 2.41 (both IT/BT) and $Wo$ = 7.61 (IT/BT = 50%), the right bifurcation (Planes 5, 7, 9) showed increased secondary flow streaming compared to the left bifurcation (Planes 6, 8, 10), while for $Wo$ = 7.61, the right bifurcation (Planes 5, 7, 9) showed decreased secondary flow streaming compared to the left bifurcation (Planes 6, 8, 10) at IT/BT = 25%, and $Wo$ = 24.1 (both IT/BT) showed equal secondary flow streaming between the left (Planes 6, 8, 10) and right (Planes 5, 7, 9) bifurcation planes, respectively. For all the conditions, within the G2, $I_2$ at the right bifurcation showed higher or equal in Plane 9 compared to Plane 7, and at the left bifurcation, Plane 10 had higher or equal $I_2$ compared to Plane 8. A similar pattern was also observed for $I_2$ variation across the planes for all $Wo$ and IT/BT conditions. 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 than axial streaming.

For a given condition at each phase, integral parameters over the complete breathing cycle showed an increase (both $I_1$ and $I_2$) with an increase in airway generation number. This increase was mainly due to incompressible flow simulations, such that flow velocity increased in higher generations with smaller diameters than the parent airway diameter. As a result, the value of integral parameters increased from G0 to G2, which was not observed in previous subject-specific model studies [31,34,35]. For a symmetric airway model with parent airway diameter (D) = 18 mm, the expected daughter airway diameter using Murray’s law (daughter airway diameter, ($d_o$) = $D/2^{1/3}$) is 14.2 mm, which was approximately the same parent airway to the daughter airway diameter relation in Banko et al. [31,34]. However, the idealized Weibel airway model used here has a daughter airway diameter less than those estimated using Murray’s law-based minimum daughter airway diameter (12.2 mm). Therefore, the flow velocity increased in daughter airways that in turn increased the integral parameters of higher generations. Despite the parent-to-daughter airway diameter ratios in the Weibel model deviating from ratios predicted using Murray’s law, we observed identical trends of increasing axial flow streaming and decreasing lateral dispersion for all generations as those reported by Banko et al. [34].

At a given cross-section, $I_1$ and $I_2$ for $Wo$ = 7.61 and IT/BT = 50% showed nearly no variation at peak inhalation and peak exhalation across all generations, which was in disagreement with the $Wo$ = 7 findings of Banko et al. [34] ( $I_1$ was more significant at peak inhalation and $I_2$ was larger at peak exhalation). Our idealized symmetric airway geometry is likely the reason for this disagreement compared to the asymmetric subject-specific airway of Banko et al. [34]. Current results for $I_1$ and $I_2$ trends agreed with those reported by Banko et al. [34] 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 the planes at all generations (e.g., $I_{1,\mathrm{G}0}$ < $I_{1,\mathrm{G}1}$ < $I_{1,\mathrm{G}2}$), suggesting that decreased airway diameters from the previous generation possibly increased the velocity of the flow, as well as the flow strengths. $I_1$ and $I_2$ showed variations between the left- and right-side airways of the same generation for curved mouth-to-glottis orientations, such 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. [34].

4. Discussion

Particle deposition and gas exchange within the human respiratory system are intricately dependent on flow through the airway. While previous studies have identified that axial flow streaming and lateral dispersion are effective transport mechanisms in the human airway [30,31,34,35,36], the effects of varying inhalation duration (IT/BT) in conjunction with the respiratory rate (RR) on airway flow characteristics remain unknown. The majority of airway fluid dynamics studies have either focused only on inhalation, or examined both inhalation and exhalation with IT/BT = 50%. Inter-subject variability of IT/BT < 50% can occur in voluntary and involuntary breathing exercises. Further, clinical therapies, such as HFOV, allow for patient-specific tuning of IT/BT and RR. Studies of varying IT/BT alongside RR can help identify fundamental differences in axial and lateral flow patterns, potentially valuable for improving airway flow distribution in HFOV therapy. We examined the role of varying IT/BT from 25% to 50% and RR from 10 bpm to 1000 bpm using 3D CFD simulations of airflow through an idealized double bifurcation airway model.

Axial and secondary flows were observed across all $Wo$ and IT/BT conditions examined in this study (refer Movie 1–3). Several differences can be observed when comparing flow fields at $Wo$ = 24.1 to those at lower $Wo$. At the upper tracheal section (Plane 1) for IT/BT = 25%, asymmetric axial velocity distribution was observed for $Wo$ = 2.41, as opposed to a centrally concentrated distribution of axial velocity for $Wo$ = 24.1 during inhalation. The lower tracheal section (Plane 2) for IT/BT = 25% uniform axial velocity distribution for all $Wo$. Asymmetric distribution of axial velocity was observed following the first and second bifurcations at lower $Wo$. Across all $Wo$ and IT/BT conditions, $Wo$ = 2.41 and 7.61 showed the strongest axial flow streaming ($I_1$) and lateral flow dispersion at Plane 9 and Plane 10 during inhalation, while during exhalation, all $Wo$ and IT/BT = 50% showed both strongest axial flow streaming ($I_1$) and secondary flow strength ($I_2$) at Planes 7–10 (G2). The latter observation regarding $I_2$ was in partial agreement with findings of a previous MRV study of flow through an idealized double bifurcation model [36], where secondary flows during inhalation and exhalation were reported to be strengthened with increasing $Wo$. Our results showed increased secondary flow strength at the end of inhalation, with an increase in $Wo$ along with equal secondary flow strengths for planes in a given generation ($I_{2,\mathrm{Plane}5}$ = $I_{2,\mathrm{Plane}6}$, $I_{2\mathrm{Plane}7}$ = $I_{2,\mathrm{Plane}8}$ = $I_{2,\mathrm{Plane}9}$ = $I_{2,\mathrm{Plane}10}$). In terms of secondary flow fields, we observed the generation of vortices at all $Wo$ during inhalation and exhalation. These secondary flow vortices can be expected to impact particle deposition in applications such as targeted drug delivery.

Similar to the observations of Banko et al. [31,34] where $Wo$ = 7 and IT/BT = 50%, a single-sided axial vortex (see Plane 1 at Phase B in Figure 3 and Figure 7 ) was observed during peak inhalation in the trachea for $Wo$ = 2.41 across all IT/BT ratios. Banko et al. [34] 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 single-sided axial swirls at lower $Wo$ (vortex structure from Figure 11). The effect of the glottis can also be seen during inhalation, through axial flow velocity contours in all $Wo$ and IT/BT in the upper trachea.

Choi et al. [20] used large eddy simulations to examine the effect of upper airway truncation on the flow in the trachea, and found that the airways above the glottis were crucial for generating turbulence in the trachea. Our results showing the formation of vortex structures 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. [17]. Obstructions such as tongue and upper mouth geometry that were not included in our model can also enhance turbulence in the airway flow, and has been noted by Lin et al. [17]. Further studies on subject-specific airways are needed to isolate how each of these anatomical structures impact fluid dynamics and examine the importance of $Wo$ and IT/BT in generating asymmetric axial flow in the airway.

Finally, for a given peak inhalation $R{e}_{\mathrm{T}}$, a change in IT/BT ratio for a given $Wo$ affected the tidal volume. With an increase in IT/BT ratio, tidal inhalation volume increased. The exhalation tidal volume was matched with tidal inhalation volume to maintain equal volumes in both inhalation and exhalation for involuntary breathing frequencies, such that lower IT/BT ratio breathing has smaller $Re$-based exhalation velocity profiles. The lower exhalation velocity profile affected the axial flow streaming in all generations, showing lower strength of primary and secondary flows. Using varying IT/BT ratio might be beneficial with peak inhalation and exhalation of equal $Re$, such that exhalation tidal volume would be more than an inhalation tidal volume, resulting in a greater amount of gases that would have been expelled out, which is indeed necessary for some situations, such as in ARDS [41].

Modified Flow Regime Diagram

Jan et al. [3] developed a diagram to classify the flow regime within different sections of the human airway, using an order of magnitude analysis starting from the Navier–Stokes’ equations. This diagram was used to characterize the influence of fluid viscosity, unsteadiness, and convective acceleration in an idealized airway bifurcation for $Wo$ ranging from 2.3 to 21.3. Their flow map examined $W{o}^2$ versus non-dimensional stroke length $2L/D$, where L is the stroke length (i.e., axial length of travel of a fluid particle) that can be calculated as the ratio of stroke volume (SV) swept in a cycle to the tracheal cross-sectional area) and D is the tracheal diameter. For a sinusoidal waveform (as considered in this study), Jan et al. [3] showed that $Re$ is related to dimensionless stroke length ($2L/D$) and $Wo$ as $Re/W{o}^2=2L/D$. However, their classification was restricted to oscillatory (sinusoidal) breathing cycle with IT/BT = 50%.

Table 3. Flow regime parameters evaluated at the upper tracheal section (Plane 1): $Re/W{o}^2$, dimensionless stroke length $2L/D$ and modified dimensionless stroke length $(2L/D)$(BT/2 IT). Stroke length (L) was calculated as: SV/$(\pi D^2/4)$, where SV = stroke volume and D = tracheal diameter = 0.018 m. SV was obtained by integrating the time-varying volumetric flow rate at the inlet ($Q_{\mathrm{in}}$) over a breathing cycle for each $Wo$ and IT/BT condition, where $Q_{\mathrm{in}}\left(t\right)=V_{\mathrm{inlet}}\left(t\right)A_{\mathrm{in}}$ (prescribed inlet velocity profile, $V_{\mathrm{inlet}}\left(t\right)$, refer to (5) and (6)). Reynolds number ($Re$) based on D and mean flow speed in trachea ($U_{\mathrm{T}}$) was calculated using an equation defined as: $Re$ = $\frac{U_{\mathrm{T}}D}{\nu }$, where $U_{\mathrm{T}}$ (=3.4 m s⁻¹). $Re$ was maintained constant at 4200 across all test conditions. The Womersley number ($Wo$) was calculated using Equation (1). Note that breathing time (in seconds), BT = (60/RR), where RR = respiratory rate (in bpm).

RR [bpm]	$\mathbf{Wo}$	IT/BT [%]	SV [m3]	L [m]	$2\mathit{L}/\mathit{D}$	$\mathbf{Re}/{\mathbf{Wo}}^2$	($2\mathit{L}/\mathit{D}$)(BT/2 IT)
10	2.41	50	1.66 × 10−3	6.51	724	723	724
10	2.41	33	1.09 × 10−3	4.3	478	723	724
10	2.41	25	0.83 × 10−3	3.26	362	723	724
100	7.61	50	1.66 × 10−4	0.65	72	72	72
100	7.61	33	1.09 × 10−4	0.43	48	72	72
100	7.61	25	0.83 × 10−4	0.33	36	72	72
1000	24.1	50	1.66 × 10−5	0.07	7	7	7
1000	24.1	33	1.09 × 10−5	0.04	5	7	7
1000	24.1	25	0.83 × 10−5	0.03	4	7	7

shows the SV, L, $2L/D$ evaluated for this study at the upper tracheal cross-section (Plane 1) across the various conditions of $Wo$ and IT/BT, where $Re$ = 4200, based on $U_{\mathrm{T}}$ and D. It can be seen that $Re/W{o}^2$ does not equal $2L/D$ when IT/BT ≠ 50%.

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/W{o}^2$ for all $Wo$ and IT/BT conditions considered in this study (Table 3), showing the importance of including IT/BT to classify the operating flow regime accurately. 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 the Reynolds number ($R{e}_{\mathrm{L}}$) and Womersley number ($W{o}_{\mathrm{L}}$) as follows:

$$R{e}_{\mathrm{L}}=\frac{V_{\mathrm{L}}{D}_{\mathrm{L}}}{\nu }$$

(10)

$$W{o}_{\mathrm{L}}=\frac{D_{\mathrm{L}}}{2}\sqrt{\left(\frac{2\pi }{\mathrm{BT}}\right)\left(\frac{1}{\nu }\right)},$$

(11)

where $V_{\mathrm{L}}$ and $D_{\mathrm{L}}$ denote the axial velocity and in-plane airway diameter, respectively, at a given cross-sectional plane (Table S1 in the Supplementary Materials). $R{e}_{\mathrm{L}}$ and $W{o}_{\mathrm{L}}$ are interrelated via the modified dimensionless stroke length as: $R{e}_{\mathrm{L}}/W{o}_{\mathrm{L}}^2=(2L_{\mathrm{L}}/D_{\mathrm{L}})(BT/2IT)$.

The values of $R{e}_{\mathrm{L}}$ and $W{o}_{\mathrm{L}}$ for a given test condition (i.e., $Wo$, IT/BT) were averaged over the different planes included within a particular generation (see Table 2 for plane locations). The regime diagram for classifying the flow regime at different generations is shown in Figure 16, where $W{o}_{\mathrm{L}}^2$ is plotted along the x-axis and the modified dimensionless stroke length ((2$L_{\mathrm{L}}$/$D_{\mathrm{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.41, the trachea (G0) was in the turbulent zone of the regime diagram.

With increasing $Wo$ beyond 2.41, flow through the trachea was in the unsteady-convective zone, meaning that both unsteady and convective acceleration terms are important in the Navier–Stokes’ equations. With increased generation number (i.e., moving further down the airway) for either $Wo$ = 2.41 or $Wo$ = 7.61, the flow regime tends towards the viscous-convective zone so that unsteady effects are not dominant (quasisteady). At $Wo$ = 24.1, all three generations are in the weak turbulence ($Re/Wo$≈ 100–200 [3]), such that unsteady effects on the flow field are the most dominant. Our flow regime results agree with MRV studies on a subject-specific anatomical model at $Wo$ = 7 [34] and on a planar double bifurcation model at $Wo$ = 6 and $Wo$ = 12 [36]. At $Wo$ relevant to HFOV, gas exchange and particle deposition at higher airway generations are expected to be most impacted by the unsteady acceleration of the flow. By contrast, viscous forces and convective acceleration are expected to be more influential at lower $Wo$.

Investigating the 3D fluid dynamics in anatomically accurate airway models and physiologically relevant breathing patterns can be computationally demanding. Inter-subject variability in geometry can present challenges [18] in identifying salient features of the fluid flow and in developing modified scaling relations as in this study. In this regard, studies of idealized airway models with sinusoidal breathing patterns can help in providing a basic understanding of flow physics at relatively low computational cost. The current study investigated the flow features inside the airways for an idealized airway geometry, varying clinically significant parameters (IT/BT and RR). The modified dimensionless stroke length identified in this work can be useful for fluid dynamic studies of clinically relevant situations where IT/BT variation (e.g., mechanical ventilation in COPD and asthma patients [4,39,40,48] and RR variation (e.g., HFOV in neonates) are observed. Future studies using subject-specific geometries are needed to examine the clinical significance of our findings.

5. 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-glottis section and two generations. Axial and secondary flows were observed throughout the model at all conditions of $Wo$ and IT/BT. For $Wo$ = 2.41, the strong non-axisymmetric axial flow was observed in the trachea during inhalation. For $Wo$ = 24.1, centrally concentrated axial flow was observed during inhalation in the upper trachea, followed by uniformly distributed weak axial flow in the lower trachea. On account of secondary flow, lateral dispersion was found to be the dominant transport mechanism across all test conditions. For all $Wo$, increase in IT/BT ratio increased the axial flow at early inhalation. In addition, secondary flow after peak inhalation increased during small $Wo$, and decreased during large $Wo$. For all $Wo$ increased IT/BT showed increased axial and secondary flows strength during each phase of exhalation, respectively. As changing IT/BT changes the stroke length (L) of an oscillatory flow, we observed a breakdown in the theoretically expected $Re$/$W{o}^2$ = 2 L/D relation for IT/BT ≠ 50%. We developed a modified dimensionless stroke length including IT/BT to correct this discrepancy. While a 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 velocity profiles have shown to modestly impact the flow physics [7], including airway wall motion has been reported to influence axial and secondary flows compared to rigid-walled airways [54]. Finally, RANS models (such as in this study) have well-known limitations in modeling unsteady turbulent flows and comparing 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.