http://orcid.org/0000-0001-5053-9360
Figures
Abstract
A theory for turbulent particle pair diffusion in the inertial subrange [Malik NA, PLoS ONE 13(10):e0202940 (2018)] is investigated numerically using a Lagrangian diffusion model, Kinematic Simulations [Kraichnan RH, Phys. Fluids 13:22 (1970); Malik NA, PLoS ONE 12(12):e0189917 (2017)]. All predictions of the theory are observed in flow fields with generalised energy spectra of the type, E(k) ∼ k−p. Most importantly, two non-Richardson regimes are observed: for short inertial subrange of size 102 the simulations yield quasi-local regimes for the pair diffusion coefficient, ; and for asymptotically infinite inertial subrange the simulations yield non-local regimes
, with γ intermediate between the purely local scaling γl = (1 + p)/2 and the purely non-local scaling γnl = 2. For intermittent turbulence spectra, E(k) ∼ k−1.72, the simulations yield
, in agreement with the revised 1926 dataset
[Richardson LF, Proc. Roy. Soc. Lond. A 100:709 (1926); Malik NA, PLoS ONE 13(10):e0202940 (2018)]. These results lend support to the physical picture proposed in the new theory that turbulent diffusion in the inertial subrange is governed by both local and non-local diffusion transport processes.
Citation: Malik NA (2019) Turbulent particle pair diffusion: Numerical simulations. PLoS ONE 14(5): e0216207. https://doi.org/10.1371/journal.pone.0216207
Editor: Roi Gurka, Coastal Carolina University, UNITED STATES
Received: September 17, 2018; Accepted: April 16, 2019; Published: May 20, 2019
Copyright: © 2019 Nadeem A. Malik. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the manuscript and its Supporting Information files.
Funding: NAM received funding from SABIC (Saudi Arabian Basic Industries Corporation) through grant number SB101011. The funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: This work was funded by Saudi Arabian Basic Industries Corporation (SABIC) through project number SB101011. This does not alter the author’s adherence to PLOS ONE policies on sharing data and materials.
1 Introduction
Turbulent transport and mixing play an essential role in many natural and industrial processes [1–9], where concentration fluctuations, which is related to the pair separation, often play a critical role. Most theories of turbulent particle pair diffusion assume Richardson’s locality hypothesis [10, 11]. However, a new theory for turbulent particle pair diffusion based on the physical picture that both local and non-local diffusional processes govern the pair diffusion process has been proposed by the author in [12].
The main prediction of the new theory in [12] is the existence of two non-Richardson regimes in inertial subranges with generalised energy spectra E(k) ∼ k−p: a quasi-local (i.e. approximately local) regime at moderate inertial subrange size, Rk ≪ ∞, where ; and a non-local regime at asymptotically infinite inertial subrange size, Rk → ∞, where
and < (1 + p)/2 < γ(p) < 2, and 1 < p ≤ 3.
The power laws, γ(p), cannot be determined from the theory alone. For this, and for investigating other predictions of the theory, ideally we would need experiments or Direct Numerical Simulations which resolves all the time and length scales in the turbulence. At the current time the capabilities of both experiments and DNS are far from being able to examine large inertial subranges. Nevertheless, we can make some progress through the use of diffusion models.
Here, the aim is to investigate the predictions of the new theory numerically using Kinematic Simulations (KS) which is a Lagrangian diffusion model and it has the advantage that it can generate very large inertial subranges which is essential to test the new theory. KS has been used extensively in the past for pair diffusion studies (see Section 3 below). KS can be indicative of the statistical scaling laws in the inertial subrange.
Dissipation range pair diffusion, where the initial particle separation is such that l0 ≪ η, is outside the scope of the present work. A theory for dissipation range release would likely yield different scaling laws because the ensemble of particle pairs would attain a distribution of separations and enter the inertial subrange at different times after release. The subsequent inertial range diffusion would in effect be an average over different virtual release times. Furthermore, in the present KS we do not pose a dissipation range energy spectrum which is essential for investigating dissipation scale diffusion.
In Section 2 we summaries the new theory developed in [12], and in Section 3 we describe the KS method. In Section 4.1 we describe the simulation results for large inertial subranges of size Rk = 106. In Section 4.2 the results for small inertial subranges Rk = 102 are described. In Section 4.3 the simulation results for the transition in the scaling laws in K as the size of the inertial subrange increases from small to very large values are described. In Section 5, cases for very small initial separations l0 ≪ η in order to isolate and expose the non-local diffusional process is investigated. In Section 6 we estimate the numerical errors in the results. We discuss and draw conclusions from the results, and speculate its significane to the general theory of turbulence in Section 7.
2 Summary of the new non-local theory
In order to characterize the pair diffusion process, Richardson assumed a scale dependent pair diffusion coefficient (turbulent diffusivity), because convective gusts of winds increase the pair separation at different rates depending on the separation [13–15]. In 1926, from observational data of turbulent pair diffusion coefficients collected from different sources, he assumed an approximate constant power law fit to the data, K(l) ∼ l4/3, [10]. This is equivalent to 〈l2〉 ∼ t3 [11, 16], and is often referred to as the Richardson-Obukov t3-regime. l(t) is the pair separation at time t and the angled brackets is the ensemble average over particle pairs. Note that the assumed 4/3-scaling is consistent with Kolmogorov turbulence K41 theory, see [12].
However, to date the validity of Richardson’s scaling law has not been established. The general consensus among scientists in the field at the current time is that the collection of observational data, experimental data, and Direct Numerical Simulation, suggests a convergence towards a Richard-Obukov locality scaling, but as noted by Salazar and Collins, “.. there has not been an experiment that has unequivocally confirmed R-O scaling over a broad-enough range of time and with sufficient accuracy” [17]
Furthermore, in [12] it was noted that one of the data-points used in Richardson comes from molecular diffusion studies and should be dis-regarded. The remaining data-points are sound, coming from geophysical turbulence settings containing extended inertial subranges. The line of best fit to this improved dataset displays an unequivocal non-local scaling, K ∼ l1.564; see Fig 1 in [12].
This indicates that non-local diffusional processes cannot be ignored a priori in a general theory of turbulent pair diffusion, which motivates the development of the present local-non-local theory.
2.1 The new theory
The problem is to determine the pair diffusion coefficient (diffusivity), K = 〈l · v〉, of an ensemble of pairs of fluid particles in a field of homogeneous turbulence with an energy density spectrum, E(k) containing an generalised inertial subrange, E(k) ∼ k−p, k1 ≤ k ≤ kη, for 1 < p ≤ 3, and such that E(k) → 0 as k → 0. The particles in a pair are located at x1(t) and x2(t) at time t, the pair displacement vector is l(t) = x2(t) − x1(t), and the pair separation is . The initial separation at some earlier time, t0, is denoted by l0 = |x2(t0) − x1(t0)|. The turbulent velocity field is, u(x, t), and the particle velocities at time t are, respectively, u1(t) = u(x1(t), t) and u2(t) = u(x2(t), t), and the pair relative velocity is v(t) = u2(t) − u1(t).
We assume point source release, which in practical terms means that the initial pair separation must be close to the Kolmogorov length scale, l0 ≈ η. Without loss of generality, it will also be assumed that, t0 = 0.
We define the size of the inertial subrange Rk to be,
(1)
where the inertial subrange part of the energy spectrum E(k) is defined in the wavenumber range k1 ≤ k ≤ kη.
We follow the usual convention and evaluate K at typical values of, l, namely at . Thus, the locality scaling is replaced by,
.
The theory is developed through a novel mathematical approach of decomposing the pair relative velocity, v(l) = dl/dt, as a Fourier integral. Assuming homogeneity, the ensemble average of the scalar product of l with v yields a Fourier decomposition of the diffusion coefficient itself,
(2)
where the integration is over the range of turbulent wavenumbers, and A(k) is the Fourier coefficient of the Eulerian velocity field u(x).
Let kl = 1/l be the pair separation wavenumber, so that
(3)
Note that kl(t) changes with time.
The new theory assumes the existence of two broadly independent diffusional processes within the inertial subrange that contribute to the pair diffusion process as a whole. The two transport processess are correlated locally and non-locally in space, respectively. This can also be expressed in terms of equivalent wavenumbers; thus, each transport process acts from its own range of wavenumbers relative to kl, labeled l (local) and nl (non-local):
- l: a local diffusion process operates at wavenumbers that are local to kl, say in the range
and such that k ≈ kl, and |k · l| ≈ 1. Within this range of wavenumbers local scalings apply.
- nl: a non-local diffusion process operates at wavenumbers that are non-local to kl, say in the range in the range
, and such that k ≪ kl, and |k · l| ⪡ 1. Non-local scalings apply in this range of wavenumbers.
is some arbitrary wavenumber that separates the two processes, and
. (A third transport process at very small wavenumbers k ≪ kl has been shown to be negligible, [12]).
This implies that the integral in Eq (2) is the sum of two integrals over different wavenumber ranges which are defined as follows:
(4)
which we rephrase as,
(5)
We assume a generalised inverse power-law energy spectrum in the inertial subrange in the swept frame of reference,
(6)
where Ck is a constant. A length scale L is necessary for dimensional consistency. L scales with some length scale that is characteristic of the large energy scales, such as the integral length scale, or the Taylor length scale. With this spectrum, and some closure assumptions detailed in [12], Eq (4) becomes,
(7)
The locality scaling is obtained by removing the non-local integral in this equation. This yields,
(8)
and Fl < 1 is a constant. For Kolmogorov turbulence, p = 5/3, this gives,
, which recovers the Richardson’s 4/3-scaling law, [18–20], which validates the derivation.
The non-local scaling in the first integral can be obtained in a similar manner. If the upper end of the inertial subrange is assumed to scale with the large scales, k1 ∼ 1/L, then this yields,
(9) γnl(p) is the non-locality scaling and it is always equal to 2 independent of p. This corresponds to strain dominated motion with
. Snl is a constant.
The overall expression for the turbulent pair diffusion coefficient is therefore,
(10)
or simply,
(11)
(12)
(13)
To complete the analysis we need to obtain an experssion for the balance of the local and non-local diffusion terms, defined to be ratio,
(14)
The quantity
(15)
is the size of the inertial subrange relative to the pair separation in wavenumber space. The quantity C is,
(16) C is an important quantity, which we call the locality kernel because it defines the extent (or size), relative to the pair separation, where local scaling is expected to dominate.
In [12] it was shown that under fairly weak restrictions on C—essentially little more than smoothness of C as a fucntion of p and Rl—over a wide range of smooth test functions for C, in the critical range of p close to Kolmogorov’s p = 5/3 the balance of local and non-local diffusional processes is close to unity,
(17)
This indicates that the assumption that both local and non-local diffusional processses are significant is physically reasonable. (Fl was chosen to be Fl = 0.25 in Eq (14)).
We assume an extension of Richardson’s second hypothesis, that in the limit of infinite inertial subrange, for every p the diffusion coefficient is described by a single power, say γ(p), across all scales. Thus, from Eq (11) K(l, p) must be a power law scaling which is intermediate between the local and non-local scalings,
(18)
Only in some special asymptotic limits do we obtain significantly different behaviour. Firstly, as p → 1 then MK ≪ 1, and therefore Knl ≪ Kl, yielding the locality limit,
(19)
Secondly, as p → 3 then MK ≫ 1, and therefore Knl ≫ Kl, yielding the non-locality limit,
(20)
Thus, as p → 1 then γ(p) → 1; and as p → 3 then γ(p) → 2. Furthermore, γ(p) must transform smoothly between these limiting cases as p goes from 1 to 3. Globally, 1 < γ(p) ≤ 2.
Eq (19) is the main prediction of the new theory, and it is equivalent to the mean square separation scaling,
(21)
This is a non-linear relation between γ and χ, so small changes in γ could produce large changes in χ.
Some specific results can be obtained from these sclaing laws. For Kolmogorov turbulence, E ∼ k−5/3, the new theory predicts that, γ > 4/3, and χ > 3.
For turbulence with intermittency μI > 0, such that E ∼ k−(5/3+μI), the scaling is,
(22)
Under the locality assumption, in real turbulence with intermittency p = 1.72 we should obtain the scaling γl = 1.36, and χl = 3.125. Thus, the classical RO-t3 regime does not actually exist! Richardson’s 1926 dataset, Fig 1, from real geophysical turbulence (i.e. including intermittency) suggests a scaling of, , i.e.
.
For clarity only 10 of the 25 cases in Table 1 are shown: p = 1.01, 1.1,1.3, 1.5, 5/3, 1.9,2.2, 2.5, 2.9, 3.0. Solid black lines of slopes 1 and 2 are shown for comparison.
A non-local 4/3-law, , is precited for some spectrum,
, where p* < 5/3 and γ(p*) = 4/3. This is equivalent to 〈l2〉 ∼ t3, with χ(p*) = 3, which is a new non-Richardson-Obukhov t3-regime for the mean square separation.
We define Mγ(p) to be the ratio of the scaling power γ(p) to and local scaling powers γl(p),
(23) Mγ(p) is equal to 1 at both p = 1 and p = 3, and since Mγ > 1 in the range 1 < p < 3, then there must be a maximum in Mγ at some p = pm for an energy spectrum
, where 1 < pm < 3.
The sketches in Figs 13 and 14 in [12] summarise the preditions from the new theory for, respectively, the non-local regimes in the limit of asymptotically infinite inertial subranges Rk → ∞, and the quasi-local regimes in the limit of short finite inertial subranges Rk ≈ 102.
Finally, we remark that, ultra-violet corrections from the high wavenumber close to kη, and infra-red corrections from the low wavenumbers close to k1 may modifying some of the scalings law specially in the limit of small inertial subranges.
2.2 Exposing the local and non-local processes
It may be possible to isolate the individual local and non-local processes in different limts. From the expression in Eq 14 it is evident that the analysis is valid only if Rl > C. The non-local process, which is the first term in Eq (11), exists in the wavenumber range ; but suppose that the inertial subrange itself is so small that it is close to the size of the locality kernel itself, Rk ≈ C? This corresponds to taking the limit Rl/C → 1 in Eq (14), which implies that Mk → 0 which is a locality dominated limit. Thus quasi-local diffusion regimes can be recovered in the non-Richardson limit of small but finite inertial subrange because of the absence of non-local processes,
(24)
As we progressively increase the size of the inertial subrange we would expect to see a smooth transition from the locality regime at moderate inertial subrange to the non-locality regime at very large (effectively infinite) subrange, as illustrated in the sketch in Fig 15 in [12].
We can also ‘turn off’ the local process in Eq (11) by simply removing the ‘local’ part of the specturm—this is equivalent to taking a very small initial separation l0 ≪ η. Then there is a spectral gap between the k0 = 1/l0 and kη ≪ k0 where E(k) = 0. If this gap is large enough then the spectrum between k1 ≤ k ≤ kη will in efffect be non-local to the pair separation process so long as σl(t) ≪ η, and we should then observe pure strain dominated pair separation
(25)
for all p, which is equivalent to exponential growth in time, σl(t) ∼ l0 exp(St), where S depends upon the form of E(k).
In the next sections we examine the predictions of the new theory using Kinematic Simulations.
3 Numerical simulations
3.1 Particle trajectories using KS
In this study, the Lagrangian diffusion model Kinematic Simulations (KS) was used to obtain the statistics of particle pair diffusion. In KS one specifies the second order Eulerian structure function through the power spectrum, [14, 21]. KS is useful here because it can generate very large inertial subranges sufficient to test pair diffusion scaling laws.
KS generates turbulent-like non-Markovian particle trajectories by releasing particles in flow fields which are prescribed as sums of energy-weighted random Fourier modes. By construction, the velocity fields are incompressible and the energy is distributed among the different modes by a prescribed Eulerian energy spectrum, E(k). The essential idea behind KS is that the flow structures in it—eddying, straining, and streaming zones—are similar to those observed in turbulent flows, although not precisely the same, and this is sufficient to generate turbulent-like particle trajectories [14].
KS has been used to examine single particle diffusion [22], and pair diffusion [14, 20, 23–25]. KS has also been used in studies of turbulent diffusion of inertial particles [26–30]. Meneguz & Reeks found that the statistics of the inertial particle segregation in KS generated flow fields for statistically homogeneous isotropic flow fields are similar to those generated by DNS. KS has also been used as a sub-grid scale model for small scale turbulence [31].
3.2 The KS velocity fields and energy spectra
An individual Eulerian turbulent flow field realization in KS is generated as a truncated Fourier series,
(26)
where Nk is the number of representative wavenumbers, typically hundreds for very long spectral ranges, Rk ≫ 1.
is a random unit vector;
and
. The coefficients An and Bn are chosen such that their orientations are randomly distributed in space and uncorrelated with any other Fourier coefficient or wavenumber, and their amplitudes are determined by
, where E(k) is the energy spectrum in some wavenumber range k1 ≤ k ≤ kη. The angled brackets 〈⋅〉 denotes the ensemble average over space and over many random flow fields. The associated frequencies are proportional to the eddy-turnover frequencies,
. There is some freedom in the choice of λ, so long as 0 ≤ λ < 1. The construction in Eq (26) ensures that the Fourier coefficients are normal to their wavevector which automatically ensures incompressibility of each flow realization, ∇ ⋅ u = 0. The flow field ensemble generated in this manner is statistically homogeneous, isotropic, and stationary.
Unlike some other Lagrangian methods, by generating entire kinematic flow fields in which particles are tracked KS does not suffer from the crossing-trajectories error which is caused when two fluid particles occupy the same location at the same time in violation of incompressibility; but because KS flow fields are incompressible by construction this error is eliminated.
The flow at a point in a KS flow field is irregular because of the presence of flow structures such as vortex tubes and a probe will experience irregular alternations between high levels of fluctuations and low levels of fluctuations. However, there is no dynamical energy transfer between different scales of motion so this type of ‘intermittency’ is at a formal level different to real turbulence, [14]. Nevertheless, it is of considerable interest to investigate how Lagrangian statistical scalings change as we adjust the energy spectrum E(k) ∼ k−p such that it mimics intermittent-like spectra with p ≠ 5/3. This is especially important because KS pair diffusion statistics, including the flatness factor, have been found to be in close agreement with DNS at low Reynolds numbers [25].
The energy spectrum E(k) can be chosen freely within a finite range of scales, even a piecewise continuous spectrum, or an isolated single mode are possible. To incorporate the effect of large scale sweeping of the inertial scales by the energy containing scales, the simulations are carried out in the swept frame of reference by setting E(k) = 0 in the largest scales, for k < k1, and an inverse power spectrum in the inertial subrange as discussed in [32],
(27)
where Ck is a constant. The Kolmogorov micro-scale is η = 2π/kη. L is some large outer length scale (such as the integral length scale, or a Taylor length scale). The inertial subranges sizes is defined by Eq (1).
A particle trajectory, xL(t), is obtained by integrating the Lagrangian velocity, uL(t), in time,
(28)
Pairs of trajectories are harvested from a large ensemble of flow realizations and pair statistics are then obtained from it for analysis.
The spectrum in (27) is normalized such that the total energy density contained in any flow realization is 3u′/2, where u′ is the rms turbulent velocity fluctuations in each direction.
In the current simulations, k1 = 1, L = 1, Ck = 1.5 (Kolmogorov constant) and u′ = 1. Then this yields,
(29) p = 1 is a singular limit which is not consider here. With Eq (29), vη = (εη)1/3 is the velocity micro-scale, and tη = ε−1/3η2/3 is the time micro-scale.
The distribution of the wavenumbers is geometric, kn = rn−1 k1, and the common ratio is r = (kη/k1)1/(Nk−1). The increment in wavenumber-space of the n’th wavenumber is . The frequency factor is chosen to be, λ = 0.5, which is typical in many KS studies. A choice of λ < 1 does not affect the scaling in the pair diffusivity, even frozen turbulence, λ = 0, has been found to yield the same scaling, which has been attributed to the open streamline topology of streamlines in 3D flows, [20].
The particle trajectories xL(t) were obtained by integrating Eq (28) using a 4th order Adams-Bashforth predictor-corrector method (4th order Runga-Kutta gives identical results). Thomson & Devenish [33] used a variable time step Δt that was small compared to the turnover time of an eddy of the size of the instantaneous pair separation, but larger than the turbulence micro-scale tη. While this speeds up the turnaround time of the calculations, here we want to avoid any extra assumptions so that unambiguous conclusions can be drawn from the results. Therefore, in all of the current simulations a very small but fixed time step, Δt ≪ tη has been used. This has the further advantage of avoiding any smoothening of the particle trajectories that is necessary when using variable time steps.
Eight pairs of a particles were released in each flow realization, placed far enough apart for each pair to be independent. It is crucial to run over a large number of different flow realizations, otherwise the statistics will not converge. Typically the ensemble was around 4000 flow realizations, yielding a harvest of 32, 000 particle pair trajectories per simulation. A simulation was run for about one large timescale, T = 2π/k1, which required around 106 time steps for Rk = 106, and about 103 time steps for Rk = 101. In most of the simulations the initial separation was l0 ≈ η; but for p → 3 the energy in the small scales even in the inertial subrange is so small that the particles take a long time to move apart significantly, so in order to accelerate the simulations for these cases we took l0 ≫ η.
4 Simulation results
4.1 Infinite inertial subrange Rk → ∞ (Re → ∞)
In this first set of simulations we investigate the theory for pair diffusion in asymptotically infinite inertial subrange, (infinite Reynolds number), and for this we take Rk = 106.
The spectra in Eq (27) were taken as input to the KS simulations. It is important to simulate cases over the whole range of energy spectra 1 < p ≤ 3, in order to examine the new theory comprehensively; 25 cases of p were selected in the range, 1 < p ≤ 3. The case p = 1 is singular, but p can be taken close to this limit; the smallest value of p chosen was, 1.01 (Table 1).
The pair diffusion coefficient is, , and the locality scaling is, γl(p) = (1 + p)/2. The mean square separation is, 〈l2〉 = tχ(p), where χ(p) = 2/(2 − γ(p)), and the locality scaling is, χl(p) = 2/(2 − γl(p)). The ratio of the power scalings is, Mγ(p) = γ(p)/γl(p).
The results obtained from the simulations for Rk = 106 are summarized in Table 1. This table shows γ(p) from the simulations in column 2; and the locality scalings γl(p) = (1 + p)/2 in column 3. For the mean square separation laws, 〈l2〉 ∼ tχ, the χ(p) = 2/(2 − γ(p)) is shown in column 4; and the locality scalings in column 5. The ratio of scaling powers,
(30)
is shown in column 6.
Log-log plots of the turbulent pair diffusion coefficient K/ηvη against the separation σl/η for different energy spectra E(k) ∼ k−p display clear power law scalings of the form,
(31)
over wide ranges of the separation inside the inertial subrange, Fig 1. The γ(p)’s are the slopes of the plots in Fig 1; these are more easily observed by plotting the compensated diffusion coefficient
against the separation, Fig 2. The γ(p)’s are obtained as the powers that give horizontal lines over the longest range of separation, a procedure that determines γ(p) to within 1% error for most p, except near the singular limit, p = 1, where the errors are around 6%, (see Section 4).
For clarity some of the plots have been spaced out vertically, hence no scale is shown on the ordinate; this does not affect the scalings.
Fig 3 shows the plots of γ(p) (black filled circles) and γl(p) = (1 + p)/2 (dotted blue line) against p. It is observed that the scaling powers, γ(p), are in the range, (1 + p)/2 < γ(p) ≤ 2, and as p → 1, then γ(p) → 1; and as p → 3, then γ(p) → 2. Furthermore, the, γ(p)’s, display a smooth transition between the asymptotic limits at p = 1 and 3.
The black filled circles are the γ(p)’s from the simulations. The dotted blue line is the locality scaling γl(p) = (1 + p)/2. The cyan squares are the ratios, Mγ = γ(p)/γl(p) (right hand scale). See Table 1.
The range of separation over which the scalings laws, , exist progressively reduce from both ends as p → 3. Ultra-violet corrections reduces the range from below due to the diminishing energies contained in the smallest scales so that the pair separation penetrates further into the inertial subrange before it ‘forgets’ the initial separation. Infra-red corrections reduces the range from above because the long range correlations are increasingly stronger as p → 3 causing the particles in a pair to become independent at earlier times and at smaller separation.
The plot of Mγ against p in Fig 3 (right hand scale) shows a peak at pm ≈ 1.8, where Mγ(pm) ≈ 1.15. Mγ remains close to this peak value in a neighbourhood of pm, in the range 1.5 < p < 2.
Fig 4 shows the simulation results of the scaling powers, χ(p), (filled black circles), and χl(p) = 4/(3 − p) (dotted blue line), against p. The same is shown in Fig 5, but focussing in on the range 1.5 < p < 2 which covers the intermittent turbulence range which is indicated by the two red vertical lines.
The black filled circles are, χ(p)’s from the simulations. The dotted blue line is the locality scaling, χl(p) = 4/(3 − p). See Table 1.
Similar to Fig 4, but focussing in on the range of intermittency spectra. The two vertical red lines are the lower and upper bounds for intermittent turbulence spectra, respectively p = 1.70 and 1.74.
Turbulence with intermittency has spectrum E(k) ∼ k−(5/3 + μI). Three extra cases with, pI = 5/3+μI, where therefore simulated. The currently accepted value for the intermittency lies in the range, 0.025 < μI < 0.075, [34–37]; and three cases that cover this range are, pI = 1.70, 1.72, and 1.74. For these spectra, the simulations produced, (Fig 3 and Table 1),
(32)
The mid-point in the intermittency range is close to, pI = 1.72, which gives the scaling law , which is an error of about 0.5% in the scaling power compared to the Richardson’s 1926 revised data as shown in Fig 1 in [12]. For, pI = 1.70, we obtain
, and the error in the scaling power is 1.2%; for pI = 1.74, we obtain
, and the error in the scaling power is 0.4%.
The corresponding scalings for the mean square separation, 〈l2〉 ∼ tχ(p), for the three intermittent spectra are respectively, ∼ t4.396, ∼ t4.505 and ∼ t4.651, Fig 5 and Table 1.
The Kolmogorov spectrum p = 5/3 produces, , and the error in the scaling power is 2.5%. Richardson’s 4/3-law,
, is 15% in error.
The simulations show a non-Richardson 4/3-scaling, , for the the spectrum E(k) ∼ k−1.4, where p* ≈ 1.4, (Fig 3 and Table 1); this yields a new non-R-O t3-regime, 〈l2〉 ∼ t3.
4.2 Short inertial subrange Rk ≪ ∞ (Re ≪ ∞)
In this section, we investigate the new theory in [12] in the limit of short inertial subrange. According to the new theory short subranges could isolate and expose the local diffusional process.
Simulations were carried out for the same cases of p as in Section 3.1, but now for a short inertial subrange of RK = 102. The results are summarised in Table 2 which shows p, γ(p), χ(p), and Mγ(p). Figs 6 to 10 are the counterparts to Figs 1 to 5 but for Rk = 102.
The pair diffusion coefficient is, . The mean square separation is, 〈l2〉 = tχ(p), where χ(p) = 2/(2 − γ(p)). The ratio of the power scalings is, Mγ(p) = γ(p)/γl(p).
For clarity only 10 of the 25 cases in Table 2 are shown: p = 1.01, 1.1,1.3, 1.5, 5/3,1.9,2.2,2.5, 2.9,3.0. Solid black lines of slopes 1 and 2 are shown for comparison.
For clarity some of the plots have been spaced out vertically, hence no scale is shown on the ordinate; this does not affect the scalings.
The black filled circles are the γ(p)’s from the simulations. The dotted blue line is the locality scaling γl(p) = (1 + p)/2. The cyan squares are the ratios, Mγ = γ(p)/γl(p) (right hand scale). See Table 2.
The black filled circles are, χ(p)’s from the simulations. The dotted blue line is the locality scaling, χl(p) = 4/(3 − p). See Table 2.
Similar to Fig 8, but focussing in the range of real turbulence with intermittency. The two vertical red lines are the lower and upper bounds for intermittent turbulence spectra, respectively p = 1.70 and 1.74.
Fig 6 shows the log-log plots of the diffusion coefficient K/ηvη against σl/η for the same 10 selected cases of p as in Fig 1. Fig 7 shows the corresponding log-log plots of the compensated diffusion coefficient, against σl/η—this time the diffusion coefficient is compensated by the locality scaling
. The plots are close to horizontal, though not quite.
Table 2 shows the γ(p)’s obtained from the simulations, which are also plotted in Fig 8 (c.f. Fig 3). The γ(p)’s are close to the locality scalings γl(p) = (1 + p)/2, as demostrated in the plot of the ratio, Mγ against p also in Fig 8. Although they differ by about 10% close to p = 1, most of this error is due to the singularity at this point, see Section 4.
In the critical range of p close to p = 5/3, which is far from the singular limit, the maximum difference from locality is about 5% which is a real physical effect most likely due to the fact that even in this short subrange there may be some small non-local straining effect. There is also infra-red corrections close to k = k1, and ultra-violet corrections close to k = kη which will penetrate some way in to the inertial subrange from both ends thus reducing the effective range over which the scalings can be observed. In short inertial subranges these end-of-range corrections will be appear relatively greater than in large inertial subranges.
As in the latter case, as p → 3, the range of separations over which the locality scaling is observed diminishes from both ends of the domain, and this probably leads to the sub-local diffusion, γ(p) < γl(p), for p > 2.2, Fig 8.
Fig 9 shows the plots of χ(p) and the locality scalings χl(p) against p for comparison, (c.f. Fig (4)). Fig 10 is the same except zooming in on the intermittency range (c.f. Fig 5).
From the simulations, the scalings obtained for the pair diffusivity are: for p = 5/3;
for p = 1.70;
for p = 1.72;
for p = 1.74.
For the same spectra, the scalings for 〈l2〉 are respectively, ∼ t3.39 ∼ t3.449 ∼ t3.509 and ∼ t3.571, Fig 10 and Table 2. Thus, a true R-O t3-regime cannot exist in reality.
Overall, the simulation results for Rk = 102 are always close to locality, and to within 5% in the range of spectra close to real turbulent spectra—we call these quasi-local regimes.
4.3 Scaling law as Rk increases
The theory in [12] predicts a smooth transition from local to non-local regimes as Rk increases. To examine this, simulations were carried out for progressively larger inertial subranges, from Rk = 101 to 106, for two selected energy spectra, namely for Kolmogorov spectrum, E(k) ∼ k−5/3, and for an intermittent spectrum, E(k) ∼ k−1.72.
Fig 11 shows the log-log plots of the diffusion coefficient K/ηvη against σl/η for p = 5/3 for the six cases of inertial subrange size considered, Table 3. Fig 12 is similar but for p = 1.72. Both sets of plots show similar trends; for the smallest Rk = 101 we obtain γ ≈ 1.1 ≪ 4/3—clearly the subrange is too short to manifest any kind of genuine inertial range scaling.
Lines of slope 4/3 and 1.525 are shown for comparison.
Lines of slope 1.36 and 1.556 are shown for comparison.
The cases Rk = 102 are close to locality scalings, as already discussed in Section 3.2.
As Rk increases further the γ(p)’s increase asymptotically towards the limiting values, γ → 1.525 for p = 5/3, and γ → 1.556 for p = 1.72. This is clearly seen in Figs 13 and 14 which show the log-linear plots of γ against log(Rk).
5 Exposing the non-local process
In [12] it was hypothesised that the non-local process could also be isolated and exposed by taking a very small initial separation l0 ≪ η. Then the early motion should be purely strain dominated relative motion so long as σl(t) ≪ η. It was also noted that this regime should be independent of the form of E(k) and also of the size of the inertial subrange. To examine this hypothesis we therefore only need to test a few cases to prove generality.
Thus we ran simulations for three sizes of the subrange, Rk = 101, 102, and 103 for the same spectrum p = 5/3, and in each case we took l0/η = 10−3 respectively. Fig 15 shows log-log plot of K/ηvη against σl/η for these cases. We also ran simulations for different spectra p = 1.4, 5/3, and 2.0 for the same inertial subrange size Rk = 102, and the results are shown in Fig 16.
The initial separations in each case is l0/η = 10−3. A line of slope 2 is shown for comparison.
The initial separations in each case is l0/η = 10−3. A line of slope 2 is shown for comparison.
In all the above cases, the results display clear scalings which is the signature of strained motion. Importantly, these regimes are valid till about σl/η ≈ 10−2 before they start to bend away from from a slope of 2 as the inertial subrange is approached from below and inertial subrange scaling begins to have its impact. This gives us a rough estimate of the size of the locality kernel which must be of the order of C ≈ 102.
6 Estimate of numerical errors
The numerical results presented here are the most comprehensive obtained to date from KS due to the very large ensemble of particle pairs and the small time steps used. The statistical fluctuations in the results are therefore small. The γ(p)’s, which are the slopes of the plots in Fig 5, can be determined to within 1% error. An exception is close to the singular limit p = 1 where the numerical errors can be large. An accurate estimate of this error can be obtained as follows, noting first that the error level in γ(p) is identical to the error level in Mγ.
As p → 1, then Mγ → 1; but very close to this limit, Mγ ≈ 1 is still a good approximation. For the Rk = 106 case, and for p = 1.01, KS produces, Mγ ≈ 1.06 (Fig 7 and Table 1). This is an error of 6% which is small considering that it is so close to the singular limit. An error of around 1% away from p = 1 is therefore reasonable. In the limit p = 3 there is no detectible error in Mγ from KS to three decimal places (Table 2).
KS is an established method used by many researchers in turbulent diffusion studies, as noted in the earlier references in this paper. However, some researchers [23, 33, 38] have argued that KS suffers from systematic erros due to the lack of true dynamical sweep of the small inertial range scales by the much larger energy containing convective scales. However, the author has recently investigated this issue in [32] and shown that such conclusions are unfounded because they were made under the assumption of locality being true. Through an analysis of pairs of trajectories the quantitative errors in KS due to the sweeping effect were proven to be negligible provided that the simulations are carried out in a frame of reference moving with the large scale sweeping flow. The scalings obtained from KS are therefore genuine, and not errors as previously thought.
7 Discussion and conclusions
Richardson conceived his scaling law to be applicable to real turbulence, not just a mathematical curiosity. The new theory, developed in [12] and tested numerically here, generalises Richardson’s scaling arguments and is also constructed to be applicable to real turbulence. The fundamentally new concept here is that turbulent pair diffusion is the convolution of local and non-local diffusional processes; and from this idea two limiting cases of non-Richardson pair diffusion regimes has been obtained.
It is important to note that Richardson’s scaling arguments were essentially kinematic in nature, being developed before the age of flow structures and dynamical energy transfer between different scales and intermittency was established. As such, the scaling laws are broadly independent of the precise mechanisms of dynamical energy transfer between different scales, much like the Kolmogorov spectrum itself. This is true also of the new local-non-local theory for the same reasons.
This may also explain why KS works well in the investigation of statistical scaling laws. KS is a Lagrangian model for tubulent diffusion, comparable to non-Markovian stochastic diffusion models. KS cannot account for all aspects of turbulence, like actual energy transfer between different scales, large scale sweeping of smaller scales, and true intermittency. However, as the theory itself is constructed for statistical scaling laws that are not strongly dependent on the precise dynamics, this may not be a critical limitation in the present context. Furthermore, it has been shown in [32] that the sweeping error in KS is negligible if simulations are carried out in the swept frame of reference as has been done here; and at least the energy spectrum of intermittent turbulence can be readily implemented in KS. The results presented here show that in the limit of locality, i.e. small inertial subrange, KS yields the correct Richardson scaling, and in the limit of asymptotically infinite inertial subrange with intermittent-like spectrum we obtain a non-local scaling law which is remarkably close to the revised 1926 data. These results provide some degree of confidence in the fidelity of the KS results; but we will have to wait for DNS or experiments that can interrogate very large inertial subranges for a definitive answer to how realistic KS statistics is of actual turbulence.
The main mathematical hypothesis of this theory is the separation of the inertial subrange in to two wavenumber ranges, which leads to the existence of two broadly independently diffusional processes that produce local and non-local scalings individually, but together they produce scale dependent diffusion coefficients with constant power laws γ(p)—the actual value of the γ(p) is dependent upon the inertial range power law ∼ k−p and also upon the size of the inertial subrange, Rk Eq (1).
The main predictions of the theory are, firstly that in turbulence with generalized energy spectra, E(k) ∼ k−p, and for asymptotically infinite inertial range, Rk → ∞, the pair diffusion coefficient scales like, , with (1 + p)/2 < γ(p) ≤ 2, in the range 1 < p ≤ 3, which is intermediate between the purely local and purely non-local scaling power laws. Secondly, for small inertial subranges, there exists quasi-local scaling regimes, with γ(p) ≈ (1 + p)/2, because the non-local range of scales are effectively absent and the entire inertial subrange then acts locally at all separations inside the inertial subrange.
The resuts from the KS simulations reported here confirm all the predictions of the theory. Most importantly, two sets of non-Richardson pair diffusion regimes are observed: (1) non-local regimes for asymptotically infinite inertial subrange, Rk = 106, Figs (1) to (5); and (2) quasi-local regimes for short inertial subrange of Rk = 102, Figs (6) to (10). A smooth transition from small to large subranges is observed in Figs (11) to (14).
The simulation results for Rk = 106 are in good agreement with the revised geophysical data in Fig 1 in [12], . In the critical range of intermittent spectra, the scalings γ are generally within about 1% of 1.564. For E(k) ∼ k−1.72, the simulations produced,
. The equivalent power law scaling in 〈l2〉 is ∼ t4.505, Table 1.
In short inertial subranges of size Rk = 102 in the critical range of intermittent spectra, the scalings γ(p) are generally within about 5% of locality scaling laws. For intermittency observed in real turbulence, E(k) ∼ k−1.72, the simulations produce the scaling law, , and the equivalent power law scaling in 〈l2〉 is ∼ t3.509, Table 2.
An important corollory of the current work is that real turbulence with intermittency does not contain the classical R-O t3-regime, even under the locality hypothesis. This is remarkable because the R-O t3-regime has been much debated for decades and its existance taken for granted.
To date no scaling greater than t3 inside the inetial subrage has been reported in experiments or in DNS, except for the 1926 dataset—this implies that current laboratory experiments and DNS where pair diffusion has been investigated are still in or below the quasi-locality limit. Indeed, the biggest inertial subranges generated in current DNS is around Rk = 102, which is consistent with our theory and simulations which suggest that this is indeed the lower limit for quasi-local regimes to be observable.
Some results apparently showing pair diffusion scaling greater than ∼ t3 have been reported in some recent DNS, [39, 40]. However, these results are for small initial separations l0 ≤ η and appear over a short time period after release. The authors themselves note that this is probably due to the influence of the separation in the dissipation range at earlier times—that is, in a DNS which contains a genuine dissipation range, particles that have left the dissipation range continue to be affected by the dissipation range some distance in to the inertial subrange. This is likely a manifestation of the ultra-violet corrections mentioned in the main text above and also in [12]. This should not be confused with genuine inertial range scaling.
Finally, we remark that the concepts investigated here and in [12] could have a significant impact on the general theory of turbulence, evoking some important questions for the future. For example, are long-range correlations in turbulence more significant than previously thought in high Reynolds number turbulence? Could turbulent diffusion be better modelled as a bi-variant process? And, if experiments ultimately do show that the KS results are close to real turbulence statistics, does this mean that the dynamics plays a less important role for lower order Lagrangian statistics?
Acknowledgments
The author would like to thank SABIC for funding this work through project number SB101011, and the ITC Department at KFUPM for making available the High Performance Computing facility for this project. The author would also like to thank Mr. K. A. K. K. Daoud for producing the parallel version of the KS code.
References
- 1. Shraiman BI, Siggia ED. Scalar turbulence. nature 405:639–646 (2000). pmid:10864314
- 2. Vaillancourt PA, Yau MK. Review of particle-turbulence interactions and consequences for cloud physics. Bull. Am. Meteorol. Soc. 81:285–298 (2000).
- 3. Boradwell JF, Breidenthal RE. A simple model of mixing and chemical reaction in a turbulent shear layer. J. Fluid Mech., 125:397–410 (1988).
- 4.
Peters N, Williams FA. 22nd Proc. Symp. on Combustion. (Seattle), 163-199 (1988)
- 5. Pope SB. Lagrangian PDF methods for turbulent flows. Annu. Rev. Fluid Mech. 26:23–63 (1994).
- 6. Huber M, McWilliams JC, Ghil A. A Climatology of Turbulent Dispersion in the Troposphere. J. Atmos. Sci. 58:22377 (2001).
- 7. Berloff PS, McWilliams JC, Bracco A. Material Transport in Oceanic Gyres. Part I: Phenomenology J. Phys. Oceanogr. 32:764–796 (2002).
- 8. Wolf MC, Voigt R, Moore PA. Spatial arrangement of odor sources modifies the temporal aspects of crayfish search strategies. J. Chem. Ecology 30(3):501–517 (2004).
- 9. Joergensen JB, Mann J, Ott S, Pecseli HL, Trulsen J. Experimental studies of occupation and transit times in turbulent flows. Phys. Fluids 17:035111 (2005).
- 10. Richardson LF. Atmospheric diffusion shown on a distance-neighbour graph. Proc. Roy. Soc. Lond. A 100:709–737 (1926).
- 11. Obukhov A. Spectral energy distribution in a turbulent flow. Izv. Akad. Xauk. SSSR. Ser. Geogr. i Geojz 5:453–466. (Translation: Ministry of Supply. p. 21 1097) (1941).
- 12. Malik NA. Turbulent Particle Pair Diffusion: A theory based on local and non-local diffusional processes. PLoS One 13(10):e0202940, (2018). https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0202940. pmid:30281611
- 13. Wilkins EM. Observations on the separation of pairs of neutral balloons and applications to atmospheric diffusion theory. J. Meteor. 17:324–327 (1958).
- 14. Fung JCH, Hunt JCR, Malik NA, Perkins RJ. Kinematic simulation of homogeneous turbulence by unsteady random Fourier modes. J. Fluid Mech. 236:281–318 (1992).
- 15. Virant M, Dracos Th. 3D PTV and its application on Lagrangian motion. Meas. Sci. Technol. 8:1539 (1997).
- 16. Batchelor GK. Diffusion in field of Homogeneous Turbulence II. The relative motion of particles. Math. Proc. Camb. Phil. Soc. 48(2):345–362 (1952).
- 17. Salazar JPL, Collins LR. Two-particle dispersion in isotropic turbulent flows. Annu. Rev. Fluid Mech. 41:405–432 (2009).
- 18. Morel P, Larchaveque M. Relative Dispersion of Constant-Level Balloons in the 200-mb General Circulation J. Atm. Sci. 31:2189–2196 (1974).
- 19.
Fung JCH, Vassilicos JC. Kinematic Simulation of Two-Particle Dispersion. In: Gavrilakis S., Machiels L., Monkewitz P.A. (eds) Advances in Turbulence VI. Fluid Mechanics and its Applications, Springer, Dordrecht, 36:617–618 (1996).
- 20.
Malik NA. Structural diffusion in 2D and 3D random flows. In: Gavrilakis S., Machiels L., Monkewitz P.A. (eds) Advances in Turbulence VI. Fluid Mechanics and its Applications, Springer, Dordrecht, 36:619–620 (1996).
- 21. Kraichnan RH. Diffusion by a random velocity field. Phys. Fluids 13:22–31 (1970).
- 22.
Turfus C, Hunt JCR. A stochastic analysis of the displacements of fluid element in inhomogeneous turbulence using Kraichnan’s method of random modes. In: Comte-Bellot G., Mathieu J. (eds) Advances in Turbulence, Springer, Berlin, Heidelberg, 191–203 (1987).
- 23. Nicolleau FCGA, Nowakowski AF. Presence of a Richardson’s regime in kinematic simulations. Phys. Rev. E 83:056317 (2011).
- 24. Fung JCH, Vassilicos JC. Two-particle dispersion in turbulent-like flows. Phys. Rev. E 57:1677 (1998).
- 25. Malik NA, Vassilicos JC. A Lagrangian model for turbulent dispersion with turbulent-like flow structure: comparison with direct numerical simulation for two-particle statistics. Phys. Fluids 11:1572–1580 (1999).
- 26. Murray S, Lightstone MF, Tullis S. Single-particle Lagrangian and structure statistics in kinematically simulated particle-laden turbulent flows. Phys. Fuids 28:033302 (2016).
- 27. Voßkuhlea M, Pumir A, Leveque E, Wilkinson M. Collision rate for suspensions at large Stokes numbers—comparing Navier–Stokes and synthetic turbulence. J. Turbulence 16(1):15–25 (2014).
- 28. Maxey MR. The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields. J. Fluid Mech. 174:441–465 (1987).
- 29. Meneguz E, Reeks MW. Statistical properties of particle segregation in homogeneous isotropic turbulence. J. Fluid Mech. 686:338–351 (2011).
- 30. Farhan M, Nicolleau FCGA, Nowakowski AF. Effect of gravity on clustering patterns and inertial particle attractors in kinematic simulations. Phys. Rev. E 91:043021 (2001).
- 31. Perkins RJ, Malik NA, Fung JCH. Cloud Dispersion Models. Appl. Sci. Res. 51:539–545 (1993). https://doi.org/10.1007/BF01082588.
- 32. Malik NA. Residual sweeping errors in turbulent particle pair diffusion in a Lagrangian diffusion model PLoS ONE 12(12):e0189917 (2017). https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0189917 pmid:29253875
- 33. Thomson DJ, Devenish BJ. Particle pair separation in kinematic simulations. J. Fluid Mech. 526:277–302 (2005).
- 34. Anselmet F, Antonia RA, Danaila L. Turbulent flows and intermittency in laboratory experiments. Planetary and Space Science, 49:1177–1191 (2001).
- 35. Tsuji Y. 2004 Intermittency effect on energy spectrum in high-Reynolds number turbulence. Phys. Fluids, 16, L43 (2004).
- 36. Meyers J, Meneveau CA. Functional form for the energy spectrum parametrizing bottleneck and intermittency effects. Phys. Fluids, 20:065109 (2008).
- 37. Tsuji Y. High-Reynolds-number experiments: the challenge of understanding universality in turbulence. Fluid Dyn. Res., 41(6):064003 (2009).
- 38. Eyink GL, Benveniste D. Suppression of particle dispersion by sweeping effects in synthetic turbulence. Phys. Rev. E 87:023011 (2013).
- 39. Buaria D, Sawford BL, Yeung PK. Characteristics of backward and forward two-particle relative dispersion in turbulence at different Reynolds numbers. Phys. Fluids, 27:105101 (2015).
- 40. Bragg AD, Ireland PJ, Collins LR. Forward and backward in time dispersion of fluid and inertial particles in isotropic turbulence. Phys. Fluids, 28:013305 (2016).