The human lung comprises about 300 million alveoli which are located on bronchioles between the 17th to 24th generations of the acinar tree, with a progressively higher population density in the deeper branches (lower acini). The alveolar size and aspect ratio change with generation number. Due to successive bifurcation, the flow velocity magnitude also decreases as the bronchiole diameter decreases from the upper to lower acini. As a result, fluid dynamic parameters such as Reynolds (Re) and Womersley (α) numbers progressively decrease with increasing generation number. In order to characterize alveolar flow patterns and inhaled particle transport during synchronous ventilation, we have conducted measurements for a range of dimensionless parameters physiologically relevant to the upper acini. Acinar airflow patterns were measured using a simplified in vitro alveolar model consisting of a single transparent elastic truncated sphere (representing the alveolus) mounted over a circular hole on the side of a rigid circular tube (representing the bronchiole). The model alveolus was capable of expanding and contracting in-phase with the oscillatory flow through the bronchiole thereby simulating synchronous ventilation. Realistic breathing conditions were achieved by exercising the model over a range of progressively varying geometric and dynamic parameters to simulate the environment within several generations of the acinar tree. Particle image velocimetry was used to measure the resulting flow patterns. Next, we used the measured flow fields to calculate particle trajectories to obtain particle transport and deposition statistics for massless and finite-size particles under the influence of flow advection and gravity. Our study shows that the geometric parameters (β and ΔV/V) primarily affect the velocity magnitudes, whereas the dynamic parameters (Re and α) distort the flow symmetry while also altering the velocity magnitudes. Consequently, the dynamic parameters have a greater influence on the particle trajectories and deposition statistics compared to the geometric parameters. The results from this study can benefit pulmonary research into the risk assessment of toxicological inhaled aerosols, and the pharmaceutical industry by providing better insight into the flow patterns and particle transport of inhalable therapeutics in the acini.

## Introduction

The acini or parenchyma within the human lungs facilitate gas exchange [(1)] between the lung and the blood by the action of about 300 million alveoli which together provide approximately 70–80 m2 of surface area. A complete understanding of pulmonary fluid mechanics and particle transport at the alveolar level is necessary to improve the assessment of the toxicological impact of harmful inhaled particulate matter such as smoke. Such information can also benefit the pharmaceutical industry by providing key insights for the design and delivery methods of inhaled therapeutics.

Macroscopically, the human lung which comprises 24 generations of bifurcations, can be divided into the upper or conducting airways, and the acini or respiratory region. The acini begin at the 17th generation and differ from the conducting airways in that the acinar airways are alveolated [(2)]. The acinar airways can be classified into three categories [(2)]: (i) Respiratory bronchioles located in the upper acini extending from the 17th to the 20th generation with one or a few alveoli attached to each bronchiole, (ii) Alveolar ducts extending from the 21st to the 23rd generation with a higher density of alveoli attached to each bronchiole, and (iii) Alveolar sacs which are clusters of alveoli located at the terminus of the acinar tree.

Although the alveolus can be approximated as a 5/6th spherical cap, the alveolar size and shape vary from upper to lower acini [(2)]. Also, the alveolar size changes with age and height [(3),4], especially for infants whose airway structure and alveoli size/shape change considerably during the first 18 months [(3),5]. Due to successive bifurcation, the flow velocity magnitude also decreases as the bronchiole diameter decreases from the upper to lower acini. As a result, the associated geometric, and dynamic parameters such as Reynolds and Womersley numbers, also vary along the acinar tree [(2)]. The varying parameters affect the local fluid mechanics which in turn cause variations in particle transport and deposition from one generation to the next. The current work focuses on the effect of geometric and dynamic parameters on fluid flow and particle transport in the upper acini (respiratory bronchioles).

In the past two decades, numerous studies focused on two-dimensional [(6,7,8,9,10,11,12)] and three-dimensional [(8),13,14,15,16,17,18,19,20] acinar models to understand alveolar fluid mechanics and particle transport using numerical [(6,7,8),13,14,15,16,17,18,19] and experimental [(6),12] approaches. Broadly, two types of acinar models were employed by these researchers: rigid-walled, and flexible-walled/oscillating. Within each type, various approximations to the acinar geometry were used: axisymmetric annular or toroidal alveoli with a central duct [(6),13], a single alveolus attached to a semi-infinite plane [(21)], and honeycomb alveoli and duct [(17)]. We have summarized the key findings of these studies in our previous manuscripts [(22),23]. In this paper, we only discuss the literature pertaining to the effect of geometric or dynamic parameters on acinar fluid mechanics and particle transport.

Karl et al. [(13)] employed an axisymmetric, rigid-walled, annular model to study the effect of alveoli aspect ratio (AR) and alveoli-to-bronchiole diameter ratio (Da/Db) on the resulting flow patterns. Their experimental and numerical results showed a single recirculating region occupying the complete alveolar volume for smaller AR, and one large and two small recirculating regions in each alveolus for larger ARs. Furthermore, they observed a streamline spanning the alveolar opening which separates alveolar and ductal fluids for all aspect ratios. A more sophisticated geometry was used by Federspiel and Fredberg [(6)] wherein they numerically modeled a rigid-walled toroidal geometry and varied Da/Db and the alveolar-opening-to-cell-length ratios (AL) to study their effect on the resulting flow patterns.

