Water Vapor Transport and Its Effects on the Deposition of Hygroscopic Droplets in a Human Upper Airway Model

The fundamentals of 3-D airflow as well as heat and water vapor transport and droplet vaporization (or hygroscopicity) are described for a human upper airway model under steady laminar-transitional-turbulent inspiratory flow conditions. Water vapor distributions from the mouth to the first four generations of the tracheobronchial tree are given in terms of relative humidity or mass fraction. The mass transfer coefficients of water vapor are correlated as a function of local flow rate and temperature-dependent diffusivity, which can be readily used for estimating the regional water loss or moisture variations in the human upper airways. Furthermore, the dynamics of hygroscopicity and deposition of isotonic saline droplets have been simulated as an example, applying the basic theory. Specifically, droplet evaporation rates and deposition pattern are analyzed and the effects of inhalation flow rates and thermodynamic air properties are discussed.

INTRODUCTION

Many types of droplets, such as saline, pharmaceutical, and fuel droplets, undergo size changes due to evaporation or condensation (hygroscopicity) when being inhaled through the respiratory tract. Especially the effects of water evaporation and condensation on aerosols in the human respiratory system have been experimentally and theoretically investigated (see Broday and Georgopoulos 2001; Eisner et al. 1990; Ferron et al. 1988, 1989; Finlay and Stapleton 1995; Gebhart et al. 1990; Hickey and Martonen 1991; Morrow 1986; among others). These studies show that the functional dependence of deposition fraction on particle size shifts towards smaller size particles with increasing hygroscopicity and relative humidity. For example, hygroscopicity can substantially increase the size and hence the airway deposition fraction of drug particles with initial dry sizes of about 0.5–2 μ m (Ferron et al. 1989). Traditionally, deposition equations, based on analytical solutions for particle diffusion, impaction and sedimentation in straight pipes, were employed in modeling deposition fractions of hygroscopic aerosols in human lungs. The temperature and relative humidity of air was assumed to be at alveolar conditions, i.e., T = 37°C and RH = 99.5%, or calculated with one-dimensional mathematical models. Experimental data on the regional deposition of hygroscopic aerosols have been hardly available until the rapid development of the single-photon emission computerized tomography (SPECT) technique as reviewed by Finlay et al. (1996). However, both theoretical and experimental studies focused on total and regional deposition of hygroscopic aerosols. Hence, further studies of effects of water vaporization or condensation on droplet transport and local deposition is most important for dose assessment of such size-changing aerosols.

Accompanied by droplet shrinkage or growth is usually the conditioning process of the inspired air, which may affect droplet transport and deposition. The investigation of heat and moisture transport in human airways is also very important in clinical studies for treating pulmonary disorders. For example, excessive loss of water from the airways can cause cilial damage, inducing pulmonary infection or asthma and altered airway mechanics (Chen and Horton 1977; Williams et al. 1996). There have been a number of reports published attempting to quantify heat and water vapor transfer in the human airways. For example, Johnson et al. (1977) obtained a convective heat and mass transfer coefficient, h _c and h _m , averaged over three bronchial generations, based on in vivo measured airway temperatures and water vapor contents of several mechanically respired mongrel dogs. However, the correlation extrapolated from the canine measurements should be carefully examined before it is interpreted for humans. Nuckols (1981) carried out a series of convective heat and water vapor transfer experiments with a cast replica of human upper airways extending from the mouth and nose to the trachea as well as with a mechanical bronchial airway model representing generations G0 to G2 after Weibel (1963). Considering inspiration only, he expressed his results in terms of Nusselt and Sherwood numbers averaged over the entire cast from mouth to trachea and the bronchial airway model. Both Johnson et al. (1977) and Nuckols (1981) did not derive inspiratory convective heat and mass transfer coefficients for each individual generation of the upper bronchial airways. Daviskas et al. (1990) developed a time-dependent model, based on a single differential equation with an analytical solution as well as experimental data, to predict the intra-airway temperature and water vapor contents as a function of air residence time at each location. Most recently, Tawhai and Hunter (2004) modeled the water vapor and heat transfer in the normal and the intubated airways employing unsteady, one-dimensional (1-D) (axial) transport equations where a power-law method was used to describe the radial distribution of temperature, water vapor, and air velocity. The airway geometries were approximated as cylindrical segments.

In summary, most mathematical models and numerical simulations employed empirically derived heat and mass transfer coefficients to estimate transfer processes between resident air and airway walls, or focused only on 1-D and 2-D water vapor and temperature distributions in the human airways (see Daviskas et al. 1990; Hanna and Scherer 1986; Ingenito et al. 1986; Naftali et al. 1998; Tawhai and Hunter 2004; among others). In our previous paper (Zhang and Kleinstreuer 2003a), 3-D heat transfer in the human upper airway was simulated. However, coupled heat and water vapor transport phenomena in the human respiratory system have not been fully analyzed.

In this paper, water vapor transport as well as vaporization and deposition of isotonic saline droplets in a human upper airway model are numerically investigated. The distributions and mass transfer coefficients of water vapor in a gas stream, droplet evaporation rates, and local depositions are reported. The effects of inhalation flow rate and thermodynamic properties of the incoming air are discussed as well. It should be noted that the methodology and qualitative results presented are also applicable to the inhalation of other types of volatile toxic or therapeutic aerosols. The assumption of rapid mixing is implemented for the discrete liquid phase, i.e., infinite conduction and diffusion are assumed within the droplet, such that liquid temperature and concentration are spatially constant but vary with time. In addition, the effects of droplet evaporation on the water vapor transport in the airways have been neglected due to the underlying assumption of low mass loading of droplets.

