Computational models are useful for understanding respiratory physiology. Crucial to such models are the boundary conditions specifying the flow conditions at truncated airway branches (terminal flow rates). However, most studies make assumptions about these values, which are difficult to obtain in vivo. We developed a computational fluid dynamics (CFD) model of airflows for steady expiration to investigate how terminal flows affect airflow patterns in respiratory airways. First, we measured in vitro airflow patterns in a physical airway model, using particle image velocimetry (PIV). The measured and computed airflow patterns agreed well, validating our CFD model. Next, we used the lobar flow fractions from a healthy or chronic obstructive pulmonary disease (COPD) subject as constraints to derive different terminal flow rates (i.e., three healthy and one COPD) and computed the corresponding airflow patterns in the same geometry. To assess airflow sensitivity to the boundary conditions, we used the correlation coefficient of the shape similarity (R) and the root-mean-square of the velocity magnitude difference (D_{rms}) between two velocity contours. Airflow patterns in the central airways were similar across healthy conditions (minimum R, 0.80) despite variations in terminal flow rates but markedly different for COPD (minimum R, 0.26; maximum D_{rms}, ten times that of healthy cases). In contrast, those in the upper airway were similar for all cases. Our findings quantify how variability in terminal and lobar flows contributes to airflow patterns in respiratory airways. They highlight the importance of using lobar flow fractions to examine physiologically relevant airflow characteristics.

## Introduction

Computational and experimental methods to examine respiratory airflows are vital tools for improving our understanding of human lung physiology. Computational fluid dynamics (CFD) models have been employed for various purposes, such as predicting particle deposition patterns for optimizing drug delivery [1–3] or for estimating health effects due to particulate matters [4–6] in idealized or patient-specific airways. In addition, CFD models have been useful in clinical studies for assessing airflow characteristics and therapeutic effects of medical treatments in specific diseases, including asthma [7,8], sleep apnea [9–11], and chronic obstructive pulmonary disease (COPD) [12,13]. Similarly, experimental studies have been used to (1) measure in vitro airflows in physical models of respiratory airways [14–16], (2) provide a means to test the validity of the results obtained from CFD simulations [17–20], and (3) show the dependency of airflow on geometry and flow asymmetry [21,22].

One of the key steps in both developing CFD models and conducting in vitro airflow measurements is the determination of boundary conditions, which specify either the flow rates or pressure values at individual inlets and outlets of the airway geometry. However, in vivo flows cannot be measured, especially for the airway branches within the lungs, and physiologically relevant boundary conditions are difficult to determine. Furthermore, because limitations in imaging resolution preclude the construction of three-dimensional (3D) geometries for all lung branches, the flows in airway branches truncated beyond the terminal branches in the 3D geometry need to be approximated. Multiscale models have previously been developed to estimate such flows. For example, impedance models estimate the impedances of airway branches beyond the truncation and assign those values as boundary conditions to the terminal branches in computing flows in the 3D geometry [23,24]. This approach also allows us to incorporate lung mechanics in small airways to compute the flow [20,25,26]. In addition, advanced custom codes for analyzing high-resolution computed tomography (HRCT) images have been developed to represent truncated airway branches more realistically [27,28], and such codes have made it possible to predict flows in the terminal branches of 3D airway geometries more accurately in patient-specific conditions [29,30]. However, such capabilities are not readily available in widely used commercial software, and several computational studies have made simplified assumptions regarding the flow conditions in the terminal branches (e.g., equal flow rates or pressure values for each) [31–34]. These assumptions, however, may not accurately represent actual lung conditions. In patient-specific cases, they can lead to errors in predicting airflows and related quantities, such as particle dosimetry [20,27].

Uncertainty regarding the flows in terminal branches also arises when conducting in vitro experiments to measure airflow. Because airways have complex structures consisting of multiple tree branches that successively bifurcate into increasingly smaller branches, perfect control of the system may not be possible during experiments, and perturbations in the system may cause uncertainties in measuring flows in the terminal branches. Such uncertainties in terminal flows are likely to influence the airflows in the airway branches of interest. Furthermore, the use of inaccurate terminal pressures or flow rates as the boundary conditions to validate computational models may lead to discrepancies between the computed and measured airflow patterns.

Given these issues, assessing the extent to which variations in terminal flow conditions affect the predicted or measured airflows is crucial in both computational and experimental studies of physiologically relevant characteristics of airflow in respiratory airways. Here, we developed computational models of an idealized airway geometry to investigate the sensitivity of airflow patterns to variation in terminal flow rates. For the same geometry, we considered four sets of terminal flow conditions: three representatives of a healthy subject and one of a subject diagnosed with COPD. To obtain physiologically relevant boundary conditions, as proposed by De Backer et al. [35], we derived the terminal flow rates by using HRCT lung images of a healthy subject and of a subject with COPD to calculate the flow fractions in the lung lobes. We conducted the study in two steps: First, we used one set of terminal flows derived from the healthy subject and examined the detailed characteristics of airflow under a steady expiratory breathing condition. To validate the computational results, we then conducted particle image velocimetry (PIV) measurements and compared them with the computed velocity contours. We then compared the observed airflow characteristics for the healthy subject with those computed for other sets of terminal flows to assess the sensitivity of airflow patterns to the boundary conditions.

## Methods

### Airway Geometry.

Figure 1 shows the idealized airway geometry, which included both upper respiratory tracts and central airways bifurcating up to the sixth generation. We used airway dimensions obtained from casting, scanning, and 3D reconstruction [36] for the upper airway, and the Horsfield model [37] for central airways, to create a 3D geometry by using computer-aided design software (autocad2012, The Autodesk, San Rafael, CA). The geometry had 17 terminals in the central airway and one outlet at the oral cavity to represent the mouth (Fig. 1). We used the same airway geometry for both in vitro and in silico studies. For the in vitro study, we used a 3D printer (Protomold, Maple Plain, MN) to manufacture the airway model.

