Effectiveness of Direct Lagrangian Tracking Models for Simulating Nanoparticle Deposition in the Upper Airways

Abstract

Direct Lagrangian particle tracking may provide an effective method for simulating the deposition of ultrafine aerosols in the upper respiratory airways that can account for finite inertia and slip correction effects. However, use of the Lagrangian approach for simulating ultrafine aerosols has been limited due to computational cost and numerical difficulties. The objective of this study is to evaluate the effectiveness of direct Lagrangian tracking methods for calculating ultrafine aerosol transport and deposition in flow fields consistent with the upper respiratory tract. Representative geometries that have been considered include a straight tubular flow field, a 90° tubular bend, and an idealized replica of the human oral airway. The Lagrangian particle tracking algorithms considered include the Fluent Brownian motion (BM) routine, a user-defined BM model, and a user-defined BM model in conjunction with a near-wall interpolation (NWI) algorithm. Lagrangian deposition results have been compared with a chemical species Eulerian model, which neglects particle inertia, and available experimental data. Results indicate that the Fluent BM routine incorrectly predicts the diffusion-driven deposition of ultrafine aerosols by up to one order of magnitude in all cases considered. For the tubular and 90° bend geometries, Lagrangian model results with a user-defined BM routine agreed well with the Eulerian model, available analytic correlations, and experimental deposition data. Considering the oral airway model, the best match to empirical deposition data over a range of particle sizes from 1 to 120 nm was provided by the Lagrangian model with user-defined BM and NWI routines. Therefore, a direct Lagrangian transport model with appropriate user-defined routines provides an effective approach to accurately predict the deposition of nanoparticles in the respiratory tract.

1. INTRODUCTION

The inhalation of ultrafine aerosols, characterized by a diameter of 100 nm or less, may present a significant health risk (Hood 2004; Kreyling et al. 2004; Maynard et al. 2004). Current studies indicate that fine (< 1 μ m) and especially ultrafine particulate matter may be more biologically active and toxic than larger size particles of the same material (Oberdorster and Utell 2002; Li et al. 2003; Kreyling et al. 2004). Common sources of these nanoparticles include ambient atmospheric reactions, radioactive decay, combustion, and industrial processes. For example, diesel exhaust particles in the nuclei mode are primarily between 5 and 50 nm, based on the particle number fraction (Kittelson 1998). Cigarette smoke particles are typically considered to have a mass mean diameter ranging from 250 to 350 nm (Keith 1982) with a mean geometric diameter of 170 nm (Ueno and Peters 1986). However, cigarette smoke particles as small as 18 nm have been observed (Bernstein 2004). Furthermore, there is an increasing concern regarding occupational exposures to manufactured nanoparticles generated in industrial processes (Hood 2004; Maynard et al. 2004).

Due to small particle sizes and low Stokes numbers, ultrafine aerosols are typically considered to deposit in the lung as a result of Brownian motion and turbulent dispersion. In the upper airways, both Brownian motion and turbulent dispersion may significantly contribute to nanoparticle deposition. In the lower airways, turbulence is usually not present and residence times are significantly higher. Therefore, Brownian motion is a primary mechanism responsible for ultrafine aerosol deposition in the alveolar region, provided that electrostatic effects are absent. Studies that have modeled the localized deposition of ultrafine aerosols in the upper respiratory airways through approximately the third respiratory generation include Balashazy and Hoffmann (1993), Broday (2004), Hofmann et al. (2003), Martonen et al. (1996), Moskal and Gradon (2002), Moskal et al. (2006), Zamankhan et al. (2006), Zhang and Kleinstreuer (2003b; 2004), and Zhang et al. (2005). Diffusional deposition of nanoparticles due to Brownian motion and convective diffusion in the lower airways has also been widely considered (Asgharian and Anjilvel 1994; Darquenne 2002). Other studies that have focused on convective dispersions associated with wall motion in the deep lung include Haber et al. (2003), Lee and Lee (2003), and Tsuda et al. (2002).

Numerical studies of fine and ultrafine aerosol deposition have typically employed either a Eulerian-Eulerian (Ingham 1975, 1991; Martonen et al. 1996; Lee and Lee 2002; Broday 2004; Shi et al. 2004; Zhang et al. 2005; Shi et al. 2006) or a Eulerian-Lagrangian (Moskal and Gradon 2002; Balashazy et al. 2003; Hofmann et al. 2003; Moskal et al. 2006; Robinson et al. 2006; Zamankhan et al. 2006) model. The Eulerian-Eulerian particle transport model treats both the continuous and discrete phases as interpenetrating fields. Based on the dilute concentrations and low Stokes numbers of inhaled nanoparticles, the Eulerian transport model typically assumes that the effects of the particles on the flow field and particle inertia can be neglected. The particles are then modeled as a dilute chemical species in a multi-component mixture model (Crowe et al. 1998). At the wall, the deposition flux is approximated by assuming either a zero concentration or a concentration dependent reaction (Longest and Kleinstreuer 2003). The primary advantage of this Eulerian model is that it can easily simulate large particle counts and low deposition rates. However, significant modifications are necessary to address inertial (Fernandez de la Mora and Rosner 1981; Friedlander 2000), slip correction, and electrostatic (Wang et al. 2002) effects.

In contrast to Eulerian simulations, the Lagrangian transport model tracks individual particles within the flow field. The primary advantage of this method is that a variety of forces such as inertia, diffusivity, electrostatic effects and near-wall terms can be considered directly (Longest et al. 2004). For respiratory airways, particle deposition is assumed to occur at initial wall contact. However, the Lagrangian model requires the simulation of a large number of particles to generate convergent deposition profiles. As a result, the Lagrangian model may become highly inefficient for cases of low particle deposition rates. Furthermore, extremely large particle counts are required to fully resolve local deposition characteristics (Longest and Kleinstreuer 2003).

The transport of ultrafine aerosols is generally considered to be governed by the Langevin equation (Fuchs 1964; Crowe et al. 1998). For the solution of particle trajectories, this equation can be directly integrated numerically with the inclusion of a randomly generated Brownian motion force term (Li and Ahmadi 1992) or a superimposed particle displacement (Asgharian and Anjilvel 1994; Balashazy 1994; Longest et al. 2004). As an alternative to direct numerical integration of the Langevin equation, Podgorski (2001) has developed a Brownian dynamics approach for the simulation of ultrafine aerosols near solid boundaries. In the particle tracking method of Podgorski (2001), the Langevin equation is integrated analytically by including the random force terms specified by Chandrasekhar (1943). The resulting expression defines a probability density function for the change in particle velocity and position due to convection and diffusion. Individual particle trajectories are calculated by applying random variants and standard deviations for velocity and position to expected values at each particle location. Moskal et al. (2006) have applied this approach for the simulation of fractal-like aerosol aggregates in a model of the human nasal cavity. An advantage of the Brownian dynamics approach reported by Podgorski (2001) is that the effects of wall boundaries can be directly included in the equations. For systems in which changes in the flow field are relatively slow in comparison with the particle response time, the Brownian dynamics approach can be simplified (Fuchs 1964; Podgorski 2001). In this scenario, the Langevin equation can be directly integrated to find particle positions as specified in other studies (Li and Ahmadi 1992; Asgharian and Anjilvel 1994; Balashazy 1994).

