A DNS Study of Aerosol Deposition in a Turbulent Square Duct Flow

A particle-laden turbulent flow through a square duct was simulated using a direct numerical solution of the Navier-Stokes equations coupled with Langrangian particle tracking. Computations of particle transport were employed to elucidate the mechanisms by which particles with varying inertia deposit to the walls of a square duct. Gravity was neglected and a one-way coupling was assumed between the particles and the fluid. The computational results demonstrate that, although the aerosol penetration through a square duct is not significantly different than through a circular pipe, there exist differences in the transport and deposition mechanisms. Most notably, the off-axis secondary flows unique to the square duct preferentially deposit higher-inertia particles closer to the corners of the duct. By contrast, the same secondary flows act to suppress the deposition of lower-inertia particles to the duct corners by efficiently transporting them back towards the duct core before deposition can occur.

INTRODUCTION

The deposition of particles suspended in a turbulent flow has been widely studied over the last fifty years, beginning with the pioneering work of Friedlander and Johnstone (1957). Historically, the primary motivation for studying deposition in internal flows was to understand the transport of particles through aerosol sampling tubing. Consequently, the bulk of the literature concerns the turbulent flow through ducts having circular or annular cross section. More recently, there has been increasing interest in the turbulent transport of particles through ducts having rectangular cross section due to the relevance of such flows to ventilation systems in buildings (Sippola and Nazaroff 2003). Understanding deposition in these more irregularly shaped ducts is critical to evaluating the penetration of outdoor air pollution to indoor environments and the recirculation of pollutants and allergens within buildings. Unfortunately, the flow complexities introduced by corners present in the duct cross section has precluded a rigorous theoretical analysis of deposition in such flows. In the present study, we seek to understand the rates and mechanisms of turbulent deposition in ducts of square cross section, focusing on the role of the turbulence-driven secondary flows that develop in the corners of the duct.

Surprisingly, the turbulent deposition of particles in circular pipes was not thoroughly understood for decades despite extensive experimental data. It was believed that the deposition rate could be fully described by the flow very close to the wall, and that the inertia-dependent response to near-wall coherent structures was completely responsible for the observed inertia dependence of the deposition rate. Friedlander and Johnstone (1957) were the first to point out the importance of inertia to the deposition of particles suspended in turbulent flows. It was postulated that deposition was driven by flow a finite distance away from the wall, where particles acquired enough momentum from turbulent eddies to fly through the viscous sublayer and impact the wall. Over the next few decades, this so-called “free-flight” deposition model was verified through experimental measurements of the input and output concentrations (C _in and C _out) of a monodisperse aerosol flowing through a pipe. From these measurements, a particle mass flux, J _w , to the wall may be defined, from which one may calculate a characteristic velocity associated with deposition:

where C is the mean aerosol mass concentration, and the deposition velocity is normalized with the friction velocity, u _τ. When experimentalists plotted k _d ⁺ as a function of particle inertia, the deposition velocity was found to increase by four orders of magnitude over a range of particle inertia spanning two orders of magnitude (Liu and Agrawal 1974). This range, later termed the “diffusion-impaction” regime, corresponds to dimensionless particle relaxation times of 0.2 < τ _p ⁺ < 20, where

and ρ _p and D _p are the particle density and diameter, and ρ and μ are the gas density and dynamic viscosity. Indeed, this observation is consistent with a free-flight model of particle deposition. Particles with relaxation times smaller than τ _p ⁺ = 0.2 deposit primarily due to turbulent diffusion, and are thus driven to the wall down a concentration gradient. Particles with relaxation times larger than τ _p ⁺ = 20 correspond to the inertia-moderated regime. These particles acquire enough momentum in the core of the duct to completely penetrate the boundary layer and deposit on the wall, and are thus unaffected by near-wall coherent structures. The importance of the near-wall coherent structures to particle deposition in the diffusion-impaction regime is undeniable. Numerous studies have demonstrated a significant increase in deposition rate with an increase in the near-wall turbulence intensity (see, for example, Fan and Ahmadi 1993 and Botto et al. 2005). However, the observed magnitude of the near-wall turbulence intensity is far too low to result in the rates of deposition reported from experiments. This was first suggested by Davies (1966), who noted that the particle's “free-flight velocity” should be equal to the local near-wall value of the root mean square velocity fluctuations. Others argued that anomalous turbulent bursts into the near-wall region were responsible for the high free-flight velocities required to match experimentally observed deposition rates.

