Numerical Simulations of Capillary Aerosol Generation: CFD Model Development and Comparisons with Experimental Data

Abstract

For a capillary aerosol generation system, the mechanisms governing droplet transport from the capillary tip through deposition in an enclosed geometry have not been previously explored. The objective of this study was to develop and validate a CFD model of transport and deposition for capillary-generated albuterol in water aerosols in a standard USP induction port used for pharmaceutical aerosol testing. Within this system, comparisons have been made between experimental measurements and numerical predictions of the jet angle, aerosol deposition in a sectioned induction port model, and size distributions of exiting particles. The CFD model employed accounts for multiscale and multicomponent flow initialized at the 57 μ m capillary tip and extending through the USP induction port with 30 L/min of co-flow air. A discrete Lagrangian particle tracking algorithm with corrections for near-wall anisotropic turbulence has been implemented to model the polydisperse particle phase including the effects of turbulent dispersion and evaporation. Results indicated good agreement between predictions of the numerical model and experimental in vitro measurements. The experimental mean (SD) total mass fraction of drug deposited in the sectioned induction port was 14.6 (1.1)%. Numerical predictions of deposited mass fraction for non-evaporating particles and evaporating droplets were 13.1% and 13.3%, respectively, resulting in relative differences of 10.3% and 8.9%. Comparisons between in vitro measurements and predictions within individual sections of the induction port resulted in relative differences as low at 0.75%. The predicted mass median diameters exiting the induction port for the particle and evaporating droplet models were 3.07 and 3.45 μ m, respectively, in comparison to an experimental value of 3.06 μ m. The numerical model developed in this study can be applied to optimize the capillary aerosol generation process and improve its delivery of aerosols to the lung.

INTRODUCTION

Capillary aerosol generation (CAG) has been demonstrated to be an effective method for creating a controllable and potentially monodisperse distribution of aerosol droplets (Hindle et al. 1998; Howell and Sweeney 1998; Hindle et al. 2004). Applications of this technology include pharmaceutical aerosol generation and delivery, toxicology testing, multiphase flow experiments, and nanoparticle coatings for precision sensors and monitors. In the CAG process, a liquid component is pumped through a micro-capillary where it is electrically heated to either a partially or a fully vaporized state. Capillary diameters are typically on the order of 100 μ m. The compressed heated vapor exits a tapered capillary tip or nozzle to create a microscale turbulent spray jet. At the capillary tip, the spray jet may exit as liquid droplets surrounded by saturated conditions, as a fully saturated vapor, or as a superheated gas. For pharmaceutical aerosol delivery, the liquid component passed through the capillary typically consists of an aqueous or non-aqueous carrier and a dissolved drug. With non-aqueous solutions, like propylene glycol, high supersaturation ratios can be achieved enhancing the formation of monodisperse droplets by homogeneous nucleation (Lesniewski and Friedlander 1998; Gupta et al. 2003). For aqueous pharmaceutical solutions, a mixture of droplets and vapor frequently exit the capillary tip, followed by air entrainment and evaporation.

Previous studies have considered various aspects of the CAG process. Hindle et al. (1998) implemented a prototype CAG model with a propylene glycol solution to demonstrate the production of a relatively monodisperse size distribution for inhaled pharmaceutical aerosols. To further control and define the aerosol size distribution generated from the CAG system, the use of reservoir chambers has been explored (Hong et al. 2002). Gupta et al. (2003) considered the effects of various non-aqueous solvents on the final aerosol size distribution for capillary-generated aerosols. The vaporization characteristics and solubility of solutes and solvents were shown to have an effect on aerosols formed by condensation (Gupta et al. 2003). The effects of thermal energy supplied by the capillary to heat and vaporize propylene glycol solutions was explored by Shen et al. (2004) where the temperature just prior to the capillary tip and the vapor jet temperature were shown to be the approximate boiling point temperature for propylene glycol under saturated exit conditions. As is often the case with aerosols designed for pharmaceutical purposes, the size range of particles or droplets that deposited between the capillary tip and the cascade impactor were largely neglected by these previous studies. However, the transport dynamics in the vicinity of the CAG are likely to have significant effects on downstream conditions.

This article begins to explore the mechanisms governing aerosol transport from the capillary tip through deposition using an enclosed geometry that is frequently employed for pharmaceutical aerosol testing. Beginning with the capillary exit, the aerosol transport process for the CAG system is characterized by droplet formation within a high Reynolds number turbulent jet resulting in a multiphase spray. Liquid droplets may be present in the efflux exiting the capillary tip and may form due to nucleation and condensation. In cases where the surrounding environment is at a relative humidity less than one, droplet evaporation may occur downstream of the capillary. For inhaled pharmaceutical aerosols, the geometry of the respiratory tract and the inhalation conditions may significantly affect drug delivery to the lung. For testing purposes, a standardized induction port is usually employed to connect pharmaceutical inhalers to the cascade impactors used to characterize their output at a constant flow rate (USP 2005). During product development, it is important to understand the processes that affect particle transport and deposition in both the respiratory tract and standardized induction port geometries.