Cohen and Asgharian (1990) have reported that inertial effects begin to influence particle deposition at Stokes numbers on the order of approximately 1 × 10^{− 5} and higher. Aerosols on the order of 100 nm in the oral airway are within this inertial range at a tracheal flow rate of 60 L/min, based on the average velocity and minimum tracheal diameter as suggested by Cheng et al. (1999). Moreover, Longest and Xi (2007) have shown that finite particle inertia may significantly influence the local deposition of particles as small as 40 nm in a double bifurcation geometry of respiratory generations G3–G5. Modifications to the standard chemical species Eulerian model are possible to account for the effects of finite particle inertia (Fernandez de la Mora and Rosner 1981; Friedlander 2000; Loth 2000). Alternatively, the Lagrangian particle trajectory approach provides a convenient and direct method to account for inertial effects and particle interaction factors. Moreover, the Lagrangian model can be directly applied to fine (100 nm–1 μ m) and coarse (> 1 μ m) size particles where inertia and gravity become increasingly important in relation to Brownian diffusion.

Lagrangian simulations of ultrafine aerosol transport are complicated by a number of factors. The very small response time of ultrafine aerosols prevents the application of standard explicit solvers such as the Runge Kutta algorithm (Ferziger and Peric 1999). Instead, an implicit solution procedure can be implemented to integrate the particle trajectory equation (Press et al. 1996). We have performed initial simulations of direct Lagrangian ultrafine aerosol tracking in an oral airway model using the commercial software Fluent 6 (Ansys, Inc.). Results of this initial study showed that predictions of the commercial code differed from the empirical data of Cheng et al. (1993, 1997a) by up to one order of magnitude and did not predict a decreasing trend in ultrafine deposition that is expected with increasing particle size. Further analysis indicated potential problems with the Brownian motion term and interpolation of fluid velocities near the boundary surface.

The objective of this study is to evaluate the effectiveness of direct Lagrangian tracking methods for calculating ultrafine aerosol transport and deposition in flow fields consistent with the upper respiratory tract. Emphasis will be placed on both regional-averaged and local deposition values. To achieve this objective, three geometries have been considered based on the availability of existing analytic and empirical deposition data for ultrafine aerosols. These representative cases include a fully developed flow field in a tubular geometry with a diameter equivalent to the fourth respiratory generation, a 90° tubular bend, and an elliptic geometry of the human oral airway. The characteristic Reynolds numbers of these systems range from 280 to 3640, which requires the consideration of laminar through turbulent flows. Analytic and experimental studies that have provided deposition data in these systems include Ingham (1975), Wang et al. (2002) and Cheng et al. (1993, 1997a). Comparisons to analytic and experimental ultrafine deposition results will be made for the chemical species Eulerian model, Fluent's ultrafine particle tracking algorithm, and user-defined functions for both Brownian motion and near-wall particle conditions. For very small nanoparticles, where inertia is largely absent, deposition results for the Eulerian and Lagrangian models should be nearly identical on a regional and local basis. While laminar and turbulent airflow has been simulated for the case of the oral airway, the effect of turbulent dispersion on Lagrangian particles has been neglected to better isolate Brownian diffusion.

2. METHODS

2.1. Flow Systems and Boundary Conditions

In order to develop and evaluate an effective Lagrangian model for ultrafine aerosol deposition in the upper airways, three idealized flow systems have been considered. These representative geometries have been selected based on the availability of analytic and experimental deposition data and are summarized in Table 1. The first model considered is fully developed laminar flow in a constant diameter axisymmetric tube (Figure 1a). The tubular diameter selected is equivalent to the fourth generation airway of an adult male based on the ICRP (1994) model (Table 1). The numbering of respiratory generations in this study is based on Weibel's convention (Weibel 1963) with the trachea representing generation G0 and the main bronchi denoted G1. For this system, inhalation tracheal flow rates of 15, 30, and 60 L/min will be considered, which correspond to sedentary, light activity and heavy activity states of exertion. The corresponding inlet flow rates and Reynolds numbers for respiratory generation G4 are provided in Table 1.

FIG. 1 Geometries considered including (a) a fully developed tubular flow field, (b) a 90° tubular bend, and (c) an idealized model of the adult oral airway. The oral airway geometry has been based on CT-imaged data and has been idealized with the use of elliptical cross-sectional segments.

TABLE 1 Dimensions and Reynolds numbers of the flow systems considered

Model	Diameter range (cm)	Radius of curvature (cm)	Inlet flow rate (L/min)	Inlet Reynolds number
Tubular fully developed flow	D = 0.45	R = 0.0	0.9375	278
			1.875	556
			3.75	1113
90° tubular bend	D = 0.46	R = 1.43	1.052	305
Oral airway model	D = 2.2–1.0	R = 3.27	15	910
			30	1820
			60	3640

Whole lung models of ultrafine aerosol deposition typically implement analytic mass transfer expressions to approximate the effect of Brownian diffusion. These analytic expressions treat each airway as a tubular segment and assume either a parabolic or blunt inlet profile. Most analytic expressions for Brownian diffusion include a dimensionless diffusion parameter defined as

where L is the length of the geometry and R is the radius. For tubular flow with a fully developed parabolic profile, Ingham (1975) suggested the following expression for deposition efficiency (DE) based on a steady Eulerian model

This expression will be used to evaluate the performance of the Eulerian and Lagrangian models applied to the fully developed tubular case.

Algebraic whole lung models often approximate inertial deposition in airway bifurcations by applying analytic correlations for particle impaction in curved tubes (Yeh and Schum 1980; Martonen 1993). In this study, a 90° tubular bend geometry (Figure 1b) has been selected based on the ultrafine deposition study of Wang et al. (2002). Flow rate and geometry conditions for this system are identical to those used in the experimental study of Wang et al. (2002), which investigated deposition losses in particle delivery systems. However, these conditions are very similar to the fully developed tubular case of respiratory generation G4 at the low flow rate. For example, the fully developed tubular flow and 90° bend diameters (0.45 vs. 0.46 cm) and inlet flow rates (0.9375 vs. 1.052 L/min) are very similar (Table 1).

The third respiratory geometry considered is an approximate model of the human oral airway. To generate this geometry, a virtual 3-D airway replica was created from CT images of a healthy adult, as described in Xi and Longest (2007). Image translation from CT format to a solid model was performed using MIMICS software (Materialise, Ann Arbor, MI). The surface geometry was then imported into ICEM (Ansys, Inc.) as an IGES file for surface modifications and meshing. Surface geometries for the oral cavity have been based on measurements of a dental impression with an approximate half mouth opening (Cheng et al. 1997b). The oral cavity and airway replica models were then assembled. As a simplification, cross-sectional segments of the resulting model have been reconstructed with ellipses of consistent hydraulic diameter and flow area (Xi and Longest 2007). Key geometric features such as the mouth entrance, larynx, and trachea have been rescaled to agree with data reported for in vitro deposition studies. The dimensions adopted in this study are from Cheng et al. (1997b), resulting in a mouth opening of 2.2 cm in diameter, a glottis of 0.87 cm² in cross-sectional area, and a trachea of 2.0 cm² in cross-sectional area (Figure 1c). In this study, nanoparticle deposition results at a tracheal flow rate of 15 L/min will be compared to the ultrafine in vitro data of Cheng et al. (1993, 1997a). Results for light and heavy activity conditions will also be considered (Table 1).