They observed greater gas mixing for larger Da/Db and smaller AL. Tsuda et al. [(11)] imparted wall motion to their previously used rigid-walled toroidal approximation [(9)] to study the effect of multiple geometric/dynamic parameters on alveolar flow patterns. They observed that AR and alveolar-to-ductal flow ratios (QA/Qd) had the most noticeable impact on the resulting flow patterns.

In contrast to the simpler axisymmetric two-dimensional model, Sznitman et al. [(16)], and Kumar et al. [(17)] employed a more complex, three-dimensional, rhythmically expanding/contracting, multigeneration numerical model to study the variation of flow patterns along the entire acinar tree. They observed a large recirculation and a saddle point in each alveolus corresponding to the upper acini. The recirculations diminish to yield radially inward/outward flow in the terminal alveoli. Tsuda et al. [(7)]. numerically modeled multiple generations of the acinar tree to observe alveolar flow patterns resulting from variations in the geometric parameters. Furthermore, Tsuda et al. [11] extended their numerical work by considering a three-dimensional model to compare the effect of AR and Da/Db between healthy human lungs and immature, developing lungs in infants [(24)].

Since the alveolar geometry, Reynolds and Womersley numbers, and the resulting flow patterns vary along the acinar tree, fluid mixing and dispersion, and particle transport may exhibit differences across acinar generations. At the same time, the orientation of a given alveolus with respect to gravity may also affect particle transport and deposition. The current work focuses on such geometric, dynamic, and gravity orientation variations in the upper acini. We consider only in-phase flow which mimics synchronous ventilation in healthy human lungs.

## Experimental Setup and Simulation Methodology

Our experimental model consisted of a single alveolus constructed from highly elastic latex film (elongation limit ∼800%) mounted over a 16.5 mm hole drilled into the on the side of a rigid Plexiglas tube representing the bronchiole (34.8 mm internal diameter, with 1.6 mm wall thickness) as shown in Fig. 1. The alveolus and bronchiole were housed in a sealed, fluid-filled, rectangular chamber which provided the required, refractive index-matched, optical access. The oscillation of the alveolus and the oscillatory flow through the bronchiole were synchronously driven by syringe pumps connected to a stepper motor stage that was computer controlled to provide sinusoidal motion. The working fluid was an 80:20 mixture of glycerol and de-ionized water with a kinematic viscosity of 55 × 10−6 m2/s and the time period of oscillation was on the order of several minutes. To ensure realistic breathing conditions in our in vitro model, three dimensionless geometric parameters were matched to the actual acinar geometry [(1)] as described in Table 1. In addition, dynamic similarity was achieved by matching model Reynolds (Re) and Womersley $(α=Dω/v)$ numbers to human alveoli. Re was calculated using the mean velocity in the bronchiole during one-half of the sinusoidal cycle, and the bronchiole diameter (D) was used as the length scale for determining both Re and α. Here, ω is the breathing frequency, and v is the kinematic viscosity of the working fluid. The physiologically relevant values for these two parameters are compared with our experimental values in Table 1. Our experimental model closely replicates the upper acini wherein the model alveolus expands/contracts in-sync with the oscillatory flow through rigid bronchiole. Additional details about our model alveolus, experimental conditions, and experimental setup are described in our previous study [(22),23]. The velocity maps for case E in Chhabra and Prasad [22] are presented as dataset 1 in this paper.

Table 1

Physiological and experimental values of geometric and dynamic parameters. The experimental values of the particle transport parameters correspond to the values used in the particle tracking simulations for dp = 1.5 μm particles

 Geometric parameters Dynamic parameters Particle transport parameters β $(DalyDbronch)$ γ $(AlveolaropeningDalv)$ $ΔVV$ (Fractional volumetric change) Re (Reynolds number) α (Womersley number) St (Stokes number) H (Gravity number) Physiological range 0.4–2 0.7 0.5 - 0.8 (≤ 0.5 in a diseased lung) <1 ∼0.15 ≤10−4 0.4–20 Experimental values 0.5, 0.65 0.7 0.25, 0.5, and 0.65 0.1, 0.2, and 0.8 0.15, 0.3, and 0.4 2 × 10−5 0.5
 Geometric parameters Dynamic parameters Particle transport parameters β $(DalyDbronch)$ γ $(AlveolaropeningDalv)$ $ΔVV$ (Fractional volumetric change) Re (Reynolds number) α (Womersley number) St (Stokes number) H (Gravity number) Physiological range 0.4–2 0.7 0.5 - 0.8 (≤ 0.5 in a diseased lung) <1 ∼0.15 ≤10−4 0.4–20 Experimental values 0.5, 0.65 0.7 0.25, 0.5, and 0.65 0.1, 0.2, and 0.8 0.15, 0.3, and 0.4 2 × 10−5 0.5
Figure 1
Figure 1

### Data Acquisition and Processing.

In the current work, we acquired PIV data to create a total of eight datasets as listed in Table 2. The datasets were designed to isolate the effect of each dimensionless parameter while keeping all other parameters fixed, and to cover a wide range of physiologically relevant geometric and dynamic parameters. However, we were limited to Re ≥ 0.1 and α ≥ 0.15 due to the extremely small velocities and long experimental durations for smaller Re and α values.

