Residual sweeping effect in turbulent particle pair diffusion in a Lagrangian diffusion model

.


Introduction
Turbulent particle pair diffusion has attained somewhat of an iconic status in the turbulence community, many researchers having addressed this topic over the decades.Nevertheless, most if not all theories of turbulent particle pair diffusion in homogeneous turbulence with extended inertial ranges have been based upon the hypothesis of locality since Richardson in 1926 [1], and Obukov in 1941 [2].
Richardson pioneered this field and introduced the idea of a scale dependent pair diffusivity as the fundamental quantity of interest in turbulent pair diffusion studies.The turbulent pair diffusivity is defined as, where l(t) is the pair displacement vector at time t, l = |l|, v(l) is the pair relative velocity, and • is the ensemble average over all particle pairs.The locality hypothesis can readily be applied to generalized power law spectra of the type, E(k) ∼ k −p , for 1 < p ≤ 3; the pair diffusivity then scales like K(l, p) ∼ σ γ l p l with γ l p = (1 + p)/2, [3], where σ 2 l = l 2 .For Kolmogorov turbulence p = 5/3, this gives the well known Richardson scaling K ∼ σ 4/3 l , which is equivalent to l 2 ∼ t 3 [2].Kinematic Simulations (KS) [4,5] has often been used to investigate turbulent pair diffusion, even though it does not yield the assumed locality scaling; for p = 5/3, KS gives γ Kol ≈ 1.53 > 4/3.For this reason, it has been widely assumed that KS must be in error.Thomson & Devenish [6] argue that the turbulent pair diffusivity must scale like, K(l(t)) ∼ S(l)τ s (l(t)), (2) where S(l) is the structure function of the turbulence velocity field and τ s (l) is an effective time scale of velocity increments.In real turbulence, assuming locality scaling for S(l) ∼ l 2/3 and for τ s (l) ∼ l 2/3 leads to Richardson's classical scaling K(l) ∼ σ 4/3 l , where we evaluate K at typical values of l, namely σ l = l 2 1/2 which is commonly assumed in these studies.
Thomson & Devenish argue that in KS because of the lack of true dynamical sweeping, the time scale must be a sweeping time scale, namely τ s ∼ l/U s which is a time scale for the large scales flow to cut through smaller local eddies, U s being the sweeping velocity scale.This leads to, K ∼ σ 5/3 l .Even when the rms turbulence velocity u is taken instead of U s , they obtained K ∼ σ 14/9 l .They concluded that whereas locality is true in real turbulence, it is not true in KS.In turbulence the large energy containing eddies carry the smaller eddies, but in KS as there is an absence of true dynamics the large scales force the fluid particles to cut through the smaller eddies in an unphysical manner, a view supported by Nicolleau & Nowakowski [7], and Eyink & Benveniste [8].
However, Thomson & Devenish's scaling argument leading to equation (2) adresses only the scaling laws for K, but does not quantified the actual errors in the diffusivity K in KS -is it large or small ?It is prudent, therefore, to seek an alternative, a more analytic, approach to address this question, which is the main concern of this work.
Here we re-examine the sweeping effect in KS with a view of quantitfying the error in the KS pair diffusivity K s compared to the physical pair diffusivity K.For this purpose, we focus upon the differences in the relative velocities along pairs of particle paths in the sweeping frame of reference.This frame of reference accounts for the physical sweeping effect of the largest energy containing scales; but a residual sweeping effect still remains due to the largest inertial range eddies sweeping the smaller inertial range eddies.Consider Fig. 1 which shows a particle pair with separation l in the inertial subrange being swept by a large scale flow.A real fluid particle pair will be swept by the physical velocity field u and will follow certain particle paths; but a KS flow will transport the pair along neighbouring particle paths due to an additional KS sweeping motion, u s , and thereby force the particles to cut through local flow structures.
The large scale physical sweeping velocity are assumed not to affect the relative motion of particles in a pair in the inertial subrange.The critical question is, are the deviations from the physical trajectories induced by KS in the pair diffusion process large or small?
In the ensuing analysis, it is the error between the KS and the physical pair diffusivities, |K s − K|, that will be calculated.We will consider generalised power law energy spectra, E(k) ∼ k −p , because the analysis must be valid for all such power spectra and this will add weight to the results and conclusions that can be drawn from this work if validated over the wholw range of p considered.
The main questions of interest are, is the KS sweeping error large or small, and in what range of separations?These questions are addressed first through a novel mathematical analysis focussing upon pairs of neighbouring particle trajectories.This is then verified against simulations using KS with very large inertial subranges.
In Section 1, we derive an expression for the error in the pair diffusivity in KS flows by analysing neighbouring trajectories in the swept frame of reference.In Section 2, the KS method is discussed and simulation results presented.In the final Section 3, we discuss the results and its implications for theory and modeling.
1 The normalized error in the pair diffusivity

