## Abstract

We have studied gas flow and particle deposition in a realistic three-dimensional (3D) model of the bronchial tree, extending from the trachea to the segmental bronchi (7th airway generation for the most distal ones) using computational fluid dynamics. The model is based on the morphometrical data of Horsfield et al. (Horsfield K, Dart G, Olson DE, Filley GF, and Cumming G. *J Appl Physiol* 31: 207–217, 1971) and on bronchoscopic and computerized tomography images, which give the spatial 3D orientation of the curved ducts. It incorporates realistic angles of successive branching planes. Steady inspiratory flow varying between 50 and 500 cm^{3}/s was simulated, as well as deposition of spherical aerosol particles (1–7 μm diameter, 1 g/cm^{3} density). Flow simulations indicated nonfully developed flows in the branches due to their relative short lengths. Velocity flow profiles in the segmental bronchi, taken one diameter downstream of the bifurcation, were distorted compared with the flow in a simple curved tube, and wide patterns of secondary flow fields were observed. Both were due to the asymmetrical 3D configuration of the bifurcating network. Viscous pressure drop in the model was compared with results obtained by Pedley et al. (Pedley TJ, Schroter RC, and Sudlow MF. *Respir Physiol* 9: 387–405, 1970), which are shown to be a good first approximation. Particle deposition increased with particle size and was minimal for ∼200 cm^{3}/s inspiratory flow, but it was highly heterogeneous for branches of the same generation.

- lung model
- aerosol
- asymmetry
- viscous pressure drop
- heterogeneous deposition

to improve the delivery of aerosolized drugs to their intended site of action in the lung and therefore increase their efficiency, it is necessary first to understand the mechanisms of aerosol transport in the pulmonary airways. However, it is experimentally impossible to measure directly the deposition of aerosols in the human lung, so computational fluid dynamics (CFD) brings a potential solution to this problem. Computing numerical solutions of the flow-governing equations (Navier-Stokes equations) and simulating the particle motion are made possible through marketed CFD software packages. Because computers are becoming increasingly more powerful, application of these numerical methods to the simulation of flow in the airways is now becoming feasible. Hence, flow characteristics (such as velocity and pressure) and particle deposition patterns that cannot be obtained experimentally can be simulated numerically. Comparison of the simulation results with experimental results whenever possible enables an evaluation of the simulation accuracy. This constitutes a major advance for the understanding of particle motion and deposition in the lungs.

Morphometric measurements of the human bronchial airways and models describing their structure have been reported by a number of investigators (11, 16, 36, 38). Some were used to mathematically predict aerosol deposition in the human respiratory tract (2, 5a, 6, 12, 21, 30, 38). Simple geometric configurations, such as a single bifurcation or a double or triple cascade of bifurcations, led to CFD simulations and gave a first insight into flow profiles and particle deposition in the upper airways (1, 3–5, 7, 9, 12, 14, 15, 22–24, 37, 41, 42).

An important part of this work involved the building of a realistic model of the human airways. All the previously proposed models (11, 16, 36, 38) lacked information concerning the angles formed by successive branching planes (i.e., the plane containing the assumed coplanar parent and daughter branches), also known as azimuthal angles Φ. The most complete study of a lung cast, taking into account the asymmetry, published by Horsfield et al. (16), gave only the length, diameter, and branching angles between the assumed coplanar parent and daughter branches. Some authors fixed the value of the azimuthal angle to be 90° (9, 20). This missing information related to the three-dimensional (3D) disposition of branches in the chest had to be found to build a realistic model of the bronchial tree.

Few approaches have been used to build geometries of a complete lung; these were based on computerized tomography (CT) scanner imaging (28, 29), or on iterative mathematical algorithms designed to appropriately fill some portion of space (17, 20). Basic CFD simulations were conducted on Perzl's simplest model (28, 29), whereas the second approach mainly led to visual improvements of the 3D lung models, leaving a gap for obtaining a computational model.

We propose a realistic 3D asymmetric computational model of the lungs (34, 35), extending from the trachea up to the segmental bronchi, based on the morphometrical data of Horsfield et al. (16) and on bronchoscopic, X-ray, and CT scanner images (32). CFD simulations were performed on this model to obtain the main flow characteristics (velocity profiles, secondary velocity patterns, and viscous pressure drops) and the particle deposition efficiency. The impact of asymmetry on the simulations results was highlighted through comparison with the available published results, all based on symmetrical models.

The commercial software package FINE/Turbo (version 6.1) of Numeca International was used to create the model and simulate flow and particle transport. This package allowed the achievement of a complete CFD simulation in four steps: *1*) creation of the geometric model delimiting the portion of space where flow will be studied, *2*) discretization of this space with a mesh, i.e., transformation of the geometric model into a so-called computational model, *3*) numerical resolution of the flow-governing equations with particles if needed, and *4*) analysis and visualization of the results.