The vapor jet exiting the capillary tip can be assumed to have a very high Reynolds number that is on the order of Re _d ∼ 10,000 for the CAG system. At much lower Reynolds number conditions, helical instabilities and other coherent structures have been observed in jet flow (O'Neill et al. 2004). For Reynolds numbers of approximately 1,000 and above, turbulent jet flow is characterized by axisymmetric instabilities propagating from the jet tip (Crow and Champagne 1971). Becker and Massaro (1968) described the flow exiting a turbulent jet with a Reynolds number greater than 10,000 as having no coherent vortex formation at the jet tip, very small disturbance wavelengths relative to the jet diameter, and fully chaotic behavior. For this high Reynolds number condition, Becker and Massaro (1968) also report an oscillation frequency of 5,700 Hz. Similarly, Weisgraber and Liepmann (1998) made particle image velocimetry measurements of a jet with an exit Reynolds number of 16,000 and reported no observable organized coherent structures. For the CAG microflow system, Shen et al. (2004) reported an initial jet angle between 20–30° without co-flow air. In a similar turbulent jet system, Lesniewski and Koch (1998) reported a jet angle of 26°.

Numerical simulations of turbulent jet systems are complicated by the miniscule time scales that must be resolved. Based on the observations of Becker and Massaro (1968), turbulent time scales are significantly less than an oscillation frequency on the order of 1,000 Hz for the CAG system. As a result, direct numerical simulations of this turbulent micro-jet are not feasible. However, Danaila et al. (1997) have conduced direct numerical simulations of low Reynolds number round jets. Large eddy simulations (LES) of high Reynolds number jets can directly account for some of the larger turbulent time scales and structures while modeling turbulence below a certain threshold. Andersson et al. (2005) have simulated the flow through an isothermal high Reynolds number jet using LES. However, this technique cannot capture the small time scales of turbulence for a very high Reynolds number jet over the relativity large time scales of a transient flow field problem, such as with an inhalation flow cycle ranging from 2 to 6 seconds.

Within the turbulent jet of the CAG system, liquid aerosols may form and grow due to nucleation, condensation, and coagulation resulting in a multicomponent spray. Homogeneous nucleation is a spontaneous gas-to-particle conversion that occurs in supersaturated gases (Sutugin and Fuchs 1968; Lesniewski and Friedlander 1998). Continued droplet growth may then occur through condensation onto existing particles and coagulation arising from droplet collisions (Friedlander et al. 1994). Lesniewski and Friedlander (1995) have shown that turbulent interactions can affect initial homogeneous nucleation and subsequent condensation and coagulation. Supersaturation conditions can be created in subsonic mixing jets (Lesniewski and Friedlander 1998) or in supersonic flows (Wen et al. 1988). In supersonic systems, upstream pressures are usually significantly higher than downstream conditions. The vapor saturation rate can increase significantly as the gas undergoes an adiabatic expansion (Wen et al. 1988). Liquid aerosols in the CAG system may form as a result of homogeneous nucleation induced by both turbulent jet mixing and supersonic adiabatic expansion. However, vapor supersaturation rates are typically reduced for aqueous systems resulting in a greater potential for condensation onto existing particles or impurities. The transport and dispersion of discrete particles in a turbulent jet, which may potentially lead to coagulation and size segregation, has been reviewed by Sbrizzai et al. (2004).

For capillary-generated aerosols, droplet evaporation may occur as the liquid spray mixes with ambient or co-flow air. Studies that have evaluated the condensation and evaporation of individual respiratory aerosols include Broday and Georgopoulous (2001), Martonen (1982), and Robinson and Yu (1998). Results of these analyses and others have been incorporated into one-dimensional whole-lung transport and deposition models for hygroscopic droplets (Persons et al. 1987; Ferron et al. 1988; Finlay and Stapleton 1995; Martonen and Katz 1996). A number of numerical studies have also considered the effects of droplet condensation and evaporation on particle deposition in detailed three-dimensional models of the respiratory tract. Schroeter et al. (2001) evaluated the effect of hygroscopic droplet growth on deposition in the nasal pathway. Gemci et al. (2002) applied a computational fluid dynamics (CFD) model to evaluate the spray mechanics of an aerosol field in a simple throat model. Longest and Kleinstreuer (2005) conduced CFD simulations to evaluate appropriate numerical models for multicomponent droplet evaporation in the respiratory system. The effects of near-wall proximity and shear stress on respiratory droplet evaporation were considered by Longest and Kleinstreuer (2004). Zhang et al. (2006a; 2006b) evaluated droplet transport and deposition in models of the oral airway through the third respiratory generation. Zhang et al. (2006a) reported that saline solutions with a concentration less than 10% did not have a significant impact on droplet evaporation and condensation characteristics. A recent experimental and analytical analysis by Stein and Myrdal (2006) indicated that for aerosol formation in metered dose inhalers, evaporation factors have a more significant impact on delivery efficiency in comparison with atomization effects.

In summary, a number of components related to the transport and deposition of capillary-generated aerosols have been considered, including turbulent jets, sprays, and hygroscopic droplet effects. These previous studies have implemented analytic, experimental, and numerical modeling approaches. However, the transport and deposition of aerosols and vapor from the CAG system have not been considered for respiratory aerosol conditions. The objective of this study is to develop and validate a CFD model of transport and deposition for capillary-generated aerosols in the standard 90° bend geometry used across the pharmaceutical industry. The region of interest extends from the capillary tip through a standard USP (U.S. Pharmacopeia) induction port (USP 2005). Validations of the numerical results will be based on comparisons to experimental data for sectional aerosol deposition in the induction port geometry and outlet particle size distributions. A laser particle analyzer will also be used to establish the initial aerosol size distribution. The CFD model will account for multiscale and multicomponent flow initialized at the capillary tip and extending through the USP induction port in the presence of 30 L/min co-flow air. Turbulent effects will be approximated using the low Reynolds number k-ω model over a 2-second period. A discrete Lagrangian particle tracking algorithm with corrections for near-wall anisotropic turbulence will be used to model the polydisperse particle phase including the effects of turbulent dispersion and evaporation. Based on successful comparisons of numerical and computational results, the developed CFD model can be applied to investigate the mechanisms responsible for deposition in the USP induction port geometry, as well as other realistic scenarios.

METHODS

The deposition of aqueous CAG aerosols in the USP induction port was determined using CFD modeling and in vitro experiments. Comparisons of experimental data and numerical predictions were used to test and validate the CFD model. To improve the quality of comparisons between experimental and numerical results, the standard USP induction port was divided into multiple sections.

Theoretical Flow Systems and Inlet Conditions

The aerosol generation system consisted of a heated micro-capillary positioned at the inlet of the standard USP induction port (USP 2005). As in previous studies, the contents of a micro-capillary were pumped and heated simultaneously, producing a partially or fully vaporized spray (Gupta et al. 2003). In this study, the capillary delivers a solution of 0.5% w/v albuterol (as sulfate) in water, at a mass flow rate of 25 mg/s over a 2-second pulse. Co-flow air at a relative humidity of 30% was pulled into the USP induction port around the capillary exit at a flow rate of 30 L/min. The tip of the capillary was 57 μ m in diameter, which induced a sonic exit velocity, high pressure and significant turbulence. The aqueous mixture exiting the capillary tip is composed of supersaturated water vapor and liquid droplets. The exiting droplets may grow due to heterogeneous nucleation in this supersaturated environment. However, rapid mixing with the co-flow air may result in relative humidity values less than one downstream of the capillary tip. For humidity levels below saturation conditions, the liquid droplets will begin to evaporate, which may reduce the probability for induction port deposition.

Inlet conditions for the specific CAG system considered in this study are illustrated in Figure 1. Regions of interest include the (1) co-flow air, (2) capillary nozzle, and (3) downstream particle size distribution (Figure 1). The co-flow air was assumed to enter at standard atmospheric temperature and pressure (25°C and 101.325 kPa), with a volumetric flow rate of 30 L/min resulting in an inlet Reynolds number of 1258. To determine conditions at the capillary nozzle exit, isentropic flow through the nozzle was assumed. Based on previous experiments, the mixture exiting the nozzle was assumed to be at the boiling point temperature of the solvent for atmospheric pressure conditions (Shen et al. 2004), which for water is 100°C. Considering the very large mass flow rate moving through the 57 μ m capillary tip, it was determined that the outlet flow was choked so that further increases in driving pressure did not increase the velocity beyond sonic conditions. Based on an ideal gas approximation, the sonic exit velocity for wet steam is c = 443.1 m/s (Figure 1). Conservation of mass was applied to estimate the average density of the aqueous solution at the nozzle exit (ρ = 22.1 kg/m³). The thermodynamic pressure at the nozzle exit was then determined from the ideal gas law and found to be 3.81 MPa. This very high static pressure at the nozzle exit experiences a sudden expansion as it adjusts to atmospheric conditions, which can further accelerate the flow and induce significant shear stresses. In Figure 1, Location 3 denotes a position 1.5 cm downstream from the capillary tip where relative humidity values approach unity and the aerosol droplets are likely well-formed and relatively stable. A laser diffraction particle size analyzer was used to classify the polydisperse size distribution at this location and the transport and deposition of these droplets was then modeled in the USP induction port geometry and compared with experimental deposition results.

FIG. 1 Inlet conditions for the CAG system including (1) co-flow air, (2) the capillary nozzle, and (3) downstream particle size distribution.

The flow field of the CAG system creates a number of challenges for numerical simulations. The region of interest for CFD analysis in this study extends from the capillary tip through the USP induction port. This region contains multiple geometric scales from a truly microflow system at the capillary tip (57 μ m) to the internal diameter of the USP induction port (31.8 mm). Therefore, the computational geometry must resolve the flow across a dimensional change of three orders of magnitude. Flow exits the nozzle at sonic velocity and a very high static pressure. The rapid expansion that occurs as the nozzle pressure adjusts to atmospheric conditions serves to accelerate the flow to supersonic speed on the microscale, which will quickly transition back to subsonic velocities. For flows at Mach numbers greater than approximately 0.3, compressible effects should be considered. The diameter-based Reynolds numbers of the co-flow air (1258) and nozzle exit (50,000) are indicative of laminar, transitional and fully turbulent conditions. Supersonic, compressible and highly sheared flows all increase the complexity of the turbulence generation and transition process. Within this environment, the combination of an aqueous solution with co-flow air creates a multicomponent mixture with concentration dependent properties. Liquid droplet aerosols are also present in the system, resulting in a multiphase flow. These liquid aerosols can grow or shrink in the multicomponent environment of air and water vapor by condensation, evaporation, coagulation, and droplet breakup. Furthermore, the presence of albuterol sulfate in the liquid aerosol may influence the droplet evaporation and condensation behavior.

Sectioned Induction Port Geometry

The standard USP induction port geometry considered in this study is shown in Figure 2. To improve the quality of comparisons between experimental and numerical results, the USP induction port has been divided into four sections, as indicated in Figure 2a. Comparisons of sectional induction port deposition data was used to avoid the previously observed scenario in respiratory aerosol systems where total deposition was matched but local or regional depositions were significantly different between model predictions and experiments (Oldham 2006; Longest and Vinchurkar 2007a). Dimensions of the sectioned USP induction port system have been illustrated in Figure 2a.

FIG. 2 Exterior and interior surfaces of the sectioned USP induction port geometry with (a) dimensions and (b) offsets used for assembly.

For the in vitro deposition experiments, the sectioned USP induction port was fabricated using a rapid prototyping process. Design of the sectioned induction port was performed in SolidWorks (SolidWorks, Concord, Massachusetts) CAD software. Individual sections were created to fit together using either offsets or pins (Figure 2b). Tolerances of 0.1 mm were used to ensure that the offset sections would slide together. The sectioned model was then created using in-house rapid prototyping capabilities on a Viper SLA system (3D Systems, Valencia, CA). This rapid prototyping system employees a 100 mW solid-state laser to selectively harden Accura 50 (3D Systems) plastic resin. The physical prototype was sanded smooth to approximate the stainless steel USP induction ports that are typically used for aerosol testing. Inspection of the assembled model indicated smooth transitions between the four individual sections.

EXPERIMENTAL METHODS

The in vitro component of this study included the measurement of droplet particle size distributions at the capillary tip, aerosol deposition losses in whole (conventional stainless steel) and sectioned USP induction ports, and aerodynamic droplet size distributions determined by cascade impaction after exiting the induction port geometry. Working conditions for the experimental CAG process were nearly identical to those simulated by the computational model. The formulation consisted of a 0.5% albuterol (as sulfate) aqueous solution (Nephron Pharmaceuticals Corp., Orlando, FL), which was pumped through the capillary at a flow rate of 25 mg/sec. Aerosols were generated as a bolus over a 2-second interval using an applied energy of between 49 and 50 J. The aerosols were sampled at a volumetric flow rate of 30 L/min via the test USP induction ports (conventional stainless steel and sectioned) into a Next Generation Pharmaceutical Impactor (MSP Corp, Minneapolis, MN) (Kamiya et al. 2004).

Following aerosol generation, washings were collected from either the conventional stainless steel induction port or the four regional sections of the fabricated sectioned USP induction port to determine the mass of deposited drug. An appropriate volume of a 1:1 admixture of methanol and deionized water was used to wash the induction port. The solutions were then assayed using a validated HPLC-UV method. The mass of drug on each section, together with the total deposition within the induction port was determined. The aerodynamic droplet size distribution of the aerosol exiting the stainless steel induction port and the fabricated sectioned induction port was also determined. All experiments were performed for greater than five trials at ambient room conditions, which have been approximated in the numerical model to be 25°C and 30% relative humidity.

The initial droplet size distribution exiting the capillary tip was determined in real time with a laser diffraction system (Spraytec, Malvern Instruments Inc., Southborough, MA). To assess the size distribution in a relatively stable section of the aerosol plume, the center of the laser was positioned 1.5 cm away from the capillary tip (Figure 1). In order to account for dilution of the vapor cloud, 30 L/min of co-flow air was pulled around the capillary and into the USP induction port during laser size analysis. The results of this size determination were represented as a time average over the 2-second CAG process.

Computational Geometry

The computational geometry of the system, which includes the interior surface of the USP induction port, is shown in Figure 3. As with the in vitro experiments, the 57 μ m capillary tip was positioned at the center of the circular USP induction port inlet. For the computational geometry, expanded inlet and outlet sections were included to approximate the boundary conditions of the physical experiment. The expanded inlet section replicates an open environment between the capillary tip and the USP induction port. Airflow was allowed to enter the circular low-x side of this expansion (Figure 3). In the physical experiments, the induction port was connected directly to the Next Generation Pharmaceutical Impactor (NGI). Because the impedance and volume of the NGI may have a direct impact on upstream flow conditions in the compressible transient system of interest, an expanded outlet was attached to the computational model of the USP induction port (Figure 3). This expanded outlet is similar in diameter to the stages of the NGI and has a volume approximately equal to that of the total NGI (Copley et al. 2005). Downstream of the expanded outlet section, the computational geometry returns to the outlet diameter of the USP induction port to ensure nearly fully developed outflow with a negligible radial component, required by the numerical model.

FIG. 3 Computational geometry of the USP induction port internal surface with expanded inlet and outlet sections to approximate experimental conditions. The expanded inlet section was used to represent exposure to an open environment. The expanded outlet section was used to approximate downstream conditions.

Continuous Phase Transport Equations

The flow field conditions generated by the CAG system are compressible, transient, locally supersonic, laminar to turbulent, unsteady, and temperature dependent. Considering the complexity of this system, the Reynolds-Averaged Navier Stokes (RANS) equations with a two equation turbulence model were selected to provide an efficient initial estimate of the continuous phase flow field. Of the available turbulence models that can be applied to the RANS equations, the low Reynolds number (LRN) k-ω approximation was employed to simulate the laminar, transitional and fully turbulent flow regimes of interest. 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 2003, 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). In the RANS equation set, the two component flow of air and water vapor will be approximated using a multicomponent mixture model. This model assumes that the various species behave as interpenetrating fields with concentration dependent properties. The individual species of air and water vapor were assumed to behave as ideal gases. The quality of these assumptions will be assessed based on comparisons of numerical and experimental results for the deposition variables of interest.