Brooke et al. (1994) addressed this controversy in a direct numerical simulation of turbulent pipe flow coupled with Lagrangian particle tracking. The authors of that study note that progress in the development of deposition models beyond the free-flight model was hindered by the lack of detailed information provided by existing experiments. Consequently, they used computations to demonstrate that the enhanced wall deposition rate could be attributed to the inertia-dependent transport of particles from the core of the pipe to the near-wall region, where they tended to concentrate. The excessively high concentration of particles near the wall coupled with the lower turbulence intensity present in that region yielded the deposition rates similar to those observed for decades. In the absence of a lateral mean flow, particles migrate from the core to the walls by turbulent diffusion and, in the case of finite-inertia particles, by turbophoresis—the inertia-dependent process by which particles travel down a gradient of turbulence intensity. It is turbophoresis that is responsible for the preferential concentration of higher-inertia particles close to the wall.

A formal theory of particle deposition in a circular pipe including turbophoresis was presented by Young and Leeming (1997), who assumed a particle phase continuum, and solved the Reynolds-averaged mass and momentum conservation equations. This model represents a valuable progression from conventional free-flight models, because it describes the relative contributions of diffusive and turbophoretic transport throughout the bulk of the flow. Unfortunately, this model cannot be easily applied to deposition in the presence of off-axis secondary mean flows, such as those present in a square duct.

There exist several computational studies focusing on turbulent single phase flow through a square duct including direct numerical simulations (Gavrilakis 1992; Huser and Biringen 1993) and large eddy simulations (Madabhushi and Vanka 1991; Kajishima and Miyake 1992; Vazquez and Metais 2002). All of these demonstrate the turbulence-driven secondary flows that act to transfer fluid momentum from the duct core towards the walls. There exist only two computational studies of particle-laden turbulent flow in a square duct. The first, an LES study reported by Winkler et al. (2004), focuses on the preferential concentration of finite-inertia particles, especially near the walls of the duct. The second, a DNS study recently reported by the authors of this study (Sharma and Phares 2006), focuses on longitudinal and lateral single particle-and pair-dispersion statistics. That study demonstrated that the secondary flows act in combination with turbophoresis to transport particles between the duct core and the walls. The present study represents the first computations of particle deposition in a turbulent square duct flow. Here, we focus on the rate of deposition as a function of the location of initial particle release, and the spatial distribution of deposited particles along the duct walls.

NUMERICAL PROCEDURE

Flow Simulation

The numerical solution of the Navier-Stokes equations and subsequent flow verification have been detailed in the previous study (Sharma and Phares 2006), so only the main points are summarized here. The straight square duct, shown in Figure 1, is bounded by no-slip walls in the cross-stream directions and the flow is assumed to be periodic in the streamwise direction. The characteristic scales are the duct hydraulic diameter, h, equal to the length of each side of the square duct, and the average friction velocity, u _τ ^a , which is known from the imposed static pressure gradient, Π. With this scaling, the dimensionless equations governing fluid flow may be written

FIG. 1 Schematic of computation domain and coordinate system.

The friction Reynolds number, Re_τ = hu _τ ^a /ν, for the simulations presented here is equal to 300 corresponding to a bulk Reynolds number, Re _b ∼ 4410. All quantities presented throughout this paper are made dimensionless using the characteristic quantities, h and u _τ ^a , unless denoted with the superscript “+”, which indicates wall units. The dimensions of the square duct are thus 4π × 1 × 1 (see Figure 1), or L _x ⁺ = 3770 and L _y ⁺ = L _z ⁺ = 300, in wall units. The discretized governing equations are solved on a uniform Cartesian grid (327 × 195 × 195) totaling 12.43 million computational grid points and the grid resolution is Δ x ⁺ = 11.67 in the streamwise direction, and Δ y ⁺ = Δ z ⁺ = 1.57 in the wall-normal directions. The grid is staggered, so that the first grid point away from the wall is located at Δ y ⁺/2 = Δ z ⁺/2 = 0.79 wall units. The Kolmogorov length scale for this flow based on the volume-averaged dissipation is estimated to be a little over 2 wall units (Gavrilakis 1992). The dimensionless integration time step was Δ t = 1.5 × 10^{− 4} or equivalently Δ t ⁺ = 4.5 × 10^{− 2}. The inertial and viscous terms in the momentum equations were approximated in the interior of the flow domain by fourth-order accurate central differences, and all other spatial derivatives in the Eulerian equations of motion are at least second-order accurate. Near the boundaries, all derivatives are second-order accurate with the exception of inertial terms for which a third-order accurate upwind-biased formula was used.

