c
1998 Cambridge University Press
The role of the turbulent scales in the settling
velocity of heavy particles in homogeneous
isotropic turbulence
By C. Y. Y A N G A N D U. L E I†
Institute of Applied Mechanics, National Taiwan University, Taipei 106, Taiwan, Republic of China
(Received 18 April 1997 and in revised form 21 April 1998)
The average settling velocity of heavy particles under a body force field is studied numerically in stationary homogeneous isotropic turbulent flows generated by the direct numerical simulation and the large eddy simulation of the continuity and Navier–Stokes equations. The flow fields corresponding to different selected ranges of turbulent scales are obtained by filtering the results of the direct numerical simulation, and employed for calculating the particle motion. Wang & Maxey (1993) showed that as a consequence of the particle accumulation in the low vorticity region and the preferential sweeping phenomenon, the average settling rate in turbulence is greater than that in still fluid. In the present study, the phenomenon of particle accumulation in the low vorticity region is found to be controlled mainly by the small scales with wavenumber kω corresponding to the maximum of the dissipation (vorticity)
spectrum. However, the increase of the average settling rate, h∆vSi, also depends
strongly on the large energetic eddies. The small eddies with wavenumber greater than 2.5kω have essentially no effect on the particle accumulation and the average
settling velocity. The large eddy simulation is thus adequate for the present study provided the smallest resolved scale is greater than 1/(2.5kω). Detailed calculations
show that h∆vSi is maximized and is of order u0/10 when τp/τk ≈ 1 and vd/u0 ≈ 0.5
for Rλ = 22.6–153, where τp is the particle’s relaxation time, τk is the Kolmogorov
timescale, vd is the settling rate in still fluid, u0 is the root mean square of the velocity
fluctuation, and Rλ is the Reynolds number based on the Taylor microscale.
1. Introduction
Knowledge of the ensemble-average settling rate of small heavy particles in turbu-lent flow is important for many applications in environmental science and engineering. Examples include the settling of aerosol particles in the atmosphere and water droplets in clouds. Consider the settling of heavy spherical particles in homogeneous isotropic turbulence with zero mean flow. It was generally assumed before the 1980s that the ensemble average of the settling velocity of particles in turbulence, hvSi, equals the
settling velocity in still fluid, vd. For example, Reeks (1977) has argued thathvSi = vd
if the flow field is treated as a stochastic random field, where hvSi and vd are the
magnitudes of hvSi and vd, respectively. Here the directions of hvSi and vd are the
same as the direction of the body force field. Maxey (1987) showed thathvSi is indeed
equal to vd for zero-inertia particles through a rigorous analysis. For large-inertia
particles, the particle aerodynamical response time is much greater than any integral timescale of the flow; the turbulent fluid velocity seen by a particle appears as an un-correlated random noise and should have negligible effect on the settling velocity. For particles with finite inertia, Maxey showed thathvSi is greater than vd in a Gaussian
random flow field. The effect of the dynamics of turbulence was included by Yeh & Lei (1991), and they also found that hvSi is greater than vd according to a simulation
of the particle motion in a decaying homogeneous isotropic turbulent flow field gen-erated by large eddy simulation (LES). However, the increase of the magnitude of the average settling velocity, h∆vSi, obtained by Yeh & Lei is smaller than that obtained
by Maxey (1987). Hereh∆vSi = hvSi−vd. Recently, Wang & Maxey (1993) studied the
settling of particles in a stationary homogeneous isotropic turbulence generated by the direct numerical simulation of the continuity and Navier–Stokes equations (DNS), and found substantially greater values of h∆vSi in comparison with those obtained
by Yeh & Lei. A significant increase in h∆vSi can occur for particles with inertial
response time (τp) and vd comparable to the Kolmogorov scales of the turbulence.
Wang & Maxey pointed out that the flow field generated by LES like that in Yeh & Lei is not adequate for evaluating the present problem. The reason is that the particle accumulation in the low vorticity region (one of the key mechanisms for the increase of the settling velocity) is related to the small scales of turbulence, which are not resolved in LES. According to their results, Wang & Maxey also suggested that the Kolmogorov scaling is important for the phenomenon associated with the increase of settling velocity, which is interesting and essentially correct, but is in contrast to the traditional approach that the large-scale, energy-containing fluid motions dominates the transport of particles.
According to the direct numerical simulations of homogeneous turbulence (for examples, see Vincent & Meneguzzi 1991; Jim´enez et al. 1993), the vorticity structure of the flow is found to be organized in coherent, cylindrical or ribbon-like, vortices, which are also known as ‘worms’ according to Jim´enez et al. A statistical study by Jim´enez et al. suggests that the vorticies are simply especially intense features of the background vorticity. The radii and the lengths of the vorticies are of the order of the Kolomogorov and integral length scales (η and Lf), respectively. The spacing between
the vorticies is of order λ, the Taylor microscale (Vincent & Meneguzzi). Jim´enez et al. also found that the circulation of the intense vorticies, γ, increases as νR0.5
λ ,
where Rλ = u0λ/ν is the Reynolds number based on the Taylor microscale, and ν is
the kinematic viscosity. Here u0 is the root mean square of the velocity fluctuation, which is also the velocity scale associated with the large energetic eddies. Let US
be the velocity scale outside the intense vorticies. With γ ∼ USη and η/λ ∼ R−0.5λ ,
it follows that US ∼ u0, which indicates that the velocity scale in the low vorticity
regions between the intense vortices (‘worms’) is of order u0.
As the diameter of the particles, dp, is in general of one order less than the
Kolmogorov length scale, the motion of a particle is exposed to the action of the whole range of the length scales of turbulence. When the particles settle under gravity (or another body force field) through the flow field, the particle inertia produces a bias in each trajectory towards regions of low vorticity due to the local centrifugal effect of the tubular vortical structures (‘worms’), and thus the particles tend to accumulate in the low vorticity (or high strain rate) region. This phenomenon of particle accumulation was simulated by Squires & Eaton (1991) and Wang & Maxey (1993). The degree of particle accumulation depends on the particle inertia, and is maximized at a certain finite value in a given flow field under the action of a specified body force. An instantaneous distribution of particles in a velocity field is
shown in figure 12 of Wang & Maxey (1993). It can be observed clearly from the figure that there are essentially no particles in the cores of the vortices (high vorticity region), but many particles are collected in elongated sheets on the peripheries of local vortical structures (low vorticity region). Furthermore, in the low vorticity region, more particles are accumulated at the downward side (with the angle between the local fluid velocity and the body force, θ, less than 90◦) in comparsion with the upward side (90◦ < θ < 180◦), which is called preferential sweeping by Wang & Maxey. As the local drag on the particle is less at the downward side and more particles are there, the ensemble average of the settling velocity of particles is thus increased.
Although the fundamental framework, including the key results and the underlying physics, of the increase of the average settling velocity of particles in homogeneous stationary isotropic turbulence has been completed successfully by Maxey (1987) and Wang & Maxey (1993), the roles of different turbulent scales are still unclear and worth further investigation. The reasons are as follows:
(1) The accumulation of particles in the low vorticity region is related strongly to the structure of the tubular vortices (of diameter η approximately), or the details of the instantaneous vorticity field. This is probably the main reason that Wang & Maxey employed DNS to generate their flow field. However, the length scale corresponding to the maximum of the vorticity (dissipation) spectrum, lω, is in
general of one order greater than the Kolmogorov length scale, η (see Tennekes & Lumley 1972 for example). Also the vorticity associated with η is of one order less than that associated with lω, and the energy associated with the Kolmogorov eddies
are negligible. It is possible that the Kolmogorov eddies might have a significant effect on the accumulation of particles, and hence provide a contribution to the increase of the settling velocity. Here lω = 1/kω, with kω the wavenumber corresponding to the
maximum of the dissipation spectrum. From the numerical point of view, should one always need to resolve accurately the Kolmogorov eddies for the correct estimation of the increase of the settling velocity? If the Kolmogorov eddies are essentially not important, one may still employ the flow field generated by large eddy simulation (LES) with suitable cutoff length scales (smaller than those in Yeh & Lei 1991) to study the average settling velocity.
(2) Although Wang & Maxey found a substantial increase of the average settling velocity when vd ≈ vk (the Kolmogorov velocity scale) and τp ≈ τk (the Kolmogorov
time scale), and demonstrated the importance of the Kolmogorov scaling for the present problem, the role of the large eddies deserves further examination. As the particles are accumulated in the region of low vorticity where the velocity scale of the flow is u0, it is expected that the large energetic eddies should play a significant role on the drag and thus the increase of the average settling velocity. Note that the Kolmogorov time scale, τk, proposed by Wang & Maxey is indeed a suitable time
scale for the present problem. The reason is that τ−1k is also the scale for the vorticity corresponding to the ‘eddies’ associated with lω, which are responsible for the particle
accumulation.
The major work of the present paper is to carry out a detailed and careful study on the effects of different length scales on the average settling velocity of particles in stationary homogeneous isotropic turbulence generated by DNS and LES. The paper is organized as follows. The formulation and numerical method are described in § 2. The results for stationary turbulence are presented and discussed in§ 3. The validity of the present calculation is first checked against the result by Wang & Maxey. Effects of small and large scales are then studied in detail. The validity of using the flow field generated by LES with suitable cutoff wavenumber for studying the settling velocity
is checked against the results using DNS. LES is then employed to extend the range of Rλ (which possesses a short inertial subrange) of the flow field for the present
study. Detailed effects of parameters on h∆vSi are presented, and the scales for the present problem are discussed. Finally, the paper is concluded in§ 4.
2. Formulation and numerical method
2.1. Formulation
Consider an incompressible Newtonian fluid with density ρ and constant kinematic viscosity ν. The governing equations for the homogeneous isotropic turbulent flow field in an inertial frame, written in index form, are
∂ui ∂xi = 0, (1a) ∂ui ∂t + ∂(uiuj) ∂xj =−1 ρ ∂p ∂xi + ν ∂ 2u i ∂xj∂xj + fi, (1b)
where ui(i = 1, 2, 3) are the velocity components, t is the time, xi(i = 1, 2, 3) are the
spatial coordinates, p is the pressure, and fi is the forcing term for maintaining the
turbulence stationary. The motion of a single spherical solid particle with diameter
dp and mass mp in the turbulent gas flow field is assumed to follow
mp
dv
dt = 3πdpρνC(u− v) + mpg, (2a)
and the particle’s trajectory is evaluated by dY
dt = v, (2b)
where Y (t) and v(t) are the instantaneous displacement and velocity of the particle at time t, u = u(Y (t), t) is the instantaneous fluid velocity where the particle is located, g is the body force per unit mass of the particle, and C is a factor accounting for the inertial effect of the drag term. In the present study, we set C = 1 + 0.1315R0.82−0.05 log10Rp
p ,
which is an experimental correlation for the standard drag curve of a single sphere in a uniformly moving stream, and is valid for 0.01 6 Rp 6 20 (see Clift, Grace &
Weber 1978, p. 112). Here Rp = |u − v|dp/ν is the particle’s Reynolds number. The
settling velocity in still fluid vd = gτp/C0, where τp = mp/(3πdpρν) is the particle’s
relaxation time (inertial response time) based on the Stokes drag law, and C0 is the value of C evaluated with u = 0. The added mass effect, the Basset force, and the buoyancy force are all neglected in (2a) since the density of the particle, ρp, is much
greater than the fluid density, ρ. Here we set ρp/ρ = 1000 in the calculations. In the
drag term of (2a), u is obtained by solving the transient three-dimensional governing equations ((1a) and (1b)) using either DNS or LES. Motions of thousands of particles are simulated numerically by solving (2a) and (2b), and employed for studying the ensemble average of the particle properties.
2.2. Numerical method using DNS
The method and procedures employed here for solving the flow field are basically the same as those in Lee & Reynolds (1985), which are based on Rogallo (1977, 1981). Such procedures have also been employed and extended successfully by Yeh & Lei (1991) for simulating the particle motion in a decaying homogeneous isotropic
turbulence generated by LES. For homogeneous isotropic turbulence, the calculation domain for (1a) and (1b) is a cube with edge length L subject to the periodic boundary conditions. Here L should be sufficiently larger than the separation at which two-point spatial correlations vanish and the length scale associated with the particle motion,
vdτp (see Yeh & Lei 1991). In accordance with the periodic boundary conditions, the
velocity and pressure fields are expanded by three-dimensional discrete Fourier series with N wavenumbers in any direction. The linear spatial terms in (1a) and (1b) are evaluated easily by direct differentiation. The pseudospectral method was used for evaluating the nonlinear term ∂(uiuj)/∂xj. The aliasing problem is handled by using
the hybrid dealiasing algorithm by Lee & Reynolds (1985), which is modified from the hybrid method proposed by Rogallo (1977, 1981). The doubly and triply aliased terms are eliminated by using a spherical truncation mask set at k2 = 2
9N
2 following Patterson & Orszag (1971). The singly aliased term is handled by the random phase shift method of Rogallo (1977, 1981). The forcing term in the wavenumber space,
ˆ
fi(k, t), is evaluated similar to Eswaran & Pope (1988). To guarantee the continuity
condition, ˆ fi(k, t) = ˆ ai− ki( ˆajkj/k2), if 0 < k < kfc 0, otherwise, (3a)
where k(k1, k2, k3) is the wavenumber, kfc is the forcing radius in the wavenumber
space, and ˆai is assumed to follow
d ˆai dt + ˆ ai Tfc = SfcRi(t)Ai for 0 < k < kfc, (3b)
with the factor
Ai= exp −10 u2 i q2 − 1 3 , q2 = u21+ u22+ u23, i = 1, 2, 3.
Here Ai is employed for reducing the anisotropic effect due to the articifical forcing,
Ri(t) is a random function distributed uniformly between −0.5 and 0.5, Sfc is the
forcing amplitude, and Tfc is the forcing time scale. In the present study, we choose
kfc =
√
8, which is the same as that in Wang & Maxey (1993). The values of Sfc
and Tfc employed in the present study are listed in table 1, and will be discussed
in § 3. Finally, (1a) and (1b) become a set of ordinary differential equations in the wavenumber space, which are solved for ˆu, the Fourier transform of the velocity. The initial three-dimensional, incompressible and isotropic field is constructed using the method in Yeh & Lei (1991). An initial energy spectrum
E(k) = 8 πu 02 0L2f0 k4 (1 + k2L2 f0)3 (4)
is employed for starting the calculation provided u00 and Lf0 are specified. An inverse
transform of ˆu gives the flow field in the physical space, u, which is required for the drag term in (2a) for the particle simulation. In order to study the effect of different scales of turbulence on the particle’s settling velocity, certain modes in the wavenumber domain can be selected for calculating the particle motion using a numerical ‘Fourier space sharp-cutoff filter’ during the inverse Fourier transform. An example is illustrated as follows. Let kmax, kmin, kCH and kCLbe the magnitudes of the
highest simulated wavenumber, the lowest simulated wavenumber, the high end cutoff wavenumber of the filter, and the low end cutoff wavenumber of the filter, respectively. If kCH < kmaxand kCL= kmin, the modes for the range kCH < k < kmaxof the simulated
Case A Case B Case C Case D Case E Case F
L3 (2π)3 (2π)3 (2π)3 (2π)3 (2π)3 (2π)3
N3 323 643 963 323 963 1283
Method DNS DNS DNS LES LES LES
Sfc 180 180 180 180 180 180 Tfc 0.2 0.18 0.195 0.195 0.2 0.2 u0f0 3.0 3.0 3.0 3.0 3.0 3.0 Lf0 1.5 1.5 1.5 1.5 1.5 1.5 ν 0.15 0.08 0.03 0.03 0.0075 0.005 u0 3.35 3.14 3.16 3.10 3.10 2.95 Lf 1.62 1.55 1.60 1.65 1.23 1.08 25.2 14.4 11.8 9.68 10.5 9.7 S −0.47 −0.47 −0.51 −0.42 −0.44 −0.44 Te1 0.46 0.49 0.51 0.53 0.40 0.37 Te2 0.45 0.68 0.85 0.99 0.91 0.90 η 0.108 0.078 0.039 0.041 0.014 0.011 τk 0.077 0.075 0.05 0.055 0.027 0.023 vk 1.39 1.04 0.77 0.74 0.53 0.47 λ 1.01 0.91 0.62 0.68 0.32 0.26 Rλ 22.6 35.7 65.3 69.8 133 153 kmaxη 1.62 1.7 1.71 0.61 0.62 0.63 u0/vk 2.41 3.02 4.10 4.21 5.85 6.28 Lf/η 15 20 41 40 87 102 Lf/λ 1.60 1.70 2.58 2.44 3.84 4.15 Te1/τk 5.97 6.63 10.1 9.67 14.7 16.1 Te2/τk 5.84 9.13 17 18 33.7 39.1 Lf/u03 1.09 0.72 0.60 0.54 0.43 0.41
Table 1.Parameters and flow properties of the stationary homogenoeus isotropic turbulence (in arbitrary consistent units).
velocity field ˆui in the wavenumber domain will be omitted in the inverse Fourier
transform, which implies that the turbulent scales with size smaller than 2π/kCH will
not be allowed to affect the particle motion. A second-order Runge–Kutta scheme is employed for the time advancement for both the flow field in the wavenumber space and the particle motion in the physical space. In order to avoid the error due to the time interpolation, we advance both the flow field and the particle motion with the same time step, which satisfies the CFL condition of the flow field and should be at least of one order less than the particle’s relaxation time (see Yeh & Lei 1991). As the instantaneous location of the particle Y (t) is not at the grid point of the flow field in general, the 13-point spatial interpolation method proposed by Yeung & Pope (1988) is employed to evaluate u(Y (t), t) from its neighbouring grid values in (2a). The triplex averaging procedure proposed by Wang & Maxey (1993) is also employed for the particle calculation in order to reduce the statistical error. Thousands of particles were introduced uniformly into the flow field after the skewness of the velocity field had reached its asymptotic value. The initial velocity of a given particle is set equal to vd plus the local fluid velocity. The particles were divided into six groups of different
orientations for the body force. The six orientations, g/|g|, are (1, 0, 0), (−1, 0, 0), (0, 1, 0), (0,−1, 0), (0, 0, 1) and (0, 0, −1). Statistical properties of the particle motion were calculated after time TS (defined later) when the average settling velocity of the
consists of three steps. First, the ensemble average is obtained among particles in each group of a specified orientation of g. Next, the arithmetic average is taken among the six group of different orientations. Finally, the time average is calculated over the interval starting from TS to the end of the particle simulation (TE). Wang & Maxey
chose TS = 2Te2 in their calculation, where Te2 is the eddy turnover time (defined in
§ 3).
2.3. Numerical method using LES
The method is the same as the method using DNS in the last section, except LES is employed to generate the flow field. The major difference between LES and DNS is that the high wavenumber components of the flow field (say, k > kC) are not resolved
in the LES, but the effect of the cutoff wavenumber components (for k > kC) on the
resolved low wavenumber components are modelled through the subgrid turbulent modelling. In the present study, the velocity and pressure fields in (1a) and (1b) are decomposed into a large-scale part and a subgrid-scale part,
ui(x, t) = ui(x, t) + u00(x, t), p(x, t) = p(x, t) + p00(x, t). (5)
The large-scale fields (resolved fields), ui and p, are defined by mean of a convolution
filter (see Ferziger 1983). The sharp-cutoff filter in Fourier space is employed in the present study. The equations governing the large-scale field are obtained by filtering (1a) and (1b). They are
∂ui ∂xi = 0, (6a) ∂ui ∂t + ∂ ∂xj ( ¯uiu¯j) =− ∂P ∂xi −∂τij ∂xj + ν ∂ 2u i ∂xj∂xj + fi, (6b) where P = p/ρ +13Qkk, τij = Qij− 13Qkkδij, with Qij= uiu00j + uju00i + u00iu00j.
The subgrid-scale Reynolds stress τij is related to the large-scale field by using the
eddy viscosity model :
τij = 2νTSij = νT ∂ui ∂xj +∂uj ∂xi , (7a)
with the eddy viscosity νT represented by
νT = (CS∆f)2(2SijSij)0.5. (7b)
Here the coefficient CS is taken from McMillan & Ferziger (1979),
CS = 0.128 1 +24.5 Rsgs −1 , Rsgs= (SijSij)0.5∆2f ν . (7c)
Note that the forcing term fi in (6b) is in general not affected by the filtering process,
since the artifical forcing is applied to the low wavenumber components while the filtering process in LES is applied to the high wavenumber components. By comparing (6a) and (6b) for LES with (1a) and (1b) for DNS, we found that the equations are of the same form for both methods if ui, P and ¯uiu¯j+ τij in LES are replaced by ui,
similar to those in the last section, with the subgrid Reynolds stress term calculated in a way similar to the viscous term. Details of the calculation of the particle motion are also the same as those in§ 2.2 provided ui in (2a) is replaced by ui.
3. Results in stationary turbulence
3.1. Flow characteristics
As discussed in Ferziger (1983), the DNS/LES calculation has three stages: the initial (developing) period, the mature period, and the final period when the assumptions in the calculation are violated. The end of the initial period is characterized when the skewness reaches its asymptotic value (around−0.4 to −0.5) and kmaxη> 1. The
mature period is the time range we employ to study the turbulence characterestics. We start our particle calculation at the beginning of the mature period, evaluate the particle’s statistical properties after the effect of the artifical initial conditions is negligible and the average settling velocity of particles reaches its asymptotic constant value (see§ 3.2), and terminate the particle calculation at the end of the mature period. Six cases were simulated in the present study, and were employed to study the particle statistics. Table 1 lists the parameters for the simulation and the flow prop-erties of the mature period (in arbitrary consistent units). The parameters for the simulation include the size of the calculation domain (L3), the number of grids (N3), the method for generating the flow field, the forcing amplitude (Sfc) and the forcing
time scale (Tfc) associated with the artifical forcing term, the initial root mean square
of the velocity fluctation (u0f0) and the inital integral length scale (Lf0) for starting
the calculations in (4), and the kinematic viscosity (ν). Note that cases C and D are similar except they are simulated using different methods (and hence different grids). The flow properties of the mature period include the root mean square of the velocity fluctuation (u0), the integral length scale (Lf), the dissipation (), the skewness (S ),
the eddy turnover time scales (Te1 and Te2), the Kolmogorov length scale (η), the
Kolmogorov timescale (τk), the Kolmogorov velocity scale (vk), the Taylor microscale
for length (λ), the Reynolds number based on λ (Rλ), and the greatest dimensionless
wavenumber simulated in the calculation (kmaxη). Here
u0= (u2)1/2, L f = 13(Lf1+ Lf2+ Lf3), = 2ν Z ∞ 0 k2E(k)dk, (8a) and S = (∂u1/∂x1) 3 [(∂u1/∂x1)2]3/2 , (8b)
where the integral length scale is defined through
u2 iLfi=
Z ∞ 0
ui(r)ui(0)dr, i = 1, 2, 3, (9a)
the energy spectrum E(k) is defined as 3 2u02=
Z ∞ 0
E(k)dk, (9b)
and (·) denotes the average of the quantity (·) evaluated from the values at the grids. With u0, Lf, and S calculated from the numercial results of the flow field, one can
evaluate Te1= Lf u0, Te2= u 02/, (10a) η = ν3 1/4 , τk= ν 1/2 , vk= (ν)1/4, (10b) λ = 15ν 1/2 u0= 151/2τku0, (10c) and Rλ= u0λ ν . (10d)
Note that Te2 is the eddy turnover time scale employed in Wang & Maxey (1993).
Also listed in table 1 are the dimensionless quantities u0/vk, Lf/η, Lf/λ, Te1/τk, Te2/τk
and Lf/u03. These quantites possess monotonically varying behaviour as Rλ increases
as expected if either case C or D is omitted from the list. Cases C and D have similar values of Rλ, but are calculated using different methods. Table 1 shows that the results
for case C and D are similar, which indicates the validity of the present LES. We have tried different values of Sfc and Tfc during the calculations, and checked the
asymptotically stationary property of the time evolution of Rλ, u0and S . Note that the
values of the skewness for all cases in table 1 have reached the theoretical asymptotic value.
Figure 1 shows the energy and the dissipation spectra for the simulated cases in the present study. Note that the vorticity spectrum is related to the dissipation spectrum through a factor 2ν. In order to check our calculations, the theoretical result based on Townsend’s (1951) model and the calculation for Rλ ≈ 150 according to Vincent
& Meneguzzi (1991) are also plotted in figure 1. Note that the definition for Rλ in
Vincent & Meneguzzi is slightly different from that in the present study. The present results for the dissipation spectra agree asymptotically with Townsend’s model for large kη, and the results for large Rλ(based on LES) agree qualitatively with the DNS
calculation by Vincent & Meneguzzi. The results at low wavenumber for large Rλ are
quantitatively different from the result by Vincent & Meneguzzi, which is expected because different forcing schemes were applied for different calculations. Also note that the results for large Rλ show a short inertial subrange with slope−5/3 like that
in Vincent & Meneguzzi. The case for Rλ≈ 65 has been simulated by both the DNS
(case C, 963 grids, R
λ = 65.3) and the LES (case D, 323 grids, Rλ = 69.8) in the
present study, and the results from both methods in figure 1 agree nicely with each other. Detailed examination of the data shows that the LES resolves approximately 95% of the total energy and 75% of the total dissipation of the results using DNS. The wavenumber corresponding to the peak value of the dissipation spectrum (or the eddies with maximum vorticity), kω, decreases and tends to approach a constant value
as Rλincreases. The values for kωη are 0.172, 0.134 and 0.122 for Rλ equal to 65.3, 133
and 153, respectively, based on the present calculations. The data for Rλ= 22.6 and
35.7 in figure 1(b) show some fluctuations at low wavenumbers, which is due to the forcing associated with the simulation of stationary turbulence. The corresponding results for decaying homogeneous isotropic turbulence are smooth (see Yang 1997), and kωη = 0.26 and 0.23 for Rω= 22.6 and 35.7, respectively.
10–4 10–3 10–2 10–1 100 101 10–2 10–1 100 101 kη D /εη (b) 10–6 10–4 10–2 100 1 10 100 kLf E /( u ′2L f /2) (a) slope = –5/3
Figure 1. (a) The energy, and (b) the dissipation spectra for the flow fields.
◦
, Rλ = 22.6; , Rλ = 35.7; 4, Rλ = 65.3 (DNS, 963 grids); +, Rλ = 69.8 (LES, 323 grids);, Rλ = 133 and O,Rλ = 153. The numerical result by Vincent & Meneguzzi (1991) is represented by the thick curve
(Rλ ≈ 150, DNS, 2403 grids). The thin curve in (b) is the theoretical result based on Townsend’s
model (1951).
3.2. The transient time scale for the particle motion
According to Riley & Patterson (1974), the initial stage of the particle calculation is affected by the imposed initial conditions, and they proposed that statistics of the particle motion which are free of initial conditions can be taken after the relative velocity between the particle and the fluid has reached the maximum value. Figure 2(a) shows the time evolution of the ensemble average of the root mean square of the relative velocity component perpendicular to vd after the first two steps of the triplex
0.5 0 3 6 9 tu′/λ 〈Ω (t )〉 /Ω ′ f (b) 1.0 1.5 0 5 10 15 t /τp 〈θv (t )〉 (a) 1 2
Figure 2.Evolution of the ensemble average of (a) the root mean square of the relative velocity component, and (b) the dimensionless enstrophy for Rλ = 65.3 and vd/vk = 1. ——, τp/τk = 0.15;
— — —, τp/τk= 0.25; - - - - -, τp/τk= 1 and — - —, τp/τk= 2.5.
averaging procedures, (hθ⊥(t)i). It is found that the time scale required for hθ⊥(t)i to reach the asymptotic stationary value is of order τp, the particle’s relaxation time.
However, the time scale for evaluating the ensemble average of the settling velocity deserves further consideration. As mentioned in§ 1, the increase of the average settling velocity is a consequence of the non-uniform particle distribution due to the highly intermittent vortical structure. Thus knowledge of the transient time scale for the particle accumulation is important for the correct estimation of the increase of the average settling velocity. As we know from § 1 that the particles tend to accumulate in the region of low vorticity, a good indicator for the accumulation of particles is the value of the mean enstrophy (square of the vorticity) where the particle locates. Figure 2(b) shows the time evolution of hΩ(t)i/Ω0f for different particle parameters. Here hΩ(t)i is the mean enstrophy at the instantaneous locations of the particles after the first two steps of the triplex averaging procedures, and Ωf0 is the mean enstrophy of the flow field evaluated at the grids, which is constant for stationary homogeneous isotropic turbulence. In the absense of particle inertia, there is no particle accumulation, and hΩ(t)i/Ω0f approaches unity at large time. For particles with a given finite value of inertia, the particles tend to accumulate in the region of low vorticity according to Wang & Maxey (1993), andhΩ(t)i/Ωf0 should approach an asymptotic constant value less than unity at large time. Figure 2(b) indeed shows such a result, and the time scale for reaching the approximate constant value is of order
–2 0 0.5 1.0 1.5 t/Ttp 〈∆ vs (t )〉 /〈 ∆ vs (Ttp )〉 (b) –1 3 0 3 6 9 tu′/λ (a) 0.5 1.0 2.0 2.5 2 1 0 0 –0.5 〈∆ vs (t )〉 /vk
Figure 3.Evolution ofh∆vS(t)i using (a) λ/u0 and (b) Ttp as time scales. The curves are the same as those in figure 2. Here Ttp= 2λ/u0+ 5τp.
can be understood as follows. As the spacing between the tubular vortices (‘worms’) is of order λ according to Vincent & Meneguzzi (1991), and the velocity scale for the particle is of order u0, the velocity scale of the flow in the low vorticity region, it follows that the average time for a particle to reach the low vorticity region is of order λ/u0. According to figure 2, the transient time scale for estimating the increase of the settling velocity is thus related to both the time scale of the particle, τp, and the
time scale of the flow, λ/u0. Figure 3(a) showsh∆vS(t)i, the magnitude of the ensemble
average of the settling velocity after the first two steps of the triplex averaging procedures, and we found that there are different transient response time for different particle inertia. However, if we plot the variation ofh∆vS(t)i with a dimensionless time
normalized by 2λ/u0+ 5τp in figure 3(b), we find that the results for different particle
inertia approach an essentially constant value when t > Ttp= 2λ/u0+ 5τp. Although
the above results in figures 2 and 3 are for Rλ = 65.3, the calculations for other Rλ not
shown in this paper also confirm such findings. Thus Ttpseems to be the appropriate
time scale associated with the transient behaviour ofh∆vS(t)i at least for the ranges of
parameters in the present study. The particle statistics is evaluated only when t> Ttp
during the particle simulation. Wang & Maxey (1993) proposed to take the ensemble average of the settling velocity twice the eddy turnover time (2Te2) after the start of
the particle calculations based on their simulations. Recall that the increase of the ensemble average of the settling velocity is maximized when τp ≈ τk according to
0 1 2 vd/vk 〈∆ vs 〉/ vk (b) 0.6 3 4 0.4 0.2 τp/τk= 1 0 1 2 τp/τk 〈∆ vs 〉/ vk (a) 0.6 3 0.4 0.2 vd/vk= 1
Figure 4.The increase of the average settling velocity for different parameters. Present results:
◦
, and4 for Rλ= 22.6, 35.7 and 65.3, respectively. Results from Wang & Maxey:•
,andN forRλ= 21.2, 31 and 61.5, respectively.
and 6.18 for cases A, B, C, E and F, respectively, according to table 1. In the present study, TS = Ttp is employed for the last stage of the triplex averaging procedures (see
§ 2.2).
3.3. The average settling velocity of particles
The present particle results are validated by comparing the increase of the average settling velocity after the triplex averaging procedures,h∆vSi, with the calculations by
Wang & Maxey (1993). The number of particles for evaluating the particle statistics is 3072× 6 (recall that there are 6 orientations of the body force field for the second step of the triplex averaging procedure). There are two particle parameters for a given flow field, τp and vd(= gτp/C0), which are related to the particle inertia and the strength of the body force field, respectively. Figure 4 showsh∆vSi for different values
of τpand vd. The present results agree nicely with the calculations by Wang & Maxey,
except for some results with large values of vd. The discrepancy is due to the fact that
Rp is sufficiently large such that the results based on the Stokes drag law (C = 1 in
(2a)) in Wang & Maxey are not accurate for large vd. The qualitative behaviour that
the results obtained by using the nonlinear drag law are less than those obtained by using the linear Stokes drag law is consistent with the previous finding according to Wang & Maxey (1993).
Consider the problem in a given flow field (given Rλ). For a fixed value of vd,
0 1 2 vd/vk (b) 0.8 3 4 0.7 0.5 τp/τk= 1 0 1 2 τp/τk (a) 0.7 3 0.6 0.5 vd/vk= 1 0.6 0.9 1.0 〈Ω〉 Ω′ f 0.8 0.9 1.0 〈Ω〉 Ω′ f
Figure 5.The ensemble average of the dimensionless enstrophy at the locations of the particles for different parameters. The symbols are the same as those in figure 4.
the centrifugal effect increases and drives more particles toward the region between vortices (low vorticity region). This effect of particle accumulation together with the preferential sweeping mechanism make the value of h∆vSi increase as τp increases. However, when τp reaches a certain finite value and increases further, the relative
importance of the ‘sideways’ centrifugal effect and the ‘downward’ body force effect decreases. It follows that the degree of particle accumulation, and thush∆vSi, decrease
as τp increases. Finally, when τp is sufficiently greater than any integral time scale
of the flow, the turbulent fluid velocity seen by the particle appears as uncorrelated random noise and h∆vSi drops to zero. Figure 4(a) indicates that h∆vSi is maximized
at τp/τk ≈ 0.5–1.0 for given values of vd/vkand Rλ, andh∆vSi increases as Rλincreases.
The phenomenon of particle accumulation is illustrated in figure 5(a), which shows the variations of hΩi/Ωf0 with τp/τk at vd/vk= 1 for different values of Rλ. HerehΩi
is the ensemble average of the enstrophy at the locations of the particles based on the triplex averaging procedures. It is found that the effect of particle accumulation in the low vorticity region is also maximized at τp/τk ≈ 0.5–1.0. By comparing the detailed
results of figures 4(a) and 5(a), it is found that for a given value of vd, the increase of
the settling velocity is consistently proportional to the degree of particle accumulation. Figure 4(b) shows that for given values of Rλ and τp/τk, there also exists a maximum
value of h∆vSi at a finite value of vd. However, the location of the maximum value of
h∆vSi/vk shifts continuously rightward as Rλ increases, which implies that vk may not
0 1 2 vd/vk (b) 1.3 3 4 1.2 1.0 τp/τk= 1 0 1 2 τp/τk (a) 1.2 3 1.1 1.0 vd/vk= 1 1.1 1.4 1.5 1.3 1.4 1.5 Nr Nr
Figure 6.The parameter for illustrating the preferential sweeping phenomenon, Nr, for different parameters. The symbols are the same as those in figure 4.
in § 3.6. The phenomenon of particle accumulation corresponding to figure 4(b) is shown in figure 5(b). It is found that the effect of particle accumulation is maximized when vd≈ 0, and vanishes as vd increases. Unlike the cases for fixed vd in figures 4(a)
and 5(a), the variation ofh∆vSi for fixed τp in figure 4(b) is not directly related to the
degree of particle accumulation in figure 5(b). The reason is that the accumulation of particles in the low vorticity region is only a necessary condition for the increase of the average settling velocity. It is also required that more particles are accumulated in the ‘downward sweeping’ region in comparison with those in the ‘upward sweeping’ region (the preferential sweeping phenomenon according to Wang & Maxey 1993). Figure 6 is employed for quantifying the preferential sweeping phenomenon. Let u be the local fluid velocity at the instantaneous location of the particle. Define Nr as the
ratio of the number of particles with u· vd > 0 (‘downward sweeping’) to those with
u· vd < 0 (‘upward sweeping’). Figure 6 shows that Nr> 1 for all cases, which implies
that more particles indeed stay in the ‘downward sweeping’ region. For vd/vk= 1, the
variation of Nr with parameters in figure 6(a) is consistent with the results in figures
4(a) and 5(a), and Nr is also maximized at τp/τk ≈ 0.5–1 for different parameters.
For τp/τk= 1, figure 6(b) shows that the preferential sweeping is weak when vd/vk is
small. In fact, there is no preferential direction, and thus no preferential sweeping, when vd = 0 (zero body force). It follows that although we have a large degree of
particle accumulation in figure 5(b) near vd = 0,h∆vSi is not maximized at vd= 0 but
0 1 2 vd/vk (b) 0.4 3 0.2 τp/τk= 1 0 1 2 τp/τk (a) 3 0.6 0.4 vd/vk= 1 0.6 0.2 〈∆ vs 〉/ vk 〈∆ vs 〉/ vk
Figure 7.Effect of small turbulent scales onh∆vSi in stationary turbulence at Rλ= 65.3.
◦
, all simulated modes;, kCHη = 0.5;4, kCHη = 0.4; and +, kCHη = 0.2.6(b), the locations of the maximum value of h∆vSi/vk are consistent with those of
Nr. Also note that when vd is sufficiently large (large ‘crossing trajectory effect’) for a
fixed value of τp, the duration of a particle in an ‘eddy’ is so short that the effect of
turbulence on the particle accumulation is negligible, and h∆vSi drops to zero.
3.4. Effect of the small scales
In order to study the effect of small eddies, we set kCL = kmin and kCH < kmax.
Thus different ranges of the high wavenumber components (kCH < k < kmax) can be
omitted by setting different values of kCH. Recall that the values of kmaxη for different
cases are listed in table 1. Figures 7 and 8 show the results of the average settling velocity and the particle accumulation for different values of τp/τk, vd/vk and kCHη at
Rλ= 65.3. The agreement between the results using all turbulent scales (kmaxη = 1.71,
all simulated modes) and those using only partial turbulent scales (kCHη = 0.5,
excluding the high wavenumber modes from 0.5 to 1.71) indicates that the high wavenumber modes have essentially no effect on the particle accumulation and the particle’s settling velocity. However, the results with kCHη = 0.2 are substantially less
than the values calculated in the flow field with all the simulated modes. The results for Rλ= 22.6 and 35.7 (not shown here) are similar to those for Rλ= 65.3. The cutoff
wavenumber in the LES by Yeh & Lei (1991) for Rλ ≈ 44 is roughly 0.2/η–0.4/η,
the vorticity structure and hence the phenomenon of particle accumulation were not simulated correctly, and thus the results of the average settling velocity in Yeh & Lei
0 1 2 vd/vk (b) 0.8 3 0.7 0.5 τp/τk= 1 0 1 2 τp/τk (a) 0.7 3 0.6 0.5 vd/vk= 1 0.6 0.9 1.0 〈Ω〉 Ω′ f 0.8 0.9 1.0 〈Ω〉 Ω′ f
Figure 8.Effect of small turbulent scales on the average dimensionless enstrophy evaluated at the locations of the particles at Rλ= 65.3. The symbols are the same as those in figure 7.
were underestimated. Note that the ‘traditional’ LES employed by Yeh & Lei is still applicable for studying the particle dispersion problem.
As in Wang & Maxey (1993), the instantaneous results within a thin slice of the flow field (with thickness equal to the grid spacing) are employed for illustrating the physics. The particle velocity vectors within the volume of the slice are projected onto the centreplane of the slice. The background plots of figures 9(a), 9(c) and 9(e) are the instantaneous dimensionless enstrophy contours at the centreplane. The corresponding distributions of the instantaneous fluid velocity are shown in figures 9(b), 9(d) and 9(f). The maximum values in figures 9(a), 9(c) and 9(e) are 37.12, 18.39 and 2.90, respectively. The phenomenon of particle accumulation in the low vorticity (enstrophy) region can be observed clearly in figure 9(a), i.e. most of the starting points of the particle velocity vectors are at the low vorticity region. The preferential sweeping phenomenon is illustrated through the consistence of the particle velocity vectors in figure 9(a) and the fluid velocity vectors in figure 9(b).
Consider the results using all simulated modes in figures 9(a) and 9(b) and the results for kCHη = 0.6 in figures 9(c) and 9(d). Although the detailed values in the
high enstrophy regions are different, the essential structure of the enstrophy contours in figures 9(a) and 9(c) are similar. Also the velocity fields in figure 9(b) and 9(d) are almost the same. The discrepancy for enstrophy at the high vorticity regions is thought to contribute negligibly to the average settling velocity since particles are rare at the high vorticity regions. Detailed examination of the data shows that the case
(a ) 10 8 6 4 2 0 g (b ) 1.5 1.2 0.9 0.6 0.3 0 (c )( d ) (e )( f ) g (a ) 1.5 1.2 0.9 0.6 0.3 0 10 8 6 4 2 0 1.4 1.0 0.6 0.2 –0.2 –0.6 (b ) (c ) F IGURE 9 ( a– f ) ( left ). F or caption see f acing page. F IGURE 14 ( a– c) ( ab o ve ). F or caption see f acing page.
0 1 2 kCHlω (c) 0.8 3 τp/τk= 1.5 0.6 1.0 Σ 0.4 0.2 4 5 6 0 1 2 kCHlω (d ) 0.8 3 τp/τk= 2.5 0.6 1.0 Σ 0.4 0.2 4 5 6 0 1 2 kCHlω (a) 0.8 3 τp/τk= 0.25 0.6 1.0 Σ 0.4 0.2 4 5 6 0 1 2 kCHlω (b) 0.8 3 τp/τk= 0.75 0.6 1.0 Σ 0.4 0.2 4 5 6
Figure 10.The relative error of the increase of the average settling velocity, Σ, for different values of τp/τk, kCHlω and Rλ at vd/vk= 1. The symbols are the same as those in figure 4.
with kCHη = 0.6 resolves approximately 75% of the total enstrophy and 98% of the
total energy of the case with all simulated modes. It follows that the sites of particle accumulation and the average settling velocity for the case with all simulated modes are basically the same as those with kCHη = 0.6. On the other hand, although the
flow field with kCHη = 0.2 in figure 9(f) still possesses some essential similar features
as in figure 9(b) (resolving 75% of the total energy), the enstrophy distribution for
kCHη = 0.2 (figure 9e) is completely different from that in figure 9(a) (resolving only
25% of the total enstrophy). It follows that the particle accumulation, and thus the average settling velocity for case with kCHη = 0.2 are different from those for the case
with all simulated modes.
Denote the value of h∆vSi for the case with all simulated modes by A, and the
value ofh∆vSi for a given kCH by A∗. Define the relative error Σ as
Σ = A− A
∗
A , (11)
which depends on Rλ, τp/τk, and vd/vk. Figure 10 plots the values of Σ for some of
our typical simulated results using lω−1as the scale for kCH. Recall that lω= 1/kω, with
kω the wavenumber corresponding to the maximum of the dissipation (or vorticity)
spectrum. The results in figure 10 together with other results not shown in the present paper suggest that the relative error is negligible for kCHlω> 2.5. It follows that LES
is still applicable for generating the flow field for the simulation of h∆vSi provided
Figure 9.(a), (c) and (e) The instantaneous particle velocity vectors on contours of dimensionless enstrophy; (b), (d) and (f) the corresponding dimensionless fluid velocity vectors: (a) and (b) are the results using all the simulated modes, (c) and (d) are the results for kCHη = 0.6, and (e) and (f)
are the results for kCHη = 0.2. The maximum values of the dimensionless enstrophy for (a), (c) and
(e) are 37.12, 18.39 and 2.90, respectively. Here Rλ= 65.3.
Figure 14.(a) The dimensionless enstrophy (maximum value = 35.57), and (b) the fluid velocity vectors of the instantaneous field for KCLLf = 2.5. (c) The instantaneous particle velocity vectors
0 1 2 vd/vk (b) 0.6 3 0.4 0.2 τp/τk= 1 0 1 2 τp/τk (a) 0.6 3 0.4 0.2 vd/vk= 1 〈∆vs〉 vk 〈∆vs〉 vk
Figure 11.Comparison ofh∆vSi between the result using DNS (963 grids, Rλ= 65.3, denoted by 4) and that using LES (323grids, R
λ= 69.8, denoted by +).
that the cutoff wavenumber in LES, kC, satisfies the condition kClω > 2.5. A result
using such LES (case D), which will be called the ‘extended LES’ later, is shown in figure 11 together with the corresponding result using DNS (case C). Both results give essentially the same prediction of h∆vSi. The agreement for the accumulation
of particles for both calculations (extended LES and DNS) is also demonstrated in figure 12 via the ensemble average of the enstrophy at the locations of the particles. Note that although the cutoff wavenumber, kC, in the present extended LES is greater
than that employed in the conventional LES, the number of grids employed here (323) is still substantially less than those used in DNS (963). The maximum wavenumber resolved by the extended LES is kmaxη ≈ 0.61 (see table 1), while kmaxη = 1.71 for
DNS. The CPU time for the extended LES is about 10% of that for the DNS. Detailed examination of the data shows that about 95% of the total energy and 75% of the total enstrophy are resolved by the extended LES for the above case. The applicability of the extended LES also allows us to study h∆vSi in a flow field with
larger Rλ which processes a short inertial subrange with 1283 grids, and is important
for studying the velocity scale for h∆vSi in § 3.6.
3.5. Effect of the large scales
The effect of the large scales onh∆vSi can also be studied using the above idea by chopping off the low wavenumber components. Consider the results of the particle calculations in the flow fields with kCH = kmax and kCL > kmin. Figure 13(a) shows
0 1 2 vd/vk (b) 0.8 3 0.7 0.5 τp/τk= 1 0 1 2 τp/τk (a) 0.7 3 0.6 0.5 vd/vk= 1 0.6 0.9 1.0 〈Ω〉 Ω′ f 0.8 0.9 1.0 〈Ω〉 Ω′ f
Figure 12.Comparison of the ensemble average of the dimensionless enstrophy at the locations of the particles between the result using DNS (963 grids, R
λ= 65.3, denoted by4) and that using
LES (323 grids, R
λ= 69.8, denoted by +).
using all the simulated modes, and the result for h∆vSi using the flow field with
kCLLf = 3.25 reduces almost to zero. Consider first the case with kCLLf = 2.5. The
vorticity contours in figure 14(a) (see p. 196) is almost the same as that in figure 9(a) (with all simulated modes). In fact, 97% of the total enstrophy has been resolved for the case using kCLLf = 2.5, and the phenomenon of particle accumulation should
be reasonably estimated. Although figure 13(b) indeed shows that the accumulation of particles in the low vorticity region for kCLLf = 2.5 is similar to that for the
case with all simulated modes, there still exists a certain amount of discrepancy between both results, which implies that the low wavenumber components also have moderate effects on the particle accumulation. Figure 14(b) shows the velocity field for kCLLf = 2.5, which is similar to the velocity field for all simulated modes in figure
9(b). Detailed examination of the data shows that 84% of the energy is resolved for the case with kCLLf = 2.5. In order to understand why we underestimate h∆vSi by
almost 50% with only 16% less energy in the case with kCLLf = 2.5, we plot the
particle velocity vector on the contours for constant values of ∆q/qall in figure 14(c).
Here ∆q = qall − q2.5, which represents the energy contributed by the large scales of turbulence with wavenumber less than 2.5/Lf. Here qall and q2.5 are the energy (uiui/2) based on the calculation with all simulated modes and that with kCLLf = 2.5,
respectively. The result in figure 14(c) indicates that in general, more particles locate instantaneously in the high energy level as expected. According to the results in
0 1 2 (b) 1.0 3 0.9 0.8 0 1 2 τp/τk 〈∆ vs 〉/ vk (a) 0.6 3 0.4 0.2 τp/τk 0 –0.2 0.7 0.6 0.5 〈Ω 〉/Ω ′ f
Figure 13.Effect of large turbulent scales on (a)h∆vSi/vk, and (b)hΩi/Ω0
f when Rλ= 65.3 and
vd/vk= 1.
◦
, all simulated modes;, kCLLf = 2.5; and4, kCLLf = 3.25.figures 13 and 14, the energy of the large eddies should play a significant role on the magnitudes of the increase of the average settling velocity.
3.6. The scales for the increase of the average settling velocity
Although Wang & Maxey (1993) have demonstrated the importance of Kolmogorov scaling for the processes of particle accumulation and the average settling velocity, the velocity scale for h∆vSi deserves further study. The results in §§ 3.4 and 3.5 show
that although the particle accumulation is controlled mainly by the small scales with maximum vorticity, the average settling rate, h∆vSi, depends strongly on the large
energetic scales. The variations ofh∆vSi/vkwith τp/τk (or vd/vk) and Rλat fixed value
of vd/vk (or τp/τk) in figure 4 show that the data spread widely apart for different
Rλ. Also shown in figure 4(b) is that the peak value of the data for a given Rλ shifts
slightly rightward as Rλ increases. In order to study more the details of the effect of
Rλ onh∆vSi, two flow fields with higher Rλ (cases E and F in table 1) were generated
by the extended LES, and employed for studying the particle motion. The results in figures 15(a) and 15(c) show that the normalized data with vk (the Kolmogorov
velocity scale) as the velocity scale spread even further apart for a larger range of Rλ.
However, if the velocity scale for the large energetic eddies, u0, is employed instead of
vk, the results in figures 15(b) and 15(d) indicate that the data for different Rλ tend
to collapse together. Furthermore, although the location of the peak value of h∆vSi
0 1 2 (c) 3 0.6 0.4 〈∆ vs 〉/ vk τp/τk 0.2 vd/vk= 1 0 1 2 (d ) 3 0.20 0.15 τp/τk 0.10 vd/vk= 1 0.05 0 2.0 4.0 (a) 6.0 1.0 0.8 〈∆ vs 〉/ vk vd/vk 0.2 τp/τk= 1 0 0.5 1.0 (b) 1.5 0.2 〈∆ vs 〉/ u ′ vd/u′ 0.1 τp/τk= 1 〈∆ vs 〉/ u ′ 0.4 0.6 8.0
Figure 15.Variations ofh∆vSi with different parameters:
◦
, Rλ= 22.6;, Rλ= 35.7;4, Rλ= 65.3; , Rλ= 133; and∗, Rλ= 153. Here vk is employed as the velocity scale for (a) and (c), but u0 is thevelocity scale for (b) and (d).
to vd/vk ≈ 3.2 when Rλ = 133, the locations of the peak values for different Rλ
in figure 15(b) stay essentially at vd/u0 ≈ 0.5. Note that the vertical scale in figure
15(b) is five times greater than that in figure 15(a). Also the statistical uncertainties are quite substantial for the increase of the average settling velocity of particles (see figure 1 of Maxey 1987 for example). The small amount of scattering of the results in figure 15(b) is consistent with the statistical uncertainties (errors) of the calculations. A more elaborate study is carried out through a series of comprehensive calculations for the whole ranges of parameters. Figures 16(a), 16(c) and 16(e) show contours for constant values of h∆vSi/vk in (vd/vk, τp/τk)-planes for three different
Rλ, which employ vk as the velocity scale. The corresponding results using u0 as the
velocity scale are shown in figures 16(b), 16(d) and 16(f). The results using vk as
the velocity scale show that both the magnitude and the vertical location of the maximum value ofh∆vSi/vk increase as Rλ increases. However, the results using u0as
the velocity scale show that the maximum values ofh∆vSi/u0have essentially the same
magnitude (h∆vSi/u0 = O(10−1)) and occur at almost the same location (τp/τk ≈ 1
and vd/u0 ≈ 0.5). Thus the more appropriate velocity scale for h∆vSi seems to be u0,
the velocity scale of the large energetic eddies. There are larger discrepancies among the data in figure 15(d) for vd/vk = 1 in comparison with those in figure 15(b) for
τp/τk= 1. This is because the results for different Rλ in figure 15(d) correspond to the
data on different horizontal straight lines with different values of vd/u0in figures 16(b),
16(d) and 16(f). Here vd/u0 = 0.415, 0.331, 0.244 and 0.171 for Rλ = 22.6, 35.7, 65.3,
and 133, respectively, in figure 15(d). If the results for a fixed value of vd/u0 (say,
vd/u0= 0.5 or 1, not included in this paper) are plotted instead of those for a fixed
value of vd/vk as in figure 15(d), the data collapse nicely.
0 1 2 (e) 8 6 τp/τk 4 2 vd vk 0 1 2 ( f ) 1.0 τp/τk 0.5 vd u′ 0 1 2 (c) 5 3 2 1 vd vk 0 1 2 (d ) 1.0 0.5 vd u′ 0 1 2 (a) 4 3 2 1 vd vk 0 1 2 (b) 1.0 0.5 vd u′ 4
Figure 16.(a), (c) and (e) Contours for constant values ofh∆vSi/vk. (b), (d) and (f) Contours for constant values ofh∆vSi/u0. (a) and (b) are the results for Rλ = 35.7, (c) and (d) are the results for
Rλ= 65.3, and (e) and (f) are the results for Rλ= 133.
time scale for representing h∆vSi may be understood as follows. According to the
description of the third paragraph of § 1, there are two essential elements in the increase of the average settling velocity: the particle accumulation and the drag on the particles. The accumulation of particles in the downward sweeping, low vorticity region where particles experience less local drag is a necessary condition for the increase of the average settling velocity. The increase of the average settling velocity is indeed proportional to the degree of particle accumulation since more samples are settling faster during the averaging process as discussed in § 3.3. However, the magnitude of the increase of the settling velocity of an individual particle depends on the drag acting on the particle. Thus the drag on the particles is also a crucial factor in the increase of the average settling velocity of particles.
The particle accumulation is due to the centrifugal effect associated with the intense vortical structures, and is thus related to the vorticity field of the flow. The scales near the peak of the dissipation (vorticity) spectrum contribute most to the vorticity. The
spectrum value at the peak in figure 1(b) is two orders greater than that at kη = 1, and thus the Kolmogorov eddies can hardly play a significant role in the phenomenon of particle accumulation.
The drag on a particle is proportional to the relative velocity of the particle and the local flow field. Consider the case for small particles such that the Stokes drag law applies. If ∆v = v− vd, (2a) may be written as
d∆v dt =
u− ∆v
τp
. (12)
As the particles tend to accumulate in the low vorticity regions where the fluid velocity has scale u0 (see figure 9), it follows from (12) that u0 is an appropriate scale for representing h∆vi. More rigorously, the local fluid velocity in the wavenumber domain is expressed in terms of a finite Fourier series in the DNS (or LES) using a pseudospectral method, with each term representing the contribution from a given turbulent scale. The relative importance of each term depends on its coefficient, which is related directly to the energy spectrum. The energy of the large energetic scales is three to four orders greater than that of the Kolmogorov scales (see figure 1a), and thus the large energetic scales contribute mainly to the local fluid velocity u, or the drag on the particle. The Kolmogorov scales play essentially no role on the drag and thus the increase of the average settling velocity of particles.
Thus the increase of the average settling velocity of particles in homogeneous isotropic turbulence is a physical problem involving quite a wide range of scales rang-ing from the large energetic scales to the small scales approximately with wavenumber 2.5kω. The small scales are primarily responsible for ‘moving’ the particles to the
re-gions where the local drag on a particle is less than the average value. The reduction of the drag on a given particle depends on the local fluid velocity, which is controlled mainly by the large energetic scales. Both the small scales with maximum vorticity and the large scales with maximum energy are crucial for the increase of the average settling velocity of particles. The Kolmogorov scales play essentially no dynamical role on the problem because both their vorticity and energy are negligible. The separation between the wavenumber corresponding to the peak of the dissipation (vorticity) spectrum and that of the energy spectrum increases as Rλ increases. When Rλis
suffi-ciently large, the large scales have an insignificant effect on the particle accumulation, but the drag on every individual particle still depends strongly on the large energetic scales. Thus there are certain physical reasons for supporting the choice of u0 as the velocity scale for representing the increase of the average settling velocity.
It is interesting that the Kolmogorov time scale, τk appears as an approximate
time scale for the present problem. Note that the results in figures 15 and 16 are maximized at τp/τk≈ 1. As the small scales of turbulences with wavenumber greater
than 2.5kω are found to have a negligible effect on h∆vSi in § 3.4, the entry of τk is
thus considered not to be related directly to the Kolmogorov eddies. A plausible way that τk might play a significant role in the problem is through the phenomenon of
particle accumulation. Since the scale for the vorticity associated with the ‘eddies’ with wavenumber kω is of order u0/λ (see Tennekes & Lumley 1972), which equals
(15)−0.5τ−1k , and the scales near kωare responsible mainly for the particle accumulation,
τk thus appears as the appropriate time scale of the present problem.
The phenomenon of the increase of the average settling velocity due to the turbu-lence occurs for limited ranges of vd/u0and τp/τk according to the above result. If this
result can be extended to cases with larger values of Rλ, we may examine its effect
related to the particle’s relaxation time τp for a given body force g. Consider first the
gravitational settling of sand particles in the atmospheric boundary layer. On taking
ρp/ρ≈ 1700, ν ≈ 1.5 × 10−5 m2s−1, u0 ≈ 0.5 m s−1 and l ≈ 200 m (see Lugt 1983),
where l is the size of the large eddies, we found ≈ u03/l≈ 6.25×10−4m2s−3, η≈ 1.52 mm, τk ≈ 0.155 s and Rλ≈ 104. If vd is chosen to be 0.5u0, we found dp≈ 68 µm and
τp ≈ 0.0292 s. It follows that τp/τk ≈ 0.188, which is much less than the unity required
for the increase of the average settling velocity to be significant (see figure 16f). However, the phenomenon could be pronounced for some engineering applications. Consider for example the gravitational settling of oil drops in the combustion cham-ber of a power plant. Due to the high temperature environment, we set ρp/ρ≈ 2000
and ν ≈ 10−4 m2s−1. On taking u0≈ 1 m s−1 and l≈ 10 m, we found ≈ 0.1 m2 s−3,
η ≈ 1.78 mm, τk≈ 0.0316 s and Rλ ≈ 1.2 × 103. If vd is chosen to be 0.5u0, we found
dp ≈ 229 µm and τp ≈ 0.0583 s. It follows that τp/τk ≈ 1.85, which is near unity for
h∆vSi to be maximized. For some smaller combustion chambers, u0 ≈ 0.5 m s−1 and
l ≈ 1 m, we found Rλ ≈ 270 and τp/τk ≈ 0.96 if vd/u0 ≈ 0.5, which falls into the
parameter ranges where the increase of the average settling velocity is significant.
4. Conclusions
When heavy solid particles settle under a body force field in homogeneous isotropic turbulence, the particles tend to accumulate in the downward sweeping, low vorticity regions (where the local drag on a particle is less than the average value) due to the local centrifugal effect associated with the intense tubular vortical structures. The accumulation effect is necessary for the increase of the average settling velocity, and is controlled mainly by the small scales of order lω, which equals the inverse of the
wavenumber corresponding to the maximum of the dissipation (vorticity) spectrum. However, the increase of the average settling rate depends also on the drag on an individual particle, which is controlled mainly by the large energetic scales. The Kolmogorov scales, which are in general of one order less than lω, have essentially
no effect on either the particle accumulation (because of their relative weak vorticity) or the drag on the particles (because of their lack of energy), and thus contribute no role in the increase of the average settling rate of particles. The flow field generated by the large eddy simulation is adequate for studying the present problem, provided the cutoff wavenumber in the simulation is greater than 2.5/lω. This is an interesting
example for turbulent transport where both the large and small scales (but not the smallest (Kolmogorov) scales) of turbulence are important for the phenomenon.
Through a series of calculations for different parameters, the increase of the average settling rate in stationary turbulence is found to be maximized and of order u0/10
when vd≈ 0.5u0 and τp/τk ≈ 1 for Rλ = 22.6–153, which implies that the appropriate
velocity and time scales for the present problem are u0 and τk, respectively. Here u0
is the velocity scale associated with the large energetic eddies, which are primarily responsible for the magnitude of h∆vSi, and τ−1k is the scale for the vorticity of the
‘eddies’ with size lω, which are associated with the accumulation of the particles in
the low vorticity region. The settling of particles in homogeneous isotropic decaying turbulence was also studied by the present authors (see Yang 1997), and the behaviour of the average settling velocity in decaying turbulence is qualitatively similar to that in stationary turbulence.
The present study is supported partially by the National Science Council of the Republic of China through the grant NSC86-2212-E-002-074. The authors are also