Table 2

Values of geometric and dynamic parameters for each data set

 Dataset number Geometric parameters Dynamic parameters β $(DalyDbronch)$ γ $(AlveolaropeningDalv)$ $ΔVV$ (Fractional volumetric change) Re (Reynolds number) α (Womersley number) 1 0.65 0.7 0.5 0.1 0.3 2 0.5 0.7 0.5 0.1 0.3 3 0.5 0.7 0.5 0.8 0.3 4 0.5 0.7 0.5 0.2 0.3 5 0.5 0.7 0.5 0.2 0.4 6 0.5 0.7 0.5 0.2 0.15 7 0.65 0.7 0.25 0.1 0.3 8 0.65 0.7 0.65 0.1 0.3
 Dataset number Geometric parameters Dynamic parameters β $(DalyDbronch)$ γ $(AlveolaropeningDalv)$ $ΔVV$ (Fractional volumetric change) Re (Reynolds number) α (Womersley number) 1 0.65 0.7 0.5 0.1 0.3 2 0.5 0.7 0.5 0.1 0.3 3 0.5 0.7 0.5 0.8 0.3 4 0.5 0.7 0.5 0.2 0.3 5 0.5 0.7 0.5 0.2 0.4 6 0.5 0.7 0.5 0.2 0.15 7 0.65 0.7 0.25 0.1 0.3 8 0.65 0.7 0.65 0.1 0.3

#### Velocity Fields.

For each data set, a total of 400 equi-spaced PIV double frames were extracted and interrogated from a complete breathing cycle. The time separation (Δt) between PIV double frames was 1 s for all datasets except 5 and 6, for which the time separation was 4 s, and 0.5 s, respectively. For all datasets, image acquisition commenced after one complete breathing cycle to ensure fully-developed flow conditions. The particle images were interrogated using 64 × 64 pixel interrogation spots with 50% overlap. The resulting velocity maps were saved for each cycle. The velocity maps from 10 successive cycles were averaged to reduce noise.

#### Particle Transport and Deposition.

Particle trajectories were computed by solving the equations of motion under the combined influence of flow-induced advection and gravity. The complete particle tracking algorithm is explained in our previous study [(22),23]. A brief description of the particle-tracking algorithm is presented here. The governing equation for particle motion can be written as:
$mp∂u→p∂t=∑F→$
1
$=F→viscous+F→gravity$
2
As suggested by Haber et al. [25], the governing equation can be rewritten in dimensionless form as follows:
$St∂u→p'∂t'=H2g→+(u→'f-u→'p)$
3
where
$St=ρpdp2ωCc18μf,H=(ρp-ρf)dp2Ccg9μfDω,u→'p=u→p/ωD,u→'f=u→f/ωD,t'=ωt$
4
where $g→$ is the unit vector in the direction of gravity, g is the acceleration due to gravity, dp is the particle diameter, ρp is its material density (∼ 1000 kg/m3), $u→p$ is the particle velocity, $u→f$ is the local fluid velocity obtained by PIV measurements [(22)], ρf is the air density, μf is the air viscosity, and Cc is the Cunningham correction factor [(26)], ω is the breathing frequency $(=2π/T)$, T is the period of the breathing cycle, D is a characteristic length scale ( = alveolus diameter), St is the Stokes number, and H is the gravity number. St and H are the only two dimensionless numbers relevant to this problem. If these two numbers are matched between the in vitro model and the actual in vivo alveolus, then all of the particle dynamics will be accurately reproduced. For all dp < 5 μm particles St is small enough that we can ignore inertial impaction (See Table 1). The particle trajectory can be obtained by integrating Eq. 3; a second order Runge-Kutta method was used to calculate particle trajectories. The particle location at the (t + 1) time step can be written as:
$xpt+1=xpt+12[upt(xpt)+upt+1(xpt+1)]dt$
5
Here dt is the time step, $upt$ is particle velocity, and $xpt$ is particle location at time = t. The particle tracking algorithm described here assumes that the particle loading is small enough that the particles do not influence the fluid velocities (one-way coupling).

In the current work, we calculated trajectories for massless (dp = 0) and one finite-size particle (dp = 1.5 μm) for three gravity orientations for each data set. The massless particle transport isolates the effect of flow advection for a chosen geometric/dynamic parameter combination, whereas the finite-size particle transport reflects the combined influence of flow advection and gravity.

## Results and Discussion

In this section, we will present velocity maps, massless particle maps, and finite-size particle maps for each dataset. Next, we will discuss the effect of each parameter on the velocity fields and particle maps.