The model-building history (*steps 1* and *2*) is explained in the methods, along with the simulation process (*step 3*). A short description of the solver is given.

The results and discussion describe flow-simulation results and focus on spherical particle deposition inside the model (*step 4*). Inspiratory flows ranging from 50 to 500 cm^{3}/s have been simulated. The particle sizes were chosen, representative of standard pharmaceutical aerosols (diameter varying between 1 and 7 μm). Particle deposition was analyzed in each branching generation, and factors influencing their deposition (inspiratory flow rate and particle diameter) were studied.

## Glossary

μ_{air}

*D*- Tube diameter
*d*_{p}- Particle diameter
- Φ
- Azimuthal angle of successive branching planes
*g*- Gravitational acceleration (9.81 m/s
^{2}) - M
_{p} - Particle mass
*L*- Tube length
*L*_{d}- Entrance length given by
*Eq. 1* - ν
_{air} - Air kinematic viscosity (1.57 × 10
^{−5}m^{2}/s) - μ
_{air} - Air dynamic viscosity (ρ
_{air}× ν_{air}) - μ
_{f} - Fluid dynamic viscosity
- μ
_{p} - Particle dynamic viscosity
- ρ
_{air} - Air density (1.2 kg/m
^{3}) - ρ
_{p} - Particle density (1,000 kg/m
^{3}) - Re
- Reynolds number
- Stk
- Stokes' number
*u*- Fluid velocity
*U*_{ref}- Maximum velocity encountered in the flow
*U*_{o}- Mean velocity
*U⃗*- Fluid velocity vector
*V⃗*- Particle velocity vector

## METHODS

#### Buildup of the model.

The proposed method to build a realistic model of the airway tree extending from the trachea to the segmental bronchi was based on the following facts. *1*) The morphometrical data of Horsfield et al. (16) were accurately measured on a lung cast up to the segmental bronchi, reflecting the asymmetry of a real lung for branch length, diameter, and branching angles. However, to the best of our knowledge, they have never been completely used. *2*) Bronchoscopic images have not been exploited so far, even though they can provide a general idea of the bronchi orientation. Indeed, the camera delivers an insight view of the bronchi, thus a quite good visualization of the rotation angles of successive branches. These data are easily available up to the segmental bronchi (32) and will be used to obtain the azimuthal angle Φ of successive branching planes. *3*) CT scanner and X-ray images were used to check the general orientation of the bronchi supplying pulmonary lobes and segments (32).

The computing resources available determined the number of mesh nodes (i.e., the simulation accuracy) used to discretize the model. The computing time needed to perform a simulation increased drastically with the number of nodes. So a compromise had to be made between the geometrical extent we wanted to model, the results accuracy, and the computing resources (time and memory) we had.

Bronchoscopic images gave the rotation angles of successive branching planes (comprised between 0° and 145°), whereas CT scanner and X-ray images served as references to check the general 3D orientation of the branches.

We used the “model 2” data given by Horsfield et al. (16) relative to length, diameter, branching angles, and curvature ratio for all of the 35 branches supplying the 18 pulmonary segments of our 3D model. Branches were modeled with tubes of constant circular section. Because the information relative to the right lower lateral basal segment was not given, we assumed the same data values for this branch as for the posterior basal segment.

The division between a parent branch and the two daughter branches was not abrupt, i.e., a smooth transitional zone existed between the parent branch and its two daughters. We used the transitional zone model proposed by Horsfield et al. (16) (Fig. 1*A*); it constituted approximately the distal 20% of the parent branch length, and its shape varied continuously to obtain two distinct daughter branches with different diameters. The transitional zone section was circular at the beginning, then evolved to form an ellipse, after which an “8-shaped” curve was progressively reached to finally form two separate circles at the end of the zone, corresponding to the beginning of the daughter branches. The walls of the two daughter branches did not make a sudden angled departure from the corresponding walls of the parent tube, but they rather headed toward the direction defined by the branching angle.

We grouped our branches in basic units called “bifurcations”: a bifurcation was made with the parent branch transitional zone and the two daughter branches without their transitional zone, with all three lying in the same branching plane. Assembling the basic units generated a 3D airway model. A total of 17 different bifurcations have been designed to cope with the asymmetry of the bronchial tree.

