Abstract

Debris clouds provide an important visual signature of tornadoes and can potentially significantly affect the wind structure, damage potential, and Doppler radar measurements of tornado wind speeds. To study such issues, the dynamics of finescale debris have been added to an existing high-resolution large-eddy simulation model of tornado dynamics. A so-called “two-fluid” or “Eulerian–Eulerian” approach is employed, together with a surface layer model for lofting and depositing debris. In this paper the debris implementation is described, three critical dimensionless parameters governing tornado debris effects are identified, and sample results from a large set of simulations of tornadoes with idealized debris are presented. The results demonstrate that the accumulation of small-scale debris within the surface layer and corner flow can significantly alter the wind speeds and flow structure of the tornado vortex within a few hundred meters of the surface. They suggest that the total mass of the debris cloud can reach tens of thousands of tons. Near the surface, the debris mass loading can be well above 1, the peak mean velocities can be reduced by as much as half, and the total momentum (air plus debris) can either significantly increase or decrease. Local air and debris velocities can differ significantly and in a nontrivial fashion, thereby complicating the interpretation of Doppler radar measurements of tornado structure. Debris fluctuations, centrifuging, negative buoyancy, and angular momentum transport are all significant mechanisms for the debris effects. A negative physical feedback reduces the sensitivity of the results to changes in the parameterization of the surface debris fluxes. The realistic simulation of tornado debris clouds and surface damage tracks should prove useful in identifying the dynamics governing their observed counterparts.

1. Introduction

There are many reasons for studying tornado debris clouds. They provide visual signatures of tornadoes, giving clues to their dynamics. Debris loading can contribute significantly to tornado damage. Understanding the differences between the local velocities of debris and air is important for interpreting Doppler radar measurements of tornado wind speeds. Contaminants could be spread long distances if lofted high into the parent storm. Finally, the presence of debris may alter the fluid-dynamic structure of the tornado itself.

To address such issues, debris dynamics have been added to a high-resolution large-eddy simulation (LES) model of tornado dynamics. In this paper, results are presented with particular emphasis on exploring the effects of debris loading on tornado fluid dynamics. The debris cloud occupies only a small part of the tornado and the volume fraction of debris within the cloud is likely very small, suggesting little impact. Previous LES studies of tornadoes, however, have shown that the dynamics in the corner flow region, where the central vortex meets the surface and the largest wind speeds are expected, is highly sensitive to the properties of the near-surface inflow (Lewellen et al. 2000; Lewellen and Lewellen 2007a). The largest debris loading might naturally be expected to occur in this inflow, raising the possibility of a significant affect on the dynamics. Given the large density ratio between, for example, soil and air, large debris mass loading is possible even for small volume fractions. Such loading would result in large changes in the effective total fluid density in the near-surface and corner flow regions and would also provide additional mechanisms for angular momentum transport through the falling or outward centrifuging of debris.

To test these possibilities, this work focuses on the corner flow region and cases likely leading to large mass loading (i.e., freely available finescale debris, such as soil or sand). Exploring the effects of debris on the fluid dynamics requires a fully interacting model. For this, a “two-fluid” approach was chosen rather than a particle tracking method (as has become common practice in engineering simulations of particulate-laden flows). The model is described in section 2, including numerical issues that had to be overcome. Some details are given in the appendix. The already multidimensional parameter space governing tornado corner flow behavior increases significantly with the addition of debris. The simplifications assumed for the present study, the simulation set considered, and a discussion of the most critical dimensionless parameters are given in section 3. Results from the simulations are presented in section 4, followed by further discussion and conclusions.

Some of this work was originally presented in Lewellen et al. (2004) and Gong et al. (2006). Debris in tornadoes and related flows has also been considered from other perspectives in previous work. The structure of dust devils has been studied both observationally and numerically, partly to better understand dust devils on Mars (e.g., Toigo et al. 2003; Kanak 2005). The best observations are probably the Doppler radar measurements of Bluestein et al. (2004); they note that because the exact size and density of the particles were unknown, it was not possible to determine how much centrifuging might have affected their results. Gu et al. (2006) used LES to model the airflow in a dust devil and used Lagrangian tracking of dust particles to determine how well different size particles follow the airflow. They show that at low levels significant centrifuging of the dust should be expected near the radius of maximum tangential velocity. Dowell et al. (2005) estimated the velocity difference in a strong tornado between air and a variety of debris. They used 1D and 2D models with an initially uniform particle distribution and no feedback on the airflow. They found that particles generally move more slowly than the air and that the velocity differences depend strongly on the particle characteristics and can be quite large, particularly near the surface. Magsig and Snow (1998) tracked debris from tornadic thunderstorms to distances beyond 150 km. They used paper checks as particularly useful target debris, for which they report fall velocities of 0.5–1.5 m s−1. The smallest sand particles considered in this study have fall velocities of this order, but we model only the central tornado flow, whereas much of their reported dispersal is due to winds in the parent thunderstorm.

2. Tornado LES model with finescale debris

Aside from the addition of debris, the LES model and simulation procedures used for the present study are as described in Lewellen and Lewellen (2007a) and references therein. It is a fully 3D unsteady implementation of the incompressible Navier–Stokes equations on staggered grids, with the grid spacing stretched independently in the three spatial directions. A leapfrog scheme in time and second-order central differences in space are employed with a variable time step. The Poisson equation for pressure is solved directly on the nonuniform grid using the method of Farnell (1980) in the horizontal directions and a tridiagonal matrix inversion in the vertical. The subgrid model includes a prognostic equation for subgrid turbulence kinetic energy (q2/2) and a subgrid turbulence length scale limited by grid spacing, distance to the surface, and stability from rotational damping or temperature stratification. The surface boundary condition is no-slip with a specified aerodynamic roughness length z0. The lateral and top boundary conditions are imposed on the velocities: fixed values except for tangential components on outflow, which are taken to be zero-slope. Generally, at a given height the lateral inflow on the square boundary is specified so that it has constant values of angular momentum and horizontal convergence.

a. Debris implementation

In engineering applications, two principal methods have been considered to include particles into flow simulations: Lagrangian tracking of individual particles or groups of particles, and treating the particles in an Eulerian framework as one or more additional species of “fluid” (see, e.g., Crowe et al. 1996; Riber et al. 2006). A modest number of Lagrangian trajectories often suffices to determine the basic transport of particles by a flow. When the effects of the particles on that flow become important, a much greater number and full two-way coupling is required to represent the interaction. At least one Lagrangian parcel (representing some fixed number of particles) is required at any given time in regions where the interaction is nonnegligible even if the particle loading is small; proportionally more are required in regions of large loading. The number of two-way coupled parcels required is at least of the order of a few times larger than the total number of grid cells, which can be numerically prohibitive in a high-resolution simulation. For small enough particles, their velocity relative to the air can be diagnosed at each time step and an “equilibrium Eulerian” approach can be employed (e.g., Ferry et al. 2005); larger particles or large fluid accelerations require independent evolution equations for the particle velocity. The so-called two-fluid or “Eulerian–Eulerian” approach has been chosen here for studying the effects of the debris on the flow structure in different regimes. A few (∼200) Lagrangian trajectories have sometimes been included to provide diagnostic checks on the debris-fluid representation. Eventually the trajectory capability will be used to incorporate large debris elements (e.g., cars or farm animals). The incorporation of a single debris “fluid” is presented here; extension to multiple debris fluids is straightforward and already implemented in the current LES.