THEORY

Upper Airway Geometry

As shown in Figure 1, the present upper airway model consists of two parts, i.e., an oral airway model, including oral cavity, pharynx, larynx, and trachea, and a symmetric triple bifurcation lung airway model representing generations G0 (trachea) to G3 after Weibel (1963). The dimensions of the oral airway model were adapted from a human cast as reported by Cheng et al. (1997). Specifically, the diameter variations along the present oral airway from mouth to trachea are almost the same as those for the hydraulic diameters from the cast. Variations to the actual cast include the circular cross sections, a short mouth inlet with a diameter of 2 cm, a modified soft palate, and a strong bend. The dimensions of the four-generation airway model are similar to those given by Weibel (1963) for adults with a lung volume of 3500 ml. The airway conduit is assumed to be smooth and rigid. The minor effects of cartilaginous rings, which may appear in the upper airways, have not been considered in the present analysis.

FIG. 1 3-D views of the oral airway model and bifurcation airway model (generations G0 to G3). B1—first bifurcation, B2.1 and B2.2—second bifurcation, B3.1, B3.2, B3.3, and B3.4—third bifurcation (the dashed lines indicate the segmental boundaries).

Airflow, Heat, and Mass Transfer

To capture the air flow structures in the laminar-to-turbulent flow regimes, i.e., 0 < Re_local < 10⁴ for the human airway configuration during inhalation, the low-Reynolds-number (LRN) k-ω model of Wilcox (1998) was selected and adapted which has been demonstrated to be appropriate for such internal flows (Zhang and Kleinstreuer 2003b). All transport equations for airflow including continuity, momentum, turbulence kinetic energy, pseudo-vorticity, and energy equations are given in Zhang and Kleinstreuer (2003a), (b).

The mass transfer equation for water vapor in the airway stream can be written as

where Y is the mass fraction, u_j is the velocity, ν _T is the turbulent viscosity, σ _Y = 0.9 is the turbulence Prandtl number for Y, and [Dtilde] is the molecular diffusivity of water vapor in air. The effect of temperature change on diffusivity is expressed by the following semiempirical correlation (Vargaftik 1975):

where T is the air temperature. The water vapor concentration (Y) can also be expressed in terms of the relative humidity (RH). The relative humidity at an air-water interface is defined as the actual water vapor concentration divided by that under saturation conditions, i.e.,

where ρ _a is the density of dry air, and ρ _v,sat is the density of saturated water vapor which can be obtained from (Finlay 2001):

The associated boundary conditions for airflow as well as heat and water vapor transport include uniform velocity, temperature and water vapor concentration profiles prescribed at the mouth inlet, and zero pressure at the outlets.

Regional variations in wall temperature were employed as the thermal boundary condition for air at the air-wall interface. The wall temperatures in human airways can be described as (Daviskas et al. 1990):

where T _m is the temperature at the lips, which is close to the temperature at the end of expiration, L is the cumulative length at the middle of each airway, [Vdot] is the airflow rate in the corresponding airway, and [Vdot] _E is the ventilation rate.

The mass fraction of water vapor at the wall (Y_w) was modeled assuming variable saturation conditions at T_w, i.e.,

Considering the vapor pressure at the wall in equilibrium with the blood serum, RH _wall is assumed to be 99.5% (Daviskas et al. 1990).

The velocity, temperature and vapor concentration fields in the oral and 4-genertion airways are obtained separately. A long exit tube that extends about nine tube diameters from the end of the trachea was used in the oral airway model to eliminate the influence of the zero-pressure assumption at the outlet. However, because the actual length from the glottis to the first carina in the trachea is about 14 cm, the mean-velocity field, turbulence quantities, temperature and mass fraction distributions in the trachea at about 11 cm from the glottis, were adjusted as the inlet conditions of the bifurcating airway segment (i.e., Generations G0 to G3) (see Figure 1). The parent tube length of the bifurcation airway model is about 3 cm. Necessary adjustments included: (1) the magnitudes of velocity, turbulence quantities, temperature and water vapor concentration which were recalculated due to the slight variation in trachea diameter (D = 1.37 cm) and Weibel's lung model (D = 1.6 cm); and (2) the profiles of variables were also reconstructed for the inlet to the bronchial airways, G0 to G3, employing an inverse distance weighting method after Shepard (1968).

Transport of Micro-Size Droplets

With any given ambient concentration of non-interacting spherical droplets, a Lagrangian frame of reference for the trajectory computations of the evaporating droplets can be employed. In light of the large particle-to-air density ratio, dilute particle suspensions, negligible particle rotation and small thermophoretic forces, drag is the dominant point force. Hence, the particle trajectory equation can be written as

where u_i ^P and m_p is the velocity and mass of the droplet, respectively; ρ _g is the density of gas mixture; and C_DP is the drag force coefficient given by

Here,

and the particle Reynolds number is