The design of every bifurcation unit was done by creating different curves representing cross section through the transitional zone, similar to those described before (circle, ellipse, and 8-shaped curves) (Fig. 1*A*). A surface lying on these curves was created and connected to the curved cylinders representing the daughter branches, so the geometric surfaces of the bifurcation were obtained. The following design rules applied, ensuring the obtaining of a smooth transitional zone. *1*) The transitional zone length was maintained as close as possible to 20% of the parent branch length. If the parent branch was shorter than 2 cm, the relative length of the transitional zone was increased; otherwise, the transitional zone cross section would have varied too abruptly, which did not visually correspond to real pulmonary bifurcation models. *2*) The bell-mouthed shape of the transitional zone outer walls was made by the heading to 8-shaped curves that were smoothly connected to the circular outlet of the parent branch. *3*) Curved circular tubes represented daughter branches until the right bifurcation angle was reached. A straight tube was added afterwards if necessary to obtain the right daughter branch length. Finally, *4*) the carinal ridge was sharp; Comer et al. (9) studied the effects of rounding this carinal ridge and found that only flow in the immediate vicinity of the carinal ridge was affected.

Once the geometry of the 17 bifurcations was created, the inner space of every bifurcation was discretized with a structured butterfly mesh (Fig. 1, *B* and *C*); ∼15 mesh blocks were necessary to mesh each bifurcation model. Meshing different zones of the bifurcations with separate blocks allows a better fitting of the mesh into the geometry and improves the control of the mesh quality parameters; this is known as the multiblock meshing approach. Careful control of the mesh quality parameters (orthogonality, aspect ratio, and expansion ratio) was performed in each block: *1*) the orthogonality of 99.99% of the cells was >36°, *2*) the maximum cell aspect ratio was 34, and *3*) the cells' expansion ratio was <3. Orthogonality measures the smallest value of the angles formed by the cell walls, the aspect ratio measures the cell length over height ratio, and the expansion ratio is a measure of the variation in the size of adjacent cells. To obtain a good mesh quality, cells must be as cubic as possible. This implies a range of values for the mesh quality parameters. The number of nodes was set sufficient to ensure mesh-independent results. Simulations performed on a single bifurcation, using the same methodology as Calay et al. (7), have demonstrated that a mean value of ∼100,000 to 150,000 nodes per bifurcation was necessary to approach a mesh-independent solution (Fig. 1, *B* and *C*). Indeed, comparisons on a single bifurcation were made with the actual number of nodes per section multiplied by 8. Computational resources increased drastically, and increasing the number of mesh nodes gave <2% maximum difference in the velocity profiles. This indicated that the chosen number of mesh nodes per section was enough to approach mesh-independent results in a reasonable computing time. Once the 17 bifurcations were meshed, they were assembled in three dimensions using successive translations and rotations to finally obtain a model of the bronchial tree extending from the trachea to the segmental bronchi (Figs. 2, *A* and *B*). Different colors were assigned for the branch generations of the model (light gray for *generations 0*, *2*, *4*, and *6*; dark gray for *generations 1*, *3*, *5*, and *7*): the model of airways extends up to the seventh generation (in the right lower lobe). Table 1 lists all distal branches with names, abbreviations, and generation numbers.

The total number of computational nodes was 2,300,000, for 425 butterfly blocks. The mesh structure is shown on Fig. 3 for a coarser grid (300,000 nodes).

#### Choice of the boundary conditions.

During breathing, air enters the lungs because of a negative difference of pressure between the alveoli and the atmospheric air. So, ideally, an adequate pressure level should have to be imposed at each model exit to mimic the respiration process. These data do not exist. However, it is possible to obtain indirectly physical boundary conditions by imposing mass flow at every exit: we assumed that the lung is uniformly ventilated, so the flow in each branch is proportional to the volume it supplies, and used the data of *model 2* of Horsfield et al. (16).

We simulated flow in our model while imposing the same pressure at each outlet: we checked whether the flows obtained from the solution matched ventilation reasonably. Such an approach was taken by Andrade et al. (1) for two-dimensional simulations in a four-generation cascade of symmetric bifurcation; they found that fluid inertia generated uneven flows, but this flow pattern was speculated to match the required homogeneity of ventilation. In the present case, with a 3D asymmetric tree, such an approach was not valid because the flows obtained in the exit branches were totally different from those necessary to obtain a uniform ventilation of the lung. Therefore, the same pressures for the boundary conditions were not used.

Airway expansion of the model during breathing was neglected. For a given respiratory flow, a constant velocity profile was imposed at the inlet of the model, and the flow percentages computed by Horsfield et al. (16) were used to impose mass flow in all of the model exits, except for one where the static pressure level was set, avoiding redundancy in the boundary conditions.

Because the exit branches were relatively short, flow was never fully developed in the exit sections. To avoid any influence of the boundary conditions on the computed solution, each branch exit was artificially lengthened with a straight tube, whose length *L*_{d} ensured that the exiting flow was almost developed and was given by (1) where Re is the flow Reynolds number and *D* is the tube diameter.