Useful information can be obtained by examining the velocity maps corresponding to four time instants within the breathing cycle: (1) peak inhalation (Fig. 2), (2) beginning of transition from inhalation to exhalation (Fig. 3), (3) peak exhalation (Fig. 4), and (4) beginning of transition from exhalation to inhalation (Fig. 5). As described in Chhabra and Prasad [22], during each breathing cycle the bronchiole and alveolar fluid undergo mixing as explained below: the bronchiole fluid is advected into the alveolus during inhalation, after which it mixes with the alveolar fluid during transition from inhalation to exhalation, the mixed fluid is ejected from the alveolus during exhalation, and the fluid retained within the alveolus under-goes another period of mixing during transition from exhalation to inhalation. In the following four subsections we will study the effect of each dimensionless parameter on the velocity maps and fluid mixing. For comparison, our baseline velocity map (dataset 1) corresponds to case E of our previous manuscript [(22)]. Furthermore, we will compare particle transport for massless and dp = 1.5 μm particles in three gravity orientations with those for case E in Chhabra and Prasad [22] (dataset 1). The particle trajectories after five and 10 cycles are shown in Figs. 6,7, respectively.

Figure 2
Figure 2
Figure 3
Figure 3
Figure 4
Figure 4
Figure 5
Figure 5
Figure 6
Figure 6
Figure 7
Figure 7

### Effect of β on Flow Patterns and Particle Transport.

The effect of aspect ratio or β (ratio of alveolar diameter to bronchiole diameter) can be observed by comparing velocity maps (a) and (b) in Figs. 2 through 5 which correspond to β = 0.5 and 0.65, respectively; β = 0.65 corresponds to the baseline case. These β values match the physiological values in the upper acini [(1)] which is the focus of our study.

A direct comparison between submaps (a) and (b) reveals that the flow topology remains largely unaffected by β: vortex rings are present during transition periods, the flow is slightly asymmetric during transition from inhalation to exhalation, and almost symmetric during transition from exhalation to inhalation. Although the fractional volumetric change in the alveolus during expansion/contraction was kept constant at the baseline value of 50%, the absolute volumetric change is smaller for β = 0.5. Therefore, the velocities during peak inhalation/exhalation are smaller as well (Figs. 2a,a, 2b,b, 4a,a, and 4b,b). Likewise, the velocities associated with the primary vortex rings also decrease with β due to the smaller amplitude of the alveolar wall motion (Figs. 3a,a, 3b,b, 5a,a, and 5b,b). For the larger β value, a secondary, weaker vortex ring with the opposite rotational sense sits atop the primary vortex ring during transition periods. The smaller volume of the alveolus corresponding to the smaller β value causes the secondary vortex ring to disappear in both the transition periods (Figs. 3a,a, 3b,b, 5a,a, and 5b b). The differences in velocity maps affect particle transport for both massless and finite particles as explained below.

At the end of 5 and 10 cycles, the deposition of massless and finite-size particles for all three gravity orientations increases as β decreases from 0.65 to 0.5 as shown in Fig. 6(a1)–6(d1) and Figs. 7(a2)–7(d2). Furthermore, the number of suspended particles for β = 0.5 is greater than that for β = 0.65 after five cycle (Figs. 6(a1)–6(d1) and 6(a2)–6(d2)). However, the opposite trend is observed after 10 cycles (Figs. 7(a1)–7(d1) and 7(a2)–7(d2)) due to additional particle deposition (or rejection for $θ=π$ only) for β = 0.5 between 5 and 10 cycles. Moreover, the particle deposition is nonuniform but well spreaded on the alveolar wall for β = 0.5 (Figs. 6(a2)–6(d2), 7(a2)–7(d2)) while it is limited to a few sites on the alveolar wall for β = 0.65 (Figs. 6(a1)–6(d1), 7(a1)–7(d1)).

The larger particle deposition observed for β = 0.5 can be attributed to the absence of the secondary vortex ring during transition periods: although the secondary vortex ring enhances fluid mixing, it also keeps the particles in suspension for longer durations. In the absence of the secondary vortex ring for β = 0.5, the primary vortex ring has a stronger influence near the deep alveolar wall due to which particle deposition is higher. Furthermore, although the velocities inside the alveolus increase with β, the proportionality constant is not strictly linear (Figs. 2a,a and 2b,b through 5a,5b ab). The smaller alveolar size implies that particles need to travel a shorter distance in order to deposit, especially when the particle settling velocity due to gravity is superposed. As a result, the overall particle deposition time decreases and deposition increases for smaller β.

### Effect of ΔV/V on Flow Patterns and Particle Transport.

The subfigures (a), (g), and (h) in Figs. 2,3,4,5 display velocity maps for ΔV/V = 0.5, 0.25, and 0.65, respectively. As before, these values of fractional volumetric change in the alveolus (ΔV/V) are in good agreement with physiological values in the upper acini. Note that a change in ΔV/V may change β by a few percent which is small compared to the variation in β (23%) considered in the previous section. The flow topology during the complete breathing cycle remains similar for different ΔV/V similar to the cases for varying β (see Sec. 3). Furthermore, a pair of vortex rings was observed for all ΔV/V values during both transition periods. Also, the flow is symmetric during peak inhalation/exhalation and transition from exhalation to inhalation, and slightly asymmetric during transition from inhalation to exhalation for all ΔV/V values (Parts (a), (g), (h) in Figs. 2,3,4,5).

Despite an overall similarity between the flow patterns, the alveolar fluid exhibits smaller velocities for smaller ΔV/V at each time step during a complete breathing cycle: the inward/outward velocities during peak inhalation/exhalation are smaller and vortex rings are weaker during both transition periods (compare parts (g) with (a) and (h) in Figs. 2,3,4,5). Such a reduction in velocities is caused by the smaller quantity of bronchiole fluid entering/exiting during inhalation/exhalation periods for smaller ΔV/V.