Pre-specified inlet velocity profiles have been selected based on the analytic or experimental system considered for each case. For fully developed tubular flow and the 90° bend, laminar parabolic velocity profiles have been implemented. Considering that the inhalation of ambient air will likely result in a nearly constant velocity profile, a blunt inlet condition has been employed for the oral airway model, which is defined as

In the expression above, r is the inlet radial coordinate, u _m is the mean velocity and R is the outer radius of the inlet. This profile is similar to a constant velocity inlet, but provides a smooth transition to the no-slip wall condition. For inlet flow rates of 15, 30, and 60 L/min, the initial turbulence intensities have been assumed to be 0.5, 1, and 2%, respectively, with an inlet turbulence length scale of l = 0.07D _inlet .

Inlet particle profiles for both the Eulerian and Lagrangian models have been specified to be consistent with the local mass flow rate (Longest and Vinchurkar 2007a). For Lagrangian particles, the local inlet mass flow rate on a finite ring, _{in, ring} , is given by

where r ₁ and r ₂ define the extent of a local ring and u _in (r) is the inlet velocity profile. The mixture density ρ _m and inlet mass fraction c _in are assumed to be constants. Initial particle velocities have been assumed to match the local fluid velocities. The inlet particle profile has been generated by applying Equation (4) to a series of 200 finite rings that encompass the entire cross section. The resulting particle profile is proportional to the inlet mass flow rate of the fluid in the absence of gravity effects. For the Eulerian model, the inlet mass fraction and mixture density have also been specified as constants. The inlet mass flow rate of particles is then determined by integrating Equation (4) over the inlet area A _in .

2.2. Transport Equations

Flows in all the geometries considered are assumed to be isothermal and incompressible. Based on the inlet Reynolds numbers, laminar conditions are expected for the fully developed tubular flow and 90° bend models (Table 1). The inlet Reynolds number for the oral airway ranges from 910 to 3,640. Moreover, the maximum Reynolds number based on conditions at the minimum cross-sectional area in the oral airway model is approximately 8,000 for the highest flow rate considered. Therefore, laminar, transitional and fully turbulent conditions are expected in the oral airway. To address these multiple flow regimes, the LRN k-ω model was selected based on its ability to accurately predict pressure drop, velocity profiles and shear stress for transitional and turbulent flows (Ghalichi et al. 1998; Wilcox 1998). This model has also been demonstrated to accurately predict particle deposition profiles for transitional and turbulent flows in models of the oral airway (Zhang and Kleinstreuer 2003a; Zhang and Kleinstreuer 2004) and multiple bifurcations (Longest and Vinchurkar 2007b). Moreover, the LRN k-ω model has been shown to provide an accurate solution for laminar flow as the turbulent viscosity approaches zero (Wilcox 1998).

For laminar and turbulent flow, the Reynolds averaged equations governing the conservation of mass and momentum are (Wilcox 1998)

where u _i is the time-averaged velocity in three coordinate directions, i.e., i = 1, 2, and 3, p is the time-averaged pressure, ρ is the fluid density, and ν is the kinematic viscosity. Overbars have not been included on time-averaged quantities to simplify the equations. The turbulent viscosity ν _T is defined as ν _T = α* k/ω · For the LRN k-ω approximation, which models turbulence through the viscous sublayer, the α * parameter in the above expression for turbulent viscosity is evaluated as (Wilcox 1998):

For laminar flow, ν _T is zero and only Equations (5a) and (5b) are solved.

Transport equations governing the turbulent kinetic energy (k) and the specific dissipation rate (ω) are

In the above equations, τ _ij is the shear stress tensor; ϵ _k and ϵ_ω represents the dissipation of k and ω, respectively (Wilcox 1998).

Dilute suspensions of nanoparticles are often treated as a chemical species in a Eulerian mass transport model. This model often neglects particle inertia and the effects of the particulate phase on the flow field, i.e., one-way coupled particle transport. The mass transport relation governing the convective-diffusive motion of ultrafine aerosols in the absence of particle inertia effects can be written on a mass fraction basis as

In the above equation, c represents the mass fraction of nanoparticles, is the molecular or Brownian diffusion coefficient and Sc _T is the turbulent Schmidt number taken to be 0.9. Assuming dilute concentrations of spherical particles, the Stokes-Einstein equation has been used to determine the diffusion coefficients for various sized particles and can be expressed as

where k _B = 1.38× 10^{− 16} cm²g/s is the Boltzmann constant in cgs units. The Cunningham correction factor has been computed using the expression of Allen and Raabe (1985)

where λ is the mean free path of air, assumed to be 65 nm. The above expression is reported to be valid for all particle sizes (Hinds 1999). To approximate particle deposition on the wall, the boundary condition for the Eulerian transport model is assumed to be c_wall = 0.

One-way coupled trajectories of monodisperse ultrafine particles ranging in diameter (d _p ) from 1 to approximately 100 nm have been calculated on a Lagrangian basis by directly integrating an appropriate form of the particle trajectory equation. Aerosols in this size range have very low Stokes numbers (St = ρ _p d _p ² C _c U/18μ D ≪ 1), where U is the mean fluid velocity and D is a characteristic diameter of the system. Other characteristics of the aerosols considered include a particle density ρ _p = 1.00 g/cm³, a density ratio α = ρ/ρ _p ≈ 10^{− 3}, and a particle Reynolds number Re_p = ρ | u−v | d _p /μ ≪ 1. The appropriate equations for spherical particle motion under these conditions can be expressed

and

In the above equations, v _i and u _i are the components of the particle and time-averaged local fluid velocity, respectively. The characteristic time required for particles to respond to changes in the flow field, or the particle response time, is τ _p = C _c ρ _p d _p ²/18 μ · The acceleration term is often neglected for aerosols due to small values of the density ratio (α ≪ 1). However, it has been retrained here to emphasize the significance of fluid element acceleration in biofluid flows (Longest et al. 2004). The drag factor f, which represents the ratio of the drag coefficient C _D to Stokes drag, is assumed to be one for submicron aerosols. The effects of gravity, lift and near-wall boundaries have not been included.

Neglecting the acceleration term (α = 0) and setting f = 1 reduces Equation (11a) to the general Langevin equation. Solution of Equation (11a) or the more general Langevin equation can be accomplished by the Brownian dynamics approach (Podgorski 2001), solution of the corresponding Fokker-Plank equation (Peters 1999), or by direct numerical integration (Li and Ahmadi 1992). In this study, it will be assumed that the particle response time (τ _p ) is much smaller than temporal changes in the flow field and the associated particle time-step (Δ t). Based on this assumption, Equation (11a) can be directly integrated, which results in a Lagrangian particle trajectory. In this approach, the effects of random Brownian motion can be incorporated as either a random force (Li and Ahmadi 1992) or displacement (Asgharian and Anjilvel 1994; Balashazy 1994; Longest et al. 2004). To remain consistent with the general form of the Langevin equation, Brownian motion effects will be included as a force term in this study.