The low Reynolds number turbulent flow was simulated starting from perturbed laminar flow initial conditions. The perturbed flow field was allowed to evolve until the flow became fully developed and statistically stationary. This numerical procedure was verified through comparisons with previous DNS and LES results on similar geometries (see Sharma and Phares 2006 for details).

Particle Tracking

Although a number of forces may act on a small particle suspended within a turbulent flow (most notably the Stokes drag force, gravity, and the Saffman lift force), only the Stokes drag is considered here. This approach is similar to that presented by Brooke et al. (1994) for particle tracking in a circular pipe. Gravity is neglected in order to better resolve the role of the secondary flows on deposition. The shear-induced Saffman lift force (Saffman 1965) is also neglected because it only assumes non-trivial magnitudes in the viscous sublayer (McLaughlin 1989). Wang et al. (1997) reported that neglecting the lift force resulted in only a slight reduction in the deposition rate computed from DNS simulations.

The dimensionless force balance equation takes the form

where the dimensionless particle relaxation time, τ _p , is defined

and is the particle-fluid density ratio ( = ρ _p /ρ _f ) and d _p is the dimensionless particle diameter. The dimensionless particle relaxation time is equal to the particle Stokes number, based on the mean friction velocity and the duct hydraulic diameter. The dimensionless particle relaxation time is commonly expressed in wall units, where τ _p ⁺ = τ _p Re _τ. In the present study, τ _p ⁺ was varied between 0 and 50, where τ _p ⁺ = 0 corresponds to a passive tracer. The density ratio for the finite-inertia particles was chosen to be = 763, which corresponds to olive oil droplets in air at 25°C. For a chosen value of τ _p ⁺, the density ratio determines the particle diameter. In each case, the selected density ratios yield particle diameters that are smaller than the estimated Kolmogorov length scale of 2 wall units.

The governing equations of particle motion are integrated using a fourth-order accurate Runge-Kutta method and a small fixed time step equal to the Eulerian flow advancement time step. Eulerian fluid velocities at intermediate points in the computational domain are determined using a tricubic interpolation scheme (Rovelstad et al. 1994). Periodicity was assumed in the streamwise direction for particles as well, allowing particle tracking over larger displacements in the x-direction. A particle was considered to deposit when it reached a distance from the wall equal to or less than its radius determined by Equation (5).

Two types of initial spatial particle distributions were used for the computations, to which we will refer as case A and case B. Case A features the release of 64,000 particles randomly distributed throughout the flow domain. Case B features the release of 120,428 particles uniformly distributed throughout planes located at a specified height, y ₀ ⁺, above each of the four duct walls (thus totaling 481,712 particles per simulation). In each case, the initial particle velocity was set equal to the Eulerian fluid velocity at that point and the total simulation time was 675 time wall units. Details of all of the simulations are listed in Table 1 and Table 2 for case A and case B, respectively.

TABLE 1 Particle transport simulation parameters for Case A—particles released at random locations throughout the flow domain

Relaxation time τ p ^+	Density ratio	Particle diameter d p	Total number of particles N p
0	—	—	64000
5	763	0.34	64000
10	763	0.49	64000
15	763	0.59	64000
20	763	0.67	64000
30	763	0.84	64000
50	763	1.09	64000

TABLE 2 Particle transport simulation parameters for Case B—particles released in planes located at specified heights above the walls

Relaxation time τ p ^+	Density ratio	Particle diameter d p ^+	Total number of particles N p
0	—	—	481712
5	763	0.34	481712
15	763	0.59	481712
30	763	0.84	481712
Initial release height, y 0 ^+ = 3, 4, 10, 20, 30, 40, 50, 60, 75, 100 and 130.

RESULTS

Deposition Rate

In accordance with previous studies of particle deposition, we have computed the deposition velocity, k _d ⁺, for a release of particles randomly distributed throughout the duct (case A). For monodisperse particles, the dimensionless deposition velocity may be written