Parts (a1)–(d1), (a7)–(d7), and (a8)–(d8) in Figs. 6,7 show particle maps for ΔV/V = 0.5, 0.25, and 0.65, respectively. The massless particle maps (parts (a1), (a7), (a8) in Figs. 6,7) show that particle deposition is highest for ΔV/V = 0.5 after 5 and 10 cycles. In fact, most of the massless particles are ejected out of the alveolus for the other two ΔV/V values (Figs. 6(a7), 6(a8) and 7(a7), 7(a8)). The transport of dp = 1.5 μm particles reveals that particle deposition is highest for ΔV/V = 0.5 in each gravity orientation. In contrast, the fraction of suspended particles is highest for ΔV/V = 0.65 for θ = 0 orientation (compare parts (b8) with (b1) and (b7) in Figs. 6,7), and for ΔV/V = 0.5 for the other two gravity orientations (compare parts (c1), (d1), with (c7), (d7), and (c8), (d8) in Figs. 6,7).

The greater particle deposition for ΔV/V = 0.5 can be attributed to the moderate expansion/contraction which encourages more particles to enter into the alveolus during the first cycle. A small fraction of these particles leaves the alveolus in at most four cycles while others either remain suspended or deposit on the alveolar wall. In contrast, the small velocities associated with ΔV/V = 0.25 are unable to advect a reasonable portion of particles during the first cycle. Consequently, the fraction of suspended particles is smaller and deposition is low. By this argument, the deposition is expected to be higher for ΔV/V = 0.65 but in fact, the deposition is small again which can be explained as follows: during transition periods the vortex rings enhance particle transport and deposition into the alveolus. For ΔV/V = 0.65, the recirculation region is small because the alveolus is small during transition from exhalation to inhalation (Fig. 5h h). Therefore, the particle transport is smaller compared to the baseline case (ΔV/V = 0.5) and lowers particle deposition.

The higher particle suspension for ΔV/V = 0.65 for $θ=π$ orientation can be explained as follows: the gravity reinforces large inward velocities during peak inhalation while counteracting the large outward velocities during peak exhalation. Hence, during a complete breathing cycle, a net inward displacement develops but its magnitude is small. Therefore, more particles remain suspended for a longer duration. Furthermore, the penetration of particles is slow due to the small net inward displacement.

### Effect of Reynolds Number on Flow Patterns and Particle Transport.

Three sets of velocity maps for Re = 0.1, 0.2 and 0.8 are shown in parts (b), (d), and (c), respectively, in Figs. 2,3,4,5. Although the flow topology does not change as the Reynolds number is increased from 0.1 through 0.8, the flow inside the alveolus becomes asymmetric throughout the breathing cycle. The largest inward/outward velocities, which were located at the center of the alveolar opening, display a noticeable shift in the bronchiole flow direction (compare (b) and (c) with (d) in Figs. 2,4). Furthermore, the vortex ring which was slightly asymmetric in the baseline case, now displays a significantly larger circulation region in the left hemisphere compared to the right hemisphere during both the transition periods which can be observed by comparing (b) and (c) in Figs. 3, and 5. In addition, the saddle point drifts towards the right during both the transition periods (compare (b) and (c) in Figs. 3, and 5).

The formation of an asymmetric vortex ring during transition from inhalation to exhalation (Fig. 3d,d) is expected due to history effects from the peak exhalation period, as explained in our previous work [(22)]. During peak inhalation the shear layer formation is promoted by the alveolar expansion. Consequently, the rapid bronchiole flow (dataset 3, Re = 0.8) causes a tilt in the peak velocities towards the bronchiole flow direction (right side, Fig. 2d,d). In contrast, the influence of the bronchiole fluid motion is suppressed by large outward velocities during peak exhalation. Therefore, the bronchiole fluid is unable to alter the flow inside the alveolus. On the other hand, the flow asymmetry developed during transition from inhalation to exhalation continues to survive until peak exhalation due to a lack of adequate damping (Fig. 5d,d). Furthermore, since the peak exhalation is unable to damp out the asymmetry developed during transition from inhalation to exhalation, this asymmetry continues through transition from exhalation to inhalation as observed in Fig. 5d d.

Particle maps for Re = 0.1, 0.8, and 0.2 are shown in parts (a2)–(d2), (a3)–(d3), and (a4)–(d4), respectively, in Figs. 6,7. A direct comparison between (a2)–(d2) through (a4)–(d4) in Figs. 6,7 reveals that particle deposition is higher for Re = 0.2 for massless as well as dp = 1.5 μm particles in all gravity orientations. In contrast, Re = 0.1 has a greater number of suspended particles after 5 and 10 cycles ((a2)–(d2) in Figs. 6,7). Particle deposition decreases dramatically as Re is increased to 0.8 for massless as well as dp = 1.5 μm particles and does not show much gravity influence until five cycles (Fig. 6(a3)–6(d3)). However, the particle deposition for Re = 0.8 increases slightly in 10 cycles for $θ=π$ gravity orientation but remains smaller compared to the other two Re values (compare (b3) with (b2) and (b4) in Fig. 7). Furthermore, the $θ=π$ gravity orientation shows a larger fraction of suspended particles compared to the other two gravity orientations after 5 and 10 cycles (in Figs. 6(b3) and 7(b3)).