The numerical timestep error
In the swept frame of reference, the relative motion of a particle pair in the inertial subrange is unaffected by the sweeping action itself.This can be mimicked in KS by setting E(k) = 0 for k < k 1 .However, there still remains a residual sweeping caused by the largest of the inertial scales sweeping the scales local to the pair separation.Consider an ensemble of particle pairs released in a field of homogeneous turbulence at time t = 0 with some small initial separation l 0 .At some time t later, the ensemble average of the separation is assumed to be well inside the inertial subrange and the relative motions are independent of l 0 [9].
Consider the particles in one of these pairs, labeled 1 and 2, as shown in Fig. 1.The particle locations are x 1 (t) and x 2 (t) respectively at time t; and the pair displacement is l(t) = x 2 − x 1 , and l(t) = |x 2 − x 1 |.The physical flow is u(x, t).All quantities are assumed at time t unless otherwise stated.
At time t the additional (or residual) KS sweeping flow u s (x, t) is 'swiched on' -this is not to be confused with the total KS velocity field which is (u + u s )(x, t), see Fig. 1.
The physical flow u(x, t) transports the particles to x 1 (t * ) and x 2 (t * ) respectively at the next time step t * = t + dt; while the KS flow (u + u s )(x, t) transports the particles to x s 1 (t * ) and x s 2 (t * ) respectively.Note that l s = x s 2 − x s 1 , and l s = |x s 2 − x s 1 |.The superscript * will refer to quantities at time t * , e.g.l * = l(t + dt).The superscript s will refer to quantities related to the KS residual sweeping, e.g.l s (t * ) = l s (t + dt).
The following quantities are defined:

is the total KS relative velocity
Simplifying the notation as much as possible, e.g.u 2 = u(x 2 , t), and u * 2 = u(x 2 , t + dt), yields Figure 1.Schematic diagram illustrating the system discussed in the text.The locations of two nearby particles, labelled 1 and 2, are located at x 1 (t) and x 2 (t) respectively, with turbulence velocities u 1 (t) and u 2 (t), at time t; their separation is l(t) = |x 2 − x 1 |.They are transported with velocities (u 1 + u s 1 )(t) and (u 2 + u s 2 )(t) respectively to the new locations x s 1 (t) and x s 2 (t) at the next time step t + dt, as shown.
ṽ(l s * ) is calculated at the new KS swept particle locations.Using Taylor expansions wherever necessary, assuming that the velocity fields are at least twice differentiable in space and at least once in time, The pair diffusivity at time t * is, K * = l * • v(l * ) -we ignore constants of proportionality, like 2, because we are interested only in the power scalings in this work.The KS equivalent is K s * = l s * • ṽ(l s * ) .Using equations ( 5) and ( 6) and ignoring terms of order dt 2 and higher, The timestep error between the KS and physical diffusivities for a given timestep dt is, (7), yields The time scale of the sweeping T s is much larger than the local time scale of the pair separation.Hence, in the last four terms l(t) is replaced by l(t * ) without affecting their magnitudes or scalings (the associated errors are ∼ O(dt 2 ) which is neglected).
All the terms in equation ( 8) are now evaluated at the same time t * , so without loss of generality t * is replaced by t and the superscritp ' * ' is dropped.The subscript ' 1 ' is also dropped because of homogeneity.Equation ( 8) now simplifies to, It is reasonable to assume that the (residual) sweeping flow field u s , which is caused by the largest of the inertial range eddies, is close to uniform across small distances, and therefore the relative velocities across local scales l that it induces is small.However, gradients of the relative velocity v(l) itself can be large.The magnitude u s = |u s | is assumed large compared to v(l) = |v(l)|, and also as compared to v s (l) = |v s (l)|.u s scales differently to v(l).
v(l) also scales differently to v s (l), the former being governed by inertial range turbulence scaling, and the latter by differences in the residual sweeping velocity across a small distance l.This can be seen clearly in the limit of uniform (parallel) sweeping flow, where the v(l) is unaffected, but v s (l) = 0 and all the terms on the right hand side in equation ( 9) are zero except for the second term.This indicates that the second term in equation ( 9) makes the dominant contribution to the error.
Consider generalized energy spectra of the form E(k) = ε 2/3 L 5/3−p k −p , for k 1 ≤ k ≤ k η and for 1 < p ≤ 3, and with k η /k 1 1.In the swept frame of reference, E(k) = 0 for k < k 1 .The rate of energy dissipation is ε ∼ U 3 /L, where U is the velocity scale in the energy containing scales.
The previous discussion implies that |∇v s | U L and therefore, The energy in turbulent inertial scales local to l is, v 2 (l) ∼ E(1/l)/l, and therefore, It is usual to assume the scaling l ∼ σ l as previously mentioned.Then, the second term in equation ( 9) is given by, All the other terms in equation ( 9), labeled respectively E 1 , E 3 , ..., E 7 , scale proportional to u s or are much smaller, and this leads to the following estimates, All of these ratios are small for σ l /L < 1 and 1 < p ≤ 3.This is true even in the expression for E 1 /E 2 because the factor L/U dt is subdued by the very small factor in the brackets.It is reasonable to conclude that the 2nd term in equation ( 9) is dominant and therefore To estimate E 2 itself, an estimate for u s (l) is needed.The inertial subrange contains only about 1% of the the total energy in the turbulence across the entire wavenumber range, 0 < k < ∞, such as in a von Karman spectrum E vk , i.e.E ks ≈ E vk /100.This means that the root mean square turbulent velocity in the inertial range is approximately, u ks ≈ U/10, where The wavenumbers that contribute to the sweeping of particle pairs at separation σ l are in the range k 1 ≤ k < k l , where k l ∼ 1/σ l .In the KS sweeping frame of reference, the larger inertial scales are the sweeping scales, and the energy in these scales is approximately, and using ε ∼ (u ks ) 3 /L ∼ (U/10) 3 /L, this gives Using this in equation (), the residual error between the physical and the KS pair diffusivities in the inertial subrange of pair separations in the swept frame of reference, per unit timestep (retaining E K to represent this quantity), is C k is the constant of proportionality, which can depend upon p.
For p = 5/3, the residual error per unit timestep is, but is nearly zero close to σ l /L = 1.In this limit, p → 1, the pair diffusion is strongly local and is not affected by long range sweeping.
For p = 3, the residual error per unit timestep is negligibly small for σ l /L < 1, In this limit, nearly all the energy is contained in the largest scales and inertial subrange scaling is no longer applicable.

The integrated error
Equations ( 17) provides a way of estimating an upper bound for the integrated residual error.Fig. 2 shows the log-log plots of the factor f (x) = √ 1 − x p−1 x (p−1)/2 for selected powers in the range 0 < p ≤ 3. The range over which √ 1 − x p−1 ≈ 0 is very short and close to x = 1; but over the rest of the range √ 1 − x p−1 ≈ 1.Hence, to a good approximation, f (x) < x (p−1)/2 for all x ≤ 1, and using this in equation ( 17) with x = σ l /L yields, The integrated residual error, E I K , over a period of time is, Assuming the pair separation scaling σ 2 l ∼ t χp , where t is the time and for some χ p > 0, yields, If the pair diffusivity scales like, K ∼ σ γp l , for some γ p > 0 then χ p = 2/(2 − γ p ) is an exact relation.
The most important quantity is the normalised integrated residual error with respect to the pair diffusivity, e I K .Using the above expression for χ p , and replacing the scalig with L by scaling with η, and absorping all ensuing constants in to a new constant A k , leads to For strict locality scaling γ p = (1 + p)/2, and this becomes Since γ p − 1 > 0, e I K decreases with increasing pair separation for all p > 1.For p = 5/3, we have γ p = 4/3 and we obtain, Thus, as σ l /η → ∞ then e I K decreases; and as σ l /η → 0 then e I K increases.Even for non-local scaling, assuming that γ p does not deviate too far from the local scaling, the above order of magnitude for e I K is still approximately true.In fact, in KS we know, Section 2.3, that γ p is slightly greater than locality so the errors will be slightly smaller than in eqution (25).
The magnitude of these errors will also depend upon A k ; but A k cannot be determined from theoretical considerations alone, although it is reasonable to assume that A k is small enough for the e I K to be negligible towards larger separations as σ l /η → ∞ -this assumption will be justified by numerical simulations.Even so this does not guarantee that the error will remain small at smaller separations, since e I K → ∞ as σ l /η → 0. However, between these two asymptotic limits there must exist an intermediate range of separations, between 1 < σ l /η < ∞, where the errors remain negligible.The crucial question is, how wide is this range of scales?If this intermediate range is short then the KS errors will be significant at almost all separations.However, if it is long then the KS errors will be negligble inside this intermediate range where true inertial subrange scaling can be expected to be manifested.
To determine the extent of this intermediate range of scales, if it exists at all, simulations with KS must be performed with very large inertial subranges.