The basic equations of the two-fluid (sometimes called “dusty gas”) model are (Marble 1970):

 
formula

These respectively represent air momentum, debris momentum, air continuity, and debris continuity. The air is treated as incompressible with constant density ρ and velocity components ui and the debris fluid (variable density ρd and velocity components udi) is treated as pressureless (i.e., particle collisions are ignored). The latter is a good approximation as long as the debris volume fraction is much less than one. Given the large ratio of the density of the debris itself (σ) to that of air for debris such as soil, this condition is satisfied even for mass loading (Dρd/ρ) well above one. Near a surface, particle collisions can significantly affect particle transport even at relatively low volume fractions (e.g., Yamamoto et al. 2001); in our model, these effects are lumped into the surface treatment discussed later. The debris falls because of gravity, and momentum is exchanged between the two fluids through the drag force Fd, which depends on the local velocity difference between the two fluids, and the particle Reynolds number Repρd|udu|/μ; that is,

 
formula
 
formula
 
formula

Here, μ is the viscosity of air and d is the particle diameter. The form of the drag coefficient CD is appropriate for spherical particles with 0 < Rep < 2 × 105 (White 1991). Lift forces, which can arise from nonspherical shapes or particle rotation, are not considered. The debris to air density ratio is assumed to be large (σ/ρ ≫ 1) so buoyant effects of displaced air are also negligible.

Outside of the very-near-surface layer (discussed below), the importance of the subgrid turbulence (encoded in τdij and Kd) is at most secondary in the tornado flow as simulated here. Large accelerations driving differences between the air and debris motions appear already here in the mean swirling flow; furthermore, a sufficiently fine grid is employed to resolve the most important scales of turbulence in the corner flow (Lewellen et al. 1997), leaving the subgrid model with much less responsibility. Accordingly, for simplicity, the subgrid contribution to the debris evolution is treated here analogously to the treatment for the fluid, assuming down-gradient diffusion of the momentum components and using the same subgrid-scale turbulent eddy viscosity for both debris and air. The subgrid diffusion of debris loading is also assumed to be downgradient with the same turbulent diffusion coefficient as used in diffusing the subgrid turbulent kinetic energy. The effect of particles damping small-scale turbulence (e.g., Druzhinin and Elghobashi 1999) is ignored in the subgrid here; its relative importance among particle effects is greater in resolved flows not characterized by large accelerations (e.g., channel flow).

b. Numerical implementation

For simplicity, the equations governing the debris fluid have been implemented numerically in the same way as the primary fluid in the existing LES model, to the extent possible. The new 3D field variables are the debris mass loading D and the three components of debris momentum 𝒰, 𝒱, and 𝒲, with 𝒰 ≡ Dud, etc. The equations were discretized so that these variables are exactly conserved at the finite difference level. Numerical stability issues arise that are not present in the discretization of the Navier–Stokes equations alone, primarily from three sources: the often dramatic differences between debris velocity and debris momentum fields (due to the strongly varying D field), the absence of a pressure force for debris, and the inclusion of the drag force. The most important of these issues are discussed in the appendix. Limits were imposed on the debris fluxes (to avoid negative D), on the debris velocities for very small D, on the time step (to satisfy Courant conditions for both air and debris velocities), and on τυ relative to the time step. To a large degree, the drag term in the debris momentum equation plays a roll analogous to the pressure gradient term for the air, but with different numerical stability issues involved.

c. Surface fluxes

To complete the debris implementation, fluxes of debris and debris momentum at the surface must be specified. These will depend on the air and debris flow just above the surface as well as on surface properties. Because we are interested in studying large debris loading, we concentrate on surfaces with loosely bound smaller particles (e.g., sand). The motion of particles in a surface layer with friction velocity u* is largely characterized by Shield’s parameter ρu2*/(σdg), essentially a dimensionless ratio of the aerodynamic force on a particle to its weight. For values between ∼10−2 and ∼1, particles bounce along the surface in a process known as saltation (Owen 1964), which is responsible, for example, for sand dune transport. More relevant here is the less well understood range above this, where the particles enter into suspension in the turbulent surface layer. Although there is considerable uncertainty, for straight-line winds the sand flux in this regime appears to increase as a modest power of u* (∼ u* or ∼ u3/2*; see, e.g., Batt et al. 1999). Accordingly, we model the vertical flux of D at the surface as a sum of three terms:

 
formula

where the subscript 1 indicates values at the first grid level above the surface and the flux component is zero if the respective inequality condition is not satisfied. Here, 𝒲a0 represents debris lofting by turbulent airflow; q has been chosen in place of u* to capture more of the behavior of the complex tornado surface layer (for a straight-line wind the subgrid model predicts u* ≈ 0.4q); 𝒲d0, of analogous form, represents the possibility of the debris fluid knocking loose additional debris from the surface. The debris friction velocity u*d is modeled analogously to its air counterpart in terms of the debris velocity, von Kármán’s constant k, and the surface roughness length z0, u*d = k|ud(z1)|/ln(z1/z0). The last term, 𝒲0, represents falling debris impinging on the surface with some fraction (0 < cb < 1) sticking and the remainder bouncing. A vertical pressure gradient induced source term has not as yet been included, although it might play a significant role for some tornadoes passing over some surfaces (e.g., in lifting strips of road pavement).

Corresponding terms are included for the modeled vertical flux of vertical debris momentum,

 
formula

The first terms are the respective debris fluxes multiplied by a corresponding velocity scale. The last is the vertical momentum transfer to the surface from falling particles, with some fraction allowed to bounce. The vertical flux of horizontal debris momentum is modeled with two terms, the first in direct analogy with the frictional loss of air momentum to the surface and the second due to the flux of falling particles,

 
formula

and similarly for 0. No term is included for the horizontal momentum of newly ejected particles at the surface because the drag force will transfer this momentum in the layer just above.

The modeling in (5)(7) is very simplified. In principle, the coefficients should depend on properties of the surface and particles (e.g., shape, moisture, adhesion, etc.). For the present work no attempt is made to correlate the boundary conditions with any specific natural surface, no limit is placed on the amount of debris available, and no distinction is made between the uptake of new particles and ones that may have been lofted previously and redeposited.

d. Limitations and testing