The effect of Brownian motion on the particle trajectories is included as a separate force per unit mass term at each time-step. The amplitude of the Brownian force is of the form (Li and Ahmadi 1992)

where ς _i is a zero mean variant from a Gaussian probability density function, S _o is a spectral intensity function directly related to the diffusion coefficient, and Δ t is the time-step for particle integration. This expression can be re-written to highlight the diffusion coefficient as

where m _d is the mass of the particle. This expression for Brownian diffusion is theoretically valid for both laminar and turbulent flows.

For a number of particles diffusing over time, the root-mean-square displacement distance in one dimension (x _rms ) should equal the standard deviation of the Gaussian probability density function σ_{1 − D} (Hinds 1999). That is,

where N is the number of particles sampled and t is the total time available for diffusion. In two-dimensions, the radial distance traveled by a particle normal to the direction of axial flow can be approximated as

2.3. Deposition Factors

Regional deposition factors account for total deposition within the geometries of interest. For the Lagrangian model, deposition factors are based directly on the number of discrete particles that deposit. The deposition efficiency (DE) for region i can be expressed as

In this study, the fully developed tubular flow, 90° bend, and oral airway geometries will each be considered as unit regions of interest.

For the Eulerian model, regional deposition efficiency is based on the ratio of the total mass deposition rate to the inlet mass flow rate,

Assuming constant property flow, the total mass deposition rate on a wall can be expressed as

where the summation is performed over region of interest i, and n is the wall normal coordinate pointing out of the geometry. The inlet mass flow rate of particles ( _in ) results from the integration of Equation (4) over the inlet area A _in .

In order to capture local deposition profiles at the sub-branch level for Lagrangian particles, a deposition enhancement factor (DEF), similar to the enhancement factor suggested by Balashazy et al. (1999), for local region j can be defined as

where the summation is performed over the i-region of interest. In the above expression, A _crit,j is a local critical area with an assumed diameter of 500 μ m or approximately 50 lung epithelial cells in length (Lumb 2000). This prescribed area is similar to the localized areas considered by other researchers (Balashazy et al. 1999; Balashazy et al. 2003; Zhang et al. 2005) and localized in vitro experiments (Oldham et al. 2000). For all deposition factors considered, the subscript i will be used to denote a regional area such as the oral airway model. The subscript j will be used to denote a local area such as A _crit,j . Therefore, DEF _j is the local microdosimetry factor in region A _crit,j . Sampling locations are taken to be nodal points.

For Eulerian model predictions of local deposition values, the local mass flow rate of particles at the surface is

The local DEF factors can then be expressed as

In this Eulerian expression for deposition enhancement, the local region A _{crit ,j} is identical to the value used in the Lagrangian approach, which enables direct comparisons between the results of the two models.

2.4. Numerical Methods

To solve the governing mass and momentum conservation equations in each of the cases considered, the CFD package Fluent 6.2.16 (Ansys, Inc.) has been employed. User-supplied Fortran and C programs have been implemented for the calculation of initial particle profiles, particle deposition factors, inlet and wall mass flow rates, grid convergence, and post-processing. All transport equations were discretized to be at least second order accurate in space. Meshes were constructed in either Gambit 2.2 (Ansys, Inc.) or ICEM (Ansys, Inc.) and consisted entirely of hexahedral control volumes. Convergence of the flow field solution was assumed when the global mass residual had been reduced from its original value by five orders of magnitude and when the residual-reduction-rates for both mass and momentum were sufficiently small. To ensure that a converged solution had been reached, residual and reduction-rate factors were decreased by an order of magnitude and the results were compared. The stricter convergence criteria produced a negligible effect on both velocity and particle deposition fields. To improve accuracy, all calculations were performed in double precision. Details of grid convergence for a similar bifurcating model system was reported in Longest and Vinchurkar (2007a). In addition, a grid convergence study was performed for the Lagrangian and Eulerian particle transport and deposition models. The criterion for mesh independence was a root-mean-square relative error of deposition values on the order of 1% for all particle sizes. Meshes that satisfied this convergence criterion for the fully developed tubular flow, 90° bend and oral airway geometry contained 120,000, 130,000, and 450,000 control volumes, respectively. These meshes have been employed for all reported results. Further details of mesh independence testing are described below.

Particle trajectories were calculated within the steady flow fields of interest as a post-processing step. The extremely small particle response time of ultrafine aerosols results in a very stiff trajectory equation that is difficult to solve numerically. As a result, standard explicit and predictor corrector methods such as the Runge Kutta algorithm are highly unstable when applied to nanoparticles. In this study, an implicit Eulerian routine has been applied to solve the particle trajectory equation (Ferziger and Peric 1999). The implicit Euler routine is stable for relatively large time-step-sizes; however, the primary limitation of this routine is that it is second order accurate. It is noted that method order determines the rate at which solution errors are reduced as the time-step is decreased. Therefore, highly accurate results can be achieved with a second order method provided that a sufficiently small time-step is employed.

To actively select an appropriate time-step throughout the solution procedure, an adaptive step-size control algorithm can be implemented (Press et al. 1996; Longest et al. 2004). Further details of this type of accuracy control routine have been described in Longest et al. (2004). A problem that arose with the application of this accuracy control routine in the Fluent 6 software was that random Brownian displacements were often larger than the error tolerance. As a result, Brownian motion was artificially damped. To avoid this problem, adaptive step size control was not applied for any of the Brownian motion models considered. Instead, solution time intervals were based on a constant number of time-steps per control volume. As shown in the results, it was found that 10 integration steps of the particle trajectory equation per control volume produced a time-step independent solution. The associated range of time-step sizes varied from approximately 1 × 10^{− 6} to 1 × 10^{− 5} s for particles from 1 to 100 nm.

In order to determine the number of particles required to produce convergent local deposition profiles, groups of 10,000 particles were tested. The number of groups tested for convergence in a geometry was increased until deposition efficiency values changed by approximately 1% and maximum local values changed by approximately 5%. A constant number of particles was then simulated for each geometry based on this convergence criterion. Due to varying deposition rates, the particle counts required to generate convergent deposition efficiencies for the tubular, 90° bend and oral airway models were 10,000, 30,000, and 60,000, respectively. Sample calculations for particle convergence in the oral airway are provided below. For the local DEF factor, 180,000 particles were used for the tubular and 90° bend geometries and 300,000 particles were implemented for the oral airway model.

2.5. Near-Wall Interpolation (NWI)

Evaluation of the particle trajectory equation requires that the mean fluid velocity u _i be determined at the particle location for each time-step. On a computational grid, determining the fluid velocity at the particle location requires spatial interpolation from control-volume centers or nodal values. Based on preliminary studies of deposition in the oral airway model using Fluent 6 with a LRN k-ω turbulence model, fluid velocities in wall adjacent control volumes were observed to maintain values approximately consistent with the control-volume center and not approach zero at the wall. For this study, a user-defined function has been developed that linearly interpolates the fluid velocity in near-wall control volumes from nodal values. Velocity values at the wall have been taken to be zero. Therefore, implementation of this user-defined near-wall interpolation (NWI) routine provides a linear estimate of fluid velocity at the particle location that goes to zero at wall boundaries. The interpolation of variables within control volumes from nodal values has been described in detail by Pepper and Heinrich (1992).