The full RANS equations for compressible multicomponent turbulent flow governing the conservation of mass and momentum are (Wilcox 1998)

In the above equations, 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 mixture density, and μ is the absolute mixture viscosity. Summation is assumed for like indices and the Kronecker delta function is represented as δ _ij (Kundu and Cohen 2004). Over-bars for all time-averaged quantities have been excluded from Equation (1) for clarity. 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 expression for turbulent viscosity is evaluated as (Wilcox 1998)

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

In the above equations, g _k and g _ω represent the generation of k and ω; ϵ _k and ϵ_ω represent the dissipation of k and ω, respectively (Wilcox 1998).

The transport of water vapor is governed by a convective-diffusive mass transfer relation (Bird et al. 1960)

In the above expression, Y _v is the mass fraction of water vapor, _v is the binary diffusion coefficient of water vapor in air, and Sc _T is the turbulent Schmidt number, which is taken to be Sc _T = 0.9. The water vapor source term S _v can be used to account for the increase in continuous phase water vapor mass from evaporating droplets. In the mass transport relation, the transport of thermal energy due to diffusion was excluded based on Lewis numbers close to one for both air and water vapor. For the two species considered, the mass fraction of air was evaluated as Y _a = 1.0-Y _v .

For incompressible temperature dependent flow, the thermal energy equation is typically solved. However, compressible temperature dependent flow can be better approximated by the solution of the combined thermal and mechanical energy equation

In this composite energy statement, the total energy e is represented as

The total enthalpy (h) is the mass fraction weighted sum of the constituent enthalpies for each of the species (h _s ). On the right-hand side of Equation (5a), the first term represents conductive transport while the second accounts for energy transport due to species diffusion. In the diffusion term, the summation is performed over the species index, s, for air and water vapor.

Properties of the multicomponent mixture are based on available empirical data and ideal gas assumptions. The variable density of the air and water vapor mixture is calculated using the multicomponent ideal gas law

where R _u is the universal gas constant and MW _s is the molecular weight of each species. Considering other transport properties, the specific heat of the multicomponent gas mixture has been computed from a mass fraction mixing law. Mixture properties for thermal conductivity (κ) and viscosity (μ) have been determined using the ideal gas mixture law from kinetic theory (Vincenti and Kruger 1986). The binary diffusion coefficient of water vapor in air has been calculated from (Vargaftik 1975)

The temperature dependent saturation pressure of water vapor is determined from the Antoine equation (Green 1997)

which is considered to be more accurate than the Clausius-Clapeyron relation across the broad range of temperatures considered. Use of the empirical property expressions defined in Equations (7) and (8) has been limited to the original temperature ranges for which these relations were determined (Vargaftik 1975; Green 1997). Limiting values were employed for temperatures beyond the original experimental conditions, which are necessary near the capillary tip.

The relative humidity of the mixture entering the induction port is dependent on the local temperature and mixture density. The local relative humidity will affect water vapor evaporation and condensation on the surface of droplets and on the walls of the USP induction port. Relative humidity of the ideal gas mixture can be expressed

where R _v is the gas constant of water vapor and ρ is the mixture density calculated from Equation [6]. The saturation vapor pressure P _v,sat is computed using Equation [8].

To determine the surface deposition of water vapor on the walls of the USP induction port, a constant wall temperature of 25°C has been assumed to approximate experimental conditions. Relative humidity at near-wall locations has been calculated for control volume centers adjacent to wall surfaces. If the near-wall relative humidity is less than one, then no water vapor deposition occurs and a zero mass flux boundary condition is applied. If the near-wall relative humidity is greater than one, then condensation of the supersaturated vapor is assumed. In this scenario, saturation conditions are evaluated at the wall surface. The mass fraction of saturated water vapor for the 25°C wall is Y _v,wall = 2.62 × 10^{− 2}. The mass flux from the near-wall location to the wall surface is calculated as

Discrete Phase Transport Equations