In defining a physical problem, modeling it with equations, and numerically implementing them, approximations are involved at each step. The finite difference implementation of (1)(7) was checked in sample test cases (e.g., advecting a cube-shaped debris cloud in a fixed wind) and with different temporal and spatial resolution to insure conservation of debris mass and momentum, numerical stability, and accuracy. Assessing the modeled equations themselves is more difficult. The chief physical limitation of the two-fluid model is that all the debris within a grid cell is assumed to possess the same velocity; in reality there will be some spectrum. The approximation should work well except where debris and air velocities differ greatly, and even then perform reasonably if the debris within the grid cell all had similar recent histories. The model performs poorly in regions where two jets of debris collide, for example. This limitation can be lessened, at significant numerical cost, by representing the debris via multiple fluids to allow different velocities for different populations within a grid cell. Each debris fluid contributes to the computational cost approximately as much as the fluid mechanics without debris. For test purposes this technique was used for tornado simulations with two identical debris species, with the total surface source distributed between them in different ways. The analogous results from the two-fluid model were found to be in good agreement with these “three-fluid model” tests, supporting the applicability of the two-fluid approximation for tornado-type flows. Comparisons were also made between two-fluid model results and the local particle velocities found along sample Lagrangian particle trajectories, again with good statistical agreement between the two (Gong 2006).

There are a wide range of surface conditions and debris that real tornadoes can encounter. For a first assessment of debris effects on tornadoes it is enough to fall somewhere within the realistic range. Comparison of the model results with the Batt et al. (1999) wind tunnel data of the erosion of sand beds under high (30–100 m s−1) straight-line winds showed at least rough agreement for this geometry (Gong 2006); there are no systematic measurements in high-speed swirling environments to compare with that we are aware of. Fortunately, as discussed in section 4d below, there is a negative feedback mechanism that greatly reduces the sensitivity to details of the surface parameterization.

3. Parameter space and simulation set

Many physical parameters affect tornado behavior and the corner flow in particular. Significant simplification is required to allow any systematic study. In the absence of debris, Lewellen et al. (2000) argued that a minimum of three parameters was required to capture the basic corner flow dynamics. For this they chose a “far field” angular momentum level, Γ, outside of the surface–corner–core flow; a characteristic core radius well above the corner flow, rc; and a depleted angular momentum flux, ϒ, flowing through the surface–corner– core regions. Choosing Γ and rc to nondimensionalize the results, one can then classify the basic behavior in terms of a single dimensionless parameter, the corner-flow swirl ratio ScrcΓ2/ϒ. The leading effects of physical parameters such as horizontal convergence, core conditions aloft, or surface roughness were argued to manifest themselves in large part by the changes they forced in rc or ϒ. Air density is assumed constant and, in the no-debris case, buoyant forcing is ignored so g and ρ drop out of the equation set; the flow is at a very high Reynolds number so μ drops out as well. As in Lewellen et al. (2000), we choose a definition for rc in terms of the peak mean upper-core swirl velocity (rc ≡ Γ/Vc); either rc or Vc may be chosen as the more basic physical parameter.

When debris is included, ρ, μ, and g no longer drop out and debris and surface properties become important as well. With the idealizations employed (spherical particles of a single size and density, uniform debris availability on the surface, etc.), the addition of debris involves five new primary physical parameters: ρ, μ, g, σ, and d. To a good approximation this set may be reduced further. For Rep exceeding a few hundred, CD is slowly varying and asymptotically constant. In this limit, μ again drops out of the equation set and σ and d appear only in the combination σd. To nondimensionalize the results we choose length, time, and mass scales rc, r2c, and ρr3c. The remaining parameter space can then be described in terms of three dimensionless parameters. For ease of physical interpretation these are chosen as

 
formula
 
formula
 
formula

As described in Lewellen et al. (2000), Sc provides a measure of tornado corner flow “type” without debris. The quantity Aa is a ratio of a characteristic radial acceleration in the flow to the gravitational one, which provides some measure of the relative importance of debris centrifuging, and Aυ is a ratio of a characteristic flow velocity to the terminal velocity of a particle in free fall, wt. This can be interpreted as a measure of the relative ease of debris lofting in the flow. For convenience in applying (8)(10) as governing parameters for the tornado debris study below, they are evaluated for the quasi-steady tornado before the introduction of debris.

In this work, different debris is introduced into a single sample tornado; Sc and Aa for the no-debris case are held fixed and Aυ varied. The no-debris case chosen as a starting point is an intense quasi-steady medium swirl (Sc ≈ 3.3) tornado. A preliminary examination of the effects of varying tornado type and intensity (and thereby Sc and Aa) as well as test cases varying μ and g are considered in Gong et al. (2006) and Gong (2006) and will be treated more extensively elsewhere.

The principal simulations considered in this work are listed in Table 1. The run names give some indication of the primary variation considered in each, including: no debris (SN); debris diameter (S2, S1, S0.5, S0.2); density ratio (SR4, SR8); debris pick-up threshold ct (STh.2, STh1); other surface parameters (SS1, SS2, SS3, SS4); doubled horizontal grid resolution (S0.5C, S1C); and translation (T1). All were computed in (2 km)3 domains with constant horizontal convergence (ac = 0.0125 s−1) and angular momentum (Γ = 104 m2 s−1) specified on the lateral boundaries, except in the bottom 10 m where Γ dropped linearly to zero. A constant vertical velocity was specified at domain top and a surface roughness length z0 = 0.02 m was employed. Except where noted, all simulations had finest resolutions of 1 m vertical and 4 m horizontal. Because the default choice σ/ρ = 2000 is within 10% of the density ratio of quartz to air at the surface, we refer to this debris as “sand.”

Table 1.

Conditions for principal simulations along with some selected summary output as discussed in the text. Entries include debris particle diameter d; particle terminal velocity in free fall wt; lofting parameter Aυ from Eq. (10); time scale for debris cloud formation tD; total debris cloud mass ∫D; peak mean swirl velocity overall (Vmax) and in the upper core (Vc); debris cloud interior cone angle αD, total height Ht, and center-of-mass height HD and radius RD. Except as noted, all simulations had debris density ratio σ/ρ = 2000, surface parameters ca = 0.12 (debris flux coefficient), cd = 0 (no added flux from debris collisions), ct = 0.5 (debris pick-up threshold), cu = cb = 1.0 (no debris bounce from surface), and were nontranslating.

Conditions for principal simulations along with some selected summary output as discussed in the text. Entries include debris particle diameter d; particle terminal velocity in free fall wt; lofting parameter Aυ from Eq. (10); time scale for debris cloud formation tD; total debris cloud mass ∫D; peak mean swirl velocity overall (Vmax) and in the upper core (Vc); debris cloud interior cone angle αD, total height Ht, and center-of-mass height HD and radius RD. Except as noted, all simulations had debris density ratio σ/ρ = 2000, surface parameters ca = 0.12 (debris flux coefficient), cd = 0 (no added flux from debris collisions), ct = 0.5 (debris pick-up threshold), cu = cb = 1.0 (no debris bounce from surface), and were nontranslating.
Conditions for principal simulations along with some selected summary output as discussed in the text. Entries include debris particle diameter d; particle terminal velocity in free fall wt; lofting parameter Aυ from Eq. (10); time scale for debris cloud formation tD; total debris cloud mass ∫D; peak mean swirl velocity overall (Vmax) and in the upper core (Vc); debris cloud interior cone angle αD, total height Ht, and center-of-mass height HD and radius RD. Except as noted, all simulations had debris density ratio σ/ρ = 2000, surface parameters ca = 0.12 (debris flux coefficient), cd = 0 (no added flux from debris collisions), ct = 0.5 (debris pick-up threshold), cu = cb = 1.0 (no debris bounce from surface), and were nontranslating.