2.6. Comparison of Particle Transport Models

To evaluate ultrafine aerosol transport and deposition in the three geometries selected for this study, four particle transport models have been employed. For each model, Fluent 6 has been used to evaluate the flow field as described above. The first particle transport model considered is the chemical species Eulerian approach described by Equation (8). Results of this model will be compared to the Lagrangian results for both regional and local depositions. The second model considered employees the Lagrangian particle tracking algorithm of Fluent 6 with the Brownian motion (BM) option turned on and without additional modification. As an alternative, the third model considered replaces the Brownian motion routine of Fluent 6 with a user-defined expression of Equation (13) to model the diffusion process. The final Lagrangian model implements user-defined functions for both the Brownian motion and near-wall velocity interpolation, as described above.

Flow field simulations for the fully developed tubular and 90° bend geometries have been conducted for laminar conditions, i.e., ν _T = 0, based on the inlet Reynolds numbers (Table 1). To approximate transitional and turbulent conditions in the oral airway model, the turbulent form of the flow field equations has been employed. However, Lagrangian particles in the oral airway model have been simulated without the effect of turbulent dispersion. The resulting mean flow particle tracking approach has been implemented to isolate the effects of Brownian motion on ultrafine aerosols. Moreover, current turbulent dispersion models may be highly inaccurate for the two-equation isotropic turbulence assumptions employed in this study (Matida et al. 2004). For the Eulerian particle transport model in the oral airway geometry, particle deposition results have been determined with and without the turbulent dispersion term. Exclusion of the turbulent dispersion term from the Eulerian model allows for direct comparisons to the Lagrangian mean flow tracking deposition results.

2.7. Model Testing and Sensitivity Analysis

Grid convergence testing for the oral airway geometry based on the deposition of 5 and 100 nm aerosols is illustrated in Figure 2. Grids containing approximately 150,000–450,000 control volumes have been considered in conjunction with the Eulerian and user-defined BM and NWI models. The number of control volumes for the oral airway has been rounded to the nearest 25,000 cells. For both 5 and 100 nm particles, increasing the number of control volumes up to 350,000 is observed to decrease the overall deposition efficiency (Figure 2). Further increasing the number of control volumes beyond 350,000 appears to have a negligible, less than 1%, effect on the total deposition efficiency. Therefore, the oral airway mesh considered for this study to provide grid independent deposition results contained approximately 450,000 control volumes. In this model, the height of the first near-wall control volume was 0.05 mm resulting in a maximum y⁺ range of 1.69 to 3.41 for inlet flow rates varying from 15 to 60 L/min. Differences in Eulerian and user-defined BM model deposition values observed in Figure 2 are discussed in the results section.

FIG. 2 Evaluation of grid convergence for overall deposition efficiency in the oral airway model at 30 L/min with (a) 5 nm and (b) 100 nm particles.

Sensitivity of the total deposition efficiency to the number of particles simulated in the oral airway is illustrated in Figure 3. The deposition of 5 nm particles at an inhalation flow rate of 30 L/min has been simulated using the user-defined BM model with NWI. Due to the random nature of BM particle deposition, significant fluctuations in deposition efficiency are observed for the simulation of 10,000 particles over multiple runs (Figure 3a). As the number of simulated particles is increased, fluctuations in overall deposition efficiency decrease. The standard deviations of deposition for groups of particles ranging from 10,000–60,000 over ten different trials are shown in Figure 3b. As the number of particles is increased, the standard deviations of deposition are observed to decrease significantly. From 50,000 to 60,000 particles, the change in the mean deposition efficiency is observed to be below approximately 1%. Results were similar for other particle sizes tested. As a result, 60,000 particles have been simulated to determine overall deposition efficiencies in the oral airway model.

FIG. 3 Effect of particle count on overall deposition efficiency in the oral airway geometry at 30 L/min for 5 nm aerosols. (a) Fluctuations in deposition efficiency over ten runs and (b) standard deviation of fluctuations in deposition efficiency with mean values. Mean values of deposition efficiency for 50,000 and 60,000 particles differ by less than 1%.

3. RESULTS

3.1. Numerical Estimate of Particle Diffusion

An effective model of Brownian motion should produce particle displacements that are consistent with the standard deviation of the Gaussian distribution function. That is, the root-mean-square of particle displacement should satisfy Equation (14). To evaluate the effectiveness of the Fluent BM routine in comparison to the user-defined model employed in this study, the diffusion of 5 and 40 nm particles has been considered in the fully developed tubular model at an inlet flow rate of 1.875 L/min. For each particle size considered, groups of 1,000 particles have been released from the center of the geometry inlet. Because the flow field is fully developed, radial dispersion in this geometry can occur only as a result of Brownian motion. The root-mean square of radial displacement in two-dimensions has been calculated and compared to the expected standard deviation value over time, as illustrated in Figure 4a. For both 5 and 40 nm aerosols, the Fluent BM model under predicts the expected standard deviation value by approximately one order of magnitude. In contrast, the user-defined BM model provides very close agreement to the expected standard deviation value for both 5 and 40 nm aerosols over the particle residence time considered. As a result, deposition predictions of the Fluent model are expected to be significantly less than predictions of the user-defined model and the analytic deposition solution for this system.

FIG. 4 Comparisons of predicted radial root-mean-square displacements (r_rms) over time to theoretical standard deviation values (σ). (a) Effect of particle transport models on predicted r _rms for 5 and 40 nm particles. The Fluent BM model under predicts the expected standard deviation value by up to one order of magnitude. (b) Effect of the number of particle integration steps per control volume on predicted r_rms values.

The force term describing random Brownian motion, Equation (13), contains the time-step used for particle integration. As a result, the Brownian motion force may be influenced by the size of the time-step selected. The time-step used in calculating particle position is inversely related to the specified number of steps per control volume. The effect of varying the number of integration steps on the radial displacement of 5 nm particles is illustrated in Figure 4b. The system considered is fully developed flow in the tubular geometry with a single particle release point at the centerline, as applied in Figure 4a. Solutions with 1, 10, and 20 integration steps per control volume have been considered, which resulted in mean particle time-steps of 4.2 × 10^{− 5}, 5 × 10^{− 6}, and 1.25 × 10^{− 6} s, respectively. The number of steps per control volume appears to significantly affect the diffusion results for the Fluent BM model. Increasing the number of steps, which decreases the time-step size, is observed to increase the error between the Fluent predictions and the theoretical standard deviation value. This result is unexpected considering that decreasing the time-step size will typically decrease the error. For the user-defined BM model, some minor variations are observed between the predicted radial displacement and the theoretical standard deviation for the lower limit of one step per control volume. However, increasing the number of steps to approximately 10, which decreases the time-step size, provides an accurate match to the theoretical standard deviation value. Further increasing the number of steps, or decreasing the time-step, was not observed to significantly improve agreement with the theoretical value for the user-defined BM model (Figure 4b).