A polydisperse distribution of droplet diameters is expected for the aqueous CAG system potentially ranging from 100 nm through 60 μ m. The Stokes number (St = ρ _p d _p ² C _c U/18μ D ≪ 1) of these droplets based on USP induction port inlet conditions can range from 1.73 × 10^{− 6} to 6.0 × 10^{− 3}. Other characteristics of the aerosols considered included an assumed droplet density ρ _p = 1.00 g/cm³, a density ratio α = ρρ _p ≈ 10^{− 3}, and a droplet or particle Reynolds number (Re _p = ρ| < eqid25 > |d _p /μ) ranging from Re _p ≪ 1 to Re _p < 100. In the vicinity of the capillary tip, a low liquid to water vapor mass concentration has been assumed, i.e., c = 5%. The resulting momentum coupling parameter (c/(1 + St)) between the discrete and continuous phases for a mean droplet diameter of 3 μ m is 0.0499 (Crowe et al. 1998). In regions far from the capillary tip, the maximum momentum coupling parameter for 3 μ m particles is 0.043. As a result, the effect of droplet inertia on the continuous phase has been neglected in this study, i.e., one-way coupled multiphase flow.

The transport and deposition of smaller droplets are considered to be the result of inertial and diffusional effects. Droplets greater than approximately 1 μ m will deposit primarily by impaction and sedimentation. Turbulent fluctuations will also significantly affect the deposition of aerosols throughout the size range considered. To address this broad range of deposition mechanisms, a Lagrangian particle tracking method was employed (Longest et al. 2004). The Lagrangian transport equations for particles ranging from 100 nm through 60 μ m can be expressed

Here v _i and u _i are the components of the particle and local fluid velocity, respectively, and g _i denotes gravity. The characteristic time required for a particle to respond to changes in fluid motion, or the particle relaxation time, is expressed as τ _p = C _c ρ _p d _p ²/18μ, where C _c is the Cunningham correction factor for submicron aerosols. The pressure gradient or acceleration term for aerosols was neglected due to small values of the density ratio (Longest et al. 2004). The drag factor f, which represents the ratio of the drag coefficient to Stokes drag, is based on the expression of Morsi and Alexander (1972)

where the a _i coefficients are constants for smooth spherical particles over the range of particle Reynolds number considered in the current study, i.e., 0 ≤ Re _p ≤ 100. The effect of the lubrication force, or near-wall drag modifications, are expected to be reduced for the aerosol system of interest in comparison to liquid flows due to near-wall non-continuum effects (Longest et al. 2004). Therefore, this term was neglected for the simulations described here. Saffman style lift was considered based on the force term described in Longest et al. (2004).

To model the effects of turbulent fluctuations on particle trajectories, a random walk method was employed (Gosman and Ioannides 1981; Crowe et al. 1996; Loth 2000; Matida et al. 2000). This method assumes that the fluid velocity used in Equation (11) is constant during the time that a particle spends in an eddy and is taken as

To determine the fluctuating component of the instantaneous velocity, u′ _i is selected from a Gaussian distribution with a variance of 2k/3. The time that the particle spends in the eddy is the minimum of the eddy crossing time (T _c ) and the random eddy lifetime (T _e ) where

In the above equations, ℓ is the characteristic eddy length scale taken to be 0.15 k ^0.5ω and r is a uniform random number between 0 and 1. The primary limitation of this eddy interaction model is that it does not account for reduced turbulent fluctuations in the wall-normal direction, which may result in an over-prediction of deposition (Kim et al. 1987; Matida et al. 2000, 2004). To better approximate turbulent effects on particle deposition, an anisotropic turbulence model has been applied where the near-wall fluctuating velocity is calculated as (Wang and James 1999; Matida et al. 2004)

In the above equations, n is the wall-normal coordinate and u _τ is the turbulent friction velocity (Wilcox 1998). The wall-normal damping function (f _n ) is typically evaluated from the wall to a maximum n ⁺ value ranging from 10 to 100; otherwise it is assumed to be 1.0. In this study, f _n was evaluated for n ⁺ values ranging from 0 to a maximum of 50.

The effect of Brownian motion on the trajectories of submicron particles has been included as a separate force per unit mass term at each time-step. This force has been calculated as (Longest and Xi 2007)

where ς _i is a zero mean variant from a Gaussian probability density function, Δ t is the time-step for particle integration, and m _d is the mass of the particle. Assuming dilute concentrations of spherical particles, the Stokes-Einstein equation was used to determine the diffusion coefficients for various size particles as

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

The simulation of liquid aerosol transport and deposition requires specifications of particle starting locations, initial velocities, spatial profiles and size distributions. The initial starting location of the liquid droplets was assigned to be the capillary tip. Even though droplets are unstable at the tip, initializing droplets at this location was considered necessary in order to capture the inertial effects of the rapidly expanding high pressure flow field. However, the initial droplet velocity at the capillary tip is not known. This initial droplet velocity will likely be in the range of zero to the sonic condition. While the gas phase is moving at sonic speed, initializing droplets at this velocity was found to significantly overestimate their momentum and deposition. In this study, initial droplet velocity was assumed to be zero at the capillary tip. Because of the very large difference in initial droplet and flow velocity, the aerosol droplets were quickly accelerated in the direction of the local gas flow. The spatial distribution of droplet concentration at the capillary tip was assumed to be uniform. Details of specifying a spatially uniform droplet distribution have been discussed previously in Longest and Vinchurkar (2007a).

The initial size distribution for the droplets was based on experimental laser measurements made 1.5 cm downstream of the capillary tip. Experimental size estimates at this location were expected to be most accurate due to reduced droplet velocity and the existence of a wide spray plume. However, to adequately model droplet inertia, droplet conditions are needed at the capillary tip. It is expected that some droplet size change occurs between the capillary tip and the location of laser measurement due to condensation and droplet breakup. In this study, the polydisperse aerosol size distribution measured 1.5 cm downstream was assumed to be the initial size distribution at the capillary tip. The mechanisms of droplet condensation and breakup were excluded from the present model in order to avoid overestimating these effects. The validity of this initial condition will be assessed based on comparisons of numerical and experimental deposition results.

A potential difficulty in simulating the polydisperse aerosol distribution considered with a Lagrangian tracking model is the number of droplets required to resolve a broad size range. For example, for every one 40 μ m droplet simulated, approximately 2.4 × 10⁶ 0.3 μ m droplets must be included to represent an equal mass balance. To address this issue, representative groups of tracer droplets were considered to approximate the experimentally reported polydisperse distribution. The experimental droplet distribution was divided into 40 size bins and an equal number of droplets were simulated to represent the behavior of each bin. The final deposition results were then corrected to reflect the behavior of the actual polydisperse droplet distribution and to provide an accurate estimate of drug mass. Because a representative tracer size distribution has been employed in the numerical simulation of droplet trajectories, mass and momentum exchange from the droplets to the continuous phase were neglected.

The aqueous droplets considered in this study will evaporate until either a solid drug particle remains or equilibrium has been reached with the surrounding environment. The presence of the drug in the droplets should reduce the evaporation rate (Finlay 2001). However, the degree of evaporation and amount of drug remaining in the droplets at the site of laser size analysis is unknown. To address this issue, two approximations of aerosol size change have been considered. In the first simulation, both condensation and evaporation were neglected and the aerosols were assumed to behave as solid particles with the same density as water (non-evaporating particles). In the second simulation, a liquid to water vapor mass concentration of 5% was assumed at the site of laser particle size sampling. The initial concentration of albuterol in the original solution was 0.5% w/v. Assuming 5% liquid water in the initial aerosol stream resulted in droplets that contained 10% drug and 90% liquid. These droplets were assumed to evaporate as pure liquid water until 10% of the original mass remained to approximate the non-volatile drug component. After evaporation of the volatile species, the droplet was treated as a non-evaporating drug particle of unit density. In both the non-evaporating and evaporating simulations, the amount of drug in each particle or droplet was assumed to be 10% of the initial mass. Comparisons of these two approaches with experimental deposition data and exit size distribution results were used to determine the likely impact of neglecting solute effects on droplet evaporation and the need to consider this aspect in future studies.

The evaporation model considered in this study assumes that droplets evaporate when the surrounding local relative humidity value falls below one. As described above, aerosol growth due to condensation has been neglected. In this article, droplet evaporation was simulated based on the conservation of energy and mass equations with a semi-empirical rapid mixing model (RMM) approach (Longest and Kleinstreuer 2005). Conservation of energy for an immersed droplet, indicated by the subscript d, under RMM conditions can be expressed (Longest and Kleinstreuer 2005)

In the above equation, m _d is the droplet mass, C _pd is the composite liquid specific heat, q _conv is the convective heat flux, n _v is the mass flux of the evaporating chemical species at the droplet surface, and L _v is the latent specific heat of the liquid component. The integrals are performed numerically over the droplet surface area, A.

