Computational Model Upper Generations Respiratory Tract

Computational fluid dynamics (CFD) studies of airflow in a digital reference model of the 17-generation airway (bronchial tree) were accomplished using the FLUENT® computational code, based on the anatomical model by Schmidt et al. (2004). The lung model consists of 6.744•106 unstructured tetrahedral computational cells. A steady-state airflow rate of 28.3 L/min was used to simulate the transient turbulent flow regime using a Large Eddy Simulation (LES) turbulence model.

This CFD mesh represents the anatomically realistic asymmetrical branching pattern of the larger airways. It is demonstrated that the nature of the secondary vertical flows which develop in such asymmetric airways varies with the specific anatomical characteristics of the branching conduits.

Calay et al. (2002) studied the respiratory flow patterns in a single first generation bifurcation distal to the trachea and in a multiple-bifurcation model, with their three generations based on the anatomy given by Horsfield et al. (1971). They used four different grid densities comprising respectively 31,104, 79,820, 159,872 and 320,980 nodes. Lee and Lee (2002) generated a 3-D conduit network model with four generations which conformed closely to Weibel’s (1963) model in order to study aerosol bolus dispersion. The total number of cells in their 90˚ out-of-plane model ranged from 40,000 to 60,000. Liu et al. (2003) studied the 3-D inspiratory flow in a 3-generation asymmetric model from the 5th to the11th branches of Weibel’s B model. Kleinstreuer and Zhang (2003) analyzed targeted aerosol drug deposition analysis in a rigid triple-bifurcation tracheobronchial airway model. Their model represents the symmetrically bifurcating generations G3 to G6 as in Weibel’s lung model with three different hemispherical tumor models placed along the side-wall of the G5-airway. Their final mesh comprised about 360,000 cells for the triple bifurcation configuration. Nowak et al. (2003) demonstrated a four-subunit CFD simulation method for the human tracheobronchial tree. They economized on computational effort by segmenting the first 12 generations into four 3.5-generation “tranches” (G0-G3, G3-G6, G6-G9, and G9-G12), without adequately accounting for the fluctuation in boundary conditions between each tranche. Their results indicate dramatic differences in the predicted particle deposition patterns between the two models. The absence of airway curvature and surface irregularities in a Weibel-based model renders the flow fields very different from those in a real human lung (CT-based model). Cebral and Summer (2004) studied the central tracheal and bronchial airways down to 4 generations by using a virtual bronchoscopy reconstruction method. Airflow patterns resulting from airway stenoses in generations from G0 to G4 were simulated computationally. van Ertbruggen et al. (2005) studied the gas flow and particle deposition in a realistic 3D model of the bronchial tree, extending from the trachea to the segmental bronchi (7th airway generation for the most distal airways), based on the morphometrical data of Horsfield et al. (1971). Considering symmetric double bifurcation models, Longest and Vinchurkar (2007b) have recently assessed the effects of upstream transition to turbulence on the flow field and particle deposition in the generations G3-G5 of the respiratory tract. Turbulence was shown experimentally to influence the local deposition of 10 mm-diameter particles, primarily by influencing the initial velocity and particle profiles. Their results underline the importance of correct inlet conditions and the need to consider upstream effects in experimental and computational studies of the respiratory tract. In addition, Longest et al. (2006) compared flow patterns and particle deposition in both normal and childhood asthma-induced constricted models of pulmonary tree generations G3-G5 and G7-G9. Both laminar solutions and the low Reynolds number (LRN) κ-ω turbulence model were employed using a Fluent-6 software package along with their own Fortran and C programs. The authors concluded that airway constriction caused local cellular-level deposition rates to increase by one to two orders of magnitude. Asthma constriction may significantly increase branch-averaged particle deposition with larger increases in local cellular-level deposition, resulting in an aggravated health risk.