As previously noticed, the velocity maps for Re = 0.2 exhibit slightly stronger recirculation during both transition periods, compared to Re = 0.1 (compare parts (b) with (d) in Figs. 3,5). The differential between the two recirculations accumulates over multiple cycles and causes a higher particle deposition for Re = 0.2 ((a4)–(d4) in Figs. 6,7). In contrast, more particles remain suspended for Re = 0.1 due to a lower recirculation in this case ((a2)–(d2) in Figs. 6,7). As Re is increased to 0.8, the velocity field becomes significantly asymmetric during both transition periods (Figs. 3c,c and 5c,c). However, the shear layer remains separated from the vortex rings. Consequently, the shear layer carries away all the particles introduced near the alveolar opening. Hence, particle deposition is lowest for Re = 0.8 for both massless and dp = 1.5 μm particles in all three gravity orientations ((a3)–(d3) in Figs. 6,7). However, gravity causes streamline crossing for some particles for the $θ=π$ orientation due to which some particles enter the recirculation region ((b3) in Figs. 6,7). Consequently, they remain suspended after 5 and 10 cycles and eventually begin to deposit on alveolar wall after 10 cycles (Fig. 7(b3)).

### Effect of Womersley Number on Flow Patterns and Particle Transport.

Parts (d), (e), and (f) in Figs. 2,3,4,5 display velocity maps for three physiologically relevant α values of 0.3, 0.42, and 0.15, respectively. It is observed that a strongly affects the alveolar flow unlike the weaker effects caused by varying the geometric parameters. Similar to the results for the previous parametric variations, the overall flow topology remains the same for all three α values, i.e., the bronchiole fluid enters/exits during inhalation/exhalation, and the alveolar and bronchiole fluid mix during transition periods (parts (d), (e), and (f) in Figs. 2,3,4,5).

Similar to Re, the effect of α is clearly noticeable at each time step in the breathing cycle. As α is decreased from 0.42 to 0.15, the vortex rings during transition periods become asymmetric which can be observed by comparing (e) with (d) and (f) in Figs. 3,5. Furthermore, the asymmetry is significantly stronger during transition from exhalation to inhalation compared to the other transition period for α = 0.15 (compare Figs. 3f,f and 5f f). In fact, during transition from exhalation to inhalation, the left circulation suppresses the right circulation, merges with the shear layer, and occupies almost the entire alveolar volume. During peak inhalation/exhalation, a small α causes a shift in the largest inward/outward towards the bronchiole flow direction. A similar effect was also observed for large Re ( = 0.8) in the previous subsection.

The asymmetry in vortex rings during transition periods can be explained as follows: although the fractional volumetric expansion/contraction during inhalation/exhalation is the same for all α values, the rate of alveolar expansion/contraction is significantly different. In fact, α = 0.15 represents 4 and 8 times slower expansion/contraction compared to α = 0.3 and 0.42, respectively. Consequently, the velocities induced by the alveolar wall motion are significantly smaller for the smaller α value. In contrast, the bronchiole flow for all α remains unchanged (Re = 0.2). As a result, the bronchiole flow-induced velocities overwhelm velocities caused by the alveolar wall motion and the vortex rings become highly asymmetric. A similar effect causes the largest velocities to shift in the bronchiole flow direction during peak inhalation/exhalation.

Parts (a4)–(d4), (a5)–(d5), and (a6)–(d6) in Figs. 6,7 demonstrate that particle deposition and the fraction of suspended particles decrease with α for massless and finite particles in all gravity orientations. On the other hand, the deposition sites remain unaffected as α is decreased from 0.42 to 0.3 but become fewer for an even smaller α (=0.15). As previously noticed, the velocity maps and the particle deposition sites are quite similar for α = 0.42 and 0.3 ((a4)–(d4), (a5)–(d5) in Figs. 6,7). However, higher inhalation velocities in combination with higher recirculation velocities for α = 0.42 (see (d) and (e) in Figs. 2,3,4,5) push more particles into the alveolus resulting in greater particle deposition and a higher fraction of suspended particles ((a4)–(d4), (a5)–(d5) in Figs. 6,7). When α is reduced to 0.15, the overall flow topology changes (compare Figs. 3f,f and 5f,f with 3d,d and 3e,e and 5d,d and 5e,e), i.e., the secondary vortex ring and the shear layer disappear and the primary vortex ring becomes highly asymmetric. Consequently, during transition from exhalation to inhalation almost all the introduced particles are swept away by the merged shear layer-asymmetric vortex ring structure (Fig. 5f f).

In summary, our results show that the overall impact on velocity maps of dynamic parameters is stronger than geometric parameters. Consequently, the dynamic parameters alter particle trajectories more significantly than the geometric parameters.

## Implications of Our Results