Conservation of mass for an immersed droplet based on the evaporating flux can be expressed as

For a semi-empirical RMM solution, the area-averaged heat flux is evaluated from

where Nu is the Nusselt number, κ is the thermal conductivity of the gas mixture, and T _∞ is the temperature condition surrounding the droplet. Similarly, the area-averaged mass flux is

where Sh is the non-dimensional Sherwood number, ρ is the gas mixture density, and Y _v,∞ is the water vapor mass fraction surrounding the droplet. This expression includes the effect of droplet evaporation on the evaporation rate, which is referred to as the blowing velocity (Longest and Kleinstreuer 2005). For the semi-empirical RMM solution, assumptions of thermodynamic equilibrium at the droplet surface and ideal gas conditions have been made, such that the gas-phase surface concentration of water vapor can be expressed as

In the above expression, P _v,sat (T _d ) is the temperature dependent saturation pressure of water vapor, calculated from Equation (8).

The non-dimensional Nusselt and Sherwood numbers employed in Equations [17a and b] have been based on the empirically derived expression of Clift et al. (1978)

These correlations are valid for droplet Reynolds numbers up to 400 and include blowing velocity effects.

NUMERICAL METHODS

To solve the governing conservation equations, the CFD package Fluent 6.2.16 (Ansys, Inc.) was employed. User-supplied FORTRAN and C programs were used to calculate the appropriate outlet pressure, surface vapor absorption, relative humidity, near-wall properties, droplet initial conditions, polydisperse mass deposition rates, and outlet size distributions. A summary of user-defined functions (UDF) that were necessary to simulate continuous and discrete phase transport in the CAG system are provided in Table 1. All transport equations were discretized to be at least second order accurate in space. The computational mesh was constructed in Gambit 2.2 (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 four orders of magnitude for each time-step 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 droplet deposition fields. To improve accuracy and to better resolve the significant change in flow scales, all calculations were performed in double precision.

TABLE 1 User-defined functions (UDF) and codes required to simulate the transport and deposition of capillary-generated aerosols.

UDF Module	Description
outlet_pressure	Module to actively modify the outlet pressure to maintain a constant outlet volumetric flow rate of 30 L/min
surface_absorption	Module to calculate vapor phase absorption on the walls in areas where the relative humidity is greater than 1.0
rel_humidity	Module to calculate the relative humidity based on local temperature and mass fraction
near_wall_interpolation	Module to interpolate near-wall particle velocities based on nodal values in a turbulent field
near_wall_kinetic_energy	Module to calculate near-wall turbulent kinetic energy based on i) anisotropic effects and ii) interpolation near the surface
turbulent_tracking	Module to track turbulent aerosol dispersion using the anisotropic kinetic energy relation
sect_deposit_report and outlet_report_histogram	Modules to post-process results, report deposition efficiencies and outlet particle sizes in the sectioned USP throat

Construction of the 3D hexahedral computational mesh for the CAG system shown in Figure 3 required a broad range of control volume sizes. Grid resolution near the 57 μ m capillary is shown in Figure 4. Consecutive control volume sizes were allowed to vary by no more than a factor of 1.15. A sufficiently dense mesh was applied in all near-wall regions so that the maximum y⁺ value of the LRN k-ω model was maintained on the order of approximately 1 or less at the first grid point above the wall. Considering grid convergence, 3D meshes consisting of approximately 325,000, 550,000, and 700,000 cells were tested. Due to variability in the meshes of the geometry, these approximate mesh sizes have been rounded to the nearest 25,000 control volumes. Negligible variations in the parameters of interest, i.e., the velocity field, local particle deposition efficiencies and outlet particle size distributions, were observed between the two finest meshes considered. As such, the final mesh for the CAG computational geometry contained approximately 550,000 cells.

FIG. 4 Multiblock structured computational mesh near the 57 μ m capillary tip. The USP inlet diameter is 31.8 mm, resulting in a three order of magnitude change in spatial dimensions that must be resolved.

A two second pulse of aerosol generation fills the USP induction port with a stream of vapor, evaporating liquid droplets, and significant momentum in a highly time dependent way. Therefore, a transient solution was employed for numerical analysis of the CAG system. The time-step selected was intended to capture major temporal changes in the flow field as well as the variable relative humidity profile. High frequency fluctuations and turbulent eddies were approximated by the time-averaged RANS equations and the LRN k-ω model. Temporal discretization was achieved using a second-order implicit method. The continuous phase results were observed to become time-step independent for a constant step size of 0.025 s during the majority of the 2-s period. Halving this time-step resulted in a negligible change in temporal velocity fields, relative humidity and deposition results.

To model the experimental condition of 30 L/min co-flow air pulled through the USP induction port by a vacuum pump, a variable low pressure was defined at the computational geometry outlet downstream from the expanded section. To maintain a constant flow rate during the aerosol generation process, which adds significant momentum to the system, the downstream pressure was adjusted at every time-step with a UDF function (Table 1). Inclusion of the expanded outlet section was found to minimize time dependent fluctuations in flow rate. The total variation in flow rate associated with capillary momentum and air entrainment implementing the downstream pressure algorithm was approximately 5% at the outlet.

Transient particle trajectories were simulated during the two second aerosol generation period using the Fluent 6.2 particle tracking routine. Particles were injected and advanced in time with each flow field time-step. The integration scheme employed to solve Equation [11] was based on the Runge-Kutta method with a minimum of 20 integration steps in each control volume. An error control routine was also employed to actively modify the particle time-step and maintain accuracy, which is comparable to the error control routine described by Longest et al. (2004). Doubling the number of integration steps within each control volume had a negligible (less than 1%) effect on cumulative particle/droplet deposition values. Due to relatively small particle relaxation times, double precision calculations were employed for the trajectory calculations (Longest et al. 2004). In order to produce convergent deposition results within each section of the USP induction port for all particle sizes considered, groups of 4,000 particles were released at each flow field time-step, resulting in a total of 320,000 particles traced over two seconds. Considering that the particle size distribution was divided into 40 bins, there were 100 tracer particles in each size bin at each time-step. The final mass deposition results were scaled to reflect the experimentally determined initial polydisperse size distribution. Doubling the number of tracer particles considered had a negligible impact on both total and sectional deposition results.

As reported by Longest and Xi (2007), the particle tracking algorithm employed by Fluent 6.2 does not appear to interpolate near-wall particle velocities and variables in a turbulent field. This simplification results in a significant over-prediction of deposition for all particle sizes in geometries consistent with the oral airway. In this study, a UDF module was used for the particle tracking algorithm to linearly interpolate near-wall particle velocity and turbulent kinetic energy between the wall adjacent control volume center and the boundary surface (Table 1). The interpolation of variables within control volumes from nodal values has been described in detail by Pepper and Heinrich (1992) and has previously been reported by Longest and Xi (2007).

RESULTS

Continuous Field Variables

Continuous field variables near the capillary after 2 seconds of aerosol generation are shown in Figure 5. Centerline velocities in excess of 100 m/s were observed to extend approximately 1 cm away from the capillary tip (Fig. 5a). At approximately 1.5 cm from the capillary tip, centerline velocities were approximately 40 m/s. Instantaneous midplane streamlines illustrate significant entrainment of the co-flow air due to the momentum of the vapor jet. The streamlines were seen to be pulled into the jet, which may result in mixing, rapid dilution, and particle evaporation. Centerline temperatures at 1 and 1.5 cm from the capillary tip were approximately 40° and 35°C, respectively (Fig. 5b). At a centerline distance of approximately 4 cm, an equilibrium temperature of approximately 30°C was established, which is maintained through the majority of the induction port geometry at 2 s.

FIG. 5 Continuous gas phase variables at the midplane and near the capillary tip including (a) velocity magnitude with instantaneous streamlines and (b) temperature after 2 s for 30 L/min of co-flow air.

A time series of midplane and cross-sectional velocity conditions is shown in Figure 6. The computational system was initialized with a steady condition of 30 L/min co-flow air, as in the in vitro experiments. Vapor flow through the capillary was then initialized at time 0 s and continued through 2 s. At 0.5 s, recirculation zones were observed in the velocity profile near the capillary tip (Fig. 6a). In the downstream vertical section of the induction port, significant vortical flow was seen as a result of the sudden increase in momentum due to capillary activation. This vortical mixing may increase particle deposition, residence time, evaporation, and can temporarily reduce the mass flow rate through the system. At 1.0 s, recirculation had largely disappeared in the velocity profile near the capillary (Fig. 6b). As expected, a zone of recirculation was observed near the inner surface of the induction port downstream of the 90° bend. However, the significant vortical flow observed at 0.5 s was eliminated at 1.0 s resulting in a well defined flow path throughout the geometry. The velocity field observed at 1.0 s remained relatively unchanged for the remainder of the 2 s simulation.