where N _d is the number of particles that deposited over a time interval Δ t _d , and N _p is the number particles initially suspended in the flow at the beginning of that time interval. Both N _p and N _d will decrease over time as the number of suspended particles is depleted. However, the ratio remains constant after an initial transient period that is sensitive to the release conditions. This is demonstrated in Figure 2, which displays the deposition fraction, N _d /N _p , computed in 100 time intervals (Δ t _d ⁺ = 6.75) over the entire simulation duration of t ⁺ = 675. For each value of particle inertia, the deposition fraction became constant with time within 75 time wall units of the simulation. The deposition velocities, averaged over the last 600 time wall units of the simulation, are plotted in Figure 3 for τ _p ⁺ = 1, 5, 10, 15, 20, 30, and 50. Figure 3 compares the computed values with previously obtained experimental results for deposition in circular pipes compiled by Young and Leeming (1997). The dashed line denotes the data envelope within which the reported deposition data lie. We note that the DNS predictions of Marchioli et al. (2003) for deposition in a circular pipe also lie within the data envelope in the diffusion-impaction regime.

FIG. 2 Computed deposition fraction plotted versus time for case A simulation.

FIG. 3 Computed deposition velocities from case A versus dimensionless particle relaxation time. The dashed line represents the envelope around the circular pipe data compiled by Young and Leeming (1997).

The most notable feature of Figure 3 is that there does not appear to be a significant difference between deposition in a circular pipe and a square duct for the values of particle inertia considered here, all of which lie in the diffusion-impact and inertia-moderated regimes. The data point corresponding to the smallest value of particle inertia (τ _p ⁺ = 1) is slightly below the data envelope, suggesting that deposition is less likely in a square duct for low-inertia particles than in a circular duct. This is consistent with a previous observation that low-inertia particles are more likely than high-inertia particles to be carried back into the core of the square duct by the mean secondary swirling flow (Sharma and Phares 2006). By contrast, higher inertia particles are efficiently carried to the walls where they are more likely to deposit or accumulate near the wall in regions of high strain rate.

Deposition Patterns

The problem with using deposition velocity to describe particle transport to the walls of ducts having a non-circular cross-section is that local effects are ignored. Whereas deposition along the inner wall of a circular pipe is relatively uniform, the non-uniform deposition patterns along the walls of the square duct suggest deposition mechanisms that differ from a circular pipe, despite the fact that the overall penetration appears to be similar.

Figure 4 displays the deposition patterns computed along the walls of the duct for five values of particle inertia from the case A simulations. The data in Figure 4 are presented as the number of particles deposited within small bins along the wall (Δ y or Δ z), averaged along the total streamwise length of the duct and among all four walls. Thus the abscissa represents the distance along any wall in the direction normal to the streamwise direction, and y,z = 0 or 1 correspond to corners. The plot reveals that there exist definite deposition patterns, and that these patterns are dependent upon particle inertia. For the lower-inertia particles (τ _p ⁺ ≤ 15), deposition is more efficient closer to the center of the walls, whereas for higher-inertia particles (τ _p ⁺ ≥ 30), deposition is more efficient towards the corners. This suggests a preferential enhancement of deposition of higher-inertia particles by the secondary flows unique to the square duct, or a suppression of deposition of the lower-inertia particles by the secondary flows, or perhaps both. Deposition suppression of particles having τ _p ⁺ ≤ 15 is not immediately obvious from Figure 3, but is consistent with the observation by Sharma and Phares (2006) that lower-inertia particles are increasingly more likely to be carried back to the center of the duct by the mean swirling flows. Sharma and Phares (2006) noted that centrifuging by the mean swirl is not a valid explanation for this phenomenon, since the Stokes number based on the mean off-axis velocity and the size of the mean circulation is much smaller than unity for all values of particle inertia considered. One explanation for the observations is a combination of the secondary flow and turbophoretic force, which acts in the direction of decreasing turbulence intensity. In a circular pipe, the turbophoretic force peaks at a height above the wall of y ⁺ = 50 and is the driving force for deposition of particles in the diffusion-impact and inertia-moderated regimes (Young and Leeming 1997). In a square duct, one might imagine that, as particles are carried away from the duct center, the turbophoretic force preferentially increases for higher-inertia particles.

FIG. 4 Particle deposition patterns from case A averaged along the streamwise length of the duct and the four walls.

Effect of Release Height