4. Results

a. Debris cloud development and general appearance

In simulation T1, a quasi-steady tornado without debris translates at 15 m s−1 into a region with loose 1-mm-diameter “sand” on the surface. Figure 1 shows the resulting sand cloud development at different times as a check of verisimilitude. The funnel cloud is generated assuming a constant total water mixing ratio (14 g kg−1) and liquid water potential temperature (303 K) and so represents essentially a constant pressure surface. The sand cloud development, especially animated, looks much like a real debris cloud. As the simulated tornado approaches the debris field, the sand is kicked up in a thin layer and drawn radially inward, spiraling around the vortex center. Within the corner flow the sand is lofted, forming over time a turbulent cloud that grows up and outward. As sand is tossed out of the regions of strongest updraft it falls downward, some to be recirculated into the corner flow by the strong radial inflow at low levels. The sand cloud itself eventually reaches a quasi-steady state in which the pickup and deposition balance on average. Figure 2 shows a simulated surface track for this case. Reminiscent of real damage tracks, there is a complex pattern of removal and deposition of the sand. Further, there is a clear weakening of the peak surface pressure drop as the sand cloud forms, an indication of a strong effect of the sand on the flow, which is considered at length below.

Fig. 1.

Time series of 3D debris cloud for a simulated tornado translating at 15 m s−1 from right to left and entering a field of 1 mm sand. Simulation times shown are at (left) 20 s (when the center of the simulated tornado reaches the sand boundary), (middle) 40 s, and (right) 5 min.

Fig. 1.

Time series of 3D debris cloud for a simulated tornado translating at 15 m s−1 from right to left and entering a field of 1 mm sand. Simulation times shown are at (left) 20 s (when the center of the simulated tornado reaches the sand boundary), (middle) 40 s, and (right) 5 min.

Fig. 2.

Simulated 3.5-km-long surface track from the tornado of Fig. 1: (a) peak surface pressure drop encountered (contour interval = 20 hPa) and (b) net sand redistribution with contours (interval = 5 kg m−2, equivalent to ∼3-mm sand depth) showing removal and grayscale deposition (darkest shade indicating 20 kg m−2). The position of the sand boundary is clear on the left-hand side of (b).

Fig. 2.

Simulated 3.5-km-long surface track from the tornado of Fig. 1: (a) peak surface pressure drop encountered (contour interval = 20 hPa) and (b) net sand redistribution with contours (interval = 5 kg m−2, equivalent to ∼3-mm sand depth) showing removal and grayscale deposition (darkest shade indicating 20 kg m−2). The position of the sand boundary is clear on the left-hand side of (b).

To obtain better quantitative statistics in this highly turbulent flow, the rest of the simulations in Table 1 were conducted without translation. Once the particle cloud and tornado reach a quasi-steady state, the time average is then axisymmetric, permitting an azimuthal average to be taken as well. It is these azimuthal–time average fields that are used for the summary output of Table 1 and for the contour plots below. Three measures of particle cloud extent are given: the height HD and radius RD of the cloud center of mass and the height Ht below which 98% of the total debris mass lies. The interior cone angle αD is defined such that 98% of the total debris mass lies outside of it. For these measures and for the total particle cloud mass, the lowest 1 m above the surface is excluded from the integrals to distinguish lofted particles from those merely bouncing along the surface. Although the nontranslating simulations give more symmetric mean fields than their translating counterparts, the basic dynamics and magnitudes of the sand clouds in the two cases remain similar (cf. S1 and T1 in Table 1).

Rough estimates of uncertainty levels for the various quantities can be inferred from the table. For simulation S0.5, averages over two time windows are provided. The outputs are all in close agreement except for the sand cloud heights; for finer particles there tends to be a rise and fall of the upper portion of the cloud on time scales that are too long to be averaged out over the time windows used. Results are also given from two simulations with the horizontal grid spacing doubled everywhere. The fields from simulations S0.5C and S1C showed good qualitative agreement with their standard resolution counterparts, S0.5 and S1, and generally good quantitative agreement of the mean summary quantities as well, although with a reduction in the peak velocities. Several runs are included that differ only in the surface debris parameterization; these are considered further below.

The appearance of the quasi-steady particle cloud depends strongly on the particle type. Figure 3 provides an example. Not surprisingly, finer sand is lofted higher into a larger cloud and thrown radially outward at a smaller interior cone angle than larger-sized sand (i.e., it is centrifuged less). The reduction of the funnel cloud with increased sand loading in Fig. 3 provides another indication of debris lowering velocities in the corner flow.

Fig. 3.

The funnel cloud and cutaway view of the particle cloud for a simulated medium-swirl tornado with fixed boundary conditions picking up (left) 2 mm (S2) and (right) 0.5 mm (S0.5) “sand.”

Fig. 3.

The funnel cloud and cutaway view of the particle cloud for a simulated medium-swirl tornado with fixed boundary conditions picking up (left) 2 mm (S2) and (right) 0.5 mm (S0.5) “sand.”

Figure 4 shows the time development for the total mass of the debris cloud for several cases. The time for the total mass to reach an initial maximum, tD in Table 1, provides a time scale for the development of the cloud; it naturally grows with the size of the cloud itself. The inner corner flow is found to equilibrate more rapidly than the particle cloud as a whole. The total mass in the quasi-steady sand cloud can reach tens of thousands of tons. Reaching such states requires both sufficient loose sand on the surface and a sufficiently long duration to approach the quasi-steady state. For simulation S0.2, the quasi-steady state should be considered more of an idealization given the quantity of fine sand and the long time scale involved. Also, some particles were lofted out through the top of the domain at 2 km in this case, suggesting they might be transported to much greater distances in a real storm. The total flux through the top compared to the total mean upward flux within the corner flow was about 20% for S0.2 (0.2-mm sand), ∼2% for S0.5 (0.5 mm), and numerically insignificant for S1 and S2 (1 and 2 mm).

Fig. 4.

Time histories of total sand cloud mass for different medium-swirl, nontranslating tornadoes from Table 1 as labeled.

Fig. 4.

Time histories of total sand cloud mass for different medium-swirl, nontranslating tornadoes from Table 1 as labeled.

Figure 5 shows mean contours of D for simulations with different sizes of sand. The increase of the sand cloud height and radius and decrease in interior cone angle with decreasing sand diameter are all readily apparent. The mass loading is modest compared to unity over most of the sand cloud but significantly exceeds it within a swirling annulus above the surface within the corner flow. For a fixed surface parameterization, tD rises slightly faster than linearly with Aυ, and HD slightly slower (Table 1); RD rises approximately linearly with Aυ from a base value (∼170 m) that likely scales with rc. The equilibrium total particle cloud mass scales approximately with HDR2D (i.e., with cloud volume). The results of SR4 and SR8 fit in smoothly with these scalings, supporting the theory that particles in the same drag class as measured by wt give similar results even when size and density vary greatly.

Fig. 5.