FIG. 6 Contours of velocity magnitude and velocity vectors at the midplane and selected cross-sectional slice locations for (a) 0.5 s and (b) 1.0 s after the capillary was turned on. Highly vortical flow was observed downstream of the capillary tip at 0.5 s due to the sudden increase in momentum associated with start-up. At 1 s, the flow was more evenly distributed. Velocity conditions remain relatively constant after 1 s.

The time history of relative humidity in the induction port with a co-flow rate of 30 L/min is shown in Figure 7. The relative humidity of the initial system and co-flow air was 30%. At 0.5 s, centerline relative humidity values of 1.0 extended approximately 2.3 cm from the capillary tip (Fig. 7a) and values in the region of Slice 2 had not appreciably increased. At 1.0 s, relative humidity values at Slice 1 were approximately 80%, while at Slice 2 they were significantly smaller (Fig. 7b). At 1.5 s the majority of the flow field had reached a relative humidity value of approximately 80% (Fig. 7c) and water vapor had accumulated in near-wall regions of low flow like the inner and outer bends of the induction port (∼ 90% with near-wall values likely higher). At 2.0 s, relative humidity values in Slice 1 were reduced slightly (Fig. 7d), perhaps because of an increase in velocity or an increase in the flow field temperature over time; values at Slice 2 were approximately 80–100%. Based on inlet mass flow rates of water vapor and air, the theoretical relative humidity value of the induction port system should approach unity for steady state, perfectly-mixed conditions. Values at Slice 2 after 2 s are likely below this theoretical limit due to imperfect mixing and upstream wall deposition losses of water vapor where near-wall values of relative humidity exceed one.

FIG. 7 Relative humidity values at the midplane and selected cross-sectional slice locations at (a) 0.5, (b) 1.0, (c) 1.5, and (d) 2.0 s. Values in the downstream section increased from an initial value of 30% to approximately 80–90% after 2 s.

The effects of the time varying humidity fields shown in Figure 7a–d on droplet evaporation are expected to be significant. Droplets moving through the induction port system will have a mean residence time of 0.1–0.15 s. Based on these short residence times, aerosols may fully traverse each of the relative humidity fields before significant changes in the continuous variables occur. Therefore, the mean relative humidity surrounding a droplet while in the induction port may vary between approximately 30 and 80% depending on the time of droplet release. As a result, transient simulations of both continuous phase variables and discrete element transport were necessary.

Relative humidity conditions have been used to estimate the approximate angle of the vapor jet exiting the capillary tip (Fig. 8). It was assumed that the jet will be visible in experiments at a relative humidity value of 1.0, which is the criterion selected for estimating the jet angle. The velocity boundary condition at the capillary tip was specified only in the x-direction (Fig. 3). The jet angle is then allowed to develop as a result of interactions between flow inertia and expansion effects. The resulting jet angle beginning at the capillary tip is observed to be 29° in the presence of 30 L/min co-flow air (Fig. 8). This value compares well with the experimental estimate of Shen et al. (2004), who reported jet angles between 20 and 30° for the CAG. Similarly, Lesniewski and Koch (1998) reported a value of 26° for a highly turbulent round micro-jet. Some variability may occur between these numerical and experimental estimates of jet angle based on how the visible jet is defined.

FIG. 8 Angle of the capillary jet based on a relative humidity value of 1.0 initiated at the capillary tip under 30 L/min co-flow conditions. The computed angle of 29° was similar to the CAG experimental study of Shen et al. (2004).

Droplet Size Distribution and Trajectories

The initial droplet size distribution measured using the laser diffraction system at a location of 1.5 cm from the capillary tip with 30 L/min co-flow air is illustrated in Figure 9. Particle diameters have been grouped into 40 bins with midpoint size values ranging from 0.14 to 60 μ m. The laser diffraction system used experimentally measures geometric particle size and calculates an associated volume fraction size distribution. In this study, unit density aqueous aerosols were assumed allowing for a direct conversion to mass fraction values. The mass median diameter (MMD) of the initial distribution shown in Figure 9 was 3.11 μ m. The majority of the droplet mass (80.2%) occurred between 1 and 10 μ m with the remainder in droplets less than 1 μ m (12.7%) and greater than 10 μ m (7.1%).

FIG. 9 Initial particle size distribution based on mass fraction measured approximately 1.5 cm away from the capillary tip using laser diffraction. Histogram blocks are shown separated for clarity.

Sample trajectories of evaporating droplets released at the capillary tip under the steady flow conditions at 1 second are shown in Figure 10. Steady flow conditions 1 second after capillary activation were selected to provide constant representative velocity (Fig. 6b) and relatively humidity fields (Fig. 7b). Co-flow air was included at a rate of 30 L/min. Droplets in the initial size range of 0.1–1 μ m were seen to be significantly influenced by turbulent fluctuations (Fig. 10a). As expected, these droplets readily follow approximated turbulent eddies due to their low inertia. The effects of turbulent dispersion appeared to be most significant in the mixing region near the capillary tip. Furthermore, submicron particles were observed to evaporate near the turbulent jet. Droplets in the size range of 1–5 μ m were also significantly influenced by turbulence (Fig. 10b). Turbulent dispersion appeared to induce significant deposition in the proximal region of the induction port just downstream of the capillary tip. Exiting particle diameters for these 1–5 μ m droplets appeared to be approximately 3 μ m or less.

FIG. 10 Trajectories and diameters for sample droplets released from the capillary tip after 1 s grouped by initial droplet size range: (a) 0.1–1 μ m, (b) 1–5 μ m, (c) 5–10 μ m, and (d) greater than 10μ m.

As droplet size increased above 5 μ m, turbulent dispersion had less of an impact on deposition (Figs. 10c and d). For particles in the size range of 5–10 μ m, the effects of turbulent dispersion were delayed until approximately the midpoint of the horizontal section (Fig. 10c). However, dispersion driven deposition still appeared to be significant in the horizontal section. Little impaction in the 90° bend was observed for droplets initially between 5 and 10 μ m. Exiting droplet diameters in this size range were predominately 5 μ m and below. For particles greater than 10 μ m, impaction was clearly the primary mechanism for deposition (Fig. 10d). Particles significantly greater than 10 μ m followed angled trajectories induced by the jet expansion and were not visibly influenced by dispersion. The maximum angle of these trajectories was observed to be approximately 30°. Due to significant inertia, these droplets readily deposited in the initial horizontal section of the USP throat. Particles greater than 10 μ m were also observed to deposit in the 90° bend section. However, some particles with diameters in the range of 10 μ m were observed to exit the induction port (Fig. 10d).

Particle and Droplet Deposition

Experimental deposition values for the CAG system under the conditions of interest were determined in both standard and sectioned USP induction ports. For all deposition experiments, a minimum of five replicates were performed and deposition fractions were reported based on drug mass. The average mass fraction of drug deposition in the un-sectioned stainless steel induction port was 13.8% with a standard deviation of 1.3%. For the sectioned induction port created with a rapid prototyping process, the total deposited drug mass fraction was 14.6% with a standard deviation of 1.1%. Deposition in the two ports was not significantly different statistically (P = 0.33; paired t-test) showing that sectioning the USP induction port and using a resin for its construction failed to alter the deposition behavior. Experimental deposition results for each of the four individual sections of the rapid prototyped model are provided in Table 2.

TABLE 2 Sectional induction port drug mass deposition values (%) for particles and droplets

	Section 1	Section 2	Section 3	Section 4	Total
In vitro deposition ^1 (% mass)	4.04	5.11	5.10	0.30	14.6
CFD non-evaporating particles (% mass)	3.93	3.73	4.92	0.47	13.1
% Error for particles	2.72	27.0	3.52	56.7	10.3
CFD evaporating droplets (% mass)	4.01	3.97	5.01	0.29	13.3
% Error for droplets	0.75	22.3	1.76	3.33	8.90