Simulations and Results
The normalised integrated error, e I K , is scale dependent and reduces with increasing separation.The KS diffusivity is given by, It is expected that if there is an appreciable intermediate range where the errors are negligible, then the power scaling in K s must be constant and asymptotic to the limiting case where σ l /η → ∞.The extent of this intermediate range is determined by the range over which the power scaling in K s is constant.Furthermore, significant levels of the sweeping error means that the fluid particles cut through KS eddies, and therefore must be accompied by high levels of noise -the larger the relative sweeping error the larger the noise level.
Thus, where the errors are negligible it is expected that the correct power law scaling, K s ≈ K, will be observed in that part of the of the inertial subrange.
On the other hand, where the errors are significant it is expected that K s will deviate from the true power law scaling for K and also be accompanied with significant statistical noise due to fluid particles being swept through local eddies in that part of the inertial subrange.Even in this case, however, it is expected that the errors and the associated noise diminish as the pair separation increases.

Frames of reference
Comparison will be made between two cases: first, where K s is obtained from KS in the physically correct sweeping frame of reference; and second, the case where large scale random sweeping velocities are explicitly added to the flow.
The analysis for the latter case is similar to that which leads to equation (24), except that the factor of 10 in the denominator is now 1.Both of these cases can therefore be written collectively as, where C = 10 in the swept frame of reference, and C = 1 when large random scales are included in the simulations -the residual errors are an order of magnitude smaller in the swept frame of reference.
Case 1: Swept frame.Set the spectrum to be E(k) ∼ k −p in the inertial subrange, and set E(k) = 0 for k < k 1 .
Case 2: Non-swept frame.Set the spectrum to be E(k) ∼ k −p in the inertial subrange as in Case 1, and add E(k) = E 0 δ(k − k 0 ) at some low wavenumber k 0 < k 1 , and such that E 0 is the energy in the von Karman spectrum in the range 0 < k < k 1 .
A very small fixed timestep, smaller than any timescale in the system dt τ η , is used in all the simulations reported here.τ η is the Kolmogorov time scale.