Our velocity maps demonstrate larger velocity magnitudes, greater fluid mixing, and smaller particle deposition for higher values of the geometric parameters parameters employed in this study. We also observed an asymmetric velocity field with greater mixing and larger particle deposition for smaller α. These observations are useful in predicting the variation of flow and particle transport as we travel deeper into the lung acinus which is characterized by increasing values of geometric parameters (β and ΔV/V), and decreasing values of dynamic parameters (Re and α) [(27)]. Therefore, our parametric study indicates that velocity maps would show symmetric flow patterns in the upper acini, and largely asymmetric flow in the lower acini. Nonetheless, the upper acini exhibit greater flow mixing yet a smaller particle suspension due to topological differences arising from the associated geometric variations. Moreover, our particle maps show that the particle deposition time will be smaller and greater deposition will occur in the upper acinar generations compared to the lower acinar generations.

## Conclusions

This paper presents alveolar flow patterns and particle transport results for a range of geometric and dynamic parameters within the acinar region of the human lung. Specifically, we considered two geometric (β and ΔV/V), and two dynamic parameters (Reynolds and Womersley numbers) to study their effect on alveolar flow patterns and particle deposition. We considered massless and dp = 1.5 μm particles for the latter study.

The results for β show that the velocities increase with β at each time step during a breathing cycle. The effect of β is more prominent during transition periods due to smaller velocities during these periods. In addition, an alveolus with a smaller β does not show the secondary vortex ring during transition due to its smaller size. Although the velocity maps do not exhibit a significant influence of β, the particle maps indicate a significantly higher particle deposition for smaller β. Our study also demonstrates that the velocities are strongly influenced by ΔV/V. Furthermore, the fraction of suspended particles and particle deposition is strongly affected by ΔV/V. The basline value of ΔV/V = 0.5 enhances particle penetration and deposition while smaller and larger values can reduce particle suspension and deposition. Although these results vary slightly with gravity orientation, they are generally true for all particles smaller than 2.5 μm (small and midsize particles).

The results for Reynolds number show that a higher Re causes asymmetric vortex rings during transition periods and shifts the maximum inhalation/exhalation velocities towards the bronchiole flow direction. Furthermore, as Re is increased, the deposition initially increases before decreasing for all gravity orientations for massless particles. The results for a display the opposite behavior: the velocity maps become highly asymmetric at each time step during the breathing cycle for smaller α. The particle maps show that deposition increases as α is decreased from 0.42 to 0.15.

Overall, our velocity maps show that the geometric parameters primarily affect the velocity magnitudes but the dynamic parameters distort the flow symmetry while also altering the velocity magnitudes. Also, the particle maps display a greater influence of dynamic parameters compared to geometric parameters. Our results indicate greater flow mixing and particle deposition but a lower particle suspension for upper acini compared to lower acini.

## Acknowledgment

This research was supported by Philip Morris USA, Inc., and Philip Morris International.

## References

