7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Closure approximations in particle dispersion coefficients
derived from pdf kinetic equations
A. Bragg* D.C. Swailes* and R. Skartlient
Mechanical & Systems Engineering, Newcastle University, Newcastle upon Tyne NE1 7RU, UK
t Institute for Energy Technology, P.O. Box 40, N2027 Kjeller, Norway
andrew.bragg@ncl.ac.uk, d.c.swailes@ncl.ac.uk, roar.skartlien@ife.no
Keywords: disperse flows, boundary layers, inhomogeneous turbulence, pdf kinetic models
Abstract
A comparison is made between two distinct formulations of the dispersion coefficients in pdf kinetic equations. By
reference to some simple constraints imposed by the underlying particle equation of motion it is shown that the
difference between these formulations can be nontrivial, and that one form of the dispersion coefficients provides a
more consistent model than the other. A method for closing the averages appearing in these coefficients is presented.
This takes into account turbulent inhomogeneities and the influence of particlesurface interactions. The approach
makes use of Eulerian twopoint, twotime fluctuating fluid velocity correlations, and illustrations of such correlations
obtained from LES of turbulent channel flow are given. The effectiveness of the approach is assessed by reference
to a simple onedimensional system, making comparison with particle tracking simulations and the alternative,
localhomogeneous closure approximation.
Introduction
The kinetic equation first introduced by Reeks (1991)
represents an important model for the transport of dis
persed particulates in turbulent flows. In particular, this
type of pdf model provides a rational basis for the con
struction of transport equations for meanfield variables
(mass, momentum and kinetic stresses) of the disperse
phase (Reeks (1992)). The coefficients that appear in
these meanfield equations, and that describe the influ
ence of the fluid turbulence, follow directly from the dis
persion coefficients in the kinetic equation. The detailed
specification of these dispersion coefficients is therefore
critical. However, in formulating expressions for these
coefficients some level of closure modeling is required,
and it is consideration of such closures that forms the
basis of the work presented in this paper. As originally
formulated by Reeks, using the approach of Lagrangian
History Direct Interactions (LHDI), the key quantities
that require closure are the correlations
S(X', t') j(x,t) t'
Here u denotes the fluctuating fluid velocity field, and
the subscripts x, v indicate a conditional ensemble av
erage, based on those particle trajectories x' = x(t')
that satisfy x(t) = x and v(t) = (t) = v.
Slightly different expressions for the dispersion coeffi
cients were derived by Swailes & Darbyshire (1997) us
ing the FurutsuNovikov (FN) correlationsplitting for
mula (see, for example, Klyatskin (2005)). In these ex
pressions the correlations (1) are replaced by
Ri(x', t'; x, t)
\ I x,v
t' < t (2)
where Rij(x',t';x,t) = (u(x',t') uj(x,t)). Super
ficially the correlations given by (1) and (2) appear
similar. However the twopoint, twotime fluid velocity
correlations defined by the tensor R are strictly Eu
lerian and are based on all realizations of u. So that
while the ensemble in (2) is still conditioned by the
constraints x(t) x, v(t) v it contains, unlike (1),
information determined by all realizations of u. The
natural question then is whether or not this difference is
siilm.iiiI and, if so, which of the two approaches can
be considered to represent a more consistent model with
respect to the underlying particle equation of motion
used in the derivation of the kinetic equation. It should
also be noted that the question relates to Corrsin's
Ii\i'1,'lcii (Corrsin (1959)) which, in essence, invokes
the assumption that (1) and (2) can be considered
equivalent. In the next section it is shown that the
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
difference between (1) and (2) can be significant, and
can lead to crucially different predictions for important
correlations involving the dispersion coefficients. The
conclusion being that the kinetic equation formulation
involving (2) offers a more consistent model than that
using (1).
Notwithstanding the difference between (1) and (2)
some closure is necessary for these averages, and the
approximation that is normally invoked is to set x' x
and to introduce Lagrangian time scales Tij that charac
terize the rates of decorrelation with respect to t t.
These time scales also require modeling. This closure
approach has proved satisfactory in turbulent regimes
that do not exhibit strong inhomogeneities. However the
effectiveness of such closures in the presence of large
turbulent gradients is less clear and, of course, the spec
ification of suitable Lagrangian time scales Tij is less
straightforward. A stringent test case for such a closure
approach is therefore provided by particle transport in a
turbulent boundary layer. Further, in the near wall re
gion of such flows any particlesurface interactions will
influence the values of the averages defined by (1) and
(2), and there is evidence to suggest that such effects will
further compromise the accuracy of this closure method
ology (Skartlien (2007)). In view of these observations
this work considers the construction of alternative clo
sures strategies for (1) and (2) that take into account the
effect of both turbulent inhomogeneities and particle
surface interactions. The effectiveness of the closure
approach introduced here is assessed for a simple one
dimensional system in which particles experience elastic
particlesurface interactions. Results are compared with
those obtained from corresponding localhomogeneous
closure approximations, and also with exact values of
the required averages as obtained from particle tracking
simulations.
in its simplest form, can be written
a
09
vjp
8xI
(Fj + I j)p
vj
a L a, [ a"
+ r + r
dvj [axi avi j
This form is obtained from both the LHDI and the
FN approaches (Reeks (1991); Swailes & Darbyshire
(1997)). Moreover, the expressions for the dispersion
coefficients v, A, /p derived from the two approaches
appear, at least superficially, to be equivalent. Specifi
cally, they can be written
t
KJ (umt;') j (x', t'; t) di
\ On \ / x,v
t I
Ai Him.(t; t')m ( x', t'; \ dt'
0 \ XV
Pj = f Him(t; t') Rm(x', t'; x, t)) dt'
(5)
The response tensor Him is defined the same in both
approaches
Kim(t; t')
6xj(t)
6fm(X',t') 6t'
however the specification of m,j in (5) differs: The
LHDI derivation gives
Rj (x',t';x,t) = fm(x',t')fj(X,t)
while the FN derivation gives
Rj (x', t'; x, t)
(fm(X', t')f (X, t))
Kinetic Equations LHDI or FN?
The starting point for the construction of kinetic equa
tions is a stochastic model for particle motion, generally
of the form
x(t) = (t) = F(x(t), v(t)) + f(x(t), t) (3)
where F and f define the mean and fluctuating acceler
ations of a particle with position x, velocity v at time t.
The pdf p(x, v, t) defining the joint distribution of x(t)
and v(t) is then governed by a kinetic equation which,
Note that (7) and (8) are intrinsically different, in as
much that (7) is a stochastic quantity whereas (8) is de
terministic. It is this difference that leads to the reduced
forms given by (1) and (2), which result when H is taken
out of the conditional averages in (5). However, this re
duction invokes a further approximation, namely that H
is independent of particle path (i.e. not stochastic). The
spatially dependent nature of f means that H must nec
essarily be stochastic, and it is only by neglecting the
gradients in f that a deterministic approximation for H
can be formulated. This device lies at the heart of the
commonly used interpretation of H as a deterministicc)
Green's function for x(t). Such additional approxima
tions for H must, of course, also be taken into account
in any analysis of the implications arising from the dif
ferences between (1) and (2).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
As a simple illustration of instrinsic differences orginat
ing from the alternative interpretations (7), (8) for R in
(5) consider the simple firstorder scalar stochastic dif
ferential equation
1(t) f(x,t) x(t) (t),
be such that
(v) f3)
(Vi* h)
(O) x (9)
(X j + Aij)
(v*Kj + Pi)
where 7(t) is a zeromean (and Gaussian) stochastic
process correlated in time. The response function for
x(t) is given by
6x(t)
6f(x', t')6t'
h(t t')
te
Consequently x(t) = x0h(t), and then
K K (t; t')R(x' t'; xt))x
x{ x (h(t
xx (h(t
t')h(t')7(t' )(t)
P)h(PI X \ /w t X
using (7)
using (8)
(11)
These alternative forms would be equivalent if h
independent of 7. That this is not the case is evi
from (10). It is clear, therefore, that care must be tak
the more general case when interpreting R in (5). (1
that the conditionality on the averages in (11) sh
not be omitted, since the realizations of 7 (and hen
must still be restricted to those that give x(t) x)
a further illustration of the consequences of the di
ences manifest in (7) and (8), and of the treatment
as a determinstic Green's function, consider the ev
tion equations for the secondorder correlations (x*
(v*v*), where x* x (x) and v* v (v).
equation of motion (3) gives the following exact sy
d
(x (F + fj)) + KvvJ)
\' i jU
d (t4v) ( (F + (F )) + (v (Fi + f1))
while the kinetic equation (4) gives (also exact and i
pendent of the LHDI or FN interpretation)
d
d (x v 3)= (x (F + ) + Ai) + (v>v*)
d
d (v ) = (v* (F + K) + pij) + ( v(Fi + i)
In (12) and (13) F, f, v, A and p are evaluated alon
particle trajectories. Comparing (12) with (13) w
that, for consistency, the definitions of v, A and p 1
In Appendix it is shown that, based on the LHDI for
mulation, (with H interpreted as the Green's function
for x(t))
() = ((x x)f) and (pj) = ((v, v) fj)
(15)
where x(t), v(t) satisfy x = v = F(x, v). Comparing
(15) with (14), it follows that, for consistency, the LHDI
formulation should also give (when x (x))
(4xK) 0, (v<*,) 0 (16)
However the same LHDI formulation for v also implies
(appendix)
were
ident
en in
Note
would
KV KJ) ( (x,
\v[n / ^ (
ce h) Clearly there is no reason to expect that the correlations
.As given by (17) should be identically zero as required by
iffer (16), and this suggests that this LHDI based interpreta
of H tion of the coefficients (5) is not strictly consistent with
volu the original equation of motion, and that the FN inter
v*), pretation of (5) provides a more consistent pdf model.
The One benefit of the conclusion that R in (5) ought to be
stem interpreted as (8), as opposed to (7), is that this simpli
fies the formulation closures for the coefficients (5) since
these are then reduced to forms involving (2) rather than
(12) (1). The correlation tensor R, which describes single
phase Eulerian fluid statistics is presumed to be known
(can be calculated) ab initio and no further assumptions
S (e.g. Corrsin's hi ,lici,,i need be invoked.
nde
By way of illustration figure 1 shows some spatial corre
lation profiles computed from LES of a turbulent chan
nel flow at Re, = 180. Specifically, the figure shows
profiles for
+ ji)
(13)
g the
e see
must
RSc (x',t;X, t)
R22 (x, t; x, t)
where the subscript 2 indicates the wallnormal channel
direction.
afi
09x,
a f)
09,
Spatial Correlations in a turbulent boundary
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
are taken in account directly in the closure. In addition,
the influence of particlesurface interactions can be ac
commodated within the formulation of approximations
to p(x',t'x, t). In this work attention is confined to
elastic particlesurface interactions (specular reflection).
In this case the problem can be simplified by recasting
the boundary as a line of symmetry: For a 2D flow
restricted to x2 > 0, with x2 = 0 denoting the flow
boundary at which particle are reflected elastically
(v2 v2), the domain can be extended to x2 < 0 by
defining (for x2 < 0)
F(x,v) J.F(J.x,J.v), J 0
0 ) (21)
As figure 1 demonstrates, the spatial correlations within
the turbulent boundary layer (along with the other
turbulence properties) exhibit strong spatial variation.
The characteristic distance that a particle moves within
the time for which the local fluid velocity field is
correlated determines how sensitive the particle is to the
inhomogeneity of the turbulence. Therefore, smaller
particles can be very sensitive to the inhomogeneity of
the turbulence making it essential to incorporate this
effect in any closure model for (Rij (', t';x, t))x for
particle transport in a turbulent boundary layer.
In the next section we consider how knowledge of Rij
can be utilized in the construction of closure approxima
tion for conditional averages that involve evaluation of
this fluid correlation tensor along particle trajectories.
Closure Modeling
The dispersion coefficients (5) appear in the transport
equations for the particle phase meanfield variables in
velocity averaged forms, 9, A, p (see, for example,
Swailes et al. (1998)), where
and similarly for u and f. Trajectories governed by
(3) in this extended, symmetryline model (i.e. with no
boundary at x2 0 and with a suitably symmetrized
initial distribution) then generate the required trajecto
ries of the original reflecting boundary model through
the mapping (Xl, X2) (Xz,'' ). Moreover, in this
symmetryline formulation we can replace (20) with the
equivalent exact representation
Ki *x' t'; Xt)
Rij(x',t';x,t) (p(x', t'x,t)
x!
... + p(J x', t'x, t)) dx'
(22)
with the interpretation that p(x', t'x, t) now represents
the conditional spatial distribution in the symmetryline
model. The advantage of this representation is that it
is now unnecessary to consider particlesurface interac
tions at X2 = 0: The influence of such interactions on
(20) is recovered via the use of (22). The distribution
p(x', t'lx, t) is approximated using
p(x', t'lx,t) = p(x', t' tx, 0)
S= Jp dv
V
etc. (19)
These forms of the coefficients involve the condi
tioned averages (as determined from the FN approach)
(Rij(x', t'; x, t))x. By definition these can written as
(Ri (x', t'; x, t)
\ / x
Rij(x', t'; x, t) p(x', t'lx, t) dx'
(20)
The aim is to close the conditional averages
(Rij(x',t';x,t))x by constructing approximations
for the conditional spatial distribution p(x', t' x, t) and
then evaluating explicitly the corresponding righthand
side of (20). In this way inhomogeneities in the
underlying fluid turbulence, as characterized by Rij,
where p(x, v, s x) gives the distribution of x(s), v(s)
(s > 0) subject to imposed 'initial' (s 0) conditions;
x(0) x, and pdf y(v, ulx) defining the joint distri
bution of v0 v(0) and uO u(x, t). The form of
the kinetic equation defining this pdf will be analogous
to (4), although the correlation between v and uO (as
defined by 4) will modify the form of the response ten
sor for x (this is discussed in Hyland et al. (1999)). For
the closure of the dispersion tensors it is only necessary
to consider the behaviour of particle trajectories x over
times intervals 0 < s < t t' of the order of the correla
tion times inherent in R and this prompts several differ
ent closure strategies. The simplest approach would be
to introduce the approximation f(x(s), t s) z f(x, t).
Figure 1:
layer
1 = p dv,
P J
V
p (x', v, t t'l x) dv
This is essentially the approach adopted by Oesterl6 &
Zaichik (2006). An extension of this idea, and the one
that is implemented here, is to model the Lagrangian
fluctuations f(x(s), t s) f(s; x) as a prescribed
stochastic process, with
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
(2005))
32 ((K/3a 0(0 UJja i)
=7 (/ ) + 6i 1+3) (u2 )
a+ (( + P) )
(27)
(f(x, t) f(xt))
x exp [a s s'l]
This, of course, presumes that f is statistically sta
tionary in time and, more problematically, introduces
the same Lagrangian decorrelation rates a = (x)
which, if known, would provide the direct, local
homogeneous closure for (2). However, by introducing
these ratescales into the modeling approach developed
here it becomes possible to treat these parameters as
variables, and to identify their precise values via an
iterative process. In addition, it is questionable whether
direct inclusion of the same Lagrangian scales into
the localhomogeneous closure approximation would
necessarily produce the same level of resolution as
that offered by the use of (20) and (22) which take
explicit account of spatial variations in the Rij and
of the influence of boundary effects. Note that the
simpler model f(x(s),t s) z f(x,t), which has
proved satisfactory in homogeneous systems (Oesterl6
& Zaichik (2006)), is a special case obtained by setting
a 0.
With model (24) for f (x (s), t s) (and with F linear in x
and v) the kinetic equation admits an exact closed form
Gaussian solution p. The distribution p, as defined by
(23), is then also Gaussian, with mean and covariances
that are easily evaluated (Swailes & Darbyshire (1997,
1999)). These parameters, defining the distribution of p,
will depend on the model for y, and the simplest first
approximation to this can be obtained from asymptotic
(t oc) expressions: Specifically, if
(x) ( vv vu )
C(x) = ouv ouu
denotes the covariances of (v un) then take Huu
')(u" u (u))and
21 +;
( 1 O)i Pn,
With a Stokes drag law for particle motion in a simple
linear shear, d(ui)/oxj = 7Yi1j2, we obtain (Reeks
1 D Test Case
A simple 1D model serves as a useful test case to evalu
ate the above approach for constructing closure approx
imations. The particle equation of motion is based on a
simple drag law;
(t) = i(t) = 3(u(x(t,) t)) (28)
with trajectories confined to the interval 0 < x < X
by reflecting boundaries at x 0 and x X. The
fluctuating velocity field u(x, t) is modelled in the form
u(x, t) = (x)w(t) where y(x) is a deterministic 'pro
filing' factor and the velocity w(t) a zeromean Gaussian
process with (w(t')w(t)) = exp[,t t']. The
twopoint, twotime Eulerian correlation function is then
R(x', t';x, t) 32(x')q(x)2 exp[,,., I t]
(29)
Although artificial this closed form expression for R
does provide a basis for evaluating closure integrals
of the type (22) and, via y, also enables large gradi
ents in the profile of the meansquare velocity (u2)
r y2(x): A simple model for y is
y(x) A (1
exp [(kx)2])
The coefficient A is just the normalization factor to
give y(X) 1 With k > X the profile for (u2)
then exhibits large gradients in the neighborhood of
x = 1/ 2k < X. Applying the approach outlined in
the previous section leads to the following approxima
tion for the conditional spatial distribution p(x', t' x, t)
(treating x 0 as a line of symmetry)
p(x',t'x,t) (27) 1 exp [ (x' x)2
(31)
Here 0, 0 6(s;x), s = t t' > 0, takes the form
(appendix)
B (s; x) = G O, ,, 0,, + G,,,, (32)
where O (vv0), O .".,") and O, = ', "").
Details of the functions G (s; x) are presented in ap
pendix 2. The relation of 08 and O,, to 0u (u2)
follows from the onedimensional form of (26), (27).
Specifically
3 +/ u
0 c%
VU a3 O (33)
Ov o p
Sf(s'; x)f(s; )
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Using these results the closure approximation for the
conditional average (R(x', t'; x, t))x (x < X) becomes
J R(x', t'; x, t) ((x', t'x, t) + p(x', t'x, )) dx'
x'>0
S32(a exp[a s] (x) (s; x)
(34)
where
+0s; x) p t
IF(S;x) j 1 \)xTp(x',t'x,t) dx'
1 F 2]
1 exp (kx/a)
a
a(s;x) = 1 + 2k2 (s;x)
is the closure approximation for (y(x'))x, x < X. For
this model the Lagrangian rate scale a(x) is given by
a1 1
/32/2(x)\
R(x', t'; x, t) xdt'
t)
j V j' (X)) ep[
P0(x
the results presented here have been generated by simply
taking a = w. Further, the results are presented for the
normalized form of the model, scaled such that a 1
and X 1.
To evaluate how results from the closure model proposed
here differ from those obtained from the corresponding
localhomogeneous approach we first compare values of
the two approximations to (R(x', t';x, t))x: For s > 0
we have, respectively,
R(x',t';xt) exp[ ] 7p.(x)(s;xx)
C(s; x)
R(x',t'; ,t))x
/32 ,2 exp[as02 (x)
Co(s; x)
Figures 2 to 4 illustrate the differences between these as
measured by Q(s; x) C/Co, x > 0, for St 0.5,
4 and 10. With a = aw this measure reduces to Q
'I(s;x)/y(x). The results were produced with k 10
and ac 10 1, and the figures show Q in the region
0 < x < 0.2 over the time interval 0 < s < 2.
St=0.5
a (t t')]dt'
(36)
so that, replacing ((x'))x in (36) with t(s;x), pro
vides an implicit equation for a. It is anticipated that
a will only differ markedly from aw (which is spatially
uniform in this simple model) in the nearwall region 
due to the influence of particlesurface interactions and
the gradients in y.
Results
From (36) it follows that if ( (x'))x y (x) for t t' ~
(aw)1 then a(x) ~ a,. This, of course will be the
case for large Stokes number (St a/3) particles.
With decreasing St this approximation becomes less ap
propriate and would certainly be expected to produce
unsatisfactory results if used in the formulation of lo
cal homogeneous closure approximations. Of course,
for St < 1 the usual model for a would be the fluid
Lagrangian rate scale (often used for St ~ 1 as well). It
is of interest, therefore, to test whether or not the same
level of approximation performs any better when applied
to the closure approach developed here, and whether the
two closure strategies would then produce significantly
different values for the dispersion coefficients appearing
in the meanfield transport equations. With this in mind
0.05
0.2 0
Figure 2: Comparison of closure model predictions:
Q C/Co
St=4
0.05
.0.1
0.15
0.2 0
Figure 3: Comparison of closure model predictions:
Q C/Co
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and the results were produced with k = 10, = 10 1
and show Q in the region 0 < x < 0.2 over the time
interval 0 < s < 2. The results are shown Figures 5 to
7.
St=0.5
Figure 4: Comparison
Q C/Co
of closure model predictions:
As expected, for St > 10 the two models produce
essentially the same approximations over this timespan.
However there are clear and icgnilik.iim differences for
St = 4 and less: Very close to the wall (x < 0.05) it is
seen that Q > 1. This is a consequence, in part, of the
fact that the closure model used here takes into account
the influence of the elastic collisions at x 0 which will
result in (x')x > x in this nearwall region; particles
in this region will, on average, sample fluctuating
velocities significantly larger than those at x. Further
away from the wall (0.05 < x < 0.2) it is observed
that Q attains values significantly less that 1. This can
be attributed to the gradients (actually the curvature) in
y which again bias the values of u sampled along the
trajectories of particles coincident with x at a given
time. The deviation between the two models is seen to
occur more rapidly in the case smaller Stokes number
particles (St =0.5) and it is clear that it is the inclusion
of temporal variations in F that capture the coupled
influences of particlesurface interactions and spatial
variations in (u2).
To investigate these effects further a particle tracking
simulation of trajectories governed by (28) was per
formed in order to compare results for (R(x', t'; x, t))x
generated by particle tracking and the closure model.
For s > 0 we have, respectively
Figure 5: Comparison of closure model and particle
tracking results for: Q = C/Cpr
St=4
1
0.8
Q 06
0.4
0.2
0.
2
0.15
0 0
Figure 6: Comparison of closure model and particle
tracking results for: Q = C/Cpr
St=10
R(x', t'; x,t))
( ',t';Xt))X
32T2 exp[a)s]u(x)'(s;x)
C(s; x)
'2,o exp[as]y(x) ((x'))x
CpT (S; X)
The particle tracking and closure model results are com
pared by the measure Q(s; x) C/CpT, x > 0, for
St = 0.5, 4 and 10. Again a = aw in the closure model
Figure 7: Comparison of closure model and particle
tracking results for: Q = C/OCp
St=10
^^^LSfZPI
0.05
x 0.1
0.15
0.2 0
0.15
0.1
x
0.05
C'0?
o C' :
2' .1
CI 
0
2
0 0
0.05
The results show that in the nearwall region
(0 < x < 0.2), where particlesurface interactions
and large gradients in y((x) are important, C can
provide satisfactory approximations to (R(x', '; x, t))x.
In the region 0.05 < x < 0.2, where there are large
gradients in y(x) the closure model predictions are in
good agreement with the particle tracking results across
the range of Stokes numbers tested. This demonstrates
that the closure model is able to describe the effect of
the spatial variations in u on (R(x', t'; x, t))x.
However very near the wall (x < 0.05) there is clearly
a difference. In this region, the particle tracking results
showed a very strong spatial variation in a(x), varying
from a, by a considerable degree for all of the Stokes
numbers tested. In the closure model, a was set equal to
a, which, in light of the particle tracking data, is clearly
not appropriate. Also, the closure model prediction very
close to the wall is worse for smaller Stokes numbers
(for which a deviates the most from a,) which suggests
that setting a equal to a, in the closure model is, to
some extent, a cause of the poor prediction in the very
near wall region. For future work, when using the
closure model for particle dispersion in a 3D turbulent
boundary layer, a 1 could be approximated by the
fluid Lagrangian timescale, for which data is readily
available. It is expected that this will produce more
satisfactory results close to the wall because the fluid
Lagrangian timescale is qualitatively more appropriate
than the Eulerian timescale as an approximation to a 1.
A second explanation for the poor performance of
the closure model in the very near wall region is the
specification of q(v, u x). In this work ~, is specified
by using a locally homogeneous approximation (26).
The implication of this is that 0, is always predicted
to be less than 0~, which, for larger particles, is not
true in the very near wall region. The effect of this
in the closure model, for larger particles, would be to
underpredict Ox(s;x). Clearly then more appropriate
specifications for 8, must be developed in order to
improve the closure model, and this will be considered
in future work.
A third possible cause of the discrepancy between the
simulation and closure model results is the approxima
tion that p(x')/p(x) 1 (see appendix 2). In the very
near wall region where there are steep gradients in p(x)
this approximation is clearly not valid. As a next step,
this will be addressed by solving the equation for p(x)
which can be obtained from the the continuum equations
for the particle phase (see Skartlien (2007)). With this
solution, p(x')/p(x) may be more appropriately speci
fied and implemented in the closure model.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Conclusions
The results from this preliminary study are encouraging
and suggest that, with further development, the method
ology outlined here has the potential to provide good ap
proximations to (R(x', t'; x, t))x for particle transport in
a regime exhibiting strong spatial variations in the flow
properties. The intention is to next test the model on a
more complex system, based on LES data, in which the
Eulerian fluid timescales exhibit strong spatial varia
tion. A key issue that needs to be addressed is the formu
lation of more precise estimates for the joint onepoint
pdf q(v, ulx), going beyond the simple Gaussian forms
suggested here. Also, it will be of interest to see what
scope there is to extend the approach and incorporate
other types of particlesurface interactions; the obvious
next case being particle adhesion. Clearly this presents
further challenges in the formulation of suitable approx
imations to p(x', t' x, t).
Acknowledgements
This work was performed as part of the FACE centre; a
research cooperation between IFE, NTNU and SINTEF.
The authors gratefully acknowledge funding through this
centre from the Research Council of Norway and by
the following industrial partners: StatoilHydro ASA,
Norske ConocoPhillips AS, Vetco Gray Scandinavia
AS, Scandpower Petroleum Technology AS, FMC, CD
adapco, ENI Norge AS, Shell Technology Norway AS.
Appendix 1
To test the equivalence of the two systems (12) and
(13) it is necessary to formulate expressions for the La
grangian averages (Aij), (cj) and (xzKj) appearing in
(14). In this appendix expressions are derived for these
averages when the forms of v, A and p are those given
by the LHDI interpretation, (5) and (7). The analysis fol
lows the approach given in Hyland et al. (1999); Reeks
(1991, 1992) and treats 7 as the Green's function for
the particle equation of motion (3). This, in essence, ne
glects gradients in f and requires F to be linear in x and
v. With these caveats, we may write
Xi(t)
vi(t)
t
x(t) + JHm(t;t')fm (x',t') dt'
t
o
v(t) + J im(; t')fm (x', t') dt'
o
where x, v satisfy x = v
that, in general, (f(x, t))
(39)
F(x, v). It should be noted
= 0 = (f(x,t)) 0, and
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
consequently, in general, x (x), and v 0 (v). Now,
with the LHDI interpretation of A
(AifAj.(,v,t))
where
J Aij (x, v, t) p(x, v, t) dx dv
X V
x v
x v 0dv dx
...p(x, v, t) dv dx
t) )X dt'
g(s) = /3 l(1 e )
Then Ox(s; x) = (x'x') is given by
Ox(s;x) g2(s)(vov) + .g(s sl)g(s
00
...(s) f(s2))ds ds2 + 2g(s) g(s s
0
.(f ))ds'
t
=r ix O I [ ( t; t') f{x',t')fj (X, t)6(x
x v 0
...6(v v)) dt'dv dx
t
= Hi t;') m(x',t'x)f ,t) dt'
0
K (x, 3 xi)0
From (44)
f(s11)f(s2))
/32 (uoo)exp
asl 21] (46)
and, using the FurutsuNovikov formula,
(40) (fs')v0 = Pu(s')vo) = uovo) exp [ as] (47)
In a similar fashion the LHDI forms of p and v together
with (39) imply that
{Pij(x, v,t))
((K
S* (41)
xm (m,
Xi (Xm
Appendix 2
Using Bayes' theorem we may write
p(x')
p(x', t'lx, t) = p(x, tx', t') (42)
where p(x') represents the particle phase concentration
at x at t'. The implication of this is that the backward
in time conditional spatial distribution (p(x',t'lx, t))
may be modeled via the forward in time conditional
spatial distribution (p(x, tlx', t')). In this work, as a first
approximation, the ratio p(x')/p(x) is set to unity.
The model for p(x, tlx',t') is constructed using (28)
with 'initial' conditions x(0) x (fixed) and v(0) = vo
a random variable correlated with u(x,t) .",..
Since (in this 1D model) vo and u are both zeromean
variables it follows that (x(s))x = x and then that
x'(s) = (s) x g(s)v0 + g(s s')f(s') ds'
Jo
Substituting these into (45) gives
where (vov), = ".'), 0 = "") and
G, g 2(s)
G 2 (r(1 r2)s (1 + r)g(rs)
r2(1 r2)
+ r(1 + r)(2r 1)g(s) + rg((1 + r)s) (49)
r 2(1 + r)g(2s))
2
Gu g((1 ) (g(rs) rg(s)
r(1 r)
with r = a/3 and T = 31. The formula for the de
generate case r 1 are not presented here, but results
based on expressions (49) converge to those for this spe
cial case as r 1.
References
S Corrsin. Progress report on some turbulent diffusion
research. Adv. Geophys., 6:161184, 1959.
KE Hyland, S McKee, and MW Reeks. Derivation of
a pdf kinetic equation for the transport of particles in
turbulent flows. J. Phys. A: Math. Gen., 32:61696190,
1999.
u(x, t s)
Ox(s;x) = GOv + G0, + Gv,, O
i)fj (xt)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
VI Klyatskin. Dynamics Of Stochastic Systems. Elsevier,
2005.
B Oesterl6 and LI Zaichik. Time scales for predicting
dispersion of arbitrarydensity particles in isotropic tur
bulence. Int. J. Multiphase Flow, 32:838849, 2006.
MW Reeks. On a kinetic equation for the transport of
particles in turbulent flows. Phys. Fluids, 3:446456,
1991.
MW Reeks. On the continuum equations for dispersed
particles in nonuniform flows. Phys. Fluids, 4:1290
1303, 1992.
MW Reeks. On the probability density function equa
tion for particle dispersion in a uniform shear flow.
J. Fluid Mech., 522:263302, 2005.
R Skartlien. Kinetic modeling of particles in stratified
flow Evaluation of dispersion tensors in inhomoge
neous turbulence. Int. J. Multiphase Flow, 33:1006
1022, 2007.
D Swailes and K Darbyshire. A generalized Fokker
Planck equation for particle transport in random media.
Physica A, 242:3848, 1997.
D Swailes and K Darbyshire. Probabilistic models for
particle and scalar transport in fluctuating flows: An
evaluation of simple closure approximations. Physica
A, 262:307327, 1999.
DC Swailes, YA Sergeev, and A Parker. Chapman
Enskog closure approximation in the kinetic theory of
dilute turbulent gasparticulate suspensions. Physica A,
254:517547, 1998.