As a result of observations made for Figure 4b, it appears that BM predictions for the Fluent model are time-step dependent and become progressively less accurate as the time-step is decreased. In contrast, radial displacement predictions for the user-defined BM model are observed to approach the theoretical value as the time-step is reduced. Accurate predictions of radial displacement are observed for the user-defined model for 10 steps per control volume, which results in a mean time-step size of approximately 5 × 10^{− 6} s.

3.2. Deposition in the Fully Developed Tubular Flow Geometry

Considering the tubular flow system, deposition efficiency results based on the four particle transport models considered are illustrated in Figure 5. Results of the numerical models have been compared to the analytic deposition correlation of Ingham (1975), i.e., Equation (2), for fully developed tubular flow. The Eulerian model closely matches the expected analytical result for all particle sizes considered and tracheal flow rates of 15, 30, and 60 L/min (Figure 5a). Therefore, numerical diffusion does not appear to adversely influence the Eulerian results for the ultrafine particle sizes considered, provided that a sufficiently fine grid is used near the wall. In contrast to the Eulerian model, the Fluent BM routine significantly under predicts diffusional deposition at all flow rates considered (Figure 5b). The underestimation of the Fluent BM routine ranges from approximately one order of magnitude at 5 nm to approximately 100% for 100 nm particles. The user-defined BM routine is shown to provide a significant improvement to the Fluent model (Figure 5c). Results for the user-defined BM model appear to be in close agreement with the analytic solution for all particle sizes and flow rates considered. Minor variations are observed at the highest flow rate for 100 nm aerosols (Figure 5c). Including the user-defined NWI routine (Figure 5d) does not appear to provide a significant improvement to the Lagrangian simulation in comparison to the user-defined BM result (Figure 5c) for the case of fully developed tubular flow.

FIG. 5 Deposition efficiency values in the tubular fully developed flow field for the (a) chemical species Eulerian, (b) Fluent BM, (c) user-defined BM, and (d) user-defined BM and NWI models in comparison to the analytic results of Ingham (1975).

Local deposition results for 5 nm aerosols in the fully developed tubular flow system are shown in Figure 6. Deposition results based on the Eulerian and Lagrangian user-defined BM model appear similar in pattern and magnitude (Figure 6a vs. Figure 6b). Both the Eulerian and Lagrangian models show a maximum DEF value at the inlet due to the specification of a constant particle concentration at this location. Local DEF values are then gradually reduced in the downstream direction and approach a mean of approximately one. Maximum DEF values are also similar between the Eulerian (DEF_max = 11.5) and Lagrangian (DEF_max = 9.8) models. These comparisons indicate that the Lagrangian model with the user-defined BM routine is adequately capturing local diffusional deposition in comparison to the Eulerian base case. However, differences are apparent due to the underlying mechanism of each approach. Results for the Eulerian model appear continuous and vary at specific locations in the axial direction. In contrast, the Lagrangian model results are somewhat discontinuous even for the relatively high count of 180,000 particles employed in this simulation. Nevertheless, the approximate maximum value and trend of the local DEF variable appear to be effectively captured by the Lagrangian model with the inclusion of the user-defined BM routine.

FIG. 6 Local DEF values in the tubular flow geometry for the (a) chemical species Eulerian and (b) Lagrangian user-defined BM models. Predictions of DEF values appear similar between the two models considered.

3.3. Deposition in the 90° Bend Geometry

Total deposition results for the four-particle transport models considered applied to the 90° bend geometry are shown in Figure 7. These numerical results have been compared to the experimental data of Wang et al. (2002) for particle sizes ranging from 5–12 nm. The original deposition results of Wang et al. (2002) have been converted from penetration with respect to tubular flow to deposition efficiency based on the deposition expression of Gormley and Kennedy (1949). The Eulerian model predictions are observed to be within the uncertainty bounds of the experimental data for most particle sizes considered (Figure 7). However, the Eulerian result under predicts the mean experimental data for all particle sizes. In contrast to the Eulerian model, the Fluent BM routine significantly under predicts the experimental data (Figure 7). The Fluent BM results are lower than the experimental values by a factor of approximately six. The Lagrangian results are significantly improved with the user-defined BM model (Figure 7). Results of this model are in close agreement with the mean experimental values for all particles sizes except for 12 nm. Still, the user-defined BM result at 12 nm is within the bounds of experimental uncertainty. Interestingly, the Lagrangian user-defined BM results are noticeably higher than the Eulerian model predictions. Differences in the Eulerian and Lagrangian models may potentially arise from inertia and slip correction effects. For this flow system, the inertia of 12 nm and smaller aerosols is expected to be negligible. However, the Lagrangian model does include slip correction effects in the drag term of individual particles. The resulting difference between particle and fluid velocities may be responsible for the enhanced deposition observed for Lagrangian particles and the associated better match to experimental results. Considering the fourth particle transport model, the inclusion of user-defined NWI does not appear to significantly impact deposition for the 90° bend geometry.

FIG. 7 Comparison of model deposition efficiency results with the empirical deposition data of Wang et al. (2002) for the 90° bend geometry. Of the models considered, the user-defined Lagrangian simulations provide the best agreement with the experimental data.

Local deposition results in the 90° bend geometry for 5 nm aerosols are illustrated in Figure 8. The Lagrangian results shown in Figure 8b are based on the user-defined BM model. As with the fully developed tubular flow cases, local deposition results in terms of the DEF variable appear similar between the Lagrangian and Eulerian models. Both models predict an expected mean DEF value of approximately one upstream and downstream of the bend. Enhanced particle-wall interaction in the region of the bend is characterized by elevated DEF values. The maximum DEF factors are similar between the Eulerian (DEF_max = 4.4) and Lagrangian (DEF_max = 6.2) cases. Still, differences between these two models are observed. Based on the use of discrete particle elements, Lagrangian contours of DEF are less continuous in comparison to the Eulerian result. As expected from the regional deposition results in the 90° bend (Figure 7), the maximum DEF of the Lagrangian model is 30% higher than the Eulerian model prediction. This difference may also be a result of particle slip effects, which are captured in the Lagrangian simulations but neglected by the Eulerian model. Despite these differences, the Lagrangian user-defined BM model appears to adequately capture local deposition characteristics in comparison with the Eulerian approximation.

FIG. 8 Local DEF values in the 90° bend geometry for the (a) chemical species Eulerian and (b) Lagrangian user-defined BM models. Predictions of DEF values appear similar between the two models considered.

3.4. Deposition in the Oral Airway Geometry