In order to more definitively determine the relative effects of the secondary flows and turbophoresis on deposition, we employ the results from the case B simulations, which involve a release of particles at specified heights above the duct walls. Figure 5 displays the number of particles that deposited per wall during the entire simulation (675 time wall units) as a function of the initial release height and for three different values of particle inertia (τ _p ⁺ = 5, 15, and 30). The data corresponding to τ _p ⁺ = 15 and 30 exhibit a peak in the number deposited at an initial release height of y ₀ ⁺ = 40. This is not an artifact of the finite simulation time, because particles became distributed throughout the duct cross section well before the end of the simulation, as evidenced by the data corresponding to τ _p ⁺ = 5 which stays relatively flat with y ₀ ⁺. Note that the maximum deposition for the particles having τ _p ⁺ = 5 occurs at the lowest release height of y ₀ ⁺ = 3. This may be explained in terms of a free-flight mechanism. Initially, each particle is given a velocity equal to that of the fluid at its location. For free-flight to the wall to occur, a particle must be within its stopping distance to the wall, which for a particle having τ _p ⁺ = 5 and an initial velocity equal to 1 (normalized with u _τ) is y _stop ⁺ = 5.

FIG. 5 Number of particles deposited per wall versus initial release height.

Certainly the free-flight mechanism cannot solely explain the maxima observed at y ₀ ⁺ = 40 for the higher inertia particles. Rather, the spatial distribution initially centered at y ₀ ⁺ = 40 gradually broadens, while at the same time becoming skewed towards wall due to turbophoresis. This behavior is evident from Figure 6, which displays the concentration, C _N , normalized with the average concentration in the duct, as a function of height above the wall for particles initially released at y ₀ ⁺ = 50 at four different snapshots during the simulation between t ⁺ = 45 and t ⁺ = 135. The top plot (Figure 6a) corresponds to passive tracers, and the bottom (Figure 6b) corresponds to finite inertia particles (τ _p ⁺ = 30). For passive tracers, the maximum in the distribution remains close to y ⁺ = 50, whereas for the finite inertia particles, the maximum gradually approaches the wall between t ⁺ = 45 and t ⁺ = 135.

FIG. 6 Plots of instantaneous normalized particle concentration at five times during the simulation for a) passive tracers and b) finite-inertia particles. In both cases, the initial release height is y _o ⁺ = 50, and the concentration, C _N , is normalized with the average concentration throughout the duct cross section.

In order to better visualize this effect for different initial release heights, we have tracked the maximum in the spatial distribution as a function time. These data are presented in Figure 7a for passive tracers, and in Figure 7b for particles having τ _p ⁺ = 30. Each plotted data set corresponds to a different initial release height, y ₀ ⁺, which may be determined by the value at t ⁺ = 0. For passive tracers, the location of the maximum in the distribution stays relatively constant, until it begins to change chaotically for each subsequent time step, indicating that the particles are uniformly distributed throughout the duct cross section. As expected, the particles become uniformly distributed sooner for particles released farther from the duct walls. By contrast, the maximum in the spatial distribution of particles having τ _p ⁺ = 30 gradually approaches the wall between t ⁺ = 45 and t ⁺ = 135 for y ₀ ⁺ = 20 and y ₀ ⁺ = 50. This behavior occurs for release heights up to y ₀ ⁺ = 60, beyond which particles become more uniformly distributed until preferential concentration of the particles in the sublayer quickly shifts the location of the maximum to the near-wall region.

FIG. 7 Location of the maximum concentration as a function of time for y _o ⁺ = 20, 50, 75, and 100 for a) τ _p ⁺ = 0, b) τ _p ⁺ = 30.

The efficient transport of finite-inertia particles towards the duct wall for initial release heights lower than y ₀ ⁺ = 60 may be attributed to either turbophoresis or transport by the lateral mean flow. In order to distinguish between these two modes of transport, we note that, for the values of particle inertia considered here, turbophoresis is inertia-dependent, whereas transport by the mean flow is not. Figure 8, which tracks the location of the maximum as a function of time for three values of particle inertia, suggests that the particles migrate towards the wall embedded within the mean flow. However, unlike passive tracers, which are carried to the wall and then back to the core, finite-inertia particles tend to deposit or accumulate near the wall at an inertia-dependent rate. This inertia-dependent effect is evident from Figure 8 beyond roughly t ⁺ = 100.