In Equation (6), u_i is the instantaneous fluid velocity with u_i = [ubar]_i +u′_i, where [ubar]_i is the time-averaged or bulk velocity of the fluid, and u′_i are its fluctuating components. Turbulence is assumed to consist of a collection of randomly directed eddies; hence, an eddy-interaction model (EIM) with near-wall correction is used to simulate the particle trajectories, and the fluctuating velocities u′_i are obtained by (Gosman and Ioannides 1981):

where ξ _i are random numbers with zero-mean, variance of one, and Gaussian distribution. Due to the assumption of turbulence isotropy, the fluctuating velocities normal to the wall calculated with Equation (10) may be higher than the actual values (Kim et al. 1987), which may overpredict the particle deposition in some cases (Matida et al. 2004). As proposed by Matida et al. (2004), a near-wall correction can be used to simulate the near-wall particle trajectories, i.e., the component of fluctuating velocity normal to the wall u′ _n can be expressed as

with

where f_v is a damping function component normal to the wall considering the anisotropy of turbulence near the wall (Wang et al. 1999), and y⁺ = u_τ y/ν is the inner variable or sub-layer-scaled distance with u_τ being the friction velocity (u_τ ≡ ρ _g , where τ _w is the wall shear stress). Usually, Equation (12) is used for y ⁺ < 10.0, while f_v = 1 elsewhere.

At the mouth inlet, presently uniform droplet distributions were prescribed. The effect of randomized non-uniform droplet inlet distributions on the deposition efficiency has been discussed by Zhang and Kleinstreuer (2001). The initial droplet velocities were set equal to that of the fluid, and one-way coupling was assumed between the air and droplet flow fields because the maximum mass loading ratio (mass of particles/mass of fluid) is below 10^{− 5} in the present analysis. Droplet diameters, positions, velocities and temperatures in the cross section O-O (see Figure 1) in the lower trachea were set or adjusted as the inlet particle conditions of the bifurcating airway segment. The adjustment process for particle position and velocity is similar to that for air flow described in the previous subsection.

Droplet Vaporization

In light of the complex liquid-phase aerosols, this analysis assumes that the droplet temperature and composition are maintained spatially uniform but evolve with time and the droplet remains spherical during vaporization (Clift et al. 1978). This is reasonable because of the small Biot number (Bi ≪ 1) for micro-size droplets. Considering the convective heat and mass transfer over the whole surface of the spherical droplets, the change in droplet mass can be expressed as (Longest and Kleinstreuer 2004)

where m_p is the droplet mass, n_s is the mass flux of each evaporating chemical species at the droplet surface, which can be calculated as (Benaissa et al. 2002; Berlemont et al. 1995)

where Sh is the Sherwood number, Sh = (1 +Re _p Sc)^1/3 max [1,Re _p ^0.077]; and Y_s,surf and Y_s,∞ are the gas-phase mass fractions of each component on the droplet surface and far from the droplet, respectively. Specifically, for water or saline droplets,

and

where ρ _v,surf and ρ _v,∞ are the water vapor density on the droplet surface and in the nearby surrounding air stream, respectively; ρ _a is the density of dry air. For saline droplets, the dissolved Na⁺ and Cl⁻ ions may reduce the water vapor concentration at an air-water interface. In other words, the relative humidity at the air-water interface on the saline droplet surface (RH_d) is not equal to one. It can be determined experimentally or described by an empirical equation as (Finlay 2001):

where x_s and x_w is the molar concentration of solute and water, respectively; and i is sometimes called the Van't Hoff factor and is determined experimentally. For dilute NaCl solutions, the measured value of the Van't Hoff factor is i = 1.85. Considering physiological saline droplets (0.9% NaCl), the calculated RH_d is about 0.995 according to Equation (17), which indicates that this solution is isotonic. Then, the actual water vapor concentration in Equation (16) on the surface of the saline droplet can be determined by

which is very close to that on the surface of a pure water droplet. More realistically, continuous evaporation of water may increase the droplet solute (NaCl) concentration; hence, actual activity data should be used to account for the solute effect instead of using Raoult's law (i.e, Equation (17)) for dilute solutions. As measured by Peng et al. (2000), a hysteresis exists for RH_d(C_s) (where C_s = m_s/(m_s+m_w), and m_s, m_w are the mass of solute and water, respectively) during the processes of crystallization and deliquescence. Specifically, when RH_d decreases from a high value (e.g., 0.995) and Y_{s, surf} is larger than Y _s,∞, an aqueous NaCl droplet loses water and its concentration increase gradually until an RH_d of about 0.45 is reached when the droplet crystallizes to form a solid particle. When the ambient RH (RH_∞) is lower than about 0.76 (i.e., the RH_d that is at equilibrium with a saturated solution with C_s = 0.267), the solid NaCl particles remain dry; otherwise, the dry particles absorb water and form droplets (i.e., deliquescence).

In the present study, the relationship between RH _d and the solute concentration C_s for the unsaturated and saturated droplet is expressed by the formula derived by Cinkotai (1971), i.e.,

The RH_d(C_s) for a supersaturated droplet is correlated based on the measured data of Peng et al. (2000) as:

Equation (20) is only applicable for the vaporization process. For the deliquescence (hygroscopic) process,

The density of saline droplets also vary with the evaporation or condensation of water, which can be expressed as:

where ρ _w ,ρ _s are the density of pure water and pure NaCl particle, respectively. ρ _s = 2.165 kg/m³ is used in this study.

The coupled heat transfer equation for liquid droplets reads,