Regional-averaged deposition efficiency results in the oral airway model at flow rates of 15, 30, and 60 L/min are illustrated in Figure 9. Particle diameters ranging from 1 to 120 nm have been considered. Numerical results for 15 L/min have been compared to the available in vitro data of Cheng et al. (1993, 1997a) (Figure 9a). These experimental results were based on inhalation flow rates of 10 and 20 L/min in realistic oral airway casts. Theoretically, the numerical results for 15 L/min should fall between the two experimental flow rates considered. For an inlet flow rate of 15 L/min (Figure 9a), the Eulerian model with turbulent dispersion appears to under predict the in vitro results for particles on the order of 1 nm. This under prediction appears to be consistent for particles less than approximately 10 nm. Results for the Fluent BM model under predict the deposition of 1 nm aerosols by a factor of approximately 3 and over predict the deposition of 100 nm particles by one order of magnitude. Implementing the user-defined BM routine significantly improves agreement between numerical and experimental results for 1 nm; however, this model significantly over predicts the experimental results beyond approximately 10 nm in the oral airway model. This over prediction of experimental results for particles larger than 10 nm, which was observed with both BM models, was found to be a result of particle velocity interpolation near wall boundaries. The trajectory calculation in Fluent 6 does not appear to interpolate particle velocity between cell centers and wall boundaries in a manner that approaches zero at the wall for turbulent flow fields. As a result, non-zero wall normal velocities were over-estimated at particle locations and deposition was over predicted. Inclusion of both the user-defined BM and NWI routines significantly improved deposition results from 1 through 120 nm particles. By implementing both of these user-defined Lagrangian models, numerical predictions at 15 L/min are now within the bounds of experimental uncertainty for most data points considered (Figure 9a).

Comparisons of Eulerian and Lagrangian model results for the oral airway geometry at flow rates of 30 and 60 L/min are provided in Figure 9b and c. Differences between the Eulerian model results with and without dispersion are discernable for the case of 30 L/min and increase significantly for 60 L/min. As expected, the Fluent BM model predicts little variation in deposition efficiency for the ultrafine aerosols considered at both flow rates. Inclusion of the user-defined BM model significantly increases the deposition of 1 nm particles providing a deposition efficiency curve that is similar to the Eulerian results for smaller particle sizes. The best agreement with the Eulerian results is observed for the user-defined BM and NWI Lagrangian model, which predicts the expected exponential decrease in deposition with increasing particle size. Still, differences between the Eulerian and user-defined Lagrangian model are apparent at 30 and 60 L/min. As with the case of 15 L/min, the Lagrangian model predicts significantly elevated deposition in comparison with the Eulerian model.

For the Lagrangian results reported in Figure 9, difficulty was observed in achieving grid convergence for the Fluent BM and user-defined BM deposition results without the NWI routine. This was because fluid velocities at the particle location did not approach zero near the wall. By reducing the near-wall control volume size beyond the limit required to achieve flow field and Eulerian transport grid convergence, near-wall particle velocities could be reduced. Therefore, successive grid refinement showed that the results without NWI could approach the NWI case. However, this is clearly an ineffective method to achieve grid convergence of the particle deposition results. A more effective alternative is to achieve grid convergence of the velocity field and then implement the NWI routine for particle velocity calculations. The results reported in Figure 9 are based on a grid that was tested for convergence of the velocity field, Eulerian model and user-defined BM and NWI routines (Figure 2). To illustrate the effectiveness of the NWI routine, this grid was then used for the Fluent BM and user-defined BM results.

FIG. 9 Model predictions of deposition efficiency in the oral airway geometry for particles ranging from 1 to 120 nm and tracheal flow rates of (a) 15, (b) 30, and (c) 60 L/min. For an inhalation flow rate of 15 L/min, results of the Lagrangian model with user-defined BM and NWI routines provide the best match to the available in vitro deposition data of Cheng et al. (1993, 1997a).

Local values of deposition based on Eulerian and Lagrangian simulations of 5 nm aerosols at 30 L/min are provided in Figure 10. For this comparison, both models have excluded the effects of turbulent dispersion and the Lagrangian result is based on the user-defined BM and NWI routine. The pattern of deposition enhancement appears very similar between the Eulerian and Lagrangian models. Both models display elevated deposition concentrations on the dorsal side of the throat where convective diffusion is enhanced due to flow curvature. A hot spot in deposition is also observed for both models on the lateral wall just distal to the glottis, or the region of minimum cross-sectional area. The primary difference in Eulerian and Lagrangian simulations is the maximum DEF value observed. The maximum DEF value for the Lagrangian model of 17.9 is significantly higher than the Eulerian maximum of 4.8. While this difference is significant, it is not unexpected considering the differences observed for Eulerian and Lagrangian deposition efficiency in the oral airway model (Figure 9). Furthermore, the Lagrangian contours of DEF appear somewhat discontinuous compared to the Eulerian model result.

FIG. 10 Local DEF values in the oral airway geometry for the (a) chemical species Eulerian and (b) Lagrangian user-defined BM and NWI models. Contour patterns of DEF values appear similar between the models; however, significant differences in the magnitudes of DEF values are observed.

4. DISCUSSION

This study has considered the effectiveness of direct Lagrangian particle tracking algorithms for calculating ultrafine aerosol transport and deposition in flow fields consistent with the upper respiratory tract. Lagrangian deposition results have been compared with a chemical species Eulerian model, which neglects particle inertia, and available experimental data. The Lagrangian particle tracking algorithms considered include the Fluent 6 Brownian motion (BM) routine, a user-defined BM model, and a user-defined BM model in conjunction with a near-wall interpolation (NWI) algorithm. Ultrafine particle deposition has been assessed based on regional-averaged deposition efficiency and local DEF values.

Results indicate that the Fluent BM model inaccurately predicts the diffusional deposition of ultrafine aerosols by up to one order of magnitude in all cases considered. In contrast, the Lagrangian model with user-defined BM and NWI routines was observed to give the best predictions of regional deposition results for all three geometries in comparison to analytic correlations and experimental data. For the tubular and 90° bend geometries, good agreement with analytic and empirical data was observed for both the Eulerian and user-defined Lagrangian models. Furthermore, both Eulerian and user-defined BM-NWI models predicted similar local deposition patterns for the relatively simple tubular and 90° bend models. Considering the complex oral airway geometry, the user-defined BM and NWI routines were required to adequately match the available in vitro deposition data of Cheng et al. (1993, 1997a) for particle sizes ranging from 1 to 120 nm. Moreover, significant differences between the Eulerian and Lagrangian model predictions of regional and local deposition were observed in the oral airway.

A significant finding of this study is that an appropriate Lagrangian model was observed to better predict available ultrafine experimental deposition data on a regional basis in both a 90° bend and an oral airway geometry in comparison to the Eulerian approach. Inclusions of user-defined BM and NWI routines were necessary to supplement the existing Lagrangian model. Surprisingly, differences between the Eulerian and Lagrangian models were observed for particle sizes as small as 1 nm in the oral airway geometry. Discrepancies between these approaches may arise as a result of inherent differences in the models, the effects of particle inertia, and the slip correction factor. The fact that excellent agreement was observed between the Lagrangian and Eulerian models applied to fully developed tubular flow indicates that these two particle transport approaches can provide nearly identical results. Therefore, differences between the model results are most likely associated with finite particle inertia and slip correction effects.

As discussed in Longest and Xi (2007), the Lagrangian model accounts for finite particle inertia at all particle sizes. In contrast, the chemical species Eulerian model employed in this study neglects particle inertia. Longest and Xi (2007) showed that the effects of particle inertia on area-averaged deposition efficiency could be neglected for particle Stokes numbers below 5 × 10^{− 5}. In the oral airway model, Stokes numbers for particle diameters up to 120 nm are below this inertial particle limit at a flow rate of 30 L/min. Particles on the order of 100 nm begin to exceed the inertial particle limit at a tracheal flow rate of 60 L/min. Based on these observations, particle inertia is expected to play a minor role in ultrafine aerosol deposition in the oral airway model and is not likely responsible for the observed significant differences in Eulerian and Lagrangian predictions of deposition efficiency for smaller particles. However, the finite inertia of nanoparticles may significantly influence local deposition patterns for particles as small as 40 nm (Longest and Xi 2007).