Because of the curved shape of the branches and the asymmetry of the tree, secondary flows did persist in some diameters downward of the bifurcating zone. Imposing a constant mass flow or a constant pressure level too close to the bifurcation zone can therefore have unwanted effects on the solution and lead to simulation convergence problems. This is why outlet boundary conditions were imposed sufficiently far from the physical exits to ensure that the flow was unaffected by the boundary conditions. The artificially lengthened zones only served during flow simulations; they were not considered for the particle deposition study.

Several flow-simulation tests were performed on the tree for different combinations of boundary conditions; mass flow was imposed at every exit, except for one where the pressure level was set, but the exit position of the pressure level varied for each combination. No significant differences were observed concerning the convergence history, flow patterns, or deposition results. We therefore choose to impose a constant static pressure level at the exit of the left lower anterior medial basal segmental bronchi (exit LS8).

The flow velocity was zero at the air-solid boundary inside the model.

#### Flow simulation process.

The solver Euranus of Fine 6.1 (13) uses the finite volume method to solve the Reynolds Averaged Navier-Stokes equations. A central explicit second-order scheme was used for spatial discretization, whereas a fourth-order explicit Runge-Kutta scheme was used for time discretization. At the highest inspiratory flow considered in this study, i.e., 500 cm^{3}/s, the maximum Reynolds number encountered in the model was 2,500. The flow regime can therefore be considered as laminar or lightly transitional, but certainly not turbulent (the critical Reynolds number averages 8,000 in that case). Indeed, we performed both laminar and turbulent simulations (with the models available in the software Fine/Turbo, i.e., *k*-ε model of turbulence) to determine which model (laminar or turbulent) was best suited in our case. We compared the flow profiles and particle deposition percentages obtained in both cases and found no significant differences. Using a turbulence model is therefore not necessary in the present case because it did not bring additional information. Previously, Zhang and Kleinstreuer (40) have used low-Reynolds *k*-ε and *k*-ω turbulence models to simulate flow in constricted tubes. They concluded that the low-Reynolds *k*-ω turbulence model was better at reproducing the behavior of laminar and transitional flow. However, even though the low-Reynolds version of the *k*-ε and *k*-ω models can reproduce a transition from laminar to turbulent regime, the physics of transition is not really modeled in these turbulence models. Therefore, the transition prediction cannot be considered as reliable. In our simulations, it is considered that the flow in the lung branches is mostly laminar because the Reynolds number is globally below the transitional value.

Additional methods such as preconditioning and multigrid were used to accelerate the convergence of the simulations: we imposed a residual decrease of four orders of magnitude on the coarse grids, and six orders of magnitude on the finest grid. The residual decrease was monitored in each block to ensure adequate convergence. The mass flow error was stabilized below 0.5%. Depending on the inspiratory flow, 400–700 iterations were needed to stabilize the residuals on the finest grid.

For the particle-tracking process, a Lagrangian approach was used. The physical particle cloud was represented by a finite number of computational parcels. Interactions between particles were neglected; this remains valid as long as the particle flow is dilute (i.e., the volume fraction of the particles is <10%, which is the case here). One-way coupling simulations were performed, i.e., particles had no influence on the flow.