Axisymmetric time averages of debris loading shown on the radial–vertical plane for quasi-steady simulations with (a) 2 mm (S2), (b) 0.5 mm (S0.5), and (c) 0.2 mm sand (S0.2). Contour intervals are 0.05 kg kg−1. Inset views are 3× blowups of the inner corner (r, z ≤ 50 m) with contour interval 1 kg kg−1.

Fig. 5.

Axisymmetric time averages of debris loading shown on the radial–vertical plane for quasi-steady simulations with (a) 2 mm (S2), (b) 0.5 mm (S0.5), and (c) 0.2 mm sand (S0.2). Contour intervals are 0.05 kg kg−1. Inset views are 3× blowups of the inner corner (r, z ≤ 50 m) with contour interval 1 kg kg−1.

b. Effects of different sized sand on tornado wind fields

Figure 6 and Table 1 show that peak mean swirl velocity drops as total particle loading increases in the simulated tornado by ∼20% in S2 (2-mm sand), by a third in S0.5 (0.5-mm sand), and is almost cut in half in S0.2 (0.2-mm sand). At the same time, the location of that peak swirl velocity is shifted progressively outward and upward. Radial and vertical air velocities in the corner flow are also reduced (cf. Figs. 7 and 8). Instantaneous 3D velocity fields for the different cases show the same trends, although the velocity reduction is not as great for S0.2 as the reduction of mean values in Table 1 implies: as will be discussed more below, the wander of the low-level central vortex, which is minimal in other cases, becomes more pronounced, so that ∼l/4 of the reduction in Vmax from SN to S0.2 is attributable to azimuthal averaging over this wander.

Fig. 6.

Axisymmetric time average contours of swirl velocity (m s−1) and momentum (kg m−2 s−1) on the radial–vertical plane for simulations as labeled. (top) Air velocity for simulations with different sand sizes: (a) no debris, (b) 2 mm, (c) 0.5 mm, and (d) 0.2 mm. (bottom) Total momentum and debris velocity pairs for (e), (f) 2 mm and (g), (h) 0.2 mm.

Fig. 6.

Axisymmetric time average contours of swirl velocity (m s−1) and momentum (kg m−2 s−1) on the radial–vertical plane for simulations as labeled. (top) Air velocity for simulations with different sand sizes: (a) no debris, (b) 2 mm, (c) 0.5 mm, and (d) 0.2 mm. (bottom) Total momentum and debris velocity pairs for (e), (f) 2 mm and (g), (h) 0.2 mm.

Fig. 7.

Axisymmetric time average contours of vertical velocity (m s−1) on the radial–vertical plane for different simulations as labeled: (a) no debris air velocity, (b) air velocity, (c) debris velocity, and (d) total momentum (kg m−2 s−1) for 0.5 mm sand; (e) air velocity and (f) debris velocity for 2-mm sand; and (g) air velocity and (h) debris velocity for 0.2-mm sand. The heavy lines in (c), (f), and (h) are contours of 1% mass loading to provide some indication of the time-averaged sand cloud boundary.

Fig. 7.

Axisymmetric time average contours of vertical velocity (m s−1) on the radial–vertical plane for different simulations as labeled: (a) no debris air velocity, (b) air velocity, (c) debris velocity, and (d) total momentum (kg m−2 s−1) for 0.5 mm sand; (e) air velocity and (f) debris velocity for 2-mm sand; and (g) air velocity and (h) debris velocity for 0.2-mm sand. The heavy lines in (c), (f), and (h) are contours of 1% mass loading to provide some indication of the time-averaged sand cloud boundary.

Fig. 8.

As in Fig. 7, but for radial velocity (m s−1) and without the heavy lines in (c), (f), and (h).

Fig. 8.

As in Fig. 7, but for radial velocity (m s−1) and without the heavy lines in (c), (f), and (h).

The largest effects of particles on the airflow are felt within the corner flow and particle cloud; however, these effects are communicated to the vortex above via core pressure gradients, reducing the peak mean upper core swirl velocities (Vc in Table 1) by more modest levels. In comparing simulations with and without debris to assess the debris effects, it should be kept in mind that the possible response of flow regions above the simulated domain have not themselves been simulated; the magnitude of the changes in the upper-core conditions produced by debris is thus affected by the choice of top boundary condition, with the degree of uncertainty so introduced increasing with the change in Vc that is found.

Although adding sand generally reduces the peak mean velocities near the surface, this does not necessarily mean that the tornado’s damage potential is reduced. Near the surface the total swirl momentum can be dominated by the sand component, as seen in Figs. 6e–h. The peak mean total swirl momentum is dramatically increased, by a factor of ∼2.5, for S2 relative to the no-debris case. For airflow alone, the kinetic energy, closely related to the flux of momentum, likely provides a better measure of damage potential. In simulation S2, the peak total kinetic energy, including that of the sand, is slightly reduced relative to that in SN. Nonetheless, given the impulsive impact of coarse sand (essentially sandblasting), this likely would result in greater damage to most structures.

c. Debris mechanisms at work

The debris affects the airflow through several mechanisms: local momentum transfer, secondary spin up, negative buoyancy, centrifuging, angular momentum transport, and other effects of debris inertia. The first is the transfer of momentum from the air to the debris. As debris is lofted and brought up to speed by the drag force the total momentum is locally conserved, so the air velocity reduces as D increases. A simple case to consider is an initially cylindrically symmetric purely swirling flow υi(r), with fine dust (initially at rest) added below some height zd and inside a radius rd with loading D. In the dust region, the drag force will quickly equilibrate the air and debris velocities to υ(r < rd, z < zd) = υi(r)/(1 + D), preserving the total momentum. As a result, however, the cyclostrophic pressure gradient balancing the inertia of the combined air and dust flow, ∂p/∂r = (1 + D)υ2/r, will differ above and below zd {being υ2i/r above and υ2i/[r(1 + D)] below}. This vertical pressure gradient drives a core updraft tending to spin up (down) the flow below (above) zd. If the dust region is small compared to that above, the adjustment will primarily be a spinup below zd to equilibrate with conditions above. A new equilibrium can be reached given υ(z < zd) ∼ υi/1 + D. In this simple case, introducing debris reduces the velocity and increases the total momentum, both by a factor on the order of ∼1 + D.

Mass loading can also change the velocity structure through negative buoyancy, producing downdrafts that were not present before (Fig. 7). The mean downdrafts are annular, surrounding the core updraft; the instantaneous downdraft is often comprised of loose spiral bands. The vertical extent is directly correlated with the depth of the particle cloud. The downdrafts produced are larger for larger Aυ because of the increased drag force and extent of the particle cloud.

Just as the particles fall due to the gravitational force, they are thrown outward in response to the centrifugal force. The effects are largest in the near surface inflow into the corner flow, effectively increasing Sc [which can be interpreted as the ratio of a swirl velocity to a flow-through velocity within the surface–corner–core flow (Lewellen et al. 2000)]. In the near surface inflow, the introduction of debris reduces the radial air velocity component more than it reduces the swirl component. Both are reduced by losing some momentum to the debris, but in addition the radial inflow is opposed by the centrifugal force, and the near-surface radial pressure gradient—which drives the radial inflow—is reduced by the debris as well. The effective increase in Sc is one way to understand the qualitative changes in corner flow structure with increasing sand loading (e.g., the outward shift of the swirl velocity peak in Fig. 6 and updraft annulus in Fig. 7).