More realistic anatomical models have been published by Kitaoka et al., 1999; Tawhai et al., 2000 and 2004; Spencer et al., 2001; Tgavalekos et al., 2003 ; Sera et al., 2003 ; Schmidt et al., 2004, using both mathematical algorithms and new experimental imaging techniques. Kitaoka et al. (1999) introduced a 3D-model of the human airway tree down to the terminal bronchioles which by a deterministic algorithm that incorporated duct branching and space division. Tawhai et al. (2000 and 2004) developed a 3D tree-growing algorithm specific to a given host geometry derived from Magnetic Resonance Imaging (MRI). Spencer et al. (2001) developed a dynamic surface modeling technique based on data from idealized models to construct 3D computer simulations of tubular pulmonary airway structures within lungs extending from the trachea (G0) to the alveoli (G23). Sera et al. (2003) developed a two-step method to visualize small airways in detail by staining excised rodent pulmonary airways with a radio-opaque solution and then visualizing the tissue with a cone-beam microfocal X-ray Computed Tomography (CT) system. Tgavalekos et al. (2003) advanced the 3D airway tree model of Kitaoka et al. (1999) to predict pulmonary function on the basis of airway structure, particularly when constriction patterns are imposed heterogeneously on the pulmonary tree in specific anatomic locations and compared their model predictions with ventilation images obtained from Positron Emission Tomography (PET) and measurements of dynamic mechanical pulmonary function. The current study includes these upper generations as well as both central and lower tracheobronchial airways within the new 17-generation model.

2. Current 17-Generation Anatomical Model

We have based our present CFD analysis on the best published anatomically explicit human lung model available to date; i.e. the 17-generation anatomical model developed by Schmidt et al. (2004). Fig. 1 represents the anatomy of the human lung of an adult male, free of pathological alterations. Fig. 1a portrays a realistic 3D surface representation of the segmented volume of the bronchial tree. The conduit model reconstructed from an abstracted and adapted topological graph, as shown in Fig. 1b, contains 1453 bronchi up to the 17th Horsfield order. To the authors’ knowledge, the computational study described herein represents the first published CFD simulation of the airflow in a digital reference model up to a maximum 17 generations of the total 26 generations of the human tracheobronchial airways with 1453 bronchi. The computer-generated symmetric human lung model of Spencer et al. (2001) for the 23-generation tree lacks the required anatomical detail and accuracy for the current computational investigation and is not as accurate as the Schmidt et al. (2004) anatomical airway model.

The mouth and glottis have not been included in the present 17-generation pulmonary model. It is unlikely that this exclusion will have a direct effect upon the downstream flow as there is a gradual recovery of pressure by the distal end of the trachea (Ma and Lutchen, 2006). The pressure difference, averaged over the luminal cross section from the proximal to the distal end of the trachea, is insignificant in comparison to the pressure difference between the upper and lower airways as indicated by Ma and Lutchen (2006), and the viscous pressure drop in the trachea was found to be negligible. Thus the trachea appears to have negligible influence on total lung resistance.

3. CFD Simulation

Our characterization of the anatomy of airways must be sufficiently realistic as we intend to extend the present model to design more efficient delivery of inhaled medications and to understand better the health effects of inhaled pollutants. Particle deposition patterns in the branching tracheobronchial network are highly influenced by the detailed branching structures of the human lung, which are important for optimizing aerosol therapy protocols.

In this study, the incompressible airflow in a comprehensive digital reference model was computed using the commercial CFD code FLUENT® (2004). FLUENT employs a finite-volume method to solve the Navier-Stokes and continuity equations on an arbitrarily shaped flow domain with appropriate boundary conditions and incorporates various two-parameter models (k- and k-) turbulence models and LES (Large Eddy Simulation). Our computations were performed on two different supercomputing platforms: 1) a RackSaver® cluster at the National Supercomputing Center for Energy and Environment NSCEE at UNLV and 2) a tera-scale machine at the Pittsburgh Supercomputing Center. The steady-state solution of the flow field was assumed to have converged when the residuals reduced to less than 10-4. A typical run time took approximately 50 hours on the NSCEE machine.

Computations were performed at 28.3 L/min for a quasi steady-state volumetric adult inhalation flow rate and as well as at 6, 12, 18, 24, and 30 L/min for the study of static pressure variation.

3.1 Human Bronchial Tree Model and Mesh Generation The computational mesh is generated using GAMBIT® software, after transporting into it the abstracted topological graph data from Schmidt et al. (2004). The tube-like surface representation shown in Fig. 1b still retains the exact anatomical asymmetry and branching angles based on a computationally affordable mesh size consisting of 6.744106 unstructured tetrahedral cells as shown in Fig. 2. The mesh generation file is 660MB and takes one hour to develop on the terascale machine of PSC. Longest and Vinchurkar (2007a) applied a Grid Convergence Index (GCI) to four mesh styles of a double bifurcation model and established grid convergence criteria for both the flow field and particle deposition results using Fluent 6 software. The number of cells needed to resolve fully the geometry of this 17-generation computational model is currently not practical. Due to current hardware limitations, one can only use approximately 6 million control volumes. Full grid convergence will necessarily be left to a future work.