### Experimental Setup: Flow System.

We designed an experimental system to obtain flow-field measurements in the hollow cast of a central airway model representative of a healthy individual. Figure 1(a) schematically shows the flow loop and its major components: the physical model of the airway, pumping system, filtration system, and temperature-control system. The apparatus was equipped to pump bulk fluid (see below for composition) at high pressures and to measure flow rates, pressures, and temperatures in situ. To circulate the flow, we used two pumps: an extreme-head turbine pump (McMaster-Carr, Elmhurst, IL; pump 1 in Fig. 1(a)) and a viking gear pump (Siewert, Rochester, NY; pump 2 in Fig. 1(a)). Pumps 1 and 2, respectively, pushed and pulled the flow through the airway model and reservoir. We controlled pump 1 by a VS1MX Microdrive (Baldor, Fort Smith, AR; driver 1 in Fig. 1(a)) and pump 2 by a VLT 5000 Aqua drive (Danfoss, Baltimore, MD; driver 2 in Fig. 1(a)).

To provide the boundary conditions for CFD simulations, we monitored the flow rates at each of the 17 terminals of the central airway. We used a clamp-on flow transducer probe (BioPro TT brand, model BCT 3/16 × 1/16 A, em-tec GmbH, Finning, Germany), which had an operational range of 0–6 l/min and an accuracy of ±3%. To generate the output signal, we connected the probe to a signal transducer linked to a flow meter board (OEM, em-tec, Finning, Germany) and an oscilloscope (54615B, Agilent Technologies, Santa Clara, CA) to average the signal and give a final readout of the probe.

A temperature-control system properly maintained the Reynolds number by preventing the fluid from overheating and causing changes in bulk fluid viscosity as it was pumped through the system. The physical loop system, especially the airway model, could also react poorly to excess heat. Therefore, we monitored the fluid temperature by using a K-type thermocouple submerged in the fluid at the base of, and immediately below, the fluid reservoir. The thermocouple (Omega, Norwalk, CT) was connected to a multimeter (HHM290, Omega, Norwalk, CT) outside of the reservoir, which displayed the temperature readout with a resolution of 0.1 °C. To cool the fluid, we used a water-cooled heat exchanger that contained a copper tube wrapped around the manifold through which the fluid flowed. During the experiment, we continually monitored and adjusted the cooling loop.

### Experimental Setup: Particle Image Velocimetry.