where c_pd is the liquid specific heat; d_p is the droplet diameter; k_g is the thermal conductivity of gas mixture; T_p and T_a are the temperature of droplet and surrounding air, respectively; γ _s is the latent heat of vaporization; and Nu is the Nusselt number, Nu = (1 + Re _p Pr)^1/3 max[1, Re _p ^0.077].

In the present analysis, the saline droplet diameter is always larger than 0.5 μ m with a initial d_p ≥ 3 μm and C_s = 0.9%. Hence, Brownian motion of particles with such large diameters in the human upper airways can be neglected (see Zhang and Kleinstreuer 2004). The Kelvin and non-continuum effects are minor for droplets with diameters d_p ≥ 0.5 μm; however, they are easily incorporated into the simulations for smaller droplets as described in the following section. The Kelvin effect for small droplets can be considered by introducing a factor K (Finlay 2001; Adamson 1990), i.e.,

where σ is the surface tension at the droplet surface (which is about 0.07 Nm^{− 1} for water at 303 K), M is the molar mass of the vapor molecules, and R_u is the universal gas constant. Combining Equations (24) and (28), we have:

The non-continuum effects (Knudsen number correction) can be considered by including a correction factor C_m and C_T in the vaporization and heat transfer equations for water or saline droplets, i.e., Equations (13) and (23) are replaced by (Vesala et al. 1997; Finlay 2001):

where

Here, Kn is the Knudsen number, defined as Kn = 2λ/d_p and λ is the mean free path of the gas surrounding the droplet; α _m and α _T are the mass and thermal accommodation coefficients, respectively. Values of α _m = α _T = 1 were used in the present study (Ferron et al. 1988; Hinds 1993).

Whatever is being released by the droplets (see Equations (26) and (27)) becomes a heat and vapor source for the air, possibly affecting the air temperature and vapor concentration (cf. Finlay and Stapleton 1995). However, these source terms are ignored in the present study in light of the low mass loading ratio of droplets (i.e., mass of droplets/ mass of inhaled water vapor < 5 × 10^{− 4}).

Deposition Parameters

The regional deposition of micro-droplets in human airways can be quantified in terms of the deposition fraction (DF) in a specific region (e.g., oral airway, bifurcation airways); it is defined as:

The local deposition patterns of micro-droplets can be quantified in terms of a deposition enhancement factor (DEF). As proposed by Balashazy et al. (1999, 2003), the deposition enhancement factor is defined as the ratio of local to average deposition densities, where deposition densities are computed as the number of deposited droplets in a surface area divided with the size of that surface area. The mathematical expression for DEF is

where A_i is the area of local wall cell (i), n is the number of wall cells in a specifc airway region, and DF _i is the local deposition fraction in local wall cell (i) which is given as

Clearly, the presence of high DEF-values indicates non-uniform deposition patterns, including “hot spots”, i.e., excessive local concentrations.

NUMERICAL METHOD

The numerical solutions of the governing equations describing fluid velocity, temperature and water vapor concentration fields were carried out with a user-enhanced finite-volume based program, i.e., CFX4.4 from ANSYS, Inc., previously AEA Technology (2001). The numerical program uses a structured, multiblock, body-fitted coordinate discretization scheme (Zhang and Kleinstreuer 2002). A Quadratic Upwind (QUICK) differencing scheme, which is third-order accurate in space, was used to model the advective terms of the transport equations. The sets of linearized and discretized equations for all variables were solved using the Block Stone's method.

Once the airflow velocity, temperature and vapor concentration fields were obtained, the droplet transport and evaporation equations were solved with our off-line F90 code with parallelized algorithms (Longest et al. 2004). For the calculation of droplet trajectories, geometry, velocity and turbulence data at all control-volume vertices are first extracted from the CFX solution and written to arrays. A second-order improved Euler predictor-corrector method (Longest et al. 2004) is then used for the integration of the droplet motion equations (see Equation (7)) including turbulent dispersion effects with near-wall correction (see Equations (11–13)) as well as the size-change of droplets (see Equations (24 and 25)). After each iteration for each droplet, the information about position, time, diameter, temperature and three components of the velocity was obtained. Droplet deposition occurs when its center comes within a radius from the wall. The near-wall effects, migration, or resuspension have been currently ignored (Longest and Kleinstreuer 2004). The number of particles at the inlet, n ≥ 10,000, was increased until the deposition fractions became independent of the number of particles simulated.

The computational mesh was generated with CFX Build4. The near-wall region required a very dense mesh. Specifically, the thickness of the near-wall cells was chosen to fully contain the viscous sub-layers and to resolve any geometric features present there. The mesh topology was determined by refining the mesh until grid independence of the flow field solution and droplet deposition fractions was achieved. The computations were performed on an IBM p690 workstation with multiple 1.3 GHZ CPUs. The steady-state solution of the flow field was assumed to be converged when the dimensionless mass residual reached (Total Mass Residual)/(Mass Flow Rate) <10^{− 3}. The convergence of other variables was monitored as well. Typical run times for the fluid flow and mass transfer simulations was approximately 160–500 CPU hours for the oral airway model and 70 CPU hours for the four-generation model under steady inhalation conditions. Utilizing the converged flow field solution, the micro-droplet trajectory and vaporization simulations required approximately 80 to 150 CPU hours for each case simulated.