Gravity and drag force (based on Stokes' law) were the only forces acting on a particle (*Eq. 2*). Brownian motion was neglected compared with other intrinsic motions because particle size was 1 μm minimum. The forces exerted on a particle were therefore (2) where F⃗ is the resulting force acting on a particle expressed as the sum of the drag force F⃗_{drag} and the gravitational force F⃗_{grav}. The *d*_{p} is the particle diameter, *V⃗* is the particle velocity vector, *U⃗* is the fluid velocity vector, μ_{air} is the air dynamic viscosity, M_{p} is the particle's mass, and *g* is the gravitational acceleration.

Both uniform and parabolic distributions of particle entering the model with the same velocity as the fluid were simulated; 50 particle probes were injected in each model inlet cell, and a total of 100,000 particles were thus tracked in the model. The number of physical particles associated with each particle probe was recorded in a file, along with its final status (i.e., exited or stuck on a wall), so deposition values could be computed. A particle hitting a wall sticks on it; this is a reasonable hypothesis because airway walls are covered with mucus.

Numerical simulations were performed on Compaq Digital Alpha Servers, GS140 (CPU: 700 MHz) and GS160 (CPU: 1 GHz), running under Compaq Tru64 Unix v.5.1A. In these conditions, each flow simulation took an average of 2 days on one processor and 2 Gb of random access memory. Each particle-tracking simulation, which constitutes mainly postprocessing, took ∼3 h.

The software was validated from cases for which the solution is analytically known (i.e., Poiseuille flow) and from previously published data (7–9). For example, we compared the velocity profiles simulated in a model of the trachea and left and right main bronchus with those of Calay et al. (7): the dimensions of the model and the Reynolds number in the trachea (Re = 1,592) were identical. The velocity profiles were taken in planes perpendicular to the branches, every half-diameter from the bifurcation, in two directions: one perpendicular to the bifurcation plane (anteroposterior), and the other one in the same plane as the bifurcation. Results obtained by simulation in the left main bronchus were compared in Fig. 4 with the results of Calay et al.: the flow profiles were very similar.

In Figs. 5*A* and 6*A*, velocity plots presented a peak close to the inner wall of the bifurcation. However, in most of the cases, a second peak is present; this is due to the more complex flow configurations created by successively rotated bifurcations. In Figs. 5*B* and 6*B*, velocity plots did not necessarily present the characteristic M-shape observed in a curved tube, for the same reasons.

Multiple different secondary flow patterns appeared because of the branches' shortness and nonplanar disposition. One, two, three, or four secondary velocity vortexes could be observed close to the bifurcation zone (Fig. 7). Small vortexes tended to disappear as the flow moved away from the bifurcating zone.

#### Viscous pressure drop study.

Pedley et al. (25–27) published an empirical formula (*Eq. 3*) giving the viscous pressure drop ΔP_{v} in bifurcating branches, based solely on velocity profiles measured in the circular section of the branches. Their study was based on a symmetrical bifurcation model. (3) where *L* is the tube length, *U*_{o} is the mean velocity, μ_{f} is the fluid dynamic viscosity, and the constant γ equals 0.327.

We computed the viscous pressure drop in all of the branches of our model based on simulation data; indeed, the viscous pressure drop equals the total pressure drop [i.e., energy losses in the extended Bernoulli equation (39)]. Because total pressure (P_{tot}) was not uniform in any sections (S), we computed mean mass averaged total pressure values as follows: (4)

For the 500 cm^{3}/s inspiratory flow simulation, we compared the simulated results with the results obtained with the empirical equation of Pedley (*Eq. 3*) applied to our model and to Weibel's symmetrical model (Fig. 8*A*). We did not neglect the effect of the bifurcation region while computing the simulated total pressure drop: the entire branch length was considered (i.e., tubular section plus bifurcation region). Compared with the results obtained by simulation, Pedley's formula overestimated viscous pressure drop in most of the cases. The constant factor γ was on average too large, as noticed by Comer et al. (9) and Pedley et al. (27): the pressure drop obtained with Pedley's formula was on average 33% higher than that of the simulated formula. We obtained very disparate pressure drop values for branches of the same generation because of the 3D-configuration, different diameter, angle, and mass flow in these branches (large standard deviation in Fig. 8*A*). The application of Pedley's formula on Weibel's symmetrical model data generated differences compared with the simulations. The calculated viscous pressure drop in Weibel's model was on average 29% higher than in Horsfield's model and 61% higher than the simulated results. We deduced new values of the factor γ on the basis of the simulation results for each generation. The new value of γ is obtained by equaling Pedley's formula with the simulated results, the only unknown being the γ factor. Mean value of γ was computed for each generation, and the results were plotted in Fig. 8*B*.

The total laminar flow resistance of the model was theoretically computed on the basis of an idealized tree made of straight tubes: the laminar flow resistance is then 4,363 Pa·s·m^{−3}. The resistance of the proposed model obtained through simulation (inspiratory flow rate of 500 cm^{3}/s) was 3.5 times higher, at 15,322 Pa·s·m^{−3}.

#### Particle deposition study.

Deposition of spherical particles (density = 1 g/cm^{3}, diameter ranging between 1 and 7 μm) was simulated in the model for inspiratory mass flows ranging between 50 cm^{3}/s (tracheal Re = 250) and 500 cm^{3}/s (tracheal Re = 2,500). The shape of the entrance particle profile was uniform or parabolic, so deposition results could be compared for both cases. This corresponded to a deposition study in the shallow region of the lung (i.e., approximately a lung depth between 50 ml and 100 ml).

Deposition increased with particle size and was generally higher with a parabolic profile of particles entering the model (Fig. 9*A*). Deposition varied with mass flow and for particles larger than 1 μm had a minimum for ∼200 cm^{3}/s (Fig. 9*A*). For 1-μm particles, deposition simply increased with mass flow. In our simulations, deposition was important at very low and high inspiratory mass flows. In Fig. 9*B*, we compared our deposition results with the experimental results of Kim et al. (18) obtained for a volumetric lung depth of 50–100 ml (this corresponded to the average lung depth covered by the model). Our results were obtained at inspiration, whereas Kim's results were obtained after a respiration cycle.

Deposition in branches of the same generation was compared (see Fig. 2 for visual representation of the model generations). The deposition was highly heterogeneous, as can be seen in Fig. 10, where deposition percentage (raw deposition data for each branches as well as mean deposition percentage ± SD) in *generations 1* to *7* was plotted for particles of 2, 4, and 6 μm. For *generation 1*, simulation results for an inspiratory flow of 500 cm^{3}/s showed that deposition in the left main bronchus was three to five times more important than in the right main bronchus (2.08 vs. 0.65% for 7-μm particles).

Complete heterogeneous deposition percentages were observed, the standard deviation being very large compared with the mean value. This occurred at every mass flow for every particle diameter.

The deposition values were highly disparate: deposition varied by a factor of 10 in branches of the same generation. Higher deposition values were found in branches whose branching angles were important. For example, deposition in the right lower apical segment (RS6, 4th generation branch), whose branching angle is close to 60°, was more than 10 times bigger than any other deposition values in the same order generation branches.

## RESULTS

#### Flow study.

We present the simulated velocity profiles for all the segmental bronchi (inspiratory flow of 500 cm^{3}/s), taken in a plane one diameter downstream of every bifurcation. These profiles were grouped in Fig. 5 for the right segmental bronchi and Fig. 6 for the left segmental bronchi. Normalized velocity (vertical axis) was plotted against normalized arc length (horizontal axis), so the general form of the velocity profiles could be compared, independently from the velocity magnitude. Velocity plots in two directions are presented: Figs. 5*A* and 6*A* present the velocity plot along the branch symmetry plane, oriented from the outer side to the inner side of the branch, whereas Figs. 5*B* and 6*B* present the velocity plot along the branch normal plane (also known as anteroposterior axis).

## DISCUSSION

#### Flow study.

The flow velocity profiles in a simple curved tube (Dean flow) present two major features (27). Note that the exterior (interior) of the bend corresponds to the inside (outside) of a bifurcation. *1*) Two secondary motion vortexes develop. These are apparent in planes perpendicular to the tube, the vortexes being directed toward the exterior of the bend and returning to the interior near the walls. *2*) Velocity profiles are distorted from Poiseuille flow and present a peak near the exterior side of the bend. Velocity plots in the tube symmetry plane present a peak near the outside wall, whereas plots along a perpendicular axis are M-shaped. Indeed, the velocity peak near the outside wall extends on the anterior and posterior walls, creating double-peak velocity profiles in the anteroposterior direction. Analogies can be found for the flow in the model's branches, because bifurcations are primarily made of curved ducts. However, the complete asymmetry of the model induces very asymmetric flow patterns (Figs. 5 and 6), much more complex than the flow in a curved pipe described previously (27).

