|Home | About | Journals | Submit | Contact Us | Français|
Accurate modeling of air flow and aerosol transport in the alveolated airways is essential for quantitative predictions of pulmonary aerosol deposition. However, experimental validation of such modeling studies has been scarce. The objective of this study is to validate CFD predictions of flow field and particle trajectory with experiments within a scaled-up model of alveolated airways. Steady flow (Re = 0.13) of silicone oil was captured by particle image velocimetry (PIV), and the trajectories of 0.5 mm and 1.2 mm spherical iron beads (representing 0.7 to 14.6 μm aerosol in vivo) were obtained by particle tracking velocimetry (PTV). At twelve selected cross sections, the velocity profiles obtained by CFD matched well with those by PIV (within 1.7% on average). The CFD predicted trajectories also matched well with PTV experiments. These results showed that air flow and aerosol transport in models of human alveolated airways can be simulated by CFD techniques with reasonable accuracy.
Numerous experimental (Heyder et al., 1975; Stahlhofen et al., 1980; Heyder et al., 1986; Darquenne et al., 1998; Darquenne et al., 1999; Darquenne and Prisk G.K., 2004; Heyder, 2004) and modeling (Yeh and Schum, 1980; Darquenne and Paiva, 1994; Darquenne and Paiva, 1996; Fleming et al., 2000; Darquenne, 2001; Asgharian et al., 2006; Harrington et al., 2006) investigations have been performed in the past decades to estimate the deposition of inhaled aerosols in the human respiratory tract. While current in vivo experimental techniques can only reliably provide total or regional deposition data, numerical simulations can often provide much more detailed information in both space and time. Computational fluid dynamics (CFD) simulation is a powerful and often indispensible tool for modeling airflow and aerosol transport in realistic and complicated airway geometry. However, one of the most important aspects of this type of simulations is that the methods and results need to be validated before such predictions can be relied on.
There were several validation studies in which numerically computed airflow in the human upper and large airways were compared with experimental measurements. To validate CFD simulations of air flow in a rigid-walled tracheobronchial airway model from trachea up to generation seven, de Rochefort et al. (2007) performed in vitro validation studies by using hyperpolarized 3He magnetic resonance phase-contrast velocimetry. Good agreement between computed and measured velocity profiles at selected cross sections was observed. Some CFD-based predictions of particle deposition profiles in the upper and large airways have also been qualitatively validated by comparing regional deposition statistics (Oldham et al., 2000), or local maps of deposition (Zhang et al., 2002). Quantitative validation of local deposition patterns has however been scarce. In a recent study by Longest and Vinchurkar (2007) the CFD predicted deposition results for 10μm particles in a model of three successive airway generations (3 to 5) were compared to in vitro data (Oldham et al., 2000) and good agreement was obtained.
There have also been some experimental validation studies of air flow in the human peripheral airways. Tsuda et al. (1995) numerically simulated air flow in an axisymmetric model of the pulmonary acinus under rhythmical expansion and contraction. This numerical study was later validated experimentally (Tippe and Tsuda, 1999) in a large-scaled in vitro model using particle image velocimetry (PIV). Qualitative agreements on the flow field were found between CFD predictions and experimental data. Both the simulations and experiments revealed highly complex flow structures in compliant alveolar structures even at very low Reynolds numbers. Karl et al. (2004) studied the low Reynolds number flow in a large scale axisymmetric rigid-walled model of a single alveolated duct. The PIV-measured and CFD-simulated flow features qualitatively agreed with each other. This study showed the influence of model geometrical factors (e.g. alveolar cavity aspect ratio) on overall flow features in the periphery of the lung. More recently van Ertbruggen et al. (2008) reported quantitative comparison of CFD simulations and experimental PIV results for fluid flow in a rigid-walled scaled-up model of an alveolated bend. At seven selected cross sections the velocity profiles were compared and good agreement was found. However, the model geometry was limited to one bend and did not include any bifurcations.
To the authors’ knowledge, there have been no previous published studies which seek to validate simulated individual particle trajectory in the context of respiratory aerosol transport. In this study, a scaled-up model of human alveolated airways was built. The model included three generations and two bifurcations and contained some essential geometrical features of the pulmonary airways. Experimental measurements of fluid flow and particle transport were obtained by PIV and particle tracking velocimetry (PTV), respectively. These experimental data were compared to CFD predictions under identical conditions. This study provides the first direct quantitative validation of numerical studies of airway flow and individual particle trajectory in three dimensional bifurcated and alveolated airways.
An in vitro scaled-up model of three generations of alveolated airways of the human lung was constructed (Fig. 1) at the von Karman Institute (VKI) using the same casting techniques as used in a previous study (van Ertbruggen et al., 2008). The cast was made of silicone whose refractive index (1.43) matches that of the carrier fluid (silicone oil, density = 970 kg/m3, dynamic viscosity = 1 Pa s). The inner diameter of all five airways was 20 mm, and the diameter of the alveolar rings was 45 mm. These dimensions correspond to an in vivo lumen diameter of 267 μm and an outer diameter of 606 μm. To facilitate experimental measurements the model size was about 75 times larger than the average dimensions of intraacinar airways based on morphometric measurements (Haefeli-Bleuer and Weibel, 1988) scaled to a lung volume of 3.5 liter (Harrington et al., 2006). The two bifurcations were symmetric with bifurcation angle of 45° for both daughter airways, as used in a previous study (Yeh and Schum, 1980). As compared to the true anatomy the model had a few modifications for ease of model fabrication and optical access for experimental measurements. Circular rings were used instead of individual alveoli since it is very difficult to fabricate a spherical cavity with transparent wall and perform optical measurements in the spherical cavity. The flow conditions were selected to represent air flow and aerosol transport in the 20th–22th generations of the Weibel model (Weibel, 1963) of human airways by matching the dynamic (Reynolds number) similarity between the two systems. The model was oriented such that gravity acted in the symmetry plane and the gravitational vector coincided with the entrance flow velocity vector. Due to symmetry conditions both in geometry and in gravity force, distal branches of one of the daughter branches of the first bifurcation were not included. One of the daughter branches of the second bifurcation was parallel to the gravity direction, while the other one was in the horizontal direction.
The physiological flow facility used in this study was described in detail in a previous study by our group (van Ertbruggen et al., 2008). Briefly, steady flow (2.04 ml/s) through the model was generated by gravity of the fluid. Based on the inlet lumen diameter and the mean velocity at the inlet, the Reynolds number was 0.13 corresponding to an in vivo mouth flow rate of about 30 l/min. The flow rate at each of the three outlets was controlled by a valve and measured by a specially designed flowmeter.
For PIV measurements, the fluid (silicone oil) was seeded with iron particles with a nominal (physical) diameter of 20 μm. These particles have a relaxation time of 1.5e-7 second and can be assumed to be neutrally buoyant in the silicone oil. At the same time, the particles had a surface area large enough to scatter sufficient light to provide a clear tracer pattern. A continuous Innova 70C Argon laser supplied the coherent light source and a standard digital camera captured the scattered light intensities. The field of view in the PIV measurements covered an area of around 82×61 mm2 (~ 0.06 mm/pixel). Image recordings were analyzed within a recursive structure yielding final interrogation windows of 1.2×1.2 mm2 with a vector spacing of 0.6 mm. The images were recorded at a rate of 10 Hz and had a typical resolution of 1372×1024 pixel. The images were processed using a cross-correlation algorithm WiDIM (Scarano and Riethmuller, 2000) developed at VKI.
To simulate aerosol transport, spherical iron (density = 7800 kg/m3) beads with nominal diameter of 0.5 mm and 1.2 mm were used. No apparent anomalies in sphericity were noticed for these beads. By matching the particle Reynolds number (see Eq. 6 below) and the velocity ratio (particle settling velocity to mean flow velocity) between the in vitro and in vivo systems, the corresponding in vivo aerosol sizes are 6.1 μm and 14.6 μm, respectively. By matching the Stokes number, the in vivo sizes are 1.7 μm and 0.7 μm, respectively (see details in discussion). The beads were released by means of an electromagnet placed sufficiently far (about 40 mm) upstream the model inlet. In the experiments, the beads were released in three separate groups and tracked separately. In each group, if multiple beads were released the beads were at the same height and released simultaneously. The first group included four 1.2 mm beads, the second group included one 0.5 mm bead, and the third group included two 0.5 mm beads. Beads in the first and second groups were released in the symmetry plane, and beads in the third group were released from locations outside of the symmetry plane.
The trajectories of the beads were measured with time resolved 3D PTV by using two cameras, one in front of the model and the second at the bottom of the model (Fig. 1). The two cameras were identical to that used for the PIV measurements, and recorded images simultaneously at a frequency of 9 Hz. Both cameras were stationary with a field of view of approximately 140×140 mm2. With an objective of 35mm, the focal plane of camera 2 had a depth of around 10 cm. As the particle image centroids were retrieved from a Gaussian fitting, the small blur created due to the particle movement perpendicular to the thick focal plane was negligible. The recorded images followed a background subtraction process and particle trajectory was obtained by the four-frame method (Malik et al., 1993). A search area was built around each particle of the first frame. All particle images of the second frame located within the area were considered to be possible partners. Based on the assumption that the particles would travel between successive frames the same distance on the same path, search areas in the third and fourth frame were located at the linear extrapolated positions from the two previous frames. After the four frame method to initialize the tracking procedure, all position predictions were therefore made from a linear extrapolation of the previous particle positions (Theunissen et al., 2006; Ruwet et al., 2008).
A virtual (CFD) model of the in vitro airway model was generated in STAR-Design (ver. 4.02, CD-adapco). The 3D coordinate system (x-y-z) that was used to build the CFD model is shown in Fig. 1. The symmetry plane of the model was in the y–z plane. A polyhedral volume mesh was generated within the same software. Two prism layers were created near the wall boundary. A plot of the polyhedral mesh at the symmetry plane of the model is shown in Fig. 2 for the grid size of 674k cells. To evaluate mesh independence, results for grid sizes of 674k and 294k cells were compared. The maximum velocity at twelve selected cross sections (see Fig. 3) differed by about 0.5%, while the particle trajectories for the seven beads were almost identical for the two grid sizes. Therefore, the results were considered to have reached mesh independence and those corresponding to grid size of 674k cells were used.
A finite-volume-based flow solver, STAR-CD (ver. 4.02, CD-adapco), was used for both fluid flow and particle tracking simulations. The same flow conditions (steady laminar flow) and material property (silicone oil) as in the experiments were used in the fluid flow simulations. Assuming steady laminar flow of incompressible Newtonian fluid, the governing flow equations are:
where u is fluid velocity vector, p is static fluid pressure, ρis fluid density, μ is dynamic viscosity, and g is gravitational acceleration.
A flat velocity profile of 6.47 mm/s was prescribed at the model entrance, and outflow split ratios as measured in experiments were used at the model outlets. The tube without alveolar ring had a flow split ratio of 0.504 to the total flow at the inlet, and for the two outlets of the second bifurcation, the horizontal and vertical tubes had flow split ratios of 0.252 and 0.244, respectively. Gravity acted downward as in the experiments. A convergence criterion of 10−5 for residuals of mass and momentum balance equations was used for controlling the number of iterations. The second-order accurate scheme MARS (Monotone Advection and Reconstruction Scheme) in the STAR-CD solver was selected for spatial discretisation. The SIMPLE (Semi-Implicit Method for Pressure-linked Equations) algorithm was used for solving the discretised algebraic finite-volume equations.
Particle motion in carrier fluid is affected by various forces such as viscous drag force, gravity force, added mass force (virtual mass force), Brownian force, and pressure force. In this study, Brownian force was ignored due to the large bead size. The added mass force was considered in a few simulations and was found to have negligible influence on particle trajectory. Thus in this study only steady viscous drag force, gravity force and pressure force were considered. The governing equation of motion for a particle of mass m is, therefore:
where up is particle velocity vector, ρp is particle density, FD is drag force given by:
where Ap is cross sectional area of the particle, CD is drag coefficient. For the low Reynolds number flow in the current study the drag coefficient is:
where Rep is particle Reynolds number defined as:
where dp is particle diameter. In Eq. 3, p is fluid pressure and includes any hydrostatic component. Thus, the buoyancy force is included in the pressure force.
The Lagrangian governing equations (Eq. 3) were discretised into first order Euler-implicit form, and solved, in conjunction with the governing equations for carrier fluid (Eq. 1 and 2), by the SIMPLE algorithm. The particle time step was selected to be the smaller of particle momentum relaxation time scale and one-third of the average time for a particle to cross a characteristic local cell dimension.
To model aerosol particles (iron beads in experiments), the motion of spherical inert particles with nominal diameters of 0.5 mm and 1.2 mm were simulated in a Lagrangian coordinate system. The interactions between different particles and the effects of particles on flow field were ignored. In the PTV experiments, the beads were released outside of the airway model. In CFD simulations, the beads were initially released into the flow field at the first point on the PTV trajectory. However, since there were optical reflections at the wall surface, the determination of distance from the wall has an uncertainty level of around 0.5 mm. These inaccuracies in PTV measurements might have contributed to some uncertainty in the positions of these first points on the trajectories. To simulate such effects, particles were released from three hundred more locations around the initial location. These release points were roughly uniformly distributed within a circle of radius 0.5 mm and had the same vertical height as the initial location. The best-fit trajectory was found by visual comparison of the trajectories, along with the root-mean-square (RMS) distance between the experimental and predicted curves. In computing the RMS distance, cubic spline curves have been used to represent the PTV and CFD trajectories to ensure comparison at the same independent variables and in the same interval. Particles were assumed to have zero initial velocity in all directions. In one simulation, particles were also tracked assuming an initial axial velocity comparable to the settling velocity (5.4 mm/s) of the large bead (1.2 mm). It was found that the initial axial velocity of the beads had negligible effect on particle trajectory. This supported the use of release locations inside the CFD airway model to simulate particles released outside of the in vitro model.
For examining the flow features and for validation purposes, velocity field was checked in the symmetry plane (Fig. 3) and at twelve cross sections (Fig. 4). The flow features observed in both the PIV experiments and the CFD simulations were typical of low-Reynolds number flow, and were consistent with those seen in previous studies (Tsuda et al., 1994; Darquenne and Paiva, 1996; Karl et al., 2004; Harrington et al., 2006). In each generation of the model, the flow became fully developed within a very short distance from the entrance, such that most of the flow field can be regarded as fully developed (Figures 3 and and4).4). Velocity profiles in the bifurcation areas deviated from a parabolic shape (see section 7 in Fig. 4) due to expanding cross-sections and change in flow directions. Flow fields were characterized by a curvilinear separation streamline at the alveolar openings, indicating little convective exchange between the alveoli and the central duct lumen. In each alveolar cavity, slow recirculation zones could be identified. Velocities inside these cavities were about two orders of magnitude smaller than the mean lumen velocity.
The contour of velocity magnitude in the symmetry plane is shown in Fig. 3. The basic features were consistent with those obtained by PIV. Velocity profiles at twelve selected cross sections in the model were compared between the PIV measurements and the CFD predictions. The twelve cross sections are defined in Fig. 3 and the velocity profiles are shown in Fig. 4. In Fig. 4, the left side of the velocity profile corresponds to the left side of an observer going downstream with the flow. Overall, the CFD computed velocity profiles agreed well with those measured in the experiments. Compared with the PIV data, the maximum velocity predicted by CFD differed by a maximum of 3.9%, a minimum of 0.05%, and an average of 1.7% (Table 1). At some points, mostly near the wall (e.g. Fig. 4, sections 2, 5, 8, 12), the deviations were larger compared to those near the tube center. Correspondingly, the signal-to-noise ratios in the images at these bad-fit points were smaller than those at the good-fit points. The signal-to-noise ratio in PIV was defined as the ratio between the highest and second highest peak and was a heuristic for the measurement reliability. This suggested that the deviations were due to limitations of experimental precisions rather than inaccuracy in CFD results.
All four 1.2 mm beads that were released in the symmetry plane and the two 0.5 mm beads that were released outside the symmetry plane deposited on the airway walls after the first bifurcation, while the 0.5 mm bead that was released in the symmetry plane exited the model without deposition. Projections of these trajectories in the symmetry plane are shown in Figures 5, ,66 and and7.7. In Figures 8 and and9,9, the CFD predicted trajectories of these beads (filament-sphere lines) were compared to streamlines (solid lines) passing through the release points. After the first bifurcation, all bead trajectories started to deviate from the streamlines due to effects of gravity. The five beads initially released in the symmetry plane remained in the symmetry plane before hitting the walls (deposited) or exiting from the model. For the two beads released out of the symmetry plane, there were some wavy motions in the out-of-plane direction. However, the magnitude of such motions was comparable to the level of uncertainty in the PTV measurements. Therefore, motions perpendicular to the symmetry plane were not compared.
By setting the release locations to be the first measured trajectory points, the numerically simulated trajectories sometimes deviated appreciably from the experimental data. Better fits to the experimental trajectories could always be found by selecting certain release locations that were within a small distance (from 0.14 mm to 0.45 mm) from the initial location (Figures 5, ,66 and and7).7). This required distance was always within the level of the experimental uncertainty which was estimated to be 0.5 mm. In the CFD simulations, when a particle hit the wall of the model, it was assumed to deposit on the wall and its tracking was stopped. In the PTV experiment, when the bead hit the wall it did not stop but rather rolled slowly along the wall. This created some differences between CFD predictions and PTV measurements towards the end of the trajectories (Fig. 5).
Validation is an important part of any modeling study. In the case of simulating aerosol particle transport in the alveolated airways, such validation has been scarce in the past. The current study used a scaled-up in vitro model to validate the CFD predictions with PIV and PTV experimental data. While this is still an in vitro model, such data provided some valuable information regarding the accuracy of numerical simulations of fluid flow and particle transport in small airways and of low Reynolds number flow in general. Oldham (2006) recently summarized some of the issues related to validating CFD-derived inhaled aerosol deposition predictions, such as lack of a consensus for what is a validated CFD model, difficulties in creating the exact desired airway geometry, and difficulties in controlling and modeling all important physical factors involved in particle deposition such as particle charge, etc. Although the issues were discussed in the context of large airways, it is also applicable to cases in alveolated airways. In the current study we have restricted our efforts to comparing CFD predicted velocity profiles and particle trajectories with those measured from experiments in an identical geometrical model.
The fast re-establishment of parabolic velocity profile after a disturbance (bifurcation) is consistent with a typical low Reynolds number viscous flow, and support the use of flat velocity profile at the model entrance. Indeed, for laminar flow in a straight tube, the entrance length LE, which is a measure of the distance required by the flow to be fully developed, is expressed by (Munson et al., 2006):
where D is the lumen diameter and Re is the flow Reynolds’ number at the entrance. At the inlet of the model in this study, LE/D is 0.0078 and LE is 0.156 mm. Thus velocity profile will be fully developed 0.156 mm from the inlet. Therefore, it is unlikely that the use of any other velocity profiles such as a parabolic profile at the entrance of the model will lead to different results. The independence on initial particle velocity showed that the system had a very short response time and that the particles reached final equilibrium velocity quickly after being released into the flow field. This is consistent with the very short relaxation time (5.5e-4 second for 1.2 mm bead and 9.5e-5 second for 0.5 mm bead) of these beads in the experiments.
To obtain the velocity profiles at the twelve cross sections (Fig. 4), a uniform distribution of 80 sampling points along the cross section were used, thus the distance between sampling points was equal. To obtain the flow field values at these arbitrarily selected points, a linear interpolation based on the values at nearby CFD mesh cell centers were used. At some points the slope of the velocity-radial distance curve was larger than those at nearby points. Therefore, the equal distance of sampling points in space did not resulted in equal distance on the velocity-distance curve, creating visual discontinuities in some velocity profiles. Such discontinuities did not have any physical significance but rather were due to the discrete nature of the data.
The particle sizes used in this study were chosen based on practical experimental limitations, such as availability of materials. Nevertheless, the interpretation of these experimental data and the corresponding CFD simulations, and the extrapolation to the in vivo system are not necessarily restricted to the assumptions made in designing the experiments. Previous analysis (Theunissen et al., 2006) showed that it is hard to match both Stokes number and velocity ratio (particle settling velocity to mean flow velocity) between in vivo and in vitro systems. In the current study the in vitro particle sizes were selected based on similarity of particle Reynolds number and velocity ratio. By these criteria, the corresponding in vivo aerosol size is 14.6 μm for 1.2 mm bead and 6.1 μm for 0.5 mm bead. Aerosol particle of 14.6 μm is very unlikely to occur in the pulmonary region due to impaction and sedimentation in the conducting airways. Including this size in the current study is mainly for validation purposes, and although the current study was designed in the context of pulmonary aerosol transport, the data themselves may not necessarily be restricted to this specific field.
By considering the Stokes number the current data may also be interpreted differently. In the experimental system the Stokes number (St = ρpdp2U/(18μD), U is mean flow velocity, and D is tube diameter) is about 2.2e-4 for the 1.2 mm bead, and 3.8e-5 for the 0.5 mm bead. Using data based on the 20th generation of the Weibel model, the Stokes numbers for 14.6 μm and 6.1 μm aerosol particles are 0.02 and 3.5e-3, respectively, which are two orders of magnitude larger than those of the beads. For the same Stokes numbers, the in vivo aerosol sizes are 1.7 μm and 0.7 μm, respectively. Therefore, although the in vitro particle sizes were initially selected based on similarity of particle Reynolds number and velocity ratio, by Stokes number similarity the corresponding in vivo aerosol sizes are smaller and fall within the size of interests for respiratory aerosols. According to the Stokes-Einstein relation, the diffusion coefficient (Ddf = KbT/(3πμdp), Kb is Boltzmann’s constant, T is absolute temperature) for 0.5 and 1.2 mm beads in silicone oil were 8.8e-19 and 3.7e-19 m2s−1, respectively. Therefore, diffusion was negligible for these beads and was not considered in CFD simulations.
Kim et al. (1998) studied the motion of a freely moving sphere injected into an initially stationary or oscillating fluid for particle Reynolds number ranging between 2 and 150 and particle to fluid density ratio ranging between 5 and 200. Results showed that the relative magnitude of the forces acting on the sphere depends strongly on the particle to fluid density ratio. When the density ratio is much larger than one, the total drag force can be well approximated by the steady viscous drag force; however, when the density ratio is small (e.g. 5), the particle history force and force due to initial velocity difference between the particle and the carrier fluid become important so that the total drag force differs slightly from the steady viscous drag force in both magnitude and phase. In the current study the history force and the force due to initial velocity difference were ignored, although the particle to fluid density ratio (7800/970 = 8.04) is much smaller than that for realistic aerosol (around 1000). Since the flow Reynolds number (<0.13) and particle Reynolds number (<0.008) in the current study were significantly smaller than those in the study of Kim et al. (1998), and based on the observation that the particle trajectory was insensitive to initial particle velocity, it seems that these ignored forces did not significantly affect the particle trajectory in this study.
Almost all previous experimental validation studies of particle deposition in the large human airways compared the total or averaged deposition efficiency using a large amount of particles. In contract, the current study has tried to directly validate the trajectory of each individual particle. The current study was designed to evaluate the accuracy of using CFD to simulate airflow and aerosol transport in pulmonary region of the human lung, using an in vitro model by dynamic similarity principle. The data in this study showed that CFD simulations can provide reasonably accurate predictions of flow field and particle trajectory for predominantly viscous flows in the in vitro model. It is expected that comparable accuracy can be obtained for CFD simulations in realistic models of the human alveolated airways using similar numerical techniques.
This work was supported by NIH grant RO1 ES011177.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.