MODEL VALIDATION

The present computational fluid-particle dynamics (CFPD) model has been validated with various experimental data sets for steady and transient laminar flows in bifurcations (Comer et al. 2001; Zhang and Kleinstreuer 2002) and for laminar, transitional and turbulent flows in tubes with local obstructions (Zhang and Kleinstreuer 2003b). Similarly, micro-particle depositions in airways were successfully compared with measured deposition efficiencies and deposition patterns (Comer et al. 2000; Zhang et al. 2002a; Zhang et al. 2005).

The validation results of droplet vaporization simulations are given in Figure 2. The simulated temporal diameter and temperature variations of a water droplet is compared with experimental measurements conducted by Smolik et al. (2001), where a water droplet with an initial diameter of 1.2 mm was suspended in a channel with moving air (T_air = 297 K, u_∞ = 0.203 m/s, and RH = 35%). It can be seen that the simulated vaporization law for water droplets is in excellent agreement with the experimental data. Considering the droplet temperature, rapid cooling occurs until the droplet wet-bulb condition is reached and then remains basically constant due to the balance between evaporative heat loss and convective heat absorption from the free stream. Concerning the experimental results of Smolik et al. (2001) for a water droplet, the present simulation matches the measured droplet temperature variation very well.

FIG. 2 Comparison of diameter and temperature evolution for a water droplet with the experimental data of Smolik et al. (2001).

The simulated mean air temperatures along the respiratory tract for T_in = 26.7°C and Q_in = 30 l/min were compared with the experimental and theoretical results of McFadden et al. (1985) and Daviskas et al. (1990), respectively, as shown in Figure 3. It should be noted that airway wall temperatures used in calculations by Daviskas et al. (1990) were directly derived from the experimental data of McFadden et al. (1985), which naturally results in a closer data match. Thus, a different temperature profile can be expected for the present 3-D simulation of complicated airway geometries when compared to the 2-D calculations with straight tubes as conducted by Daviskas et al. (1990).

FIG. 3 Mean air temperatures in the present simulation compared with experimental values of McFadden et al. (1985) and calculated results of Daviskas et al. (1990) with T_in = 26.7°C and Q_in = 30 l/min.

In summary, the good agreements between experimental findings and theoretical predictions instill confidence that the present computer simulation model is sufficiently accurate to analyze laminar-to-turbulent fluid-droplet thermodynamics in three-dimensional oral and upper bronchial airways.

RESULTS AND DISCUSSION

The airflow structures and temperature distributions in the human upper airway model have been discussed by Zhang and Kleinstreuer (2003a). This paper focuses on water vapor transport as well as droplet evaporation/hygroscopicity, transport, and deposition while effects of the initial solute concentration on droplet deposition are discussed in Zhang et al. (2005b).

Water Vapor Distributions

The 3-D water vapor concentrations are similar to those of fuel vapors, which are described by Zhang and Kleinstreuer (2003b). The quantitative variations of water vapor concentrations in the upper airway model (mouth to end-of-trachea) under different inspiratory flow conditions are displayed in Figure 4 in terms of cross-sectional area-averaged relative humidity (RH). Inhalation flow rate as well as air temperature and relative humidity at the mouth inlet are three important factors influencing the variations of water vapor concentrations in human airways. With the transport of water vapor from the airway wall to the unsaturated air stream, the RH increases gradually from the mouth to the trachea and upper bronchial tree. The lower the inlet values (RH_in), the higher is the increase rate for the RH due to the larger vapor gradients and wall vapor fluxes (see Figure 4). The RH-variation is larger with higher inlet air temperature because the saturation vapor pressure at the liquid/gas interface increases with temperature. The averaged RH in the airway conduits with a lower inlet air temperature (T_in) increases axially much more rapidly than at a higher T_in because of the lower saturated vapor pressure corresponding to the lower temperature (see Figure 4). The RH after the trachea can be as high as 84% and 97% at RH_in = 20% and 80% for a low inlet temperature of T_in = 283 K. Starting with an initial air temperature of 283 K and an RH of 80%, the maximum RH-value can reach about 104% (i.e., supersaturation) in the pharynx/larynx region. Supersaturation was also reported in the theoretical work by Ferron et al. (1985) for nasal breathing at a low inlet air temperature but for a high inlet RH. The occurrence of supersaturation may be attributed to the rapid variations in temperature and vapor concentration.

FIG. 4 Variations of cross-sectional averaged relative humidity with Q_in = 30 l/min in: (a) the oral airway model; and (b) the bifurcation airway model G0–G3.

In the airway segments G0 to G3, the mass transfer patterns are strongly affected by the skewed and asymmetric velocity profiles, which indicate different RH distributions in different airway tubes. The mean, cross-sectionally averaged RH values at the outlets of generations G0 to G3 are given in Figure 4b. The RH values at the outlet of G3 vary from 81% to 96%, for 20% ≤ RH_in ≤ 80%, Q_in = 30 l/min, and T_in = 303 K at the mouth inlet. For the case of T_in = 283 K and RH_in ≥ 20%, the RH reach the equilibrium value of 99.5% after generation G1.

Calculation of the Mass Transfer Coefficient

Calculation of the mass transfer coefficient of water vapor in human airways, h_m, is helpful in quantitatively predicting the regional water loss or variations of relative humidity of the inhaled air. The mass balance for one airway unit is:

