An approach to modeling particle motion in turbulent flows—I. Homogeneous, isotropic turbulence

Atmospheric EnvironmentVol. 29, No. 3. pp. 423-436. 1995 ElsevicrScience Ltd Printed in Great Britain

Pergamen

1352-2310(94)00269-x

AN A.PPROACH TO MODELING PARTICLE MOTION IN TURBULENT FLOWS-I. HOMOGENEOUS, ISOTROPIC TURBULENCE Q. Q. LU Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ08544, U.S.A. (First received 15 April 1994 and in final form 28 July

1994)

Abstract-A model based on time-series analysis concepts is developed for generating the velocity field at particle positions using Eulerian temporal and spatial correlations functions. The particles are then moved

in Lagrangian frame in response to the velocity field and body forces. The model was examined by simulating particle motion in homogeneous, isotropic turbulence and two grid-generated turbulent flows of Snyder and Lumley [J. Fluid Me&. 48,41 (1971)], and Wells and Stock [J. Fluid Mech. 136, 31 (1983)]. The predicted asymptotic particle longitudinal and transverse dispersion coefficients are shown to agree with the formula of Csanady [J. atmos. Sci. 20, 201 (1963)]. Comparison is also made for the computed particle fluctuating velocity with the experimental data of Snyder and Lumley and of Wells and Stock. Fairly good agreement is observed. In addition, particle longitudinal dispersion, velocity correlation and corresponding integral time scales are calculated both in presence and absence of particle drift velocity so as to investigate the crossing-trajectory, continuity and particle inertia effects. The trends are in accordance with previous theoretical studies of Reeks [J. Fluid Mech. 83,529 (1977)], Pismen and Nir [J. Fluid Mech. 84, 193 (1978)] and Nir and Pismen [J. Fluid Mech. 94, 369 (1979)]. Key word index: Particle dispersion, turbulence, Eulerian temporal correlation, Eulerian spatial correlation, time series, two-phase flow.

1. INTFtODUCTION

The transportation and dispersion of particles by turbulent flows, is of central importance in a number of industrial processes and environmental problems. Consequently, a great deal of experimental, theoretical and computational work has been done on the subject and a number of discoveries have been made recently. For example, experiments in laboratories (Rashidi et al., 1990) have demonstrated that, for certain ranges of particle time constants, heavierthan-fluid particles tend to segregate forming streaky configurations in the wall region of turbulent flows. Direct numerical simulations of particle behavior in wrill turbulence by Pedinotti et al. (1992) are in agreement with these observations. This type of behavior has also been seen by several other investigators, notably by Squires and Eaton (1990) for homogeneous, isotropic turbulence. Such segregation phenomena may lead to turbulence modulation by particles, even at rather dilute concentrations, and are therefore an interesting recent finding. Considerable activity using direct numerical simulation (DNS) is now going on in this area. Direct numerical simulations (DNS) are, however, very demanding of computational resources and still impractical for most problems except at relatively low

Reynolds numbers. While large eddy simulations (LES) offer an attractive route for the future, the methodology still has not been exploited extensively for situations of practical interest and there are developments necessary for subgrid models that take into account particle motion. The need, therefore, still exists for calculatioli of particle behavior in turbulent flows using less computationally intensive methods. Such methods are constructed along two main lines (see, for example, the review by Banerjee 1990). (1) The continuous and dispersed phases are considered in an Eulerian frame. This is often used when particle concentration is high and volume or ensemble averaging is done to develop field equation for interpenetrating continua. (2) The motion of the continuous phase is considered in an Eulerian frame, whereas the dispersed phase is calculated using a Lagrangian approach, i.e. the particles are moved by the forces acting on them, including those due to the interactions with the continuous phase, and the appropriate forces are fed back in calculating the continuous phase motion. The latter approach was used successfully for certain spray combustion problems by Williams (1985), and developed further by Gosman and Ioannides (1981), Shuen et al. (1983, 1985), Ormancey and Martinon (1984), Berlemont et al. (1990) and Lu et al. 423

AE 29-3-I

424

Q. Q. LU

(1992, 1993a-c) among many others. However, calculations of the dispersion of particles in turbulent flows, even at very low concentration, is still fraught with difficulty. The crux of the problem lies in that most methods for the calculation of continuous phase motion rely on the Reynolds averaged equations of motion (the alternative being DNS or LES which are still not feasible for most problems, as discussed earlier). This only yields estimates of the mean velocity field, and at best, the Reynolds stresses. It is well known that even these estimates may be inaccurate because several constants arise in the closure models. These constants may well be influenced by the physical scales and geometry of a problem since averaging is done over all turbulence scales. Be that as it may, a fluctuating velocity field has to be constructed from this type of averaged information and then used to move the particles. As we know when a particle moves around in a turbulent flow it encounters a series of eddies along its trajectory and interacts with them. To compute the instantaneous position of the particle by solving the particle motion equation we must devise a method to construct the information reflecting the particle-eddy interaction along the particle path. Namely, we need to know the instantaneous velocity of fluid at the particle position. Knowing the instantaneous position and velocity of particles all the statistic quantities of particles can be worked out. Along these lines, Gosman and Ioannides (1981) proposed to construct the fluctuating velocity field by sampling from a Gaussian probability density function whose variance is proportional to the local turbulent kinetic energy. There are some difficulties with the model in that differences in particle longitudinal and transverse long-time diffusion coefficients (see Csanady, 1963; Reeks, 1977; Nir and Pismen, 1979) cannot be predicted. In homogeneous, isotropic and incompressible turbulence, for example, this effect arises from the well-known difference in the longitudinal and transverse spatial correlation functions for the velocity fluctuations-sometimes called the “continuity” effect. Several improvements in such approaches have been suggested by Shuen et al. (1983), Kallio and Reeks (1989) and Burnage and Moon (1990) amongst others. Though different, these methods share the same underlying idea. Essentially, at the start of the particle-eddy interaction, the velocity of an eddy is found by making a random selection from the probability density function (p.d.f.) of velocity (taken to be an isotropic Gaussian p.d.f. with standard deviation (2k/3)“‘, where k is the turbulent kinetic energy). A particle is assumed to interact with an eddy for a time which is the minimum of the eddy lifetime t, and the time At, taken for a particle to pass through the eddy. The eddy lifetime is estimated to be t, = L/(2k/3)“2, assuming that the eddy size is L=clk3i2/E, where ci is a constant and E is the mean rate of the turbulent kinetic energy dissipation. The particles are assumed to interact with an eddy as long