The above effects can contribute even for fine dust that perfectly follows the airflow (the limit Aυ → ∞). For smaller Aυ, differences between the air and debris velocities contribute further modifications. The velocity difference increases with the particle wt and with the local acceleration. Just as a free-falling particle reaches a terminal velocity, a particle in a locally constant airflow acceleration will reach (over time scales large compared to τυ) an equilibrium “slip” velocity for which the drag force per mass balances the acceleration. This provides an estimate for differences between the particle and air velocity fields:

 
formula

Thus greater |udu| is seen in S2 than in S0.2 (e.g., Figs. 6b and 6f versus Figs. 6d and 6h, Figs. 7e and 7f versus Figs. 7g and 7h, and Figs. 8e and 8f versus Figs. 8g and 8h), and greater differences in the radial velocity component than in the swirl exist for each case. The debris velocity contours given in the figures are mass weighted averages, 〈ud〉 = 〈Dud〉/〈D〉. In S2, the mean peak radial outflow velocity for the sand is more than double that of the air, whereas the sand peak radial inflow velocity is almost halved relative to that of the air (Fig. 8). Differences between air and particle velocities must be kept in mind when interpreting Doppler radar measurements of tornadoes, particularly in the near-surface and corner flow regions. Unfortunately, the properties of the particle scatterers are generally not known with any certainty in the measurements. Moreover, in practice (11) applied to the time averaged velocity fields does not provide a very good estimate for the actual velocity difference, even for very small particles, because the contributions of fluctuations are significant and because the mass weighted averaging differs significantly for ud and u. An example of the effects of fluctuations can be seen in Fig. 7f. Because of strong centrifuging, the central core is only sporadically populated by sand. Only those fluctuations with large vertical velocities populate this region, leading to the interior annulus of mean vertical debris velocity maximum in Fig. 7f and, through the drag force, to the weaker local maximum in the vertical air velocity at the same location seen in Fig. 7e.

Differences between air and debris velocities provide an additional mechanism for angular momentum transport. Centrifuging transports angular momentum outward away from the core; debris fall transports it downward. This tends to lessen velocities by inhibiting the inward radial transport of angular momentum, but for some range of particles it can allow angular momentum to collect in a near-surface debris annulus (e.g., Fig. 6e). Particles falling from above can contribute angular momentum to the near-surface inflow (again effectively increasing Sc); ones falling to the surface provide an angular momentum sink in addition to surface friction.

Finally, the inertia of the debris and nature of the drag interaction can force other changes when the mass loading is high. For D ≪ 1 the debris responds to the airflow, but for large D this dynamic is reversed: the airflow largely responds to the debris. This tends to inhibit flow structures with large local accelerations in regions with large D. In simulations of low-swirl tornadoes, we have found that sufficient mass loading tends to eliminate the configuration of an end-wall jet (transforming it essentially to a medium corner-flow swirl ratio case); in simulations of high-swirl tornadoes, sand loading can inhibit or even eliminate multiple secondary vortices (Gong et al. 2006; Gong 2006). Furthermore, the response of a swirling debris cloud to deviations from axisymmetry differs from that of the air. The airflow is more stable because it is incompressible to good approximation. Figure 9 provides an example. Although still highly turbulent on smaller scales, larger-scale fluctuations in the swirl velocity are suppressed by pressure forces in the incompressible airflow without debris (Fig. 9a). On the other hand, the sand flow is highly compressible. Centrifuging of populations leading to large deviation from axisymmetry is not directly suppressed. This leads to large instantaneous asymmetries in D (Fig. 9c) and, through the drag force, even in the air velocity (Fig. 9b). The swirling sand cloud is stabilized only indirectly via its interaction with the air through the drag force. When the sand cloud acquires enough inertia in comparison to the airflow (e.g., simulation S0.2), a greater wander or distortion of the entire central vortex can arise.

Fig. 9.

Instantaneous horizontal cross sections at 50-m height: (a) swirl velocity (m s−1) defined about the central axis from simulation SN; (b) swirl velocity from simulation S0.5; (c) debris loading from simulation S0.5. Areas shown are 600 m on a side. Contour intervals are 10 for velocities and 0.1 for debris loading with darker shading indicating greater values.

Fig. 9.

Instantaneous horizontal cross sections at 50-m height: (a) swirl velocity (m s−1) defined about the central axis from simulation SN; (b) swirl velocity from simulation S0.5; (c) debris loading from simulation S0.5. Areas shown are 600 m on a side. Contour intervals are 10 for velocities and 0.1 for debris loading with darker shading indicating greater values.

d. Sensitivity to surface treatment

Simulations SS1–SS4, STh.2, and STh1 differ from S0.5 only by choices of surface parameters. Even with the idealization of a single particle type, surfaces can physically differ, for example, in how strongly the particles adhere or how they are situated relative to surface roughness elements (e.g., vegetation). Here, the aim in altering surface parameters was to make a first assessment of sensitivity rather than to try to replicate a specific physical surface. The resulting changes are qualitatively all in the expected direction, but at a smaller level than one might expect. For example, doubling or halving the surface flux coefficient ca (simulations SS1 and SS4) increases or decreases the total mass of the sand cloud, but only by about 25%. A factor of five change in threshold coefficient ct (STh1, STh.2) changes the sand cloud mass, but by less than 20%. Adding a source term from particle collisions with the surface (coefficient cd), with or without allowing a fraction of particles hitting the surface to rebound (by altering cb, cu), increases the total debris cloud mass ∫D, but again by little. In all these cases the changes in other statistics from Table 1, as well as mean spatial distributions of D (not shown), are commensurate with the change in ∫D or smaller. There is a strong negative feedback limiting the sensitivity of the debris effects to details of the surface parameterization, at least in quasi-steady state. A change that makes the debris easier to pick up (e.g., raising ca or cd or lowering ct) raises D in the lowest layer. As a consequence, the wind velocities there fall because of the transfer of momentum to the debris. This lower velocity then picks up less debris, limiting the total loading and its effects.

5. Concluding remarks

This work represents a first attempt at assessing the effects of finescale debris on tornadoes. Much remains to be done to improve modeling of debris within tornadoes, to extend the parameter space explored, and to compare with real tornadoes. Nonetheless, the present results suffice to illustrate the importance of including debris when studying tornado dynamics: at least in some regimes the accumulation of small-scale debris within the surface layer and corner flow can significantly alter the wind speeds and flow structure of the tornado vortex within a few hundred meters of the surface. The local debris mass near the surface can significantly exceed that of the air and remain important to altitudes of several hundred meters, particularly when fine particles are available to be lofted. In cases simulated so far the addition of debris has reduced maximum wind velocities. Near the surface, however, total momentum (air plus debris) is sometimes significantly increased because of large debris loading. Given the potential for “sand blasting” effects, it is likely that incorporation of small-scale debris could either increase or decrease damage potential, depending on circumstances. Centrifuging, negative buoyancy, and angular momentum transport were all found to be important mechanisms through which debris can alter the tornado corner flow. Local air and debris velocities can differ significantly, particularly within the corner flow and near-surface inflow layer. These differences are an important consideration in interpreting Doppler radar observations within these regions.