Chang (8) experimentally measured velocity profiles in a system of planar asymmetric curved tubes. These velocity profiles gave a first idea of what possibly happened in successive curved tubes. However, the number of measured points for each profile remained limited, as well as the number of locations where profiles were measured. The planar configuration of the branches prevented obtainment of very distorted profiles; maximum velocity peaks were therefore always located in the symmetry plane of the model.

In our case, velocity plots presented a peak close to the inner wall of the bifurcation (Figs. 5*A* and 6*A*). However, in most of the cases, a second peak was present; this was due to the more complex flow configurations created by successively rotated bifurcations. In the anteroposterior planes, velocity plots did not necessarily present a characteristic M-shape, for the same reasons. Similar conclusions applied at lower mass flow.

In case of successive non-coplanar bifurcations, the velocity peaks were deflected at the bifurcation zone, creating peaks that were not initially close to the bifurcation inner walls (i.e., secondary peaks). M-shaped velocity plots along the anteroposterior direction were therefore perturbed. The ideal situation present in the curved tube could be obtained in sufficiently long tubes. However, branches were generally too short for this phenomenon to be observed.

Multiple secondary flow vortexes were observed (Fig. 7); their repartition was not necessarily symmetrical about the bifurcation plane. This was different from the ideal situation present in a curved tube. As stated by Comer et al. (9), secondary flow fields depended on the upstream flow configuration and on the bifurcation zone shape. Because many branches were short compared with their diameter, very complex flow field patterns arrived on the bifurcation ridge, generating again more complex flow fields in the daughter branches. Secondary velocities are extremely important because they govern particle deposition, causing deposition in zones where the secondary velocities (i.e., flow transverse motion) are important.

#### Viscous pressure drop study.

Pedley et al. (25–27) published an empirical formula (*Eq. 3*) giving the viscous pressure drop ΔP_{v} in bifurcating branches, based solely on velocity profiles measured in the circular section of the branches. Their study was based on a symmetrical bifurcation model. The value of the constant γ is 0.327 in their study. They neglected the viscous pressure drop in the bifurcation region, arguing that most of the viscous dissipation occurs in the thin boundary layer that develops from the flow divider (27).