as both the time and the relative velocity distance of interaction, Ax,,, satisfy the criteria: At,,< t, and 1Ax,1
425

Modeling particle motion--I

needs the information of the fluid Lagrangian correlation. Any constants that arise in the formulation are determined by studying the behavior of fluid particles, as one limit of the particle dispersion problem. The results of the proposed procedure are then tested against a number of experiments in which the particles have densities different from the fluid, without adjusting any of the constants that we determine using fluid particle behavior as a limit.

2. METHOD As we will see shortly, the method is based on calculation of particle motion in a fluid velocity field afresh at very time step. First we discuss particle motion, and then the generation of the fluid velocity. 2.1. Particle motion In the present study, only the drag and gravity terms will be included in the equation of particle motion. The Basset, added mass, lift and pressure gradient terms are omitted. This is justified when the ratio of the particle density to the fluid density approaches 1000 (see Boothroyd, 1971). This requirement is generally satisfied in gaseous systems. However, all except the Basset terms can be readily incorporated into the present model, so we are by no means limited to equations (1) and (2) below. In the present study the following equations are used: dV ppdt=

3 -4dprC,(V-U)IV-UI

+(Pp-Pfk

(1)

P

(2) where pp and pf are the particle and fluid densities, respectively, V and U are the instantaneous velocities of particle and fluid, respectively, d, is the particle diameter and g is the gravitational acceleration. The drag coefficient CD is given by C,=$-(l+C,.L5(Re)~68’)

2.2. Model Consider, without loss of generality, that the mean flow considered has one principal direction, as in several previous studies (Snyder and Lumley, 1971; Reeks; Pismen and Nir, 1979; Nir and Pismen; Wells and Stock, 1983). Establish a coordinate system moving with the mean velocity, O-X YZ. Let the X-axis coincide with the mean flow direction. At the instant t, a particle is at the position denoted by the position vector X0.. After one time step of computation, say At, it moves to a position denoted by the position vector X, (see Fig. 1 for a schematic of the coordinates, and the various parameters associated with the problem). Let us set up a local coordinate system O’-@‘Z’Q whose @‘-axis passes through the points denoted by Xv and X,, while the directions of the other two axes are orthogonal but otherwise arbitrary. In many situations, such as decaying grid turbulence, the flow fields are not entirely homogeneous, isotropic and stationary; therefore the fluctuating velocity u, may be normalized by the square root of its local variance in the expectation that the normalized fluctuating velocity field approximates a homogeneous, isotropic and stationary situation. The normalized fluctuating component in the i-direction is denoted by Wi, i.e. Wi = ui/fi

(i, = l,2,3),

where the subscripts 1, 2 and 3 represent, respectively, the directions of the relative O’, 2’ and R’ axes. The quantities with overbars indicate the ensemble average values. This artifice is later useful in dealing with some of the experiments examined here. It is clear from Fig. 1 that when the particles arrive at the position X,, the normalized fluid fluctuating velocity at the position X,. will change from Wi(X,.), to Wi(Xo,),+,,. Here the subscripts t and t + At indicate the time. If we limit ourselves to the isotropic case the normalized fluctuating velocities of the fluid at the instants t and t+ At and at the positions Xo. and X, are related by the following correlations:

= Wi(Xo,)*Wi(X,.), Fii(At)

(i = 1,2,3),

(3)

for (Re),<200,

ep

where (Re),, the particle Reynolds number, is defined as (ee),=p.

the problem reduces to estimating the fluctuating fluid velocity at the particle locations from these quantities.

IU-VII, ”

Furthermore, we do not consider the forces, due either to particle rotation (Rubsinow and Keller, 1961) or to fluid velocity gradient (Saffrnan, 1965breasonable for the problems considered here but not a fundamental limitation of the procedure. The effect of fluid turbulence on particles is accounted for but the modification of the turbulent field due to the presence is ignored as mentioned earlier. Therefore, all the mean data of the fluid field, including the mean velocity, the mean turbulent kinetic energy and the mean turbulent kinetic energy dissipation rate, are regarded as being known. They could be obtained by measurements or by use of some turbulence model for the continuous phase. To obtain the statistical properties of the particle motion, representative particles must be followed along their trajectory by integrating equations (1) and (2). To do this, it is necessary to know the instantaneous velocity of the fluid at the positions where particles are located. Since the fluid mean velocity and some average turbulence parameters are the only things known,

= W,(X,.),+,, Wi(XO’)t+&i(A~l)

(i= L&3).

(4)

Fii(At) and Gii(As,) in relations (3) and (4) are the Eulerian temporal and spatial velocity correlation functions in the coordinate system O’-@‘Zn’, respectively. In this paper repeated index does not mean summation. Applying the idea of time series (Box and Jenkins, 1976) to W,(Xc), and Wi(Xv)r+ar and Wi(X~)~+~, and Wi(X,), re-

Fig. 1. The coordinate

system established model.

for the

426

Q. Q. LU

spectively, leads to the following two relations: Wi(Xo.),+ar=aiWi(X~)t+ai w,(x,)=b,w,(x,.),+,,+pi

(i=1,2,3),

(5)

(i=1,2,3).

(6)

Here ni and bi are two constants to be decided and ai and pi are two mean zero random variables with the correlation properties: ai(t)aj(t)=O

(i#j),

(7)

Bi(t)pj(t)=o

(i Zj).

(8)

In the later reasoning we shall put a constraint on them. To simplify the matter we assume that WJX,.), and ai are statistically independent. Of course this assumption can be removed and certain correlation is adopted for Wi(Xo,), and cti. This may be a research direction in the future. Now we start to find out the parameters involved in equations (5) and (6). Multiplying equation (5) by Wi(Xo.), and then making ensemble average gives ni = Fii(At).

(9)

bi = G,,(AsJ.

(10)

G,,(As,)=exp(-As,/A,),

(13)

Gzl(As,)=G3~!As,)=exp(-AsllA2),

(14)

where ~~ is the Eulerian integral time scale and l\l and A2 are the longitudinal and transverse Eulerian integral length scales, respectively. For isotropic, homogeneous and incompressible turbulence A 1= 2A2. At high Reynolds numbers both Fii(At) and G,,(As,) can be represented over some region of At and As, as exponential functions (see Hinze, 1975; Sato and Yamamoto, 1987). This representation clearly breaks down as At and As, -0 since there is a point of inflection near the vertex. However, for At > n/u’and As 1> 1,where u’is the turbulence intensity and I the Taylor microscope, the exponential representation is reasonable. It is known that the transverse spatial correlation function has a negative loop due to the incompressibility of the fluid field, which suggests the invalidity of an exponential correlation form for the transverse spatial correlation function. However, in the following we are able to show that the present method correctly predicts such quantities as the long-time diffusion coefficient for particles with significant inertia. Consider the simplified form of equation (I):

Similarly we have

Substituting equation (5) into equation (6) to eliminate Wi(X,.),+,, and rewriting the equation we obtain WJX,) = aibi Wi(XW), + yi (i = 1,2,3),

(11)

where yi = biari+ pi. Now we would like to make an additional assumption for yi to limit its property. We assume ui(t) to be the incremental Wiener process. The reason for doing so is simply that this process is completely determined by its variance and no more higher-order moments are necessitated. From equation (11) we can see that the higher-order moments of yi are determined by the same-order moment of Wi. When the order is higher than 2 the data for the moment of Wi are quite often unavailable in experiments against which the present model will be tested. Therefore, we content ourselves to the second-order moment. In order to define the process yi we need to know its variance. Squaring equation (11) and making ensemble average yields the variance of y,: (OJi=Jm

(i=1,2,3).

dV ;=

-K(V-U)+g.

(15)

Here K is a constant, the reciprocal of the particle time constant defined later in equation (31). Taking ensemble average equation (15) gives the stationary solution V,= v=g/~, the particle terminal velocity. Now we discuss the case that the particle time constant is so significant that it does not apparently respond to the surrounding fluid fluctuation. After simple manipulation we have the following relation for the long-time particle diffusion coefficients Di(co): 1 dX;(t) Di(co) = lim - t-m2 dt

exP(-At/r,-As,/AJd(At).

(16)

0

In the presence of gravity, As 1= At V,, where V, is the magnitude of V,. Then equation (16) yields

(12)

Equation (11) is the relation we are seeking, which connects the instantaneous velocities of fluid at two consecutive particle positions. With the aid of this equation we can solve the particle motion equation for an individual particle. Here we would like to compare equation (11) with the Langevin equation for the fluid diffusion in turbulence. The Langevin equation is a relation describing the evolution of the instantaneous velocity of the particle along its trajectory and it includes the Lagrangian time scale. In the Langevin equation there is a constant coefficient in the place of aibi of equation (11). As for equation (11), it relates the fluid instantaneous velocities at two different instants and at different spatial positions. The positions can be arbitrary in space, not necessarily on the particle path. When the two spatial positions are all on the particle path, equation (11) is also a relation describing the evolution of the instantaneous particle velocity along its trajectory. But, in this instance, As, will be random along the particle trajectory and so do aibi, a fact which distinguishes equation (11) from the Langevin equation. In addition, equation (11) always contains the Eulerian temporal and spatial correlation, never the Lagrangian correlation. It is well known that equations (5) and (6) exclusively lead to exponential functions for Fii(At) and Gii(Asl). Namely, the following forms are implied for FJAt) and G,,(As,) in the present study:

D,(m)=fin~

1 dX;(t)

Zdt

(17) Obviously, in this case direction i= 1 coincides with the gravity direction. Equation (17) implies that the long-time particle diffusion coefficient is inversely proportional to the particle terminal velocity and the transverse diffusion coefficient is half the longitudinal one (the one in the gravity direction). When V, = 0, i.e. in the absence of gravity, equation (17) gives Di( xc) = ZT,. We also have the following expression for the particle fluctuating velocity correlation: +m

s

ui(al)u&m,=K7

exp(-xAt)exp

0

x (-AtIrE-As,/AJd(At).

(18)

In the presence of gravity, substituting As, = At V, gives 1

Ui(co)Ui(cO) =J K +

I/r,

+

V,/Ai

-A. ZKt4--1.

K

(19)

Results given by equations (17) and (19) are in agreement with the ones of the previous studies (see Reeks, 1977; Pismen and Nir, 1978). The above discussion indicates that

427

Modeling particle motion-I when the particle inertia is significant enough the present model gives a prediction consistent with the past work. It is noted that the relative coordinate system O’-0’E’R’ employed in Fig. 1 is local and it will change continuously in the course of the computation. In principle, the present modeling strategy is applicable to more general cases. But in these cases much more information is needed about the Eulerian temporal and spatial correlation. 2.3. The calculation procedure (1) At t=O, the particle position and velocity are given. The initial fluctuating fluid velocity components ui(i= 1,2,3) at the particle position are obtained from Gaussian variables satisfying p.d.f. with local variances 2, as in Parthasarathy and Faeth (1990). The instantaneous fluid velocity can then be obtained by adding the known mean velocity and the fluctuating part. (2) Using the given initial particle velocity and the instantaneous fluid velocity just obtained, through equations (1) and (2), the position ‘of the particle at the instant At, X,, is obtained by the Ruuge-Kutta method. Displacement AsI can then be calculated (see Fig. 1). (3) The relative coordinate system O’-O’BQ’ is established. The random terms yi having the standard deviations given by equation (12) are generated as in Vetterling et al. (1988). The fluctuating velocity of the fluid point at X, (the new particle position computed at step 2) can then be found from equation (11). (4) Let X, be the starting point, i.e. X0,, for the next time step and repeat steps 2 and 3 till the completion of the computation.

:I. VERIFICATION

3.1. Correlation functions and associated parameters In the present paper, the Eulerian time and length scales, zE and Ai, are estimated from ?L=ClQ.

(20)

h3=A2=0.5A1,

(21)

GA1

The Lagrangian from

(22) tL=p’ integral time scale zL can be found ZL=CJ/E,

(23)

where ~=~(~+~~+~) and E is the turbulent kinetic energy dissipation rate. The coefficients C1, Cz and C3 must now be determined. In this study let C2 = 0.36. This is the value obtained theoretically by Philip (1967) by using Corrsin’s conjecture expressing the Lagrangian correlation tensor as the average of the Eulerian correlation function taken with respect to the uncertain positions of the fluid point. A value, 0.359, for C2 is inferred from Reek’s study. This value falls within the domain 0.15-0.6 discussed by Hinze (1975, p. 424) in detail. C1 is also equal to 0.36, also given by Philip. Such a value is supported by the experiment of Sato and Yamamoto where it was reported that C1 varied from 0.3 to 0.6. In a.later study Kallio and Reeks used 0.4 for C1. Turning to C3, it is determined according to the relation C3 =:0.588C2 (Hinze, 1975, pp. 398 and

426). Therefore, &=0.212 which is near the values proposed by Hinze, 0.235 and 0.25. Relation (21) is a theoretical result for isotropic, incompressible turbulence. In the cases of isotropic or weakly anisotropic turbulent flows, the following relations hold for the fluctuating fluid velocity in the absolute coordinate system O-X YZ and in the relative coordinate system O’_@‘SQ’: ~=I:~+rn:~+n:~,

(24)

1--u,=l:u~+m~ujT+n:u~,

(25)

~=l$~+rn:~+n2;

3

I9

(26)

where lI,m,,nI; 12,m2,n2; l,,m,,n, are the direction cosines of the O’, E’ and a’ axes relative to the x, y, z axes, respectively. 3.2. Particle dispersion

At first we have performed a computation for particle dispersion in stationary, homogeneous, isotropic and incompressible turbulence. This is an idealized case. However, it can help us to understand some mechanisms related to particle dispersion in turbulent flows. In our study, this case will also be used to calibrate the values of the coefficients C1, C2 and C3 to see whether the computed asymptotic diffusion coefficient matches the theoretical result of Taylor (1921) in the limiting case of fluid dispersion. Another purpose of this simulation is to examine the lower and upper bounds that must be imposed on the computation time step At. In this study, the mean turbulent flow conditions are c=655(cms-‘),

TT;=O,

z;?=z=$=z=

196(cm2s-2),

rL=32(ms)

(i=1,2,3),

c=O,

(27) (28) (29)

with air at atmospheric conditions being the fluid. These mean turbulent quantities correspond to the values taken at the section x/M = 68 in the experiment of Snyder and Lumley, which will be taken up again later. The particle diffusion coefficients are defined as the time derivatives of the corresponding meansquare displacements xz(tX i.e. (30) We denote the longitudinal (parallel to gravity) and transverse (normal to gravity) diffusion coefficients by D, and D,, respectively. It is known from the study of Taylor (1921) that at long times Di(t) will become constant. In this part of the simulation, two major cases will be distinguished depending on whether the crossingtrajectory effect is present or not. As is known, particle dispersion in turbulence is influenced by particle inertia and gravity. From the study of Csanady, it is

428

Q. Q. LU

known that gravity, because of the continuity effect, leads to differences in the particle diffusion coefficient in the directions normal and parallel to gravity. At first we intend to examine the variation of At on the computed long-time particle diffusion coefficient and therefore the lower and upper bounds on At. In the computation all the statistics are for 6000 particles. Particles were emitted at x =0 with the fluid local instantaneous velocity, and the statistics were computed beginning from x = 900 (cm) so as to eliminate the affect of the initial conditions. Computations were performed for the dispersion of fluid and particles whose time constants, Tp is defined as

T&$. P

Figures 2a and b show the computed long-time particle diffusion coefficients by using At= 1,5,10 ms in the presence and absence of gravity, respectively. As the reference we calculate the fluid diffusion coefficient by

using Taylor’s result for the fluid long-time diffusion coefficient: zL;;?=6.272 (cm2s-I). For D,, and DI in the presence of gravity Csanady derived D,,=D’[l+c:V,z/Z]D,=D’[1+4C;V,/;;z]-

l/2 , Ii-7.

(32) (33)

The two formulae give 5.6 and 4.41 cm2 s- ’ for the Tp= 20 ms particles and 2.3 1 and 1.22 cm s- ’ for the Tp= 100 ms particles, respectively. The numerical results and these theoretical predictions are obviously in good agreement. In this sense we think the present set of coefficients Cl, C2 and C, to be acceptable and they will be used in the following simulations. At the section x/M =68 (the location where we chose the present mean turbulence quantities (27)(29)) in the experiment of Snyder and Lumley, the Kolmogorov time scale 7K and Taylor dissipation time scale are 9.9 and 38.2 ms, respectively. Extensive simulations indicate that following criterion is a conservative one for the calculated results to be independent of At: (34)

Fluid 0 Tp=20 (ms): transverse + Tp=20 (ms); longitudinal 0 Tp=lO+l (ms); Iransverse X Tp=lOO (ms): longitudinal A Tp=20 (ms); transverse (Csanady) ----. Tp=20 (ms); longitudinal (Csanady) ..... Tp=lOO (ms); transverse (Csanady) ....

Tp=lOO (ms); longirudinal (Csanady) -.- -

I

I

I

0.015

0.03

0.045

Timestep (s)

(b) Fluid

0

Tp=20 (ms) + Tp=lOO(ms)

0

a

I

0.015

0.03 Tiiestep

I 0.045

(s)

Fig. 2. Variation of particle diffusion coefficient as the time step At in the (a) presence and (b) absence of gravity.

In order to examine the overall effects of crossing-trajectory, continuity and inertia on particle behavior in turbulence, Figs 3 and 4 present the computed particle mean-square dispersion and the Lagrangian velocity autocorrelation when gravity is present. Table 1 gives the predicted particle diffusion coefficient, D,velocity fluctuation, u’, and the Lagrangian integral time scale, zL, of the particle velocity correlation. These results were all obtained using At = 5 ms. It is clear that the crossing-trajectory effect significantly reduced the particle dispersion rate and made the particles disperse faster in the direction of gravity. The present study does not manifest strong dependence of the correlations on the direction of gravity. It is observed that gravity exerts unequal influence on particle fluctuation and velocity correlation and the integral time scales, in the directions normal and parallel to it. These trends are consistent with previous studies of Reek (1977). For the Tp= 100 particles, (u~)~/(u;)’ = 1.51, approaching the theoretical asymptotic value of 2 for Tp being large enough (see Reeks, 1977). Figures 5 and 6 are the computed particle dispersion and velocity correlation in the absence of gravity. In this case, particle motion is statistically isotropic, so the dispersion and correlation are depicted only in one direction. Here we can see that in the absence of gravity, particles disperse faster than the fluid and the particle diffusion coefficient increases constantly with their inertia. These results were also found earlier by Reeks and by Nir and Pismen. Table 2 lists the computed particle diffusion coefficient, velocity fluctuation and integral time scales. In Section 2.2 we have shown that the results computed from the present

429

Modeling particle motion--I

-1

(4 Fluid ,Tp=ZO(ms) ----. Tp=lOO (ms) .....

za P 5 so 6

(4 Tp=loO (ms) .....

_’

I::/: ,/’

-

_A

,/

‘@ 40 P a 30

20

/

_/-

_/

,/

,/

/

,/-

___......~~

,/

,.-’

10 0

,/

_/

_.

__._......

___....

_.__....

__.....-.

___.......

__,....

0

1

2

3

4

so

5

loo

150

Time 6)

lb) Fluid Tp=213(ms) ----. Tp=lOD (ms) .....

200 250 Time (ms)

300

350

400

(b)

1

Tp=20 (ms) ----. Tp=lOO (ms) .-.-.

10 0

0

I

3

2

4

5

Time (s)

Fig. 3. Particle dispersion in the direction of (a) gravity and (b) normal to gravity in isotropic, homogeneous and incompressible turbulence in the presence of gravity.

model are consistem with the past studies for particles with significant inertia. Whether the present model is accurate or not for light particles is an open question. In this context, the extreme case is to simulate the motion of a fluid particle wandering in turbulence. For the present model we need the turbulence intensity, Eulerian temporal and spatial integral correlation scales. The output can be the particle Lagrangian correlation integral time scale and particle diffusion coefficient. Table 2 shows the predicted results for the fluid particle Lagrangian correlation integral time scale rL and the long-time fluid diffusion coefficient D. If the predictions are correct their values should be not far from 32ms (see equation (29)) and 6.272 cm2 s- ‘, the result from Taylor’s theory, respectively. The computed values listed in Table 2 mean that for fluid dispersion the present model also works well.

“0

so

100

150

200 250 Time (ms)

300

350

400

Fig. 4. Particle Lagrangian velocity correlation in the direction of (a) gravity and (b) normal to gravity in isotropic, homogeneous and incompressible turbulence in the presence of gravity.

3.3. Simulation of the experiment of Snyder and Lumley

(1971) One of the most comprehensive sets of experiments on particle motion in turbulent flows was made by Snyder and Lumley. By use of a grid system, they produced a nearly isotropic decaying turbulent air flow (air density= 1.205 x 10V3 gcmm3 and kinematic viscosity = 14.937 x lo-’ cm2 s- ‘). The particles ranged from light particles which follow the turbulent flow to heavy particles which experience both inertia and crossing-trajectory effects. Their diameters and densities are given in Table 3. The particles were injected into the turbulent flow at x/M =20. The measurements were carried beginning from x/M = 68, where x is the distance from the grid and M = 2.54 cm is the grid spacing.

430

Q. Q. LU Table 1. Predicted diffusion coet%cient, fluctuation and Lagrangian time scales in the presence of gravity D,(cm’s-‘)

D,,(cm2s-‘) Fluid T,=20 T,,=lOO

u;(cms-I)

6.17 4.78 1.25

6.17 5.27 2.21

u;(cms-‘) 13.7 10.2 3.51

13.7 10.6 4.32

r: (ms)

r,L(ms)

33.0 47.2 116

33.0 46.0 101

Table 2. Predicted diffusion coefficient, velocity fluctuation and Lagrangian time scales in the absence of gravity &cm2 s-t) 70 -

Fluid Tp=20(ms) ----.

6.17 7.28 8.91

Fluid T,=20 T,=loO

TpelOO (ms) .....

From 1

3

2

4

5

on2

y

39.4 (x/M-

Tiie (s)

Taylor’s

1

0.8

33.0 61.5 148

(38)

Y

frozen

hypothesis

and

kinetic energy dissipation

1 2 42.4(x/M-16)2+39.4(x/M-12)2 Fluid Tp=20 (ms) ----. Tp=lOO (ms) .....

the relation

kinetic energy, is

k=f(z+$+i&,

the turbulent tained as

(ms)

12)’

dk/dt = -E, where k, the turbulent defined as

Fig. 5. Particle dispersion in isotropic, homogeneous and incompressible turbulence in the absence of gravity.

7L

13.8 10.9 7.76

7= 3=7. 2 0

u’(cms-‘)

(39) rate is ob-

1.

(40)

In this simulation, the time step was At = 1 ms. This time step is much smaller than the Kolmogorov time scale tx = 9.8 ms at the first station x/M = 68 of computation. In the present case 7K increases with x. To examine the effect of the time step At, Table 4 lists the asymptotic diffusion coefficients obtained by using At= 1,4,7 and 14ms. Here the subscript 11represents the direction of the principal mean flow, i.e. the xdirection, and I is the transverse direction. 0.2 From Table 4 one cannot see any major trend in the computed results with At. This suggests again that the time step At can be less than the Kolmogorov time 0 250 300 350 400 0 50 100 150 200 scale zK, a conclusion already reached in the last Tiie (ms) section. Figure 7a compares the computed and experimental particle transverse displacement. Good agreeFig. 6. Particle Lagrangian velocity correlation in isotropic, homogeneous and incompressible turbulence in ment is observed. Figure 7b shows particle dispersion the absence of gravity. along the x-direction so as to test the effects of crossing-trajectory and continuity. It is clear that the crossing-trajectory effect reduces the particle dispersion The principal direction bf the flow was vertically rate, while continuity makes the particle dispersion in upward. The experimental mean turbulence data are the longitudinal direction faster than in the normal given by direction. Figure 8a shows the numerical and experimental particle velocity, where u’ is the mean-square z=655 ems-‘, q=O, E=O, (35) root of the particle velocity variance. The results are normalized by U * U, U representing E. Here a larm2 z= (36) ger difference is seen for the hollow glass particles. X 42.4 (x/M-16)’

Modeling particle motion--I

431

Table 3. Parameters of the particles used in the experiment of Snyder and Lumley

Diameter (pm) Density (g cme3)

Hollow glass

Corn pollen

Glass

Copper

46.5 0.26

87.0 1.00

87.0 2.50

46.5 8.90

Table 4. Asymptotic particle diffusion computed by using different At Hollow glass At (ms) 1.0 4.0 7.0 14.0

‘5-u

Corn pollen

-

Copper

D,,

DA

D,,

9~

D,

D,

D,,

D,

5.132 5.467 5.617 5.520

5.440 5.525 5.516 5.396

3.441 3.780 3.696 3.792

2.881 2.900 2.863 2.741

2.472 2.442 2.624 2.480

1.360 1.450 1.378 1.511

2.213 2.435 2.233 2.289

1.275 1.337 1.340 1.362

I

I

,

(a) -

Glass

Hollowglass Corn pollen ----’ Glass ....Copper Hollowglass 0 Cornpollen + Glsss 0

35 ,

I

,

I

I

Fluid 30 _ Hollowglass -.-.Cornpollen ----. Glass ...-. Copper 25

Hollowglass o

0

_......

x

Corn pollen +

Glass

1

(a)

??

0.2

0.1

0.3

0.4

0.5

Tiic (s)

(b) 5 -

Hollowglasr -

Corn pollen ---Glass .....

(b)

30 Fluid -

-

0

Hollowglass -.-.Corn pollen

----Glass ..-.. Copper

I

1

I

I

I

0.1

0.2

0.3

0.4

0.5

Time (s) Fig. 7. (a) Computed (lines) and experimental (symbols) transverse particle dispersion and (b) computed longitudinal partical dispersion in the experiment of Snyder and Lumley.

1

Fig, 8. (a) Computed (lines) and experimental (symbols) particle transverse fluctuating velocity and (b) computed particle longitudinal fluctuating velocity in the experiment of Snyder and Lumley.

432

Q. Q. LU Table 5. Computed particle Lagrangian integral time scale (ms)

xlM

Hollow glass

Corn pollen

711

711

34.03

13

71

31.66

44.23

7L

46.50

Glass 711

60.86

Copper 71

711

52.11

Fluid 71

60.14

7L

62.28

28.86

longitudinal and transverse direction at the section x/M =73. Here z,_ is the fluid Lagrangian integral

time scale estimated according to equation (22). In addition, Fig. 9a and b depicts the Lagrangian particle velocity correlations in the transverse and longitudinal directions, respectively. 3.4. Simulation of the experiment of Wells and Stock (1983) Wells and Stock used a grid system identical to Snyder and Lumley’s to produce a turbulent air flow, but the main direction of the flow was horizontal. Mean data for the turbulent field are

0.2

t

OL 0

. 20

40

60

80 100 120 Tii (ms)

140

160

180 200

ems-‘,

c=655

q=O,

c=O,

m*

7= * 53.224(x/M-7.053)’

7= ’

Z=f I

(41) (42)

m’

(43)

56.546 (x/M - 8.867)’

(44)

Y’

The turbulent kinetic energy dissipation rate is 1 2 + 56.546 (x/M - 8.867)’

Fig. 9. The Lagrangian correlation of particle (a) transverse and (b) longitudinal velocity in the experiment of Snyder and Lumley.

to the low sampling rate used by Snyder and Lumley in the experiment because it is expected that the velocity fluctuation of the hollow glass should not be far from that of the fluid. Figure 8b shows the computed particle fluctuating velocity in the longitudinal direction. In order to investigate the influence of particle inertia, crossing-trajectory, and continuity on particle velocity correlation, Table 5 gives the Lagrangian integral time scale of the four kinds of particles in the

This may he attributed

1’

(45)

Here E is obtained in the same way as in the experiment of Snyder and Lumley. In this experiment 5 and 57 pm glass particles were charged before the grid and an adjustable, uniform electrical field within the test section was used so as to change the resultant force acting on particles and, therefore, the particle terminal velocity. The densities of the 5 and 57 pm particles were 2.475 and 2.420 g cm - ‘, respectively. This experiment is difficult to simulate. The particles were put into flow field before the grid, while all the particle quantities were measured at x/M =30. The measured quantities given are not the relative ones, as in the experiment of Snyder and Lumley, but the absolute values relative to the initial quantities. This means that the whole turbulent field before x/M = 30 contributed to the measured quantities. However, it can be seen from equations (42) and (43) that 7.053 and 8.867 are two singular points. This means that we cannot inject particles before x = 0, like the experiment. Therefore, we must find a method to reproduce the process encountered by particles in the

433

Modeling particle motion--I experiment. The 5 /Am particles are introduced at x/M = 10. Their fluctuating velocities are equal to the fluid ones; the transverse positions are also given by a random variable satisfying a Gaussian p.d.f. with the variance being OScm’, the value extrapolated from the measured particle transverse dispersion. For the 57 pm particles we cannot follow the procedure as for the 5 pm particles because the particles are affected significantly by their inertia and gravity. Instead, we put the 57 pm particles into the flow field at x/M = 20. The initial particle fluctuating velocities are sampled from random variables satisfying Gaussian p.d.f. with the measured variance; the particle transverse positions are given by a random variable satisfying Gaussian p.d.f. with the measured particle dispersion as the variance. The time step was At=6 ms. More smaller time steps, At, were a.lso tried and they showed almost the same results as the given ones. Figures 10 and 11 show the experimental and computed particle transverse dispersion for the 5 and 57 pm particles, respectively. Here the “ K” represents the particle terminal velocity resulting from the combined action of gravity and the electrical force. It is seen that as the particle terminal velocity increases, the difference between the computed and experimental data becomes larger, which is particularly pronounced for the 57 pm particles. This may be because the computational process devised above for the 57 pm particles does not match what happened in the experiment. However, the computed results reflect the strong influence of t:he crossing-trajectory on particle dispersion. In the absence of gravity it can be observed that the 57 pm particles disperse faster than the 5 pm particles do. Figure 12a and b are the 5 pm particle fluctuating velocity in the longitudinal direction (the x-direction) and transverse direction, respectively. Figure 13a and b are the 57 pm particle fluctuating velocity in the

4 3.5

I’

/.

,’ I’

vt=o.oo----.

3

,’

vm39.7 ...-V&1,5 Vt=l21.6 vt=O.Oo 0 vt=39.7 + Vt=81.5 0

5 2.5 “E 0 2 s ‘2 B .v: 1.5 n

I’

I’

,’ I’ 0 ,’ ,’

: _... _..’ . ..’ ,,,....... ._” ,... ._,’ + ,../”

,A

._I’,,_.... .” I’ O,..‘.,,,. f..-..

vt=121~

1 0.5 0 40

20

0

60

80

100

x/M

Fig. 11. Computed (lines) and experimental (symbols) particle transverse dispersion in the experiment of Wells and Stock for the 57 pm particle.

I41 _

‘8 5 3 ? 56 El 4

I

(4

Fluid Vt=Om Vt=5.86 vt=17.06 Vt=23.65 vt=o.Oo VtJ.86 Vk17.06 Vt=23.65

..... -.-.----. 0 + 0 x

/

2

0’ 0

I

I

I

I

20

40

60

80

100

x/M

I

14

I

I

I

(b)

4-, 3.5

12 -

t 10 -

3 -

’ 7

E2.5 +k 0 6 2 ‘2 e: .z 1.5 CI

Ve5.86 .-... vt=17.06 -‘-‘Vt=23.65 ----. vt=o.oo 0 Vt=5.86 + Vm17.06 0 Vt=23.65 x

1 ??

0.5 0-1

0

8-

3 ?

20

40

01

1

I

60

80

0 100

Fluid vt=0.@J Vt=5.86 vt=17,lJ6 -.-.Vtz23.65 --__.

1

I

I

I

I

20

40

60

80

100

x/M

x/M

Fig. 10. Computed (lines) and experimental (symbols) particle transverse dispersion in the experiment of Wells and Stock for the 5 pm particle.

Fig. 12. (a) Computed (lines) and experimental (symbols) particle longitudinal dispersion and (b) computed particle transverse velocity in the experiment of Wells and Stock for the 5 pm particle.

434

Q. Q. LU

longitudinal and transverse directions, respectively. Disagreement is apparent between the predicted and experimental results for the 57 pm particles. This may be attributed to the initial velocity assigned to particles. In the computations, we put the particles into the flow field at x/M = 20 and let their initial velocity be equal to the local mean fluid velocity. This means the initial particle fluctuating velocity at x/M = 10 is zero, while the experimental data tend to infinity near

X/M = 10. Table 6 lists the 5 pm particle Lagrangian integral time scales computed at x/M = 45. Here the values of T;‘P are obtained by using the ratio of the particle Eulerian velocity correlation to its Lagrangian one when V,=O. Table 7 lists the 57 pm particle Lagrangian integral time scale. Figure 14a and b are the 5 pm particle Lagrangian velocity correlation in the longitudinal and transverse directions, respectively. Figure 15a and b shows the 57 pm particle

I

18 -

(a) 0.8

16 -

Flutd vt=o.oo ----

14 -

vw39.7 -. -. Vt=81.5 Vt=121.6 - - vt=o.OO 0 Vt=39.7 + Vt=81.5 0 kk121.6 x

8 ?

12 -

?

10 -

,’ ,I

/

,’

*’

...’ ..:

,.I I

;”

60.2

0 0

0

40

20

60

10

20

30

40 50 Time (ms)

60

70

20

30

40 50 Time (ms)

60

70

100

80

. J 80

x/M

1

2oI

(b)

18 -

0.8 16 14 z! 12 7 = 10 -

Fluid vt=o.OO Vk39.7 Vt=81.5 Vki21.6

---....’ - - -

__...-y -’

2?

)

,..

8-

,,

.’

___... ___..-

6 -

____---,

4 -

I/

,__. ..

____---

0.2

~

2 -

0

0'

0

40

20

60

80

-I

0

10

100

80

x/M

Fig. 13. Computed (lines) and experimental (symbols) particle (a) longitudinal and (b) transverse velocity in the experiment of Wells and Stock for the 57 pm particle.

Fig. 14. (a) Computed (lines) and experimental (symbols) Lagrangian longitudinal particle velocity and (b) computed Lagrangian transverse particle velocity correlation in the experiment of Wells and Stock for the 5 pm particle.

Table 6. 5 pm particle Lagrangian integral time scales computed at xfM = 45

Table 7. 57 pm particle Lagrangian integral time scales computed at x/M = 45

v(cms-‘) t (ms) ~~(ms) ~~~~ (ms)

0.0 19.5 19.9 19.6

5.86 19.7 19.9

17.06 16.7 18.7

23.65 15.2 18.4

K(cms-‘) t,,(ms) 7I (ms) 77, (ms)

0.0 45.6 46.8 37.6

39.7 30.7 32.4

81.5 22.3 25.8

121.6 21.4 21.2

435

Modeling particle motion-I I

I

I

I

I

I

(a) 0.8

vt=o.cxl

-

VW39.7 ----.

Vt=81.5 ....Vtz121.6 ....... vt=O.Oo 0

ZJ 0

20

40

60 (m”s’: Time

100

120

140

(b) vt=o.oovw39.7 Vtz81.5 Ve121.6

0

20

40

60 80 Time (ms)

----’ ..-..

100

120

140

Fig. 15. (a) Computed (lines) and experimental (symbols) Lagrangian longitudinal particle correlation and (b) computed Lagrangian transverse particle correlation in the experiment of Wells and Stock for the 57 pm particle.

Lagrangian velocity correlation in the longitudinal and transverse directions, respectively. The symbols in Figures 14a and 15a are the data inferred from the experimental Eulerian velocity correlations.

4. DISCUSSION AND CONCLUSIONS

Based on the idea of a time-series analysis (Box and Jenkins), a Lagrangian model has been established for particle dispersion in turbulent flows. The model was applied in simulating particle dispersion in stationary, homogeneous, isotropic and incompressible turbulence and in the two grid-generated turbulent flows of Snyder and Lumley
can reproduce the findings in the studies of Reeks, Pismen and Nir, and Nir and Pismen. For the experiment of Wells and Stock apparent discrepancies are found. This may be because we could not introduce particles as in the experiment. It is clear that the present model can be used for cases where more terms are added to the equation for particle motion. However, there are two important issues open for discussion about the values of the coefficients C 1, Cz and CJ. Our computational experience shows that the relations in which they occur and their magnitude were important in making predictions. In fact, there is some diversity in the values of these coefficients. Based on the Corrsin’s conjecture, Saffman derived a value of 0.4 for C2 a value of 0.359 for C2 can be inferred from the study of Reeks and Pismen and Nir obtained. In addition, from the study of Pismen and Nir a value, 0.73, can be derived for C1. However, from Pismen and Nir’s comparison with the experimental particle dispersion of Snyder and Lumley, there seems some evidence which supports a larger value for C1. In their comparison, good agreement is observed for the copper particle but for the corn pollen their method overpredicts significantly. This is understandable because Corrsin’s conjecture becomes realistic as the particle time constant increases. In fact, in this case the effect of the Lagrangian turbulence structure already partially influences the corn pollen dispersion. Since the computed result so largely exceeds the experimental values considerably, it may be inferred that their method leads to too large a Lagrangian integral time scale, i.e. the value of C’, is too large. However, these arguments need to be further clarified and checked by some independent methods. One choice is probably by use of direct numerical simulation (DNS). Without introducing any models such a method may provide information unavailable to us before, most notably for time and length scales. This, in one respect, may help to check our model and to develop new ideas regarding modeling. REFERENCES

Banerjee S. (1990) Modeling considerations for turbulent multiphase flows. Invited paper in Engineering Turbulence Modeling and Experiments (edited by Rodi W. and Ganic E. N.), pp. U-866. Elsevier, New York. Berlemont A., Desjonqueres P. and Gouesbet G. (1990) Particle Lagrangian simulation in turbulent flows. Int. J. Multiphase Flow 16, 19-34. Boothroyd R. G. (1971) Flowing Gas-Solid Suspmsions. Chapman & Hall, London, Box G. E. P. and Jenkins G. M. (1976) Time Series Analysis. Holden-Day, Oakland, CA. Burnage H. and Moon S. (1990) PrCdttermination de la dispersion de particules matirielles dans un tcoulement turbulent. C. R. Acad. Sci. Paris Serie 11 310, 1595~1600. Csanady G. T. (1963) Turbulent diffusion of heavy particles in atmosphere. J. armos. Sci. 20, 201-208. Durbin P. A. (1980) A stochastic model of two-particle dispersion and concentration fluctuations in homogeneous turbulence. J. Fluid Mech. 100, 279-302.

436

Q. Q. LU

Frenkiel F. N. (1948) Etude statistique de la turbulencefonctions spectrales et coefficients de corr&lation. Rapport Technique, ONERA No. 34. Gosman A. D. and Ioannides E. (1981) Aspects of computer simulation of liquid-fueled combustors. Presented at the AIAA 19th Aerospace Science Mtg., St. Louis, MO, Paper 81-0323. Hinze H. 0. (1975) Turbulence, 2nd Edn. McGraw-Hill, New York. Kallio G. A. and Reeks M. W. (1989) A numerical simulation of particle deposition in turbulent boundary layers. Int. J. Multiphase Flow 15, 433-446. Lu Q. Q., Fontaine J. R. and Aubertin G. (1992) Particle motion in two-dimensional confined turbulent flows. Aerosol Sci. Technol. 17, 169-185. Lu Q. Q., Fontaine J. R. and Aubertin G. (1993a) Numerical study of the particle motion in grid-generated turbulent flows. Int. J. Heat Mass Transfer 36, 79-87. Lu Q. Q., Fontaine J. R. and Aubertin G. (1993b) Particle dispersion in shear turbulent flows. Aerosol Sci. Technol. 18, 85-99. Lu Q. Q., Fontaine J. R. and Aubertin G. (1993~)A Lagrangian model for solid particles in turbulent flows. Int. J. Multiphase Flow 19, 347-367. Nir A. and Pismen L. M. (1979) The effect of a steady drift on the dispersion of a particle in turbulent fluid. J. Fluid Mech. 94, 369-381. Ormancey A. and Martinon A. (1984) Prediction of particle dispersion in turbulent flows. Physico-Chem. Hydrodynam. 5,229-240.

Parthasarathy R. N. and Feath G. M. (1990) Turbulent dispersion of particles in self generated homogeneous turbulence. J. Fluid Mech. 220, 515-537. Pedinotti S., Mariotti G. and Banerjee S. (1992) Direct numerical; simulation of particle behaviour in the wall region of turbulent flows in horizontal channels. Int. J. Multiphase Flow 18, 927-941. Philip J. R. (1967) Relation between Eulerian and Lagrangian statistics. Phys. Fluid 10, Suppl. 9, S69-71. Pismen L. M. and Nir A. (1978) On the motion of suspended particles in stationary homogeneous turbulence. J. Fluid Mech. 84, 193-206.

Rashidi M., Hestroni G. and Banerjee S. (1990) Particle-turbulence interaction in a boundary layer. Int. J. Multiphase Flow 16, 935-949. Reeks M. W. (1977) On the dispersion of small particles suspended in an isotropic turbulent fluid. J. Fluid. Mech. 83, 529-546.

Rubinow S. I. and Keller J. B. (1961) The transverse force on a spinning sphere moving in a viscous fluid. J. Fluid Mech. 11, 447-459. Saffman P. G. (1961) An approximation calculation of the Lagrangian autocorrelation coefficient for stationary homogeneous turbulence. Appl. Sci. Res. 11, 245-255. Saffman P. G. (1965) The lift on a small sphere in a slow shear flow. J. Fluid Mech. 22, 385-400. Sato Y. and Yamamoto K. (1987) Lagrangian measurement of fluid-particle motion in an isotropic turbulent field. J. Fluid Mech. 175, 183-199.

Shuen J. S., Chen L. D. and Faeth G. M. (1983) Evaluation of a stochastic model of particle dispersion in a turbulent round jet. A.1.Ch.E. J. 29, 167-170. Shuen J. S., Solomon A. S. P., Zhang Q. F. and Faeth G. M. (1985) Structure of particle-laden jets: measurements and predictions. AIAA 23, 396-404. Snyder W. H. and Lumley J. L. (1971) Some measurements of particle velocity autocorrelation functions in a turbulent flow. J. Fluid Mech. 48, 41-71. Squires K. D. and Eaton J. K. (1990) The interaction of particles with homogeneous turbulence. Report No. MD55, Dept. of Mechanical Engineering, Stanford University, CA. Taylor G. I. (1921) Diffusion by continuous movements. Proc. R. Sot. Series A 20, 196-211. Vetterling W. T., Teukolsky S. A., Press W. H. and Flannery B. P. (1988) Numerical Recipes Example Book. Cambridge University Press, Cambridge. Wells M. R. and Stock D. E. (1983) The effects of crossing trajectories on the dispersion of particles in a turbulent flow. J. Fluid Mech. 136, 31-62. Williams F. A. (1985) Combustion Theory, 2nd Edn. AddisonWesley, Menlo Park, CA.