Figure 1(b) shows the imaging setup, which consisted of a laser sheet and a camera to measure flow fields within the physical airway model. We measured the flow fields by PIV, whereby we seeded a fluid with hollow glass microspheres (Sphericel, item # 110P8, Potters Industries LLC, Malvern, PA) and used a camera system to capture images of the seed particles as they tracked the flow. The fluid was composed of water, glycerin, and sodium iodide (weight ratio, 30.6:20.4:49) to match the refractive index between the fluid and model wall and reduce distortion in the captured images due to light refraction. The fluid mixture had a viscosity of 5.5 cP and a specific gravity of 1.58. The specific gravity of the particles (1.10) resulted in a calculated settling velocity of 5 *μ*m/s [38], which is orders of magnitude lower than our measured velocities.

We used a light sheet (thickness, <1 mm), generated by a laser beam (532-nm Nd-YAG Class 4 laser, New Wave Research, Fremont, CA) and positioned such that it is parallel to the coronal plane at the center of the airway geometry to illuminate the seed particles running along the fluid flow. The camera (1.3 Megapixel 8 Hz CCD PIV Camera, TSI Inc., Shoreview, MN), which was positioned so that its line of sight was 90 deg in relation to the light sheet, captured the light scattered by the seed particles with a small time differential (8–200 *μ*s). We used a laser pulse synchronizer (model 610034, TSI Inc.) to synchronize the camera and laser and the software Insight 3G (TSI Inc.) to collect the images simultaneously. We performed cross-correlation analysis on these images and determined the flow-field velocity over a two-dimensional (2D) plane.

To process the images for PIV analysis, we used a program that employed the window deformation and iterative multigrid method, which utilized an iterative cross-correlation algorithm to measure particle displacement [39,40]. The PIV experiment allowed us to measure the velocity components in the coronal plane (i.e., *x–y* plane) created by the laser sheet.

The instantaneous velocity measurement at a single measurement point contains an error that can be decomposed into a systematic error and a random error. As discussed earlier, the particle and fluid densities were similar, and the viscosity of the fluid was relatively high. Accordingly, the settling velocity of the particles, and hence the systematic error, was negligible. The uncertainty of PIV measurements based on digital images has been studied extensively [41,42]. Assuming adequate density of seed particles, uncertainty due to random errors is a function of the interrogation window size and the diameter of an individual particle's image. A particle image diameter of 2–3 pixels and a 32 × 32 pixel interrogation window (conditions for the current experiment) lead to a random error of ∼0.1 pixels, which corresponds to ∼5% of the maximum velocity within any PIV field of view [43].

As the number of samples *N* approaches infinity, the value of the ensemble average will converge to exactly that of the true mean of the velocity as long as the flow is steady in time. For all the measurements shown, 500 instantaneous vector fields were processed to determine the mean flow properties. Therefore, the resulting uncertainty of the mean velocity is less than 0.1 pixel, which corresponds to <3% of the reported maximum velocity in any reported region.

### Numerical Computation of Flows.

To simulate airflows in the airway geometry for different input conditions, we performed CFD simulations under steady exhalation conditions. We sought to investigate flow rates within the physiological range. Therefore, when matching the simulated and experimental data, we set the branch inflow conditions so that the Reynolds number at the trachea was 4000 (i.e., total flow rate of air = 770 ml/s), which corresponds to the mean breathing flow rate during moderate activity [44]. Because the flow rates we examined here were in the turbulent flow regime, we adopted the shear stress transport *k–ω* turbulence model for incompressible steady-state flow [45]. At all boundaries, we used a turbulent intensity value of 5% and a turbulent viscosity ratio of 10.

We solved the flow equations in the finite-volume-based solver fluent (version 15.0.7, ANSYS, Canonsburg, PA), using the SIMPLEC method for pressure–velocity coupling and a second-order scheme for the convective terms. To ensure that near-wall meshes were sufficiently resolved, we used ICEM CFD (version 14.0, ANSYS) to adopt a hybrid-meshing scheme that smoothly transitioned from fine to coarse grids away from the wall. To verify that the results did not depend on grid size, we varied the mesh density, using either the original (12,423,894 cells) or refined meshes (19,858,046 cells) and ensured that the relative differences between velocity magnitude values along the centerline across the trachea (50 mm above the bifurcation) were less than 2%. We used the refined mesh data to generate the results presented in this paper. We ran the simulations in parallel, using 160 2.8-GHz Intel Xeon quad-core Nehalem processors at the U.S. Department of Defense (DoD) Supercomputing Resource Center, located at the U.S. Army Research Laboratory in Adelphi, MD. The typical run time was ∼8 h for each simulation.

### Patient-Specific Lobar Flow Fractions.

To experimentally and computationally examine the flow patterns in the idealized airway geometry, we initially determined the boundary conditions (i.e., sets of flow rates in the terminal branches) for different lung conditions. Figure 2(a) schematically shows the steps by which we accomplished this here. First, we collected thin slice computed tomography (CT) images of the whole lung at two breath-holding conditions (full inspiration and full expiration) following a previously described protocol (SOMATOM Definition Flash scanner (Siemens, Erlangen, Germany) [27]. By segmenting the five lobes of the lung, we calculated the change in volume of each lobe from inspiration to expiration to obtain the flow distribution among the five lung lobes (hereinafter, lobar flow fractions). Subsequently, we specified the flow rates in the terminal branches by using the lobar flow fractions as constraints for both healthy and diseased (i.e., COPD) lung conditions.

The lobar flow fractions were calculated from two subjects, one healthy and one with COPD, at the University of Virginia (Charlottesville, VA). The healthy subject was an asymptomatic 72-year-old female with no history of smoking and normal results on chest X-ray and pulmonary function tests [forced expiratory volume in 1 s (FEV_{1}) = 102%, ratio of FEV_{1} to forced vital capacity (FEV_{1}/FVC) = 0.94]. The COPD subject was a 51-year-old female smoker (30 packs/year) with an FEV_{1} of 49% and FEV_{1}/FVC of 0.57 on the pulmonary function test. The protocols for subject recruitment and study procedures were approved by Institutional Review Boards at the University of Virginia and the U.S. Army Medical Research and Materiel Command Human Research Protection Office, Fort Detrick, MD. All subjects gave written informed consent.

*Q*, in lobe

_{k}*k*, by the volume differential between the images obtained at the two breath-hold conditions,

*ΔV*, as a percentage of the sum of the lung volume differentials, $\u2211k=15\Delta Vk$

_{k}The lobar flow fractions for the healthy subject and COPD subject thus obtained are shown in Fig. 2(b), denoted as healthy III and COPD, respectively. We then used these lobar flow fractions as constraints to define the corresponding boundary conditions in experimental (in vitro) and computational (in silico) studies, with the aim of validating our CFD model and assessing how changes in the boundary conditions would affect airflow patterns in the central and upper airways.

We also used two other boundary conditions in this study. In the in vitro work, we experimentally mimicked airflow patterns during steady expiration for a healthy lung by manually adjusting the flows to the physical model of the flow setup (Fig. 1(a)) to match the lobar flow fractions for the healthy subject. We set the flow fractions at individual branches to the asymptotic values obtained after the flow reached a steady-state with all control valves at each outlet kept open. Then, to characterize the details of the airflow patterns in the physical model (Fig. 1(b)), we used the same flow setup and conducted PIV measurements on two separate days (day 1 and day 2). The repeated measurements allowed us to sample a sufficient number of PIV images at multiple branches. The flow rates in the individual terminal branches varied across days (Table 1), even though the lobar flow fractions on both days were similar to those of healthy III. (The intraclass correlation coefficient among the three healthy conditions was 0.8, indicating excellent agreement [46].) Accordingly, we designated these two sets of terminal flow rates as separate boundary conditions for a healthy lung (denoted correspondingly as healthy I and II in Fig. 2).

Consequently, in the in silico study, we ran simulations to assess the sensitivity of airflow patterns to four sets of terminal flow rates, three derived from one healthy lung and the other from one diseased lung: healthy I, used on day 1 of the PIV experiment; healthy II, used on day 2 of the PIV experiment; healthy III, directly derived from the lobar flow rates of the healthy subject, assuming equal distribution of flows among the terminal branches within a lobe; and COPD, directly derived from the lobar flow rates of the COPD subject under the same assumption. Even distribution of flow within a lobe is one of the simplest assumptions that can be made for branches beyond the truncation, without involving multiscale models [23,24,26,30]. Figures 2(b) and 2(c), respectively, show the lobar and terminal flow rates (as fractions of the total flow rate), which we used as boundary conditions for the CFD simulation.

### Measure of Airflow Similarity.

*R*) and the root-mean-square of differences (

*D*

_{rms}) for velocity contours at axial planes. Because we computed the airflows in unstructured grids for the CFD simulations, we rediscretized the cross-sectional plane into rectangular lattice points (

*x*,

_{i}*y*) by using matlab (R2016a, MathWorks, Natick, MA). We assigned the corresponding velocity magnitude (

_{j}*v*,

_{i}*v*) as matrix element (

_{j}*i*,

*j*) and defined

*R*as follows:

where $A\xaf$ and $B\xaf$ denote the mean velocity values of matrices **A** and **B**, respectively.

## Results

### Comparison of In Vitro and In Silico Airflows in the Idealized Airway.

We found broad agreement between the PIV-measured and computed velocity contours, albeit with slight regional variations in velocity magnitudes, for both healthy I and II. Because the velocity contours simulated for these conditions were indistinguishable, we compared the measured velocity contours (Fig. 3(a)) only with the simulated contours for healthy I (Fig. 3(b)).

For both measured and computed results, we observed a high-velocity stream along branches *a*, *c*, and *e* of the right lung, which subsequently continued into branch *i*. In both results, the high-velocity domain in these branches was skewed medially (i.e., in the positive and negative *x* directions for the right- and left-lung branches, respectively; Figs. 3(a) and 3(b1)) from the outset, with a peak velocity of 4.8 m/s at branch *a*. This flow asymmetry in branch *a* likely occurred because 76% of the airflow came from terminals 6–8, where the momentum was in the medial direction. Further downstream, the peak velocity in branch *c* and that in branch *e* increased relative to that in branch *a*, and the flow remained skewed in the same direction.

Similarly, we found general agreement between the measured and computed velocity contours for the left-lung branches. For both results, the peak velocity values in branches *f* and *h* were less than those in branches *e* and *c* of the corresponding generations in the right lung (Figs. 3(a) and 3(b1)). In contrast, the stream of flow in branch *g* had a peak velocity higher than that in branch *f* or *h*.

The measured and computed velocity contours in branches *e* and *i* did not agree as well as those in other branches. Although the peak velocity in branch *e* was asymmetric for both experimental and computational results, its measured magnitude was 19% lower than the computed magnitude (Figs. 3(a) and 3(b1)). In branch *i*, although the computational and experimental results both showed a high-velocity stream that likely continued from branch *e*, the computed peak velocity was skewed toward the positive *x* direction at the entrance of branch *i* (Fig. 3(b1)), unlike the measured peak velocity, which was skewed slightly in the opposite direction (Fig. 3(a)).

### Computed Airflow Characteristics in the Idealized Airway.

Given the broad agreement between the experimental and computational results, we next examined the flow characteristics in the idealized airway in greater detail. We computationally analyzed the axial velocity contours and the directions of secondary flows (arrows) in the axial plane (Figs. 3(b1) and 3(b2)-b). In the right-lung branches, the axial velocity contours changed progressively along branches *a*, *c*, and *e*, with the added flow momentum from branches *b* and *d* in the medial direction increasing the skewness of the high-velocity domain and the Dean numbers. The Dean numbers were 430 at branch *a*, 1200 at branch *c*, and 1700 at branch *e*. In branch *a*, the shape of the axial velocity contour was elliptical, and the secondary flow was directed toward the positive *x* direction. In branches *c* and *e*, the contours were crescent-shaped, and their peak velocity values increased, respectively, by 17% and 8% relative to the value in branch *a*. In addition, we observed two pairs of counter-rotating vortices in the axial plane of branch *c*, as reported in other airway geometries [20,31,47]. Similar vortical structures were present in the axial planes of other branches (*f*, *h*, and *i*).

The computational results revealed flow separation and accompanying secondary flow structures in curved branches *b*, *d*, *f*, and *h*, whose Dean numbers were 490, 560, 520, and 820, respectively. In branch *b*, the flow was almost symmetric in the axial plane near the entrance of the branch, with a high-velocity domain at the center (Figs. 3(b1) and 3(b2)-b). However, further downstream, the high-velocity domain was skewed toward the outer curve of the branch, accompanied by flow recirculation on the opposite side (Fig. 3(b1), red arrow). In branch *d*, the flow was separated along a broad area at the inner curve of the entrance from terminal 1 (Fig. 3(b1), black arrow). Unlike the flow in branch *b*, the flow in terminal 1 maintained a straight path until it entered branch *d*. However, it then turned direction sharply by 40 deg and separated at the inner curve of branch *d*. This flow separation was mitigated downstream after flows from terminals 2 and 3 merged in branch *d*, although the high-velocity domain became highly skewed toward the negative *y* direction and crescent-shaped in the axial plane (Fig. 3(b2)-d). Interestingly, in branch *f*, the flow recirculated at its outer curve (Fig. 3(b1), blue arrow at branch *f*). The high-velocity domain in branch *f* was skewed medially at its entrance because of the imbalance in the flow rates upstream. The high-velocity stream started from the medial side, where 73% of the flow merged from terminals 10–12 and 27% entered from terminal 13. Although, the curvature of branch *f* directed the flow slightly toward the positive *x* direction, its effect was not sufficient to separate the flow at its inner curve as in branches *b* and *d*. Instead, flow recirculated at the outer curve near the junction to branch *g*.

In branch *h*, we observed flow separation along its inner curve (Fig. 3(b1), gray arrow at branch *h*). After the two streams of flow from branches *f* and *g* converged, the flow in branch *h* formed two pairs of counter-rotating vertical structures with a peak velocity at the center of the axial plane (Fig. 3(b2)). This presumably occurred because the amounts of flow from branches *f* and *g* were equal; hence, the flow from these branches likely exerted comparable forces upon one another when they converged in branch *h*. Further downstream, the high-velocity domain was displaced toward its outer curve owing to the influence of the curvature.

### Airflow Sensitivity to Variation in Terminal Flow Rates.

To examine how variation in the terminal flow rates (i.e., boundary conditions) affects flow patterns relative to those for healthy I (used on day 1), we computed flows for three different boundary conditions—healthy II (used on day 2), healthy III (directly derived from the healthy subject), and COPD (directly derived from the diseased subject; Fig. 2(c))—using the same geometry and same total flow rate (770 ml/s). The terminal flow rates for healthy II and III were different from those for healthy I, whereas the lobar flow rates were similar for all three conditions (Fig. 2(b)). The terminal flow rates for COPD, which were obtained from the lobar flow rates estimated for the COPD subject, differed from those for healthy I.

Figure 4 shows the velocity contours computed for the four conditions in the central airways (Fig. 4(a); coronal plane) and upper airway (Fig. 4(b); sagittal plane). In the central airways, including airways from the trachea and bronchial trees, the computed velocity magnitudes for healthy II and III were similar to those for healthy I. In contrast, the computed velocity magnitudes for COPD were different from those for all other conditions. For the healthy conditions, high-velocity flows commonly present along the right-lung branches *a*, *c*, and *e* propagated to branch *i*. In contrast, the computed high-velocity flows for COPD originated from terminals 1 and 2, merging at branch *d* and continuing in branches *e* and *i*.

In the left-lung branches, the velocity profiles for healthy II and III were similar to those for healthy I. In contrast, the velocity profile for COPD differed from those for the other three conditions. The peak velocity values in branches *h* and *f* decreased by 26% and 55%, respectively, and increased by 17% in branch *g* relative to healthy I.

The changes in the velocity profiles for COPD were accompanied by flow separation and recirculation in curved branches. For COPD, we did not observe the flow recirculation commonly present at the inner curve of branch *b* for the healthy conditions (Fig. 4(a), red arrow at branch *b*). Similarly, for COPD, flow did not separate at the outer bend of branch *f*, where it recirculated for the healthy conditions. Instead, in branch *h*, flow separation occurred in a domain larger than that in the healthy conditions, even though the effect of branch curvature was expected to decrease because of the reduced flow rate. The enlarged flow separation may be due to the imbalance of upstream flows in branches *f* and *g*. For COPD, the jetlike stream of flow in branch *g* had a peak velocity 20% higher than that for healthy I. Because branches *g* and *h* have the same bending direction, the jetlike stream in branch *g* for COPD continued to turn in the same direction in branch *h*; thus, the flow was continually subjected to centrifugal force and underwent separation in a larger domain of the inner bend in branch *h*. For the healthy conditions, this effect was counteracted by the flow from branch *f*, which underwent bending in the opposite direction. However, for COPD, the flow rate in branch *f* was 38% of that for healthy I (Fig. 2(c)); hence, the flow from branch *f* had less of a counteracting force on that from branch *g*.

The flow differences observed in branches *e* and *h* for COPD propagated to branch *i*. Because the amount of flow entering from branch *e* was 20% higher than that from branch *h*, the stream along branch *i* was pushed slightly more toward the positive *x* direction. However, this difference vanished in the upper respiratory tract. For all simulated conditions, we observed a strong jet in the larynx that accompanied flow recirculation in the expansion near the glottis. The jet flow continued to expand in the oral airway in all cases (Fig. 4(b)).

Figure 5 shows the computed axial velocity contours for the four conditions in the axial plane at branches *a*–*i* (Fig. 5(a)) and the metrics (Figs. 5(b) and 5(c))—the correlation coefficient, *R* (Eq. (2)), and the root-mean-square of differences in velocity, *D*_{rms} (Eq. (3))—which reflect the differences in these contours. We used *R* to measure the similarity in the shapes of the velocity contours. *R* is 0 if the shapes of the contours have no similarity and 1 if they are identical. In computing *R*, we excluded low-velocity domains (i.e., less than 50% of the maximal value) from the contours because they consistently resided near the wall owing to the no-slip boundary condition. *R* values were not highly sensitive to the threshold values used to exclude the lower velocity domains from the contours: for example, at plane *h*, *R* increased by 5% when the threshold was lowered from 50% to 30%. However, *R* alone was insufficient to distinguish flow differences if the contours had similar shapes but different magnitudes, as would be expected from the insensitivity of *R* to affine scaling. Therefore, we used *D*_{rms} to highlight differences in velocity magnitude between contours. For example, in comparing healthy conditions and COPD for branch *a* (Fig. 5(a), top), *R* was 0.97 in all cases (Fig. 5(b), leftmost points) because their velocity contours showed patterns similar to those of healthy I. However, the *D*_{rms} values for branch *a* were 0.23, 0.15, and 1.05 m/s for healthy II, healthy III, and COPD, respectively (Fig. 5(c), leftmost points).

The *R* value for COPD was distinctly smaller than the values for healthy II and III (Fig. 5(b)). The difference between the healthy conditions and COPD was largest in branch *h*, where *R* for COPD was 0.26—almost three times less than those for healthy II (*R* = 0.89) and III (*R* = 0.96). Similarly, in branch *c*, COPD showed an *R* value twice as small as that of healthy II. In contrast, *R* values for healthy II and III ranged between 0.80 and 0.97, respectively, in all branches. They were larger than that for COPD in all other branches except branch *a*.

The *D*_{rms} values for healthy II and III were similar in all branches, whereas those for COPD were much larger in most branches (Fig. 5(c)). For example, healthy II and III showed a maximal difference of 0.45 m/s in branch *b*. In the same branch, *D*_{rms} for COPD was 1.73 m/s larger than that for healthy III. The difference between the healthy conditions and COPD was maximal in branch *d*, where the *D*_{rms} value for COPD was almost 10 and 12 times the values for healthy II and III, respectively.

## Discussion

We used an idealized airway geometry to measure and compute the flow characteristics for a healthy lung. We further examined how such characteristics were affected by changes in the boundary conditions, which defined the flow distribution among individual terminal branches. We found that variations in the terminal flow rates among three conditions (healthy I, II, and III) based on a healthy volunteer had little impact on flow patterns in both the central and upper airways, whereas a notable difference in lobar flow rates in a diseased condition (i.e., COPD) markedly affected airflow patterns in the central but not upper airways.

### Potential Sources of Differences Between In Vitro and In Silico Flow Fields.

In most branches, we found broad agreement between the experimentally obtained and computed velocity contours. Consistent with the computed results, the experiment successfully captured the asymmetry of the flow along the right- and left-lung branches. In addition, with the exception of branch *e*, the peak velocity magnitudes measured by PIV differed from the computed results by less than 6% in each of the branches.

In branch *e*, the measured peak velocity magnitude was 19% smaller than the computed result. This discrepancy increased downstream in branch *i*. Compared with the computed result, the measured peak velocity was 21% smaller and located slightly toward the right side of the lung (negative *x* direction). The mismatches between the PIV measurements and computed results could be attributed to uncertainties in our experimental methods, numerical methods, or both. At locations where a mismatch was observed, the flow was in the transitional regime with a Reynolds number above 3000. Although direct numerical simulations can provide an accurate description of the flow, this was not feasible here because the temporal and spatial resolution required to sufficiently compute turbulent fluctuations would have been too high given the complexity of the airway geometry and flow structure. In this study, we employed the Reynolds-averaged Navier–Stokes equations with a two-equation *k*–*ω* turbulence model. Although this model was established to describe flows in the turbulence regime [48], it may not have adequately captured turbulent velocity fluctuations [49]. Specifically, a discrepancy could have occurred between the trajectories of the microspheres and those of the computed flows because the microspheres did not show perfectly neutral buoyancy. If so, then the PIV-measured velocity profile should be broader than the computed results. Indeed, the discrepancy between the measured and computed results increased in large airways where the Reynolds numbers were higher than those in other branches, whereas in smaller airways (e.g., branches *a* and *c*) the measured and computed results showed minimal differences in peak velocity values. Such a reduction in the peak velocity magnitude measured in branch *e* would likely reduce the momentum of flow directed toward the positive *x* direction relative to the computed results and shift the peak velocity toward the negative *x* direction in branch *i*.

### Secondary Flows in the Idealized Geometry for a Healthy Condition.

We also observed flow separation and recirculation in airway branches *b*, *d*, *f*, and *h.* For any particular branch, these effects can be attributed to three factors: its geometrical characteristics, the Reynolds number for the flow within it, and the balance of flows upstream of it. The curvature of these branches contributed to the development of a pressure gradient in the transverse direction, which was large enough to separate the flow and promote recirculation. Although the curvature ratio (i.e., the ratio of the curvature radius to branch diameter) ranged between 4 and 8 for these branches—which is larger than the ratio (1.5) at which flow separation is likely to occur owing to the adverse pressure gradient near the inner wall of the bend [50,51]—previous studies have shown flow separation in bent conduits with a curvature ratio between 1 and 5, similar to the range examined here (4–8), with sufficiently large Reynolds numbers [20,47,52].

In our lung geometry, flow recirculation required a Reynolds number above 1500. The Reynolds numbers in branches *b*, *d*, *f*, and *h* ranged between 1500 and 3000 and were higher than those in smaller branches. For example, we observed recirculation of flow at the inner curve of branch *b*, whose curvature ratio was 8. However, flow separation did not occur in the branch where the flows from terminals 14 and 15 merged, even though the branch had a curvature ratio of 1.2. This is presumably because the Reynolds number of the flow in this branch was 870, which was less than that for the flow in branch *b* (i.e., 1600).

In addition to a high Reynolds number and a small branch curvature ratio, the upstream flow profile is an important factor in the development of flow separation. For example, the flow from branch *e* maintained a straight streamline until it curved to enter branch *i*. At this junction, the curvature ratio was 3 (less than that of branch *b*). In addition, the Reynolds number was 3100 (higher than that of branch *b*). Nonetheless, we did not observe flow separation. This is presumably because the centrifugal force due to bending in branch *e* was counteracted by the flow from the neighboring branch *h*; of the flow from branch *h*, 49% merged at this junction with a momentum opposite to that from branch *e* along the *x* direction. Furthermore, the upstream flow profile changed the site of recirculation in the branch. In branch *f*, flow recirculated at the outer curve of the branch but not at the inner curve, and a stream of high-velocity flow from upstream branches entered the branch highly skewed toward its inner curve. Although the curvature ratio was 6.3 (less than that for branch *b*) and the Reynolds number was 2100 (higher than that in branch *b*), the flow did not separate at its inner curve because the high-velocity flow from upstream remained skewed toward the inner curve. In contrast, the area of the branch increased toward the outer curves while a jetlike high-velocity stream remained skewed on the other side. Subsequently, the flow expanded at the outer curve and underwent recirculation.

These considerations suggest that upstream flows are important factors in determining the flow patterns in central airways in addition to the local airway geometries. In the Insensitivity of Airflow Patterns to Uncertainty in Terminal Flows for Healthy Conditions section, we discuss the extent to which variations in upstream flows lead to changes in flow patterns further downstream.

### Insensitivity of Airflow Patterns to Uncertainty in Terminal Flows for Healthy Conditions.

We examined how upstream flows in the terminal branches affected the flow patterns by comparing those computed for healthy I with those for three other conditions. Our results showed that the flow patterns examined were insensitive to differences in individual terminal flow rates if their lobar flows were similar. The flow fractions in individual lung lobes for the healthy conditions were similar to those estimated from CT images of a healthy subject (Fig. 2(b)). Although their flow rates at individual terminals varied by up to 40% for healthy II (Fig. 2(c), terminal 16) and 140% for healthy III (Fig. 2(c), terminal 3) relative to those of healthy I, their flow patterns in major branches *a*–*i* showed common characteristics, especially in terms of the asymmetry of flow and the development of flow separation.

In particular, the differences in the terminal flow rates between healthy I and II represent the range of uncertainty in measuring and controlling the flow rates in this experiment. However, their flow patterns were similar in both axial and coronal planes; the correlation coefficients were larger than 0.89 (in branches *c* and *h*), and the root-mean-squares of velocity magnitude differences were less than 0.57 m/s (the value in branch *b*) in all branches. This suggests that the experimental uncertainty in the terminal flows, defined by the variation observed between healthy I and II, do not meaningfully alter the velocity patterns.

As described in the Methods section, in healthy III, flow rates were assumed for simplicity to be equal among the terminal branches within the same lung lobe when only the lobar fractions were known. A comparison of the flows for healthy III with those for other healthy conditions revealed no notable differences in the flow patterns. In addition to the similarity in the locations of major flow streams in the coronal plane, the shapes of the velocity contours and the velocity magnitudes at each axial plane were similar across the healthy conditions (Fig. 5(a)). The correlation coefficients for flows in the axial plane were greater than 0.80 (branch *c*), and the root-mean-squares of the velocity magnitude differences were less than 1.02 m/s (branch *b*). These results suggest that velocity fields may not be sensitive to flow rates in the terminal here for the healthy conditions.

### Differences in Flow Characteristics Between Healthy and Chronic Obstructive Pulmonary Disease Conditions.

To compare the flows for a diseased condition with those for the healthy conditions, we used the same geometry to investigate the flow characteristics for a COPD condition. In modeling the COPD condition, we analyzed CT images of a COPD subject and used the obtained flow fractions for individual lung lobes. Notably, the COPD subject examined in this study showed flow fractions distinctively different from those for the healthy subject. The flow fractions for the healthy subject were largely uniform among the lung lobes: right upper (21%), right middle (12%), right lower (20%), left upper (25%), and left lower (22%). In contrast, volume changes between inhalation and exhalation conditions in the COPD lung occurred mostly in the right upper (43%) and left upper (29%) lobes and only sparingly in the other lobes (middle right (5%), lower right (13%), and lower left (10%)). The agreement between these lobar flow fractions was poor (intraclass correlation coefficient = 0.3) [46].

The heterogeneity observed in the lobar flow fractions for COPD was manifested as an altered flow pattern in central airways. In particular, the main streams of high-velocity flow for COPD were from the upper lobar branches, with lower velocity magnitudes relative to healthy I in the remaining branches (Fig. 4(a)). Such changes deformed the shapes of the velocity contours in the axial plane, whose similarities to those of healthy I were quantified by correlation coefficients. In comparison, COPD showed noticeably smaller correlation coefficients in all branches. Interestingly, the greatest deformation in the contour shapes did not occur at the branches immediately downstream of the terminal branches, as might have been expected. Instead, the smallest correlation coefficients, 0.26 and 0.61, were observed at branches *h* and *e*, respectively. This implies that changes of flow patterns in upstream branches are augmented in downstream branches under expiratory flow conditions. In addition, our results showed that the area of flow separation changed in COPD, presumably because of the altered flow distribution among the upstream branches. For COPD, flow recirculation increased in area relative to the healthy conditions in branch *h*, whereas it disappeared in branches *b* and *f* where the flow recirculated locally.

In contrast to the flow in the central airways, the flow in the upper respiratory tracts was insensitive to the terminal flow conditions. The difference in flow patterns between healthy and COPD conditions propagated downstream to the tracheal region (i.e., branch *i*), although it was not as pronounced as in other branches. However, in the upper respiratory tracts for all conditions, a jet stream developed at the constriction in the larynx, accompanied by a large recirculation around the jet due to expansion and a large pressure drop, which are commonly observed for expiratory flow [52–54] and often associated with pathological conditions such as obstructive sleep apnea [9,55,56]. Consequently, the flow patterns downstream of the larynx were predominantly determined by the jet stream.

Our results demonstrate the importance of using lobar flow fractions relevant to the disease condition to study the flow characteristics in central airways. In particular, compared with healthy conditions, COPD showed a different velocity profile in terms of flow asymmetry or the location and size of flow recirculation as well as a different axial velocity profile. Such changes can produce errors in predicting disease- or patient-specific respiratory physiology. For example, the efficacy of drug delivery or particle dosimetry in lungs depends highly on the velocity fields within the geometry. In particular, recirculation of flow promotes deposition efficiency in its vicinity by trapping aerosols or particles dominated by impaction and interception [57]. Therefore, neglecting to consider differences in patient-specific lobar boundary conditions may lead to inaccurate predictions of airflow characteristics and particle dosimetry. Using accurate lobar flow fractions could contribute for predicting airflow characteristics specific to disease and patient cases. Thus, there is an urgent need to develop more comprehensive computational models [58–60].

## Limitations

One limitation of this study is that we used 2D PIV instead of 3D PIV to measure velocity fields in coronal planes at multiple locations for the in vitro experiment. Although 3D PIV yields more information [61,62], 2D PIV provided sufficient information to validate the computational model and show the effect of varying the boundary conditions. The computational results and their ensuing conclusions were primarily based on streamwise velocities (within the *xy*–– plane in Fig. 3). Therefore, our measurements of both *x* and *y* components within the coronal plane would appear to be sufficient and appropriate. Future studies that investigate either mixing, or phenomena more closely related to nonstreamwise velocity, might warrant 3D PIV measurements.

The other major limitation is that we assessed the influence of flow distribution for healthy and COPD boundary conditions and made a number of simplifying assumptions to isolate the effects of boundary conditions on the flow patterns while controlling other factors that influence airflow characteristics. First, we assumed that the airway was rigid and static. A number of computational [31–33] and experimental studies [18,63,64] have made this assumption largely because accurate data on the material properties of the airway walls or airway motion in vivo are lacking. However, large airway branches, including the upper respiratory airway and trachea, may collapse under pathological conditions [65,66]. Airways, by undergoing expansion or collapse to reduce the pressure load, are likely to have altered pressure and flow velocity profiles. The extent to which velocity profiles are sensitive to airway dynamics is unclear, especially in comparison with other factors that determine airflow patterns. To address this issue, more comprehensive computational models that include accurate descriptions of airway dynamics will need to be developed.

Second, we used an idealized airway geometry for both healthy and diseased cases. However, this simplified model is likely to be inadequate in addressing specific problems. For example, the central airway branches in this model had circular cross sections lying in a single plane. Considering that flow patterns are influenced by bending and rotation in the geometry [67], those at individual branches are likely to be different from those computed for patient-specific airways. Nonetheless, if the branches immediately upstream are sufficiently long, we expect that even in a planar geometry the asymmetry and formation of secondary flows would be qualitatively similar to those in a nonplanar geometry, as long as the bifurcation angle and bending curvature of the branches remain the same.

Third, we focused on steady-state unidirectional expiratory flow rather than flows during the full breathing cycle. We previously showed that flow fields are affected by upstream but not downstream flows [31]. Expiratory flow is more adequate than inspiratory flow for assessing the effects of lower airway boundary conditions on flow patterns. We expect that the effects of unsteady inlet flow will be non-negligible [20,68] when the inlet flow rate changes rapidly, as is found in forced breathing conditions. Nonetheless, if the change in the inlet flow is gradual, the influence of unsteadiness in the inlet flow would be minimal, and velocity patterns would be similar to those for steady-state inlet flow [69].

Fourth, we derived the flow distribution among the lung lobes by using two sets of lung images measured after peak inspiration and peak expiration, because data at intermediate time points were unavailable. However, the flow distribution in the lung changes nonlinearly over time and depends on the breathing condition [20,70]. In addition, our study assessed flow sensitivity by using exemplar lung images from one healthy subject and one COPD subject, although intersubject variability is expected within a healthy or diseased group [71,72]. To develop more comprehensive models, studies that simultaneously measure both HRCT lung images and breathing flow rates at multiple time points will be required [70,71].

Finally, although the focus of our study was on isolating the effects of lobar flow fractions associated with healthy and COPD lung conditions, other physiological parameters, such as airway structure and flow rates during forced breathing, also differ between these conditions. Typically, FEV_{1} values are lower for COPD cases than for healthy cases [73]. Accordingly, the average velocity magnitudes at individual branches for the COPD case will be lower than for the healthy case, and the difference in their flow profiles, quantified in terms of *D*_{rms}, will increase. Similarly, if diseased lungs are accompanied by structural deformations, such as bending and constriction of airway branches, we would expect more heterogeneous flow patterns, induced by flow recirculation and jet streams near obstructions, compared to healthy cases [31,74–76].

## Conclusions

In this study, we performed experimental and computational studies to examine the airflow characteristics for a healthy lung condition and found broad agreement in the velocity fields in most branches. We then used our validated CFD model and investigated how variation in terminal flows affects airflow patterns in a steady expiratory condition. We derived three healthy conditions and one COPD condition from the lobar flow fractions acquired from HRCT images of one healthy lung and one COPD lung, respectively, while using the same airway geometry and expiratory flow rate for all conditions. The flow patterns for COPD, derived from the lobar flow fractions of a COPD subject, were distinctly different from those for three conditions derived from a healthy subject. However, the three healthy conditions showed similar flow patterns despite differences in their individual terminal flow rates. These results specify the range of uncertainty allowable in developing computational models or designing experimental studies to investigate airflow characteristics in airways. In addition, by revealing differences between healthy and COPD flow patterns, our results highlight the importance of employing boundary conditions based on accurate estimation of lobar flow fractions.

## Acknowledgment

High-performance computing resources were made available by the U.S. Department of Defense High Performance Computing Modernization Program. The authors thank Dr. Ching-Long Lin and Dr. Eric Hoffman for providing the HRCT imaging protocol for lobar flow fraction analysis.

The opinions and assertions contained herein are the private views of the authors and are not to be construed as official or as reflecting the views of the U.S. Army or of the U.S. DoD.

## Funding Data

This research was sponsored by the U.S. DoD Defense Health Program, managed by the Military Operational Medicine Program Area Directorate, U.S. Army Medical Research and Materiel Command, Fort Detrick, MD, and by the U.S. Army's Network Science Initiative.