While comparing the results predicted with the formula and the simulated results (Fig. 8*A*), we found that the constant factor γ was too large, as already noticed by Comer et al. (9). We analyzed the pressure drop for lower mass flow and the same conclusions applied: the constant γ was overestimated. We found that Pedley's formula overestimated viscous pressure drop in the circular sections of the ducts compared with simulated results; this was consistent with the results of Comer et al. This overestimation was reduced by considering the viscous pressure drop simulated in the bifurcation zone; in that case, the gap between the simulated results and those obtained with Pedley's formula decreased. Again, this was consistent with the results of Comer et al. However, the agreement Comer et al. found between their simulated results and the predicted results was on average better than what we found. But their model was made of three planar symmetrical bifurcations, which was closer to the geometrical hypothesis used for Pedley's empirical formula than the model we used, which was asymmetrical and nonplanar.

We computed new values of γ, which varied between 0.162 and 0.566 (Fig. 8*B*). Generally, the values of γ obtained with the simulation results were smaller than the value given by Pedley et al. (26). Comer et al. (9) computed a new value of γ for each generation of their three bifurcation models, on the basis of the simulations results, and found values around 0.2, which were similar to our results. In our model, the computed value was higher than 0.327 only for *generations 6* and *7*; however, the values of γ obtained by simulation in these generations might not be statistically representative of the mean value of γ for the generation because the model contained only two branches of *generations 6* and *7*.

The empirical formulas showed their limitations in the case of complex geometries because they were based on very simple and symmetric bifurcations (branching angle of 35° for each daughter branch). So simulations on asymmetric models allowed obtainment of more accurate information, but *Eq. 3* remained a relatively good first approximation for such a complex structure.

The total laminar resistance of the model was 3.5 times smaller than the resistance obtained through simulation, for an inspiratory flow rate of 500 cm^{3}/s. Indeed, the real tree geometry was made of curved bifurcating tubes, oriented in 3D. Therefore, the flow resistance of the model was higher than the laminar flow resistance obtained while considering a tree made of idealized straight tubes of constant section and ignoring the resistance at the bifurcation. This is consistent with the fact that the rate of dissipation of energy by viscosity is minimum in the ideal case of laminar flow in a straight tube.

#### Particle deposition study.

We simulated particle deposition in our model during inspiration. The inspiratory mass flow varied between 50 and 500 cm^{3}/s. Because the cartilage rings present in the central airways were not modeled, their potential effect on deposition was not considered; this constitutes a limitation of the actual model. A constant and a parabolic profile of particles were injected at the inlet of the model. Parabolic profile generated higher deposition percentage than the constant profile: the deposition value increased by a factor varying between 3 and 18% (Fig. 9*A*). A parabolic profile seemed closer to reality than a constant profile, because the particles had to pass first through the glottis opening before reaching the trachea. If a constant profile of particles passes through a narrowing, their repartition after that glottis narrowing is not constant anymore because the particles close to the walls will hit the glottis walls first. In their study of deposition in alveolated ducts, Tsuda et al. (33) found that the injection of a constant profile of particles generated higher deposition than a parabolic one. However, in the zone of the lung they considered, gravitational sedimentation had a strong influence on the particles' trajectories because velocities were small, whereas in the region we studied, impaction was predominant. In the upper part of the lung that we considered, a parabolic profile generated higher deposition values because more particles were present in the core of the flow than in constant profile, so more particles were trapped around the bifurcation zone of the first generations.

Deposition increased with particle size (Fig. 9*A*). This was consistent with the experimental data of Kim et al. (18), in which regional deposition of inhaled particles was assessed through serial bolus delivery method. This is logical because the mechanisms of sedimentation and inertial deposition are proportional to the square of the particle diameter. Indeed, gravitational sedimentation depends on the sedimentation velocity *V*_{s} (5) and inertial deposition depends on the particle Stokes' number Stk (6) In *Eqs. 5* and *6*, ρ_{p} is the particle density, *d*_{p} is the particle diameter, μ_{p} is the particle dynamic viscosity, and *u* is the fluid velocity.

Deposition varied with mass flow and had a minimum for ∼200 cm^{3}/s (Fig. 9*A*). Deposition data of Kim et al. (18) were obtained for a complete respiration cycle, whereas we computed deposition at inspiration only. So, the experimental deposition values of Kim et al. were therefore higher than our simulated results. If we suppose that the deposition efficiency at inspiration and expiration is identical, the simulated results should be about two times smaller than the experimental results. This was not exactly the case (Fig. 9*B*): the experimental results obtained for a complete respiratory cycle were *1*) at low mass flow (150 cm^{3}/s) more than twice those obtained with an inspiratory simulation, and *2*) at high mass flows (250 and 500 cm^{2}/s) less than twice those obtained with an inspiratory simulation. Indeed, at expiration, the role of inertia is less important than at inspiration because the bifurcating region does not face the particles' trajectories. Expiratory deposition values at high mass flow are therefore less important than in the inspiratory case. This is consistent with the findings of Kim et al. (19). However, at lower mass flows, the experimental deposition was more important; this could be explained by the fact that sedimentation probably played a key role during expiration at low mass flows, acting in opposition with the particle velocity for most of the branches in the model, thus increasing the deposition.