3.2 Numerical MethodsThe steady-state mass conservation (Eq. 1) and momentum equations (Eq. 2) were solved using a finite-volume method to compute the particle-free fluid flow field. Air was assumed here to be an incompressible, Newtonian fluid with constant density , viscosity and fluid static pressure . The continuity and steady-state Navier-Stokes equations for viscous incompressible Newtonian fluid are given in the original paper.

In order to model the turbulent airflow, we employed a LES (Large Eddy Simulation) model, in which large energy–containing eddies are computed directly, while the remaining smaller eddies are modeled using a subgrid-scale (SGS) turbulence model.

3.3 Boundary Conditions The quasi steady-state inlet velocity is taken here as 2.896m/s, based on a tracheal inlet area of 162.86mm2 and a tracheal diameter of 14.4mm. A uniform inlet velocity field boundary condition is applied at the tracheal inlet section (G0). The pressure at each outlet (total of 720 outlets) must be known a priori, reflecting the effect of the remaining truncated generations G18 through G26 at the alveolar level. An alternative approach to setting the flow division ratios would be to use the flow ratio values given by Horsfield et al. (1971). Ma and Lutchen (2006) attempted this, but the pressures computed at the thirteen outlets of the G6 large airway model were unacceptably inaccurate, since this is a pressure-driven flow. Accordingly we recommend the use of mass distributions to determine outflow boundary conditions. The conclusions of Ma and Lutchen (2006) regarding flow versus outlet pressure boundary conditions at the end of G6 do not appear relevant to our model which comprises almost the complete anatomical replica of the pulmonary tree down to the 17th generation level, at which level a constant output p = 0 boundary conditions is clearly more appropriate than at the upstream G6 level. From Ma and Lutchen (2006) it appears that turbulence does not have a significant effect on the total input impedance.

We used an inflow boundary condition at G0 and specified pressure outlet boundary conditions at all terminal bronchi exits. The total diameter-dependent mass flow rate at all bronchial exits approximates the inlet value to within 11.210-6 g/sec. No-slip boundary conditions were assumed at the tracheobronchial airway walls. Generation G0 and the first bifurcation were represented by a less dense mesh (1.063106 cells) than for the remaining generations (5.680106 cells) as the flow was more uniform at this upstream location. The computational mesh, with a total of 6.744106 cells, was divided into 1.063106 tetrahedral cells in the G0 zone and into 5.680106 unstructured tetrahedral cells in the remaining zones.

4. Results and Discussions

4.1. Tracheobronchial Model up to a Maximum 17 Generations During inhalation, the flow transitions from laminar to turbulent flow in the larynx and in the tracheal section, and continues into the first main generation. In the conducting zone, containing the trachea, bronchi, bronchioles, and the terminal bronchioles, there is no gas-blood exchange. The downstream increase in cross-sectional area causes a drop in airway pressure and flow rate, such that final air-blood exchange occurs by diffusion at the alveolar level. Our CFD calculations indicate that the static pressure drop is 50 Pa for a volumetric gas flow rate of 28.3L/min over the complete span as shown in color-coded Fig. 3. Pressure drop in the bifurcating airways plays an essential role in the respiratory process. Secondary flows in the tracheobronchial airways can be very pronounced. Pedley et al. (1971) and Pedley (1977) determined the static pressure variation to be 60 Pa along Weibel’s symmetric bronchial tree, in reasonably good agreement with the experimental measurement by Hyatt and Wilcox (1963) of 75 Pa across the entire respiratory tract. Our comprehensive 3D CFD simulation of the asymmetric bronchial tree model predicted a reasonable static pressure drop of 50 Pa for a steady-state flow rate of 28.3L/min, corresponding to normal adult breathing conditions. The viscous pressure drop increases nonlinearly from 6 Pa at 100 cm3/s to 54 Pa at 500 cm3/s, as shown in Fig. 4. The viscous pressure drop over two selected airway path-lines is shown in Fig. 5. The local pressure drops differ because of the anatomical asymmetries between right and left lungs along the two selected path-lines. The same trend was observed for a CFD simulation of 7 generations (see Fig. 8 in van Ertbruggen et al., 2005 and in the formula for airway pressure drop of Pedley, 1977). Flow partitioning in the lung is due to anatomical asymmetry of the right and left lungs and influences aerosol deposition. Fig. 6 depicts the percentage change in mass flow rate in the first three generations of the right and left lungs. Generally, the lateral daughter branch attracts much more fluid flow than does the medial branch, as was also observed by Liu et al. (2003).