Predicted deposition locations as well as comparisons between experimental and numerical deposition values on a sectional basis are shown in Figure 11. Models for non-evaporating particles (Fig. 11a) and evaporating droplets (Fig. 11b) were considered. Deposition patterns with and without droplet evaporation appeared largely similar. As observed with sample trajectories, the larger droplets considered evaporated little and deposited by impaction due to either angled trajectories associated with the jet or the effect of the 90° bend. For particles and droplets on the order of 10 μ m, deposition occurred due to impaction, which extended down the outer vertical section of the induction port. Deposition due to turbulent dispersion of particles ranging from submicron up to approximately 5 μ m was observed to be distributed somewhat uniformly among Sections 1–3. However, potentially enhanced deposition of smaller particles in Section 1 was hidden by the larger particles. The deposition locations shown in Figure 11 are based on the representative uniform distribution of tracer particles used to simulate the polydisperse aerosol of interest. The associated deposition fractions reported in Figure 11 have been scaled based on the experimentally determined initial particle distribution and the mass fraction of non-volatile drug estimated to initially be in each particle.

FIG. 11 Sectioned deposition results of a transient polydisperse aerosol released over 2 s from the capillary tip for (a) non-evaporating particles, and (b) evaporating droplets. The reported CFD mass fractions are based on the original mass of drug estimated to be in the droplets.

Comparisons of in vitro and computational deposition values in each section of the induction port have been illustrated in Figure 11 and tabulated in Table 2 with relative error values. The non-evaporating particle model closely matches the experimental results in Sections 1–3 with a relative error as low as 2.7% in Section 1. However, the experimental and predicted values in Section 4 (0.30 vs. 0.47, respectively) result in a relative error value of 56.7%. The relative error for total deposition based on the experimental and computational results without droplet evaporation is 10.3%. For the evaporating droplet model, agreement between experimental and computational results within individual sections is improved with relative errors ranging from 0.75% to 22.3%. The relative error in total deposition with the particle evaporation model is 8.9%. As a result, both aerosol models appear to adequately match the experimental deposition results with the evaporating model providing improved agreement in Section 4.

Droplets that deposit in the induction port are expected to not be available for lung delivery. Therefore, it is important to consider the size distribution of deposited aerosols. For the droplet evaporation model, the distribution of deposited sizes is shown in Figure 12 on a mass fraction basis. Approximately 20% of the deposited mass in the induction port occurs due to droplets less than 1 μ m (Fig. 12). Little deposited mass is observed for droplets in the range of 1–4 μ m. The remaining majority of deposited mass arises from droplets greater than approximately 5 μ m, as expected.

FIG. 12 Distribution of deposited mass fraction based on final droplet size. A majority of the deposited mass results from particles that are either 1 μ m and below or greater than approximately 5 μ m.

Distribution of Exiting Particles and Droplets

The mass fraction distribution of particle sizes exiting the induction port is shown in Figure 13. Experimental results were based on NGI measurements of drug mass. Equivalent particle size bins were then used to report the computational results for the particle and droplet evaporation models. For aerosols greater than approximately 1 μ m, the non-evaporating particle model matched the experimental results very well. The resulting MMD for the experimental and computational models are 3.06 and 3.07 μ m, respectively, resulting in a relative error of 0.33%. However, the non-evaporating particle model failed to predict the correct size trend for submicron aerosols. In contrast, the evaporating droplet model matched the trend in aerosol size distribution for the submicron range through approximately 10 μ m. Similarities in the particle size distribution between the experiment and droplet model include an elevated mass fraction for aerosols less than 0.5 μ m, a minimum mass fraction around 0.76 μ m, an a peak in mass fraction around 3 μ m. However, it appears that the droplet evaporation model may under predict the mass fraction of aerosols in the range of 2 μ m. The MMD for the evaporating model was 3.45 μ m, resulting in a relative error of 12.7%. A shift in the peak of the evaporating droplet profile in the direction of larger particles is likely due to a complex interaction of evaporation, deposition, and the initial polydisperse and multimodal particle size distribution. A comparison of the particle and evaporating droplet models indicates that the droplet model provides a better match to experimental exit conditions up to an aerosol diameter of approximately 1 μ m. For particle sizes greater than 1 μ m, the particle model provides a better estimate of the outlet size distribution (Fig. 13). Furthermore, aerosol size changes within the impactor due to continued evaporation and interactions with the side walls may account for some differences between the model results and experimental data.

FIG. 13 Mass fraction distribution of droplets exiting the USP induction port based on measurements made in the Next Generation Impactor (NGI) and CFD simulation with and without evaporation.

DISCUSSION

In this study, a CFD model for the transport and deposition of capillary-generated aerosols in the USP induction port was developed and tested with direct comparisons to experimental results. The numerical and experimental systems considered were nearly identical for capillary aerosol generation in a standard USP induction port geometry with 30 L/min of co-flow air at ambient conditions. The droplet size distribution exiting the capillary tip was measured using laser diffraction and implemented as the initial polydisperse distribution in the numerical model. Comparisons between numerical and experimental results were made for the jet angle near the capillary tip, the mass fraction of drug depositing in a sectioned USP throat model, and the outlet particle size distribution. Aerosols were assumed to behave as either non-evaporating particles or evaporating liquid droplets of water.

Results of this study indicate good agreement between predictions of the numerical model and experimental measurements. Based on ideal gas assumptions at the capillary inlet, the jet angle was predicted to be approximately 29°, which was consistent with earlier reports (Shen et al. 2004). Considering aerosol deposition, the experimentally measured total mass fraction in the sectioned USP induction port was 14.6% as compared to 13.1% for the particle model and 13.3% for the evaporating droplet model. The latter model only offered a small improvement in prediction for droplets depositing in Section 4 of the induction port. Predicted mass median diameters exiting the induction port for the particle and evaporating droplet models were 3.07 and 3.45, respectively, in comparison to an experimental value of 3.06. Furthermore, the exiting particle size distribution for non-evaporating particles provided the best match to experimental results for aerosols greater than 1 μ m showing that the evaporating droplet model overestimated the likely size changes for CAG aerosols passing through and exiting the USP induction port.

A significant finding of this study is that a user-enhanced CFD model can accurately predict the transport and deposition of capillary-generated aerosols in a standard 90° bend testing system. Simulation of capillary-generated aerosol transport and deposition presents a challenging numerical problem due to the breadth and complexity of the associated multiphase flow field. The continuous phase is complicated by the presence of transitional or fully turbulent flow, compressible effects, locally supersonic conditions, temperature dependencies, and the presence of multiple mixing species. Furthermore, the continuous phase includes geometric scales from the capillary tip (57 μ m) through the USP induction port diameter (19–31.8 mm), which span a three order of magnitude change in spatial dimension. The inclusion of a multiphase aerosol spray introduces the challenges of turbulent particle dispersion, evaporation and condensation of multicomponent droplets, nucleation and coagulation. Varying relative humidity values associated with the filling of the geometry with a vapor phase and the sudden increase in momentum associated with capillary activation prevent the assumption of steady flow conditions.

In order to address a number of factors associated with the CAG system, user-defined functions were developed to supplement the existing CFD model. User-enhanced modules that were necessary to successfully match the experimental results included functions for controlling the outlet pressure to remain consistent with experimental conditions, a model for the calculation of relative humidity, and evaluation of vapor deposition on the walls. The inclusion of modules that interpolated near-wall particle quantities (Longest and Xi 2007) and accounted for near-wall anisotropic effects were necessary to match the experimental deposition results. Still, a number of simplifying assumptions had to be made such as the exclusion of discrete phase effects on the continuous field. Nevertheless, good agreement was observed between the predicted results and experimental measurements for the deposition of drug mass fraction in each section of the induction port and the MMD of exiting particles.

The numerical model developed in this study for the transport and deposition of capillary-generated aerosols can be applied to improve the delivery of pharmaceutical aerosols to the lung. Based on the specific CAG conditions considered, it was found that a significant portion of submicron aerosols (approximately 20% by mass) deposited in the induction port geometry (Figure 12). This high fraction of submicron particle deposition was the result of turbulent mixing and dispersion in the jet near the capillary tip. Reducing the turbulence levels by changing the specific CAG conditions may help to avoid these losses. Furthermore, significant theoretical differences in the relative humidity fields within the induction port were observed between 1 and 2 seconds following activation of the CAG system (Fig. 7). Activation of the CAG process for 1 second may enhance drying of these aerosols and further enhance the delivery of fine droplets. Analysis of the turbulent jet near the capillary tip and the associated turbulent particle dispersion in this region can also be used to develop mouthpiece designs to mitigate unwanted deposition.