FIG. 8 Location of the maximum concentration as a function of time for τ _p ⁺ = 5, 15, and 30. In each case, the initial release height was y _o ⁺ = 50.

One would expect the effect of the secondary flows to be most prominent for the deposition of particles some distance from the center of the duct wall. For example, the deposition patterns shown in Figure 4 exhibit the off-center deposition peaks located a distance of 0.1 from the corners (i.e., at y,z = 0.1 or 0.9). In light of the discussion in this section, one would expect maximum off-center deposition for particles that are initially released near y ₀ ⁺ = 40, and a smaller effect of initial release height on deposition close to the center of the wall. This behavior is evident from the deposition patterns plotted for varying initial release heights displayed in Figure 9 for τ _p ⁺ = 30 (Figure 9a), τ _p ⁺ = 15 (Figure 9b), and τ _p ⁺ = 5 (Figure 9c). That figure demonstrates that initial release height only affects the deposition rate for higher inertia particles (τ _p ⁺ ≥ 15), and only close to the corners. The deposition rate of lower inertia particles is always greater towards the center of the wall, and appears relatively insensitive to y ₀ ⁺ throughout.

FIG. 9 Particle deposition patterns for a) τ _p ⁺ = 30, b) τ _p ⁺ = 15, and c) τ _p ⁺ = 5. The data are plotted as the number of particles deposited within small bins along the wall (Δ y or Δ z), averaged along the total streamwise length of the duct and among all four walls.

DISCUSSION

From the computations presented here, one may conclude the following about the mechanism governing the deposition of particles in a square duct. The secondary flows carry particles from the duct core to the corners. Higher-inertia particles (τ _p ⁺ > 10) deposit efficiently close to the corners, gaining momentum towards the wall from the turbulent fluctuations in the outer region of the boundary layer. Lower-inertia particles have smaller stop distances, and thus tend not to deposit close to the corners. Rather, they are carried by the mean secondary flow towards the center of the wall, and then back into the duct core. The mean secondary flow velocity has been previously shown (Gavrilakis 1992) to be higher along the walls than along the diagonals of the duct. Thus particles that do not deposit immediately in the corners may be efficiently carried away from the corners and towards the wall midpoints. Since there appears to be a critical value of roughly τ _p ⁺ = 10, above which corner deposition is more efficient, it is interesting to note that this corresponds to particles that have a stop distance equal to the size of the viscous sublayer, assuming an initial velocity equal to u _τ.

The efficient corner deposition of particles having τ _p ⁺ = 15 and 30 from an initial height of y ₀ ⁺ = 40 and the corresponding insensitivity to initial height of the deposition along the center of the wall suggest unique transient flow structures in the corner that are conducive to deposition. These are likely related to the high mean streamwise vorticity gradients present in the corners (Gavrilakis 1992).

The extension of these deposition mechanisms to higher Reynolds number flows, such as those encountered in ventilation ducts, requires that the relative contributions of turbophoresis and the mean secondary flows to particle transport remain similar. One might imagine that increasing the Reynolds number should not significantly affect the deposition rates, since the turbulence variations close to the wall would not be significantly altered (Young and Leeming 1997). However, additional work is needed to definitively address the effect of Reynolds number on the secondary flows, turbophoresis, and thus deposition.

CONCLUSIONS

Direct numerical simulation of a turbulent flow through a square duct coupled with Lagrangian particle tracking was used to investigate the mechanisms by which particles deposit to the duct walls. The non-uniform deposition patterns along the walls of the duct indicate that the secondary flows play an important role in the deposition process. Our computations suggest that they tend to transport particles from the core of the duct towards the corners of the duct cross section and that they tend to preferentially recirculate lower inertia particles back into the duct core. Consequently, the secondary flows cause an increase in deposition closer to the corners for the higher-inertia particles and a decrease for the lower-inertia particles. This results in a greater sensitivity of the deposition velocity to particle inertia, compared to a circular pipe. The deposition rate along the wall close to the centerline appears less affected by the presence of the secondary flows.

Acknowledgments

The DNS computations were performed at the Texas A&M Supercomputing Facility and on the NSF Terascale Computing System at the Pittsburgh Supercomputing Center under NSF-NPACI grant no. CTS030005P. Partial financial support from the Comprehensive Test Ban Treaty Organization (CTBTO) through the Lovelace Respiratory Research Institute under contract no. JK020160 is also gratefully acknowledged.