4.2 Airflow Velocity FieldsAirflows in the upper airways can periodically exhibit the characteristics of laminar, transitional, or fully developed turbulent flow during inhalation due to a succession of conduit cross-sectional area changes occurring between the trachea and the respiratory surfaces of the lungs.

A uniform velocity profile was imposed computationally at the inlet to the trachea with an average value of 2.896 m/s, corresponding to a steady-state inhalatory flow rate of 28.3L/min and a tracheal diameter of 14.4 mm. At normal breathing rates of 28.3L/min, the airflow in the larynx and in the tracheal section of the upper airways is slightly turbulent. Distal to the first few main upper generations, the airflow quickly transitions to laminar flow. The effects of the turbulent laryngeal jet on the gas flow field extend throughout the length of the larynx, causing recirculatory flows and turbulent eddies even beyond the tracheal bifurcations (Gemci et al., 2002, 2003). Stapleton et al. (2000) and Allen et al. (2004) have shown that the choice of a reasonable turbulence model is critical. For this reason we performed our airflow simulation with the LES turbulence model, even though it is computationally more expensive than RANS-based turbulence models.

The velocity fields in more anatomically realistic tracheobronchial airway models are more complex than those in the Weibel-based models. Airway curvature and complex anatomical shapes at junctions contribute to those complex flow patterns as noted in Novak et al. (2003). Figs. 7 and 8 depict the secondary flow patterns in the main left and right bronchi. In previous symmetric and asymmetric Weibel-based tracheobronchial airway models (Pedley et al., 1971; Guan and Martonen, 2000; Calay et al., 2002; Zhang et al., 2002; Liu et al., 2003; Nowak et al., 2003; Shi et al., 2004), the secondary flow in the in-plane cross-section of the bronchus was directed from the center of the inner wall towards the center of the outer wall; it then divided equally into a double vortex, also known as Dean flow. Dean flow was also observed here.

5. Conclusions

CFD evaluations of the quasi-steady incompressible airflow in the human respiratory tract have been presented here with a view to their subsequent extension to the study of the uptake of particulate matter by the tracheobronchial airways in flows varying periodically in time.

We have adapted the anatomical graph data of Schmidt et al. (2004), based on measurements derived from High-Resolution Computer Tomography (HRCT). Assumed pressure drops across successive bifurcations agree reasonably well with earlier published results using more limited anatomical models. Laminar-turbulent flow transitions have been accounted for with a LES model. The computed flow fields also reflect the development of secondary swirling flows resulting from airway curvature.

This digital reference model of the human pulmonary tree is proposed as a first step toward the development of detailed time-varying deposition models with a view to the quantitative design and assessment of therapeutically targeted drug delivery systems via the respiratory tract and to the regional evaluation of the ingestion of particulate pollutants into the respiratory tract.

Acknowledgements

We express our thanks to Dr. Andres Kriete at the Coriell Institute for Medical Research, Camden, NJ for providing us with the abstracted topological graph data. The access to the supercomputers at NSCEE of UNLV and Pittsburgh Supercomputing Center (Grant #CTS010015P) is also gratefully acknowledged.

Fig. 5 Variation of the pressure drop along two selected path lines for the flow rate of 28.3 L/min.

Fig. 6 Percentage variations of the mass flow rate portioning in the first main bifurcation and its distal daughter branches. Total mass flow rate is 0.6125 g/s (corresponding 500 cm3/s or 30 L/min)

Fig. 7 The velocity vector (right frame) and velocity magnitude (left frame) are shown at an in-plane cross-section at the left main bronchus in the top frame. Pictured here is an oblique plane view of the first left bronchus in a direction perpendicular to the flow.

Fig. 8 The velocity vector (right frame) and velocity magnitude (left frame) are shown at an in-plane cross-section at the right main bronchus in the top frame. Pictured here is an oblique plane view of the first right bronchus in a direction perpendicular to the flow.

ABOUT THE AUTHOR: Richard CollinsDr. Richard Collins is a former professor of medicine, chair professor of biomedical engineering, professor of mechanical and aerospace engineering.

His current activities include research and development and litigation support for law firms nationwide, with demonstrated accuracy in research and scientific analysis in mechanical and biomedical engineering and regulatory assessments. He is qualified as an expert witness in civil and criminal courts and is internationally recognized for his engineering and medical expertise.

While every effort has been made to ensure the accuracy of this publication, it is not intended to provide legal advice as individual situations will differ and should be discussed with an expert and/or lawyer.For specific technical or legal advice on the information provided and related topics, please contact the author.