In order to better evaluate the delivery of respiratory aerosols to the lung in pharmaceutical applications, more accurate geometries of the respiratory tract may be necessary. Several researchers have reported differences in deposition between highly simplified and more realistic oral airway models (Li et al. 1996; Xi and Longest 2007b, 2007a). In contrast to the USP induction port, more physiologically realistic models of the oral airway geometry have been developed and evaluated (Stapleton et al. 2000; Zhang et al. 2005; Zhang, Chia, and Finlay 2006; Xi and Longest, 2007b; Xi and Longest, 2007a).

A significant assumption made in this study was that water vapor exited the capillary tip entirely in the vapor phase, and that droplets were then rapidly formed by nucleation and coagulation downstream. The actual partitioning of water vapor and liquid droplets would be excessively difficult to determine experimentally at the 57 μ m capillary tip. However, the assumption of 100% water vapor exiting the capillary tip affects both the momentum of the system and the reported relative humidity fields. To better determine whether the momentum of the system was accurately captured by this boundary condition, a second numerical test case was conducted in which the droplet to vapor mass fraction at the capillary tip was assumed to be 50%. For this reduced momentum system, total deposition for the droplet model was decreased by approximately 63% and deposition in Sect. 1 was reduced by a factor of 3.4. Therefore, it appears that to accurately model the free expansion process, and the associated vapor and particle momentum, the assumption of 100% water vapor at the capillary tip is reasonable. Regarding relative humidity, the one-way coupled model employed in this study neglects vapor loss due to the formation of droplets. As a result, the relative humidity fields reported in this study represent an upper maximum of possible values. Still, the computational model appears to over-predict droplet evaporation (Fig. 13).

A primary limitation of the current aerosol transport model is the approximation of droplet evaporation characteristics. Modeling the aerosol droplets as non-evaporating particles is an oversimplification of expected physical behavior. In contrast, assuming that these aerosols evaporate as pure liquid droplets neglects the effect of the drug on the evaporation process. This assumption may be valid for small droplets with low drug concentrations. However, as droplet evaporation becomes more significant, the concentration of the non-volatile drug will increase and reduce the evaporation rate. Zhang et al. (2006a) have recently reported that saline droplets are affected by the presence of the solute at concentrations of approximately 10%. Moreover, the presence of a solute drug may limit the total amount of evaporation possible in a given relative humidity environment such that an equilibrium diameter is reached before full drying has occurred.

Based on results of this study, droplet evaporation characteristics do not appear to have a significant effect on sectional deposition values. The primary deposition mechanisms for the CAG and induction port system appear to be turbulent dispersion of smaller particles and the impaction of larger particles due to either the jet angle or 90° bend. It can be concluded that the significant evaporation of smaller droplets (approximately less than 5 μ m) does not alter deposition rates due to turbulent dispersion. Droplets larger than 5 μ m were not observed to significantly evaporate (Fig. 10c). Furthermore, the deposited mass fraction is calculated as the amount of non-volatile drug that deposits. Therefore, it is reasonable that droplet evaporation does not appear to largely affect the deposition result.

In contrast to deposition, the model used for droplet evaporation did significantly influence the size distribution of aerosols exiting the USP induction port (Fig. 13). It was observed that the evaporating droplet model provided a good match to the experimental outlet size distribution for aerosols that were 1 μ m and smaller. However, the non-evaporating particle model provided a significantly better estimate of the outlet size distribution for aerosols larger than 1 μ m. Differences in the experimental results and model predictions of outlet size distribution may result from a number of assumptions in the droplet evaporation model. As described above, neglecting the effect of the drug solute on droplet evaporation may significantly overestimate evaporation. The assumption of a 5% liquid to vapor mass fraction in the initial aerosol stream was used as an estimate of predominately vapor conditions to determine the corresponding initial mass of drug solute in each droplet. However, assuming an initial liquid to vapor mass concentrations of 10 and 15% in additional simulations did not appreciably alter the droplet outlet distribution results. Perhaps the most significant simplification of the current droplet model is the omission of the evaporating droplet vapor in the surrounding relative humidity field. Finlay and Stapleton (1995) showed that a cloud of evaporating aerosols will affect local relative humidity values and continued evaporation of individual droplets for inhaled respiratory aerosols. Accounting for evaporating droplet mass in the surrounding vapor phase is expected to impede the diameter change of larger droplets more than smaller droplets. As a result, the inclusion of two-way mass coupling should significantly improve the performance of the evaporating droplet model in predicting outlet particle size distributions for aerosols greater than 1 μ m. Furthermore, future studies should also consider condensation effects for cases in which the relative humidity exceeds one.

A major difficulty associated with predicting solute effects is the lack of data regarding the vapor pressure characteristics of different drug aerosols. Raoult's law may provide a first approximation, but this relation is theoretically only valid for dilute solutions. While concentration dependent expressions have been developed for saline solutions (Cinkotai 1971) and selected pharmaceutical aerosols (Peng et al. 2000) more information of this nature needs to be made available for albuterol and other formulations.

In addition to a simplified droplet evaporation model, a number of other assumptions were made for these numerical simulations. For the continuous phase, ideal gas approximations have been used to estimate the boundary conditions at the capillary tip and to evaluate pressure and density coupling in the flow field. The inclusion of a compressibility correction factor did not significantly alter the boundary conditions implemented at the capillary tip. However, the use of ideal gas assumptions to solve the flow field will likely result in under and over estimations of variables in the microscale flow near the capillary tip. Corrections for non-continuum gas effects have not been applied to the micro-jet employed by this aerosol delivery system. Downstream impedance and resistance of the experimental setup have been approximated by the inclusion of an equivalent volume. Furthermore, the two-equation LRN k-ω turbulence model only approximates the effects of turbulent eddies and time scales. Significant assumptions made in the discrete phase model include neglecting a number of spray phenomena such as nucleation, condensation, coagulation, and breakup. Reductions in evaporation rate due to the Kelvin and Fuchs effects have been neglected. Shear stress conditions may also significantly impact droplet evaporation rates. In order to enable the simulation of a polydisperse aerosol cloud over a broad range of particle sizes, the effects of the discrete phase on the flow field have been neglected in this initial model. Even so, good agreement was observed between the numerical simulations and experimental measurements of deposition and outlet MMD. As a result, it appears that the current model adequately captures the transport phenomena most responsible for droplet transport and deposition in this CAG system. Future studies are required to evaluate the effects of additional variables on the transport and deposition of capillary-generated aerosols.

A number of recent studies have employed large eddy simulations (LES) to model flow fields and particle transport in geometries consistent with the upper respiratory tract (Breuer et al. 2006; Matida et al. 2006). Large eddy simulation provides a powerful numerical technique that directly simulates a majority of the eddies and time scales associated with turbulence while modeling the smaller scales to improve efficiency. However, application of LES is expensive from a computational perspective. Grid resolutions for LES applications are approximately four times denser than the LRN k-ω model. Moreover, time-steps to simulate large eddies in the oral airway using LES are on the order of 1 μ s (Matida et al. 2006). This time-step size may need to be reduced further to capture the turbulent effects of the high Reynolds number micro-jet considered in this application. Implementation of time-steps on the order of 1 μ s and below is currently not practical for simulating transient effects that occur over a 2-s period. In contrast, results of this study indicate that an anisotropic turbulence approximation in conjunction with a transient LRN k-ω model can adequately approximate particle transport and deposition in the complex flow field environment of the CAG system. Therefore, the anisotropic turbulence model appears to be an effective alternative to LES in cases where significant time periods (on the order of seconds) must be resolved.

In conclusion, a CFD model has been developed and tested for the transport and deposition of capillary-generated aerosols. Good agreement was found between the experimental results and numerical predictions for aerosol deposition rates in a sectioned USP induction port and the exiting aerosol MMD. The inclusion of user-defined functions was necessary to account for near-wall anisotropic turbulence and other factors specific to the CAG system. Droplet evaporation was shown to have a minor impact on particle deposition and a significant influence on particle sizes exiting the induction port. Future applications of this model include optimization of the CAG system for enhanced lung deposition. Future studies of capillary-generated aerosols are needed to evaluate factors such as the effects of drug solutes on aerosol evaporation, the influence of downstream impedance, and the effects of the droplet phase on the continuous field.

Acknowledgments

This work was supported by Chrysalis Technologies, a division of Philip Morris USA. The authors are grateful for the critical advice of Kenneth Cox of Chrysalis and the help of Jeffrey Hinkins of VCU's Department of Mechanical Engineering for the design and fabrication of the sectioned USP induction port.

Notes

¹Mass of albuterol deposited expressed as a percentage of the theoretical delivered dose from the CAG.