Another potential factor responsible for the observed difference in Eulerian and Lagrangian predictions of deposition efficiency in the oral airway model is the phenomenon of particle slip due to rarefied gas effects. The influence of slip in terms of the Cunningham correction factor influences the diffusion coefficient of both the Eulerian and Lagrangian models (Equation [9]). However, the effects of slip can also cause both particle devisions from streamlines and increases the residence times, which are not considered by the chemical species Eulerian model. Therefore, effects of slip correction on Lagrangian particle trajectories and residence times may provide a potential mechanism for the observed increase in Lagrangian particle deposition. Differences between Eulerian and Lagrangian results due to slip correction should be larger for smaller particles. In contrast, differences in model predictions due to inertia effects should increase with particle size. Deposition efficiency results for the oral airway geometry indicate larger differences between the Eulerian and Lagrangian model predictions at smaller particle sizes (Figure 9). This difference appears to decrease as particle size increases. Therefore, the elevated deposition predicted by the Lagrangian model and the associated improved agreement with experimental results may be largely due to slip correction effects. However, it is not apparent why slip correction effects did not result in differences between the model predictions for fully developed tubular flow.

Based on the results of this study, the Lagrangian model with appropriate BM and NWI routines appears to adequately model regional and local ultrafine aerosol deposition in the idealized respiratory geometries considered. Moreover, the Lagrangian model with user-defined functions provides a better match to area-averaged ultrafine experimental deposition data in the 90° bend and oral airway geometries, compared with the Eulerian approximation. An additional advantage of the Lagrangian model is that it can be extended to larger particle sizes where inertia plays an increasing role in deposition. For the larger ultrafine particles considered in this study, inertial effects are likely responsible for differences in Lagrangian and Eulerian model predictions at a flow rate of 60 L/min (Figure 9c). As indicated by Longest and Xi (2007), particle inertia will significantly influence the regional and local deposition values of fine particles ranging from 100 nm to 1 μ m in the upper airways. Having one model that can directly account for the effects of slip correction and inertia and that can be applied across a spectrum of particle sizes from ultrafine to coarse is clearly advantageous.

The primary disadvantage of the Lagrangian model employed in this study is the number of particles and the associated computational resources required to fully resolve the area-averaged and local deposition values. For the oral airway model, the simulation of 300,000 particles was required to produce a convergent maximum DEF value. However, the contour field of DEF values remained discontinuous due to the inherent discrete nature of the Lagrangian approach. In contrast, the Eulerian model provided continuous contours of particle deposition that may be more consistent with ultrafine deposition in vivo for physically relevant exposure scenarios. A potential enhancement to the existing Lagrangian model would be a probability density function approach used to account for particle diffusion along mean trajectories. Interactions of the probability field and the wall could then be used to model diffusional deposition and reduce the overall number of required Lagrangian elements.

In this study, the equation governing nanoparticle transport has been directly integrated with the inclusion of a random force to describe the BM. This approach requires the assumption that the particle response time (τ _p ) is much smaller than the characteristic time for changes in the flow field or the particle integration time-step (Δ t). As indicated by Podgorski (2001), the assumption that τ _p ≪ Δ t may not be valid as particles approach solid boundaries. A potentially more robust approach to the solution of the Langevin equation is the Brownian dynamics method proposed by Podgorski (2001) and implemented by Moskal et al. (2006). This Brownian dynamics approach allows for the inclusion of coupling effects and particle-wall interactions. A solution to the Langevin equation based on evaluation of the corresponding Fokker-Plank equation has been described by Peters (1999). However, the disadvantage of the theoretically complete and robust approaches of Peters (1999) and Podgorski (2001) is a significant increase in computational complexity. Based on comparisons to available analytic and experimental deposition results made in this study, the direct numerical integration of the particle transport equation including a Brownian motion force term appears to be adequate for general deposition predictions. Future studies are suggested that evaluate the performance of the Brownian dynamics and Fokker-Planck equation approaches in comparison with direct integration of the particle trajectory equation.

In this study, a transitional or turbulent flow field was assumed to occur in the oral airway model. However, effects of turbulent dispersion on nanoparticle trajectories have been neglected for two reasons. First, turbulent dispersion may become more significant than BM diffusion in some cases depending on the particle size and turbulent intensity. Turbulent dispersion was neglected to isolate the effects of BM, which was the focus of this study. The second reason for neglecting turbulent effects on trajectories is that particle dispersion models for isotropic conditions are generally inaccurate (Matida et al. 2004). A potentially better method for evaluating transitional and turbulent flow fields and the dispersion of nanoparticles is to perform Large Eddy Simulations (LES). Good agreement has been reported between LES predictions and micro-particle deposition experiments in a number of studies, such as the curved tube geometry considered by Breuer et al. (2006). However, the LES approach may become too computationally expensive for simulations over significant time periods, such as a respiratory inhalation waveform. A potentially more efficient method for simulating particle deposition in transient respiratory dynamics is the anisotropic turbulence approach suggested by Matida et al. (2000, 2004).

In addition to the absence of turbulent particle dispersion, other limitations of this study include the assumptions of simplified inlet and outlet conditions, the use of idealized geometries, and steady flow. Other studies have shown that geometric surface effects, such as cartilaginous rings (Martonen et al. 1994; Zhang and Finlay 2005), may significantly enhance the deposition of particles beyond diffusional effects. While correlations based on the tubular and 90° bend geometries have been used for many available whole lung models, these geometries are not very representative of realistic upper airways (Nowak et al. 2003). Moskal and Gradon (2002) and Zhang et al. (2002) have shown that transient flows may have a significant impact on deposition in vivo, especially in the interim between inhalation and exhalation. While the assumptions made in this study directly impact model realism, they were largely necessary to approximate the available experimental systems that were considered. Furthermore, the effect of turbulent dispersion on Eulerian particle transport was shown to be relatively minor in the oral airway geometries at flow rates of 15 and 30 L/min (Figure 9a and Figure 9b). These limitations must be addressed in future studies in order to develop more physically realistic models of particle deposition for fine and ultrafine aerosols.

In conclusion, a Lagrangian model of ultrafine aerosol transport appears to effectively simulate deposition in idealized geometries of the upper respiratory tract in comparison with available experimental results. Necessary modifications to the commercially available particle tracking routine include an appropriate Brownian motion model and near-wall interpolation of fluid velocities at particle positions. Future studies will be necessary to improve model realism before these results can be directly applied to make dose-response and health effects predictions. A Lagrangian based probability density function model of ultrafine aerosol transport may be an effective alternative that can potentially capture the effects of finite particle inertia, slip correction, and significantly improve simulation efficiency.

Acknowledgments

This work was sponsored by Philip Morris USA (Dr. Mohammad R. Hajaligol, Program Manager).