which yields the wall mass flux of water vapor N, so that

In the above equations, the subscripts “0” and “i” denote the properties at the inlet (parent) and the ith daughter tube of a selected airway unit, respectively. In the oral airway model, the “i” denotes the trachea and n = 1; Y_m is the mean mass fraction of the “free” airstream in this region; [mdot] is the mass flow rate determined as

For simplicity, we only used the average values of the inlet and outlet cross-sectional mass fractions to determine the regional , i.e.,

Although the relative humidity and temperature of the incoming air at the mouth inlet (i.e., RH_in and T_in) have a significant influence on the variation of the mass fraction in airway conduits, their effects on the vapor mass transfer coefficient (h_m) is quite minor (see Figure 5). This may be attributed to the fact that both the numerator and denominator in Equation (34) are almost equally dependent on the inlet mass fraction which makes h_m only a weak function of RH_in and T_in. Hence, the mass transfer coefficient is mainly a function of inspiratory flow rate and diffusion coefficient. Once h_m is known, the mass fraction at the exit of an airway unit (Y _out ) can be calculated with the following equation which is derived from a mass balance of a micro-element:

where D _m and L are the mean diameter and length of the airway segment, respectively.

FIG. 5 Effects of inlet relative humidity and air temperature on the mass transfer coefficient of water vapor (Q_in = 30 l/min).

The mass transfer coefficient can also be nondimensionalized in terms of the Sherwood number (Sh), which is an indicator of convective dispersion, i.e.,

Figure 6 depicts the regional Sherwood number (Sh) for the oral and tracheobronchial airway models vs. the product of Reynolds (Re) and Schmidt (Sc) numbers, based on different inhalation flow rates and temperature-dependent diffusion coefficients. The diameters in the trachea and the parent tubes are employed to determine the Sh- and Re-numbers for the oral airways and bifurcation units, respectively. The correlation between the Sh-number and the ReSc-group for the oral airway can be expressed as:

where r² = 0.99. The convective mass transfer coefficients of water vapor for airway generations G0 to G3 (see Figure 6) are obtained for each individual bifurcation. Affected by the non-homogeneous airflow structures and mass concentration distributions, the regional Sh-numbers are slightly different for each bifurcation. However, a best-fit correlation yielded (Figure 6):

where r² = 0.94. It can also be seen from Figure 6 that the mass transfer is stronger in the bronchial airways than in the oral airways due to the presence of carinal ridges in G0–G3. That is, enhanced water loss (mass transfer) occurs just upstream of the carinal ridges due to the large concentration gradients in these regions. However, the long flow passage makes the oral airways a significant pathway for humidification of the incoming air (see Equation (37)).

FIG. 6 Regional Sherwood number (Sh) versus the product of Re and Sc for the human upper airway model.

Droplet Evaporation Rate

The evaporation rate of droplets (or percentage of droplet evaporation) can be defined as the ratio of total mass loss of all droplets to the total droplet mass at the inlet of the airway section, i.e.,

Where n is the total number of droplets and d_p,i,out and d_p,i,in is the diameter of individual droplet at the outlet and inlet of the airway section (i.e., oral airways and G0–G3), respectively. Figure 7 shows the percentages of isotonic saline droplet evaporation in the oral airway model, where the initial NaCl-concentration (C_s) was 0.9%. Saline droplets, initially with d_p ≤ 5 μm, may evaporate more than 95% of their mass in the oral airway when RH_in ≤ 60%, Q_in ≤ 30 l/min and T_in = 303 K. Thus, droplet evaporation can be very important when simulating the deposition of volatile aerosols in the human upper airways. The percentage of droplet evaporation increases with decreasing inlet RH-values, or an increasing inlet air temperature, or a decreasing initial droplet diameter (see Figure 7). The effect of inhalation flow rate is somewhat complicated. For relatively large micro-droplets (say, d_p > 5 μm), evaporation increases with decreasing inhalation flow rates due to the longer residence time. However, many droplets with a small initial size (say, d_p < 5 μ m) may start to absorb water in the oral airway due to the rapidly increasing solute concentration; hence, the total percentage of evaporation can be somewhat lower for a lower inhalation flow rate (say, Q_in = 15 l/min). The percentage of micro-droplet evaporation in the oral airway is very small (less than 6%) with incoming air at a low inlet air temperature (T_in) and a high inlet RH (e.g., T_in = 283 K and RH_in = 95%) (see Figure 7c).

FIG. 7 Variations of evaporation rate of isotonic saline droplets in the oral airway model under different inhalation conditions: Influence of (a) inhalation flow rate; (b) inlet RH and high T_in; and (c) inlet RH and low T_in.