1.
Weibel
,
E. R.
, 1963,
Morphometry of the Human Lung,
,
New York
.
2.
Bleuer
,
B. H.
, and
Weibel
,
E. R.
, 1988, “
Morphometry of the Human Pulmonary Acinus
,”
Anat. Rec.
,
220
, pp.
401
414
.
3.
Zeltner
,
T. B.
,
,
J. H.
,
Gehr
,
P.
,
Pfenninger
,
J.
, and
Burri
,
P. H.
, 1987, “
The Postnatal Development and Growth of the Human Lung. I. Morphometry
,”
Respir. Physiol.
,
67
, pp.
247
267
.
4.
Angus
,
G. E.
, and
Thurlbeck
,
W. M.
, 1972, “
Number of Alveoli in the Human Lung
,”
J. Appl. Physiol.
,
32
, pp.
483
485
.
5.
Burri
,
P. H.
, 2006, “
Structural Aspects of Postnatal Lung Development—Alveolar Formation and Growth
,”
Biol Neonate
,
89(4)
, pp.
313
322
.
6.
Federspiel
,
W. J.
, and
Fredberg
,
J. J.
, 1988, “
Axial Dispersion in Respiratory Bronchioles and Alveolar Ducts
,”
J. Appl. Physiol.
,
64
, pp.
2614
2621
.
7.
Tsuda
,
A.
,
Henry
,
F. S.
, and
Butler
,
J. P.
, 2008, “
Gas and Aerosol Mixing in the Acinus
,”
Respir, Physiol, Neurobiol
,
163
, pp.
139
149
.
8.
Darquenne
,
C.
, and
Paiva
,
M.
, 1996, “
Two- and Three-Dimensional Simulations of Aerosol Transport and Deposition in Alveolar Zone of Human Lung
,”
J. Appl. Physiol.
,
80
, pp.
1401
1414
.
9.
Tsuda
,
A.
,
Butler
,
J. P.
, and
Fredberg
,
J. J.
, 1994, “
Effects of Alveloated Duct Structure on Aerosol Kinetics I. Diffusional Deposition in the Absence of Gravity
,”
J. Appl. Physiol.
,
76
, pp.
2497
2509
.
10.
Tsuda
,
A.
,
Butler
,
J. P.
, and
Fredberg
,
J. J.
, 1994, “
Effects of Alveloated Duct Structure on Aerosol Kinetics II. Gravitational Sedmimentation and Inertial Impaction
,”
J. Appl. Physiol.
,
76
, pp.
2510
2516
.
11.
Tsuda
,
A.
,
Henry
,
F. S.
, and
Butler
,
J. P.
, 1995, “
Chaotic Mixing of Alveolated Duct Flow in Rhythmically Expanding Pulmonary Acinus
,”
J. Appl. Physiol.
,
79
, pp.
1055
1063
.
12.
Tippe
,
A.
, and
Tsuda
,
A.
, 2000, “
Recirculating Flow in an Expanding Alveolar Model: Experimental Evidence of Flow-Induced Mixing of Eerosols in the Pulmonary Acinus
,”
J. Aerosol. Sci.
,
31
, pp.
979
986
.
13.
Karl
,
A.
,
Henry
,
F.
, and
Tsuda
,
A.
, 2004, “
Low Reynolds Number Viscous Flow in an Alveolated Duct
,”
J. Biomech. Eng.
,
126
, pp.
420
429
.
14.
van Ertbruggen
,
C.
,
Corieri
,
P.
,
Theunissen
,
R.
,
Riethmuller
,
M.
, and
Darquenne
,
C.
, 2008, “
Validation of cfd Predictions of Flow in a 3d Alveolated Bend with Experimental Data
,”
J. Biomech. Eng.
,
41
, pp.
399
405
.
15.
Sznitman
,
J.
,
Heimsch
,
F.
,
Heimsch
,
T.
,
Rusch
,
D.
, and
Rosgen
,
T.
, 2007, “
Three-Dimensional Convective Alveolar Flow Induced by Rhythmic Breathing Motion of the Pulmonary Acinus
,”
J. Biomech. Eng.
,
129
, pp.
658
665
.
16.
Sznitman
,
J.
,
Heimshch
,
T.
,
Wildhaber
,
J. H.
,
Tsuda
,
A.
, and
Rosgen
,
T.
, 2009, “
Respiratory Flow Phenomena and Gravitational Deposition in a Three-Dimensional Space-Filling Model of the Pulmonary Acinar Tree
,”
J. Biomech. Eng.
,
131
, pp.
1
15
.
17.
Kumar
,
H.
,
Tawhai
,
M. H.
,
Hoffman
,
E. A.
, and
Lin
,
C. L.
, 2009, “
The Effects of Geometry on Airflow in the Acinar Region of the Human Lung
,”
J. Biomech. Eng.
,
42
, pp.
1635
1642
.
18.
Li
,
Z.
,
Kleinstreuer
,
C.
, and
Zhang
,
Z.
, 2007, “
Particle Deposition in the Human Tracheobronchial Airways Due to Transient Inspiratory Flow Patterns
,”
J. Aerosol. Sci.
,
38
, pp.
625
644
.
19.
Li
,
Z.
,
Kleinstreuer
,
C.
, and
Zhang
,
Z.
, 2007, “
Simulation of Airflow Fields and Microparticle Deposition in Realistic Human Lung Airway Models. Part I: Airflow Patterns
,”
Eur. J. Mech. B Fluids
,
26
, pp.
632
649
.
20.
Henry
,
F. S.
,
Laine-Pearson
,
F. E.
, and
Tsuda
,
A.
, 2009, “
Hamiltonian Chaos in a Model Alveolus
,”
J. Biomech. Eng.
,
131
, p.
011006
.
21.
Haber
,
S.
,
Butler
,
J. P.
,
Brenner
,
H.
,
Emaneul
,
I.
, and
Tsuda
,
A.
, 2000, “
Shear Flow Over a Self-Similar Expanding Pulmonary Alveolus During Rhythmical Breathing
,”
J. Fluid Mec.h
,
405
, pp.
243
268
.
22.
Chhabra
,
S.
, and
,
A. K.
, 2010, “
Flow and Particle Dispersion in a Pulmonary Alveolus–Part I: Velocity Measurements and Convective Particle Transport
,”
J. Biomech. Eng.
,
132
, p.
051009
.
23.
Chhabra
,
S.
, and
,
A. K.
, 2010, “
Flow and Particle Dispersion in a Pulmonary Alveolus–Part II: Effect of Gravity on Particle Transport
,”
J. Biomech. Eng.
,
132
, p.
051010
.
24.
Tsuda
,
A.
,
Filipovic
,
N.
,
Haberthur
,
D.
,
Dickie
,
R.
,
Stampanoni
,
M.
,
Matsui
,
Y.
, and
Schittny
,
J. C.
, 2008, “
The Finite Element 3d Reconstruction of the Pulmonary Acinus Imaged by Synchrotron X-Ray Tomography
,”
J. Appl. Physiol.
,
105
, pp.
964
976
.
25.
Haber
,
S.
,
Yitzhak
,
D.
, and
Tsuda
,
A.
, 2003, “
Gravitational Deposition in a Rhythmically Expanding and Contracting Alveolus
,”
J. Appl. Physiol.
,
95
, pp.
657
671
.
26.
Cunningham
,
E.
, 1910, “
On the Velocity of Steady Fall of Spherical Particles Through Fluid Medium
,”
Proc. R. Soc. A
,
83
, pp.
357
365
.
27.
Pedley
,
T. J.
, 1977, “
Pulmonary Fluid Dynamics
,”
Ann. Rev. Fluid Mech.
,
9
, pp.
229
274
.