For this work the debris was treated as a pressureless fluid, interacting with the air through drag forces. “Eulerian–Eulerian” LES treatments of multiphase flows are still largely in their developmental stage (e.g., Riber et al. 2006). Several numerical issues had to be overcome to achieve a successful implementation for this work. Further improvements in efficiency and accuracy should prove useful for other LES applications as well.

Three dimensionless parameters governing debris dynamics in tornadoes given a single debris type were also identified. These are the corner flow swirl ratio Sc (Lewellen et al. 2000), providing a measure of the “type” of tornado; a ratio of a characteristic tornado velocity scale to the terminal velocity of the debris in free fall, Aυ, providing a measure of the ease of debris lofting; and a ratio of a characteristic radial acceleration in the core flow to gravity, Aa, providing a measure of the relative importance of debris centrifuging. This work has concentrated on effects of varying Aυ within a strong, medium-swirl, quasi-steady tornado. The results suggest that debris cloud mass, spatial dimensions, and the reduction in velocities due to debris loading all rise monotonically with Aυ.

The time scales for formation of quasi-steady particle clouds (a few times tD in Table 1) can be several minutes for fine particles. If this exceeds the time over which the underlying large-scale tornado flow is changing, then new interactions between the debris and tornado dynamics may arise beyond those considered here. For example, in adding idealized sand to one of the rapidly evolving “corner flow collapse” cases discussed in Lewellen and Lewellen (2007b), it was found that the sand cloud was largest well after the tornado had passed its most intense stage. The peak flow velocities at the most intense stage were also found to be reduced by a third.

Several extensions of this work are under investigation. Quantitative parameterization of debris effects on tornado structure and damage potential and of the difference between air and debris velocity fields are needed. The sensitivity of the particle cloud properties to the particle size indicates that a more accurate assessment of debris effects on tornadoes will require multiple debris species and a proper accounting of their relative availability. Simulations are now being conducted both with multiple debris species and with surface debris layers of limited depth. Preliminary results show the tendency of the corner flow to accumulate debris; thus, significant debris clouds can develop in translating tornadoes even with very limited debris depths available on the surface. For the present work a simple parameterization of the surface debris and momentum fluxes was used. Representing the fluxes for a specific surface remains more challenging. Fortunately, it was found that a negative feedback greatly reduces the sensitivity to details of the surface parameterization. Nonetheless, more work is warranted here, particularly for realistic soils involving a spectrum of particles. The interaction of tornadoes with vegetation and man-made structures adds a rich variety of debris possibilities as well.

The debris cloud, condensation funnel, and damage track are the primary visual signatures of tornadoes; realistically simulating their appearance may allow us to correlate the observed structure of some tornadoes and/or their damage tracks with the known structure and velocity measurements available from simulations. Producing simulated radar returns from the modeled debris clouds could aid interpretation of actual radar observations. The range of tornado types, strengths, and evolutions—and the surface conditions that they could encounter—is large; only a small subset has been sampled to date. Transient effects may prove of particular interest. Simulations of evolving tornadoes without debris sometimes show dramatic effects that do not have quasi-steady counterparts (e.g., corner flow collapse). Given that debris can dramatically change the near-surface flow, and that changes in the near-surface inflow can lead to dramatic changes in the corner flow, it is important to assess whether, for example, changing debris availability (e.g., translating into a freshly plowed field, or across a pond) could itself lead to dramatic transient intensification (or de-intensification) in some scenarios.

Acknowledgments

This work was supported by the National Science Foundation, Grants ATM236667 and ATM-0635681. We thank Paul Lewellen for developing the software for rendering figures 1 and 3 from our simulation data.

REFERENCES