As indicated in Equation (26), the diameter of each droplet, starting at the mouth inlet, decreases gradually due to mass loss by evaporation, and then may increase due to water condensation by hygroscopicity, driven by the airflow field and convective heat/mass transfer. Affected by the highly non-uniform, occasionally turbulent airflow structures in the oral airways, the trajectories for different droplets vary; hence, the changes of droplet diameters are different because of the different local airflow as well as heat and mass transfer. Figure 8 depicts the evolution of both trajectories and diameters in the oral airways, considering four different isotonic saline droplets at Q_in = 30 l/min, RH_in = 60%, and T_in = 303 K. Even when starting with the same diameters, i.e., 3 and 7 μ m, but different release positions, they undergo contrasting changes depending upon local (moist) environment. Specifically, the diameter of droplet #1 gradually decreases to about 5.5 μ m from the mouth inlet to the glottis due to vaporization; however, it may increase slightly in the trachea when it moves near the airway wall where RH_∞ is very high. Droplet #2 may decrease from 7 μ m to 2.6 μ m before it starts to grow after the trachea. Small-size droplets (say, d_p = 3 μm, i.e., droplets #3 and #4) may continuously shrink from the inlet until reaching an equilibrium size where they are neither losing nor absorbing water, and then they start to grow or shrink as they move through a high- or low-humidity environment. Clearly, the enhanced solute concentration has a significant effect on vaporization of initially small saline droplets (say, d _p ≤ 5 μ m) when compared to cases assuming constant solute concentration (i.e., solute and solvent vaporizes together) or single-component droplets (e.g., pure water droplets). With the assumptions of invariable solute concentration, almost all initially small-size droplets evaporate in the oral airways when RH_in ≤ 60%, Q_in ≤ 30 l/min and T_in = 303 K and hence very few droplets can enter the tracheobrochial airways. However, the absolute change of deposition fraction may tend to be minor in the human upper airways because inertial deposition is small for small-size microparticles. The hygroscopic growth of such droplets is larger in the lower airways and alveolar regions. The effect of enhanced solute concentration is relatively minor for large microdroplets.

FIG. 8 Trajectories and diameter evolution of selected isotonic saline droplets in the oral airway model with Q_in = 30 l/min, T_in = 303 K and RH_in = 60%.

Figure 9 shows the percentages of isotonic saline droplet evaporation in the bifurcating airway model G0 to G3. The negative values of evaporation rate in Figure 9b indicate water condensation on the droplet surface. With water evaporation or condensation and air humidification in the oral passages, droplet vaporization or hygroscopicity becomes more complicated in the upper tracheobronchial airways. Large isotonic saline droplets (say, d_p > 5 μm) may continuously vaporize with incoming air at high T_in (say, T_in = 303 K) and RH_in ≤ 80%. The effects of RH_in and d_p,in on evaporation are similar to those in the oral airways for such droplets (see Figure 9a). However, hygroscopicity is of importance for small-size droplets (d_p < 5 μm) even with a high T_in (say, T_in = 303 K and RH_in ≤ 80%; see Figure 9b) because the enhanced solution concentration due to evaporation in the oral passages greatly reduces the vapor concentration on the droplet surface (see Equation (19)). Droplet vaporization becomes minor (< 2%) for all size droplets with d_p,in ≤ 10 μm in the tracheobronchial airways G0–G3, when the incoming air temperature is low (say, T_in = 283 K) and RH_in ≥ 60% due to the air humidification in the oral passages (see Figure 9b).

FIG. 9 Evaporation rate of isotonic saline droplets in the bifurcation airway model G0-G3 under different inhalation conditions: Influence of (a) inlet RH, initial particle size and high T_in; and (b) initial particle size and inlet air temperature.

As shown in Figure 8, the monodisperse droplets entering at the mouth inlet turn polydisperse during vaporization. Figure 10 shows the cumulative distribution function of the droplet size at the outlet of the oral airways, or the inlet to the first bifurcation, for different inhalation rates and inlet air temperatures. The fraction less than a particular size can be obtained directly from the graph. The fraction of particles having diameters between two sizes can be determined from the difference of the cumulative fractions for these two sizes. The initial droplet size was a uniform distribution with d_p = 10 μm. As expected, the change of droplet size decreases with increasing flow rate due to decreasing residence times. The variation in such large droplet size is larger with lower inlet relative humidity because the water vapor pressure in the free air stream decreases with RH_in. The mean droplet diameters for initially 10 μm-droplets after the oral airways are d_p = 7.57, 8.51, 9.16 μm for Q_in = 15, 30, 60 l/min with T_in = 303 K and RH_in = 60%, respectively, and d_p = 5.99, 7.43, 9.36 μm for RH_in = 20%, 40% and 80% with Q_in = 30 l/min and T_in = 303 K. The variations in droplet size should greatly influence droplet deposition downstream. Clearly, the variation of droplet size is smaller with lower inlet air temperatures because the local RH in the conduits increases with decreasing temperature. The mean droplet diameters after the oral airways are d_p = 8.84, 9.24, 9.58, 9.84, and 9.98 μm for RH_in = 20%, 40%, 60%, 80%, and 95% with Q_in = 30 l/min, T_in = 283 K and initial droplet diameter of 10 μm, respectively.

FIG. 10 Cumulative distribution of isotonic saline droplet sizes at the inlet of the bifurcation airway model: Influence of (a) inhalation flow rate; (b) inlet RH and high T_in; and (c) inlet RH and low T_in.

The illustrations for evaporation and size change (see Figures 7, 8, 9, 10) also demonstrate that vaporization is dominant for isotonic droplets with a relatively large initial size (d_p,in > 6 μ m) in the human upper airway when compared to hygroscopicity for incoming air of RH_in ≤ 95% and T_in ≥ 283 K. For the small-size isotonic droplets (e.g., d_p,in < 6 μm), hygroscopicity may become significant after the trachea even if the inlet RH is relatively low (say, RH_in = 60%) and T_in is high (say, T_in = 303 K).