In our simulations, deposition was important at very low and high inspiratory mass flows. Indeed, at very low mass flows, air velocity was small, so sedimentation was effective for particles close to the walls, but particle inertia was low because their velocity was small (small Stk). The sedimentation effect was especially marked for large particles. At high flow rates, the inertia of the particles increased (because fluid velocity increased) and impeached them to strictly follow the streamlines. Their sedimentation velocity, which was independent of the mass flow, became negligible compared with the particle inertial motion. So deposition by impaction increased with mass flow and particle diameter. Darquenne at al. (10) studied deposition of particles with normal gravity and in microgravity. For an inspiratory mass flow of 450 cm^{3}/s, they found almost no differences of deposition for small particles (1 and 2 μm) at shallow volumes of penetration (150 ml). We compared deposition of small particles (1 and 2 μm) in our model at high inspiratory mass flow (500 cm^{3}/s) with and without gravity. The results were very similar: 1.25 and 1.38% of deposition in microgravity and 1.21 and 1.45% of deposition with gravity for 1-μm and 2-μm particles, respectively. This showed that gravitational sedimentation became negligible compared with the mechanism of inertial impaction for high inspiratory mass flow.

To summarize, when inspiratory flow increased, sedimentation velocity became negligible but particles inertia increased, leading to increasing deposition percentage. For a flow rate of ∼200 cm^{3}/s, simulations indicated that the sum of these two effects had a minimum: particles had enough speed so deposition occurred less by sedimentation but they had not enough inertia to deviate much from their trajectories either.

Deposition in the left main bronchus was three to five times more important than in the right main bronchus (2.08 vs. 0.65% for 7-μm particles) because the left main bronchus was 2.5 times longer than the right bronchus and also more curved because the branching angle of the left main bronchus was twice that of the right bronchus (73 vs. 35°). Inertial impaction and gravitational sedimentation occurred therefore more frequently in the left main bronchus.

The deposition values were highly disparate in branches of the same generation (Fig. 10). This was no surprise because the model was fully asymmetric: length, branching angles, diameters, and mass flows were very different for branches of the same generation. Because deposition percentage depends also on the size of the branches, the longer the branch, the higher the probability a particle would deposit. This effect cannot be reproduced either by simulations on symmetrical models of the lungs or by theoretical studies on simplified models.

In this model, the deposition was deterministic and depended on the sedimentation and the local flow velocity. We observed high particle deposition in the bifurcation zone and in the anterior and posterior walls of the bifurcation due to the secondary flow velocity. Deposition by impaction increased with increasing bifurcation angles and mass flows, whereas sedimentation was predominantly marked in zones where gravity was directed toward a branch wall. Therefore, deposition was highly heterogeneous in the model.

In conclusion, using CFD techniques, we simulated flow, viscous pressure drop, and particle deposition in a realistic model of the airway tree, extending from the trachea to the segmental bronchi. Velocity profiles were more asymmetrical than those obtained with symmetrical models (8). Some analogies with the flow in a curved pipe could be observed: peak velocity occurred near the interior wall of the bifurcation daughter tubes, and double-peaked velocity profiles were observed transversally. However, additional velocity peaks and secondary velocity patterns appeared and were due to the complex flow configurations generated by successive bifurcations. Viscous pressure was approximated by Pedley's formula. The value of the constant γ used by Pedley was too large, leading to overestimated pressure drop. New values of γ based on the simulations were computed for each generation and were generally smaller than the value given by Pedley. The model's asymmetry generated heterogeneous particle deposition for branches of the same generation. Deposition increased with particle size and had a minimum at ∼200 cm^{3}/s because the combined effect of sedimentation and impaction was minimum around that value. Because respiratory flows of ∼200 cm^{3}/s correspond to a minimum particle deposition in the considered upper airways; particle transfer to the peripheral region of the lung is therefore maximized, which is important for medical aerosols that have to go deep in the lung. The use of CFD and realistic models allows obtainment of simulated flow and particle deposition patterns that cannot be measured in vivo. More detailed studies of the deposition sites can therefore be made, enhancing our understanding of aerosol transport in the lung.

## GRANTS

C. van Ertbruggen is supported by a research fellowship of the Fonds pour la Formation à la Recherche dans l'Industrie et l'Agriculture, Brussels, Belgium.

## Acknowledgments

C. van Ertbruggen thanks Colinda Francke, Vincent Bouffioux, and Frédéric Wilquem (Numeca International, Brussels, Belgium) for all the advice and the technical support provided throughout this research. The authors thank Chong S. Kim of US Environmental Protection Agency for providing the numerical data used in Fig. 9*B*.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2005 the American Physiological Society