REFERENCES
Batt
,
R. G.
,
M. P.
Petach
,
S. A.
Peabody
, and
R. R.
Batt
,
1999
:
Boundary layer entrainment of sand-sized particles at high speed.
J. Fluid Mech.
,
392
,
335
360
.
Bluestein
,
H. B.
,
C. C.
Weiss
, and
A. L.
Pazmany
,
2004
:
Doppler radar observations of dust devils in Texas.
Mon. Wea. Rev.
,
132
,
209
224
.
Crowe
,
C. T.
,
T. R.
Troutt
, and
J. N.
Chung
,
1996
:
Numerical models for two-phase turbulent flows.
Annu. Rev. Fluid Mech.
,
28
,
11
43
.
Dowell
,
D. C.
,
C. R.
Alexander
,
J. M.
Wurman
, and
L. J.
Wicker
,
2005
:
Centrifuging of hydrometeors and debris in tornadoes: Radar-reflectivity patterns and wind-measurement errors.
Mon. Wea. Rev.
,
133
,
1501
1524
.
Druzhinin
,
O.
, and
S.
Elghobashi
,
1999
:
On the decay rate of isotropic turbulence laden with microparticles.
Phys. Fluids
,
11
,
602
610
.
Farnell
,
L.
,
1980
:
Solution of Poisson equations on a nonuniform grid.
J. Comput. Phys.
,
35
,
408
425
.
Ferry
,
J.
,
S. L.
Rani
, and
S.
Balachandar
,
2005
:
A locally implicit improvement of the equilibrium Eulerian method.
Int. J. Multiphase Flow
,
29
,
869
891
.
Gong
,
B.
,
2006
:
Large-eddy simulation of the effects of debris on tornado dynamics. Ph.D. thesis, West Virginia University, 187 pp
.
Gong
,
B.
,
D. C.
Lewellen
, and
W. S.
Lewellen
,
2006
:
Effects of fine-scale debris on different tornado corner flows. Preprints, 23rd Conf. on Severe Local Storms, St. Louis, MO, Amer. Meteor. Soc., P10.3. [Available online at http://ams.confex.com/ams/pdfpapers/115317.pdf.]
.
Gu
,
Z.
,
Y.
Zhao
,
Y.
Li
,
Y.
Yu
, and
X.
Feng
,
2006
:
Numerical simulation of dust lifting within dust devils—Simulation of an intense vortex.
J. Atmos. Sci.
,
63
,
2630
2641
.
Kanak
,
K. M.
,
2005
:
Numerical simulation of dust devil-scale vortices.
Quart. J. Roy. Meteor. Soc.
,
131
,
1271
1292
.
Lewellen
,
D. C.
, and
W. S.
Lewellen
,
2007a
:
Near-surface intensification of tornado vortices.
J. Atmos. Sci.
,
64
,
2176
2194
.
Lewellen
,
D. C.
, and
W. S.
Lewellen
,
2007b
:
Near-surface vortex intensification through corner flow collapse.
J. Atmos. Sci.
,
64
,
2195
2209
.
Lewellen
,
D. C.
,
W. S.
Lewellen
, and
J.
Xia
,
2000
:
The influence of a local swirl ratio on tornado intensification near the surface.
J. Atmos. Sci.
,
57
,
527
544
.
Lewellen
,
D. C.
,
B.
Gong
, and
W. S.
Lewellen
,
2004
:
Effects of debris on near-surface tornado dynamics. Preprints, 22nd Conf. on Severe Local Storms, Hyannis, MA, Amer. Meteor. Soc., 15.5. [Available online at http://ams.confex.com/ams/pdfpapers/81449.pdf.]
.
Lewellen
,
W. S.
,
D. C.
Lewellen
, and
R. I.
Sykes
,
1997
:
Large-eddy simulation of a tornado’s interaction with the surface.
J. Atmos. Sci.
,
54
,
581
605
.
Magsig
,
M. A.
, and
J. T.
Snow
,
1998
:
Long-distance debris transport by tornadic thunderstorms. Part I: The 7 May 1995 supercell thunderstorm.
Mon. Wea. Rev.
,
126
,
1430
1449
.
Marble
,
F. E.
,
1970
:
Dynamics of dusty gases.
Annu. Rev. Fluid Mech.
,
2
,
397
446
.
Owen
,
P. R.
,
1964
:
Saltation of uniform grains in air.
J. Fluid Mech.
,
20
,
225
242
.
Press
,
W. H.
,
S. A.
Teukolsky
,
W. T.
Vetterling
, and
B. P.
Flannery
,
1992
:
Numerical Recipes in Fortran 77: The Art of Scientific Computing.
2nd ed. Cambridge University Press, 992 pp
.
Riber
,
E.
,
M.
Garcia
,
V.
Moureau
,
H.
Pitsch
,
O.
Simonin
, and
T.
Poinsot
,
2006
:
Evaluation of numerical strategies for LES of two-phase reacting flows. Proc. Summer Program 2006, Stanford, CA, Center for Turbulence Research, 197–211
.
Toigo
,
A. D.
,
M. I.
Richardson
,
S. P.
Ewald
, and
P. J.
Gierasch
,
2003
:
Numerical simulation of Martian dust devils.
J. Geophys. Res.
,
108
.
5407, doi:10.1029/2002JE002002
.
White
,
F. M.
,
1991
:
Viscous Fluid Flow.
McGraw-Hill, 614 pp
.
Yamamoto
,
Y.
,
M.
Potthoff
,
T.
Tanaka
,
T.
Kajishima
, and
Y.
Tsuji
,
2001
:
Large-eddy simulation of turbulent gas–particle flow in a vertical channel: Effect of considering inter-particle collisions.
J. Fluid Mech.
,
442
,
303
334
.

APPENDIX

Some Details of the Numerical Model Debris Implementation

The numerical implementation of the debris fluid follows that of the air where possible; the main differences are summarized here. The exact conservation of the debris loading (D) and debris momentum components (𝒰i) has been maintained at the finite difference level by employing flux forms for the advection terms and imposing limits, when required, on the fluxes rather than on D and 𝒰i themselves. For simplicity, consider the 1D versions of the debris loading and momentum equations with only advection terms included:

 
formula
 
formula

which we discretize (leapfrogged in time) on staggered grids in space (D defined on grid centers with momenta defined on corresponding grid faces) as

 
formula
 
formula

The simplest discretization for the flux, ℱi ≡ (𝒰2/D)i = (𝒰i+1/2 + 𝒰i−1/2)2/4Di, is unstable for any size time step. By assuming a more general form, ℱi = [β1𝒰i+1/2𝒰i−1/2 + β2(𝒰2i+1/2 + 𝒰2i − 1/2)]/[α1Di + α2(Di+1 + Di−1)], demanding second-order accuracy, and performing a (somewhat lengthy) von Neumann stability analysis (e.g., Press et al. 1992), it was found that the form

 
formula
 
formula

gives a theoretically stable discretization, provided that κ > 1/2 and a Courant condition,

 
formula

is satisfied at all grid points. In three dimensions there are off-diagonal flux terms as well, which are more straightforwardly discretized and centered at the necessary grid points; for example,

 
formula

We choose κ = 1 and run with a variable time step. Courant conditions for air and debris velocities are checked at each grid point and time step and maintained to within an additional safety factor, with dt halved as needed or doubled as allowed.

Several complications arise in regions where D → 0. The flux in (A5) and computed debris velocities become ill defined, opening the door to additional instabilities and a requirement of an artificially small time step to satisfy (A7). These complications are treated in the momentum equations as follows: (i) the local debris velocity Courant condition is used to limit dt only if D > Dthresh; otherwise the debris velocities appearing in the flux terms are themselves limited, if necessary, to satisfy the Courant condition; (ii) the fluxes at a point are set to zero if D centered on that point falls below some numerical cut δm; and (iii) if Di < Dthresh, 𝒰i−1/2 > 0, and 𝒰i+1/2 < 0, then the discretization for the diagonal flux term is changed from (A5) to ℱi = (𝒰2i+1/2 + 𝒰2i−1/2)/2Di.

Modification of (A3) is required to avoid instabilities that can arise if D is driven negative locally. The debris fluxes in (A3) have been modified (𝒰 → 𝓊?) to limit the maximum fraction of debris within any grid cell that can be removed by a single flux component in any time step; that is,

 
formula

For the simulations in this work we have taken Dthresh = 0.02, δm = 10−5, and df = 0.25, which prove to be small (or, in the case of df , large) enough that results are unaffected by the precise value chosen but sufficient to avoid instabilities or numerically significant D < 0.

In an incompressible fluid, the pressure force insuring incompressibility serves to damp and smooth many velocity fluctuations. For the compressible debris fluid, these stabilizing effects can occur only indirectly, through the drag force. The connection between the drag and pressure forces is most clear in the limiting case of uniformly distributed fine dust: the two-fluid equations then reduce to incompressible equations for a single fluid of combined density, with the drag force simply proportional to the pressure force (Marble 1970). Not surprisingly, care is required in the drag implementation or rapidly growing instabilities can arise. Evaluating the drag terms using intermediate air and debris velocity values (i.e., a standard leapfrog approach) produced a strong oscillating instability, so values from the time step being updated were used instead. If the drag time scale τυ is too small relative to dt then the drag terms overdamp and a different oscillatory instability arises. By performing a von Neumann stability analysis on the reduced air and debris momentum equations featuring just the drag terms (and constant τυ),

 
formula

one finds that stability requires

 
formula

In regions with D → 0 (in which ud will be ill defined), τυ as computed from (3) and (4) can sometimes give unphysically small values, making (A11) difficult to satisfy. In analogous fashion to the handling of (A7), we ensure that dt is small enough to satisfy (A11) for physically reasonable particle velocities and then use (A11) to put a lower (numerical) limit on the actual τυ employed at a given grid point if necessary. Finally, to avoid a diffusion instability arising on rare occasions in regions with essentially no debris loading, the subgrid diffusion of D is turned off when D < δm.

Footnotes

Corresponding author address: D. C. Lewellen, MAE Dept., WVU, P.O. Box 6106, Morgantown, WV 26506–6106. Email: dclewellen@mail.wvu.edu