Droplet Deposition Pattern

The 3-D surface views of the local droplet deposition patterns in terms of deposition enhancement factor DEF (see Equation (31)) for isotonic saline droplets with different combinations of d_p,in, Q_in, T_in, and RH_in are shown in Figure 11. Clearly, micro-droplet deposition during inhalation is mainly due to impaction, secondary flow convection, and turbulent dispersion. Thus, they mainly deposit at stagnation points for axial particle motion, such as the tongue portion in the oral cavity, the outer bend of the pharynx/larynx, and the regions just upstream of the glottis and the straight tracheal tube. As shown in Figure 11, the maximum DEF-values are in the range of 100 to 900, which vary with the flow rate, initial particle size and inlet air parameter. This implies that the deposition patterns of micro-droplets in the oral airway are highly non-uniform, and hence a small surface area, where the maximum DEF occurs, may receive hundred times higher dosages when compared to the average value for the entire airways. Such a site of massive particle deposition is usually located in the glottis region and/or the outer bend in, or just after, the curved pharynx region for relatively large-size droplets (say, d_p = 10 μ m). Relatively small-size droplets (say, d_p = 5 μ m) may vaporize very rapidly and only deposit at the initial part of the oral airways (see Figure 11).

FIG. 11 3-D distributions of deposition enhancement factor (DEF) in the oral airway model for micro-size isotonic saline droplets.

Figure 12 depicts the distributions of DEFs in the bifurcation airway model G0 to G3, i.e., part of the bronchial tree. As expected, for micro-particles the high DEF values appear mainly around the carinal ridges due to inertial impaction, but the specific distribution of DEFs at each carina is different and varies with the inhalation parameters as well. Some micro-droplets also land outside the vicinities of the carinal ridges due to secondary flows and turbulent dispersion. However, the deposition sites of vaporizing droplets are not as broad as for non-vaporizing particles (see Zhang et al. 2005a). In the bifurcation airway model, the maximum DEF-values range from 400 to 1500 for the cases shown in Figure 12. The variations of maximum DEF-values for the given initial droplet diameter and inlet airflow conditions are complex, generating similar patterns as those for non-vaporizing particles (see Zhang et al. 2005a). The regional droplet deposition fractions (DF) (see Equation (30)) are also given in Figures 11 and 12. Clearly, the change of particle size due to evaporation greatly affects isotonic saline droplet deposition. The variation in DF is a function of initial droplet diameter, inlet temperature, RH, and flow rate, when considering droplet vaporization.

FIG. 12 3-D distributions of deposition enhancement factor (DEF) in the airways G0 to G3 for micro-size isotonic saline droplets.

CONCLUSIONS

The methodology for simulating water vapor and heat transfer as well as transport and deposition of vaporizing droplets has been discussed. As an example, the evaporation/ hygroscopicity, transport and deposition of isotonic saline droplets (NaCl ≈ 0.9%) have been simulated for the human upper airways and analyzed for different steady inhalation conditions. The results indicate the following:

1.	With the transport of water vapor from the airway wall to the unsaturated air stream, the relative humidity (RH) of the air flow increases gradually from the mouth to the trachea and upper bronchial tree as a function of inhalation flow rate, temperature and RH-value of air at the mouth inlet. The RH-values after generation G3 are in the range of 81% to 100%, when the inlet conditions are 20% ≤ RHin ≤ 80%, Qin ≤ 30 l/min, and Tin ≤ 303 K.
2.	The inlet relative humidity and air temperature have almost no effect on the mass transfer coefficients of water vapor in airway conduits. Thus, the mass transfer coefficients of vapors in the upper airway in terms of Sh-number can be correlated as a function of the product of Reynolds (Re) and Schmidt (Sc) numbers, i.e., and
3.	Saline droplets, initially with dp ≤ 5 μm, may evaporate more than 95% of their mass in the oral airways when RHin ≤ 60%, Qin ≤ 30 l/min and Tin = 303 K. Increasing the relative humidity or decreasing the temperature of the incoming air, or increasing the inspiratory flow rate, is helpful in reducing vaporization of saline droplets.
4.	Vaporization is dominant for isotonic saline droplets with a relatively large initial size (dp,in > 6 μm) in the human upper airways when compared to hygroscopicity for incoming air of RHin ≤ 95% and Tin ≥ 283 K. For small droplets (e.g., dp,in < 6 μm), hygroscopicity may become significant after the trachea even if the inlet RH is relatively low (say, RHin = 60%) and Tin is high (say, Tin = 303 K).
5.	Deposition sites of vaporizing droplets are more confined when compared with constant-diameter particles; however, the deposition enhancement factors are in the same order of magnitude (i.e., 102–103) for these two types of aerosols.

Acknowledgments

The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the Air Force Office of Scientific Research or the U.S. Environmental Protection Agency.

This effort was sponsored by the Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant number FA9550-04-1-0422 (Dr. Walt Kozumbo, Program Manager) and the US Environmental Protection Agency (Dr. C. S. Kim, Program Monitor). The U.S. Government is authorized to reproduce and distribute reprints for governmental purposes notwithstanding any copyright notation thereon. The use of both CFX software from ANSYS Inc. (Canonsburg, PA) and the IBM p690 at the High Performance Computing Center at North Carolina State University (Raleigh, NC) are gratefully acknowledged as well.