Kinematic Simulations
Kinematic Simulation [4,5] is a Lagrangian method for particle diffusion in which the velocity fileld is prescribed as a sum of energy-weighted Fourier modes.It is akin to the widely used random flight type of statistical models in which the dynamical interactions between turbulent length scales is not explcitly modeled, rather the overall effect on the statistical moments of particle diffusion is mimicked.In KS this is accomplished by specifying the energy spectrum E(k).KS continues to be used in turbulent diffusion studies for both passive and inertial particle motion, including cases with generalized power-law energy spectra of the form E(k) ∼ k −p for p > 1, Maxey [10], Turfus [11], Fung & Vassilicos [12], Malik & Vassilicos [13].Meneguz & Reeks [14] carried out a DNS of inertial particle motion, and compared it to results from KS which they found to agree well with the DNS.The turbulent difusivity itself can be computed in two ways.Directly from the forumla K(l) ∼ l • v(l) , i.e. the ensemble average of the scalar producted of v and l.But it has been found that using the equivalent formula, K(l) ∼ d l 2 /dt, i.e. the derivative of the l 2 , converges faster statistically needing a much smaller ensemble of trajectories, although the two methods give identical results for large enesmbles of particle trajectories.The latter method has been adopted here.
Lagrangian statistics are the physically meaningful output from KS.It is not correct to compare the kinematically generated flow fields directly with DNS flow fields.As such, KS is like Lagrangian methods such as Random Walk models where an individual particle trajectory has no physical meaning, but the ensemble average over many such random trajectories produces physically meaningful Lagrangian statistics.
In Case 2, large scale random sweeping were added at the low wavenumber k 0 = 1/10, with E(k) = E 0 δ(k − k 0 ).The energy in these sweeping scales, E 0 , was equal to the energy contained in the von Karman turbulence spectrum for k < 1. k 0 corresponds approximately to the location of the peak in the von Karman spectrum.
In both cases, three different power spectra were considered with, respectively, p = 1.01, 5/3, and 3.With 8 pairs released in 5000 flow realizations, the Lagrangian statistics were obtained from 40,000 particle pair trajectories.
Fig. 3 shows log-log plots of the pair diffusivity K/(ηv η ) against σ l /η, for Case 1 (red lines) and Case 2 (green lines).The energies in the two cases are different, so Case 2 plots have been shifted vertically in order to compare the two cases directly.This does not affect the scalings (the slopes) which is the main interest here.Hence the ordinate is shown without scale.
For p = 1.01, the two cases align with a constant power-law scaling, γ 1.01 ≈ 1.07, over most of the inertial subrange of scales and there is very little statistical noise, indicating that e I K (σ l ) 1 at all separations in this part of the inertial subrange in both cases.In this limit, locality is very strong, and the relative motion is unaffected by the long range sweeping.The obtained slope is indeed very close to the exact locality scaling of 1.005.
For p = 5/3 (Kolmogorov turbulence), in Case 1 (red) a clear power-law scaling is observed, γ Kol ≈ 1.53 > 4/3, and very little statistical noise in the range 1 < σ l /η < 10 5 , indicating that e I K (σ l ) 1 in this range of scales.Case 2 (green) deviates increasingly from Case 1 at small inertial separations where it is also accompanied with increasing levels of noise.Nevertheless, the agreement between the two cases for σ l /η > 10 2 is good.
For p = 3, the two cases overlap with a power-law scaling, γ 3 = 2, with almost no statistical noise.In this limit nearly all the energy is in the large scales and inertial range scaling is no longer applicable; rather uniform strained motion with the characteristic slope of 2 is obtained.

Discussion and Conclusions
All the results in Fig. 3 are consistent with the numerical analysis and the theoretical predictions in section 2. The constant power law scaling over most of the inertial subrange of separations, 1 < σ l /η < 10 5 , and the very low level of statistical noise in the swept frame of reference are especially important.(The departure from this for σ l /η < 1 observed in Fig. 3 is outside of the inertial subrange.) The KS sweeping errrors in this frame of reference is thererfore negligible for most practical purposes.It is possible that KS could produce negligible sweeping errors in even bigger intermediate ranges than reported here, but the current simulations are the maximum size of inertial subrange possible, k η /k 1 = 10 6 , with double-precision accuracy.
It is remarkable that even when large scale sweeping is included, the KS sweeping errors remain small in the range σ l /η > 10 2 .
It is also noted that some Direct Numerical Simulation (DNS) show pair diffusion which appear to display locality scaling, see [16] for example.However, the maximum inertial range obtained in DNS to date is around k η /k 1 ≈ 10 2 , which is much shorter than required to test pair diffusion scaling reliably -that requires k η /k 1 > 10 4 .Thus the current KS results cannot be compared directly with DNS at the present time.However, for low Reynolds KS has been validated against DNS for turbulent pair diffusion by Malik & Vassilicos [13]; here not only did the pair diffusion from KS closely match the DNS results with the same energy specturm, but the fourth order statistic, the kurtosis in the pair separation, also matched remarkably well.
The main contribution of this work is that it has been demonstrated that in a reference frame moving with the large energy containing scales the sweeping errors in the turbulent pair diffusion process in KS is negligible in the inertial subrange where, 1 < σ l /η < 10 5 .This is significant not only because it amends the previous theory of [6][7][8] , but it is important for turbulent diffusion modeling and applications in general, and especially for pair diffusion studies.
If the sweeping errors are negligible, why is locality scaling not oberved for p = 5/3 where KS yields γ Kol ≈ 1.53 > 4/3?There are two possible explanations.First, other errors could be present in KS, not due to the sweeping, but to as yet unknown sources; but no one has proposed any such source of error in, and so this remains very speculative.
The second possibility is that the locality hypothesis itself may be error.This goes against the widely accepted theory of Richardson [1], and no one has proposed alternative theory.On the other hand, it should be noted that the locality scaling for the pair diffusion has never been confirmed unequivocally as noted by Salazar & Collins [17] -so there is room for new thinking in this field.
These are currently the subjects of active research by the author.What can be said for certain at the present time is that the departure from locality scaling observed in KS cannot be attributed to the sweeping effect.