2.1 The numerical timestep error

The idea of sweeping scales of motion implies a clear separation of scales between the sweeping and the swept scales. In reality, although there is a broad division between large and small scales, there is no well-defined cut-off between them because the spectrum of turbulence is continuous.

However, the situation in KS is easier because all scales are already separated by construction as a Fourier series, Eq (27), so it is a trivial matter to ‘remove’ scales of motion simply by eliminating those terms in the Fourier series. In the present case, the large sweeping scales are approximated as those scales where k < k₁. Thus, in the swept frame of reference the sweeping action itself can be approximatedd in KS by setting E(k) = 0 for k < k₁. However, there still remains a residual sweeping caused by the largest of the inertial scales sweeping the scales local to the pair separation.

In the ensuing analysis we will make use of Taylor expansions applied to KS flow fields. It is important to note that this is possible because although KS is a Lagrangian method, it is not in the genre of stochastic models such as Random Walk models. In KS, we generate an ensemble of flow fields, {uⁿ}, n = 1, 2, 3, … in which each flow field uⁿ is a Fourier series with given coefficients and is therefore smooth and differentiable. The randomness in KS comes from choosing each Fourier mode and each coefficient in the Fourier expansion in each flow field uⁿ randomly from specified probability distributions, the square of the magnitudes of the coefficients being proportional to a given energy spectrum, see Section 3.2.

Consider an ensemble of particle pairs released in a field of homogeneous turbulence at time t = 0 with some small initial separation l₀. 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₀ [8].

Consider the particles in one of these pairs, labeled 1 and 2, as shown in Fig 1. The particle locations are x₁(t) and x₂(t) respectively at time t; and the pair displacement is l(t) = x₂ − x₁, and l(t) = |x₂ − x₁|. u(x, t) is the flow due to inertial range of turbulent motions that are effective in increasing the pair separation. We will refer to u(x, t) as the ‘physical’ flow. 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 ‘switched on’—this is not to be confused with the total KS velocity field which is , see Fig 1. In a real flow field, the small scale turbulence would be carried along dynamically by the larger scales, and therefore we could assume u^s(x, t) = 0. In KS, the larger scales force a particle to cut through the smaller scales, hence u^s(x, t) ≠ 0; but the spatial gradients of u^s are much smaller than the spatial gradients of u.

The flow u(x, t) transports the particles to x₁(t*) and x₂(t*) respectively at the next time step t* = t + dt; while the KS flow (u + u^s)(x, t) transports the particles to and respectively. Note that , and .

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:

1. u = u(x, t) is the ‘physical’ fluid velocity field
2. u^s = u^s(x, t) is the additional (residual) sweeping velocity field
3. v(l) = u(x₂) − u(x₁) is the ‘physical’ relative velocity
4. v^s(l) = u^s(x₂) − u^s(x₁) is the additional (residual) relative velocity
5. is the total KS velocity
6. is the total KS relative velocity

We will simplifying the notation as much as possible, e.g. u₂ = u(x₂, t), and .

In the following analysis, we consider the situation where particle pairs have been released at some earlier time, and now at time t later, the pair separation l = |l| is still small compared to some characteristic large length scale L and inside the inertial subrange, i.e. η ≪ l ≪ L; where η is the Kolmogorov scale. we will assume that inertial subrange is large enough, L/η ≫ 1, such that inertial range scaling will apply.

Furthermore, the following Taylor expansions in the KS velocity fields will be assumed, (3) (4) (5) to leading order in dt, where the increment, dt, is small compared to all other timescales in the system, governed only by the need for numerical accuracy.

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, (6)

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 . Using Eqs (5) and (6) and ignoring terms of order dt² and higher, (7)

The incremental timestep error between the KS and physical diffusivities for a given timestep dt is, E_K = |K^s* − K*|. Using the expansion in Eq (7), yields (8)

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²) which is neglected).

All the terms in Eq (8) are now evaluated at the same time t*, so without loss of generality t* is replaced by t and the superscript ‘*’ is dropped. The subscript ‘₁’ is also dropped because of homogeneity. Eq (8) now simplifies to, (9)

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 the local separation scales l that it induces is small, i.e. gradients like |(l · ∇)u^s| ≈ 0. However, gradients of the relative velocity v(l) itself can be large. The magnitude of 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 Eq (9) are zero or nearly zero compared to the second term. This indicates that the second term in Eq (9) makes the dominant contribution to the error.

Consider generalized energy spectra of the form E(k) = ε^2/3L^5/3−pk^−p, for k₁ ≤ k ≤ k_η and for 1 < p ≤ 3, and with k_η/k₁ ≫ 1 [9, 15]. In the swept frame of reference, E(k) = 0 for k < k₁. The rate of energy dissipation is ε ∼ U³/L, where U is the velocity scale in the energy containing scales.

The previous discussion implies that and therefore, (10)

The energy in turbulent inertial scales local to l is, v²(l) ∼ E(1/l)/l, and therefore, (11) (12)

It is usual to assume the scaling l ∼ σ_l as previously mentioned. Then, the second term in Eq (9) is given by, (13)

All the other terms in Eq (9), labeled respectively E₁, E₃, …, E₇, scale proportional to u^s or are much smaller, and this leads to the following estimates, (14)

All of these ratios are small for σ_l/L < 1 and 1 < p ≤ 3. This is true even in the expression for E₁/E₂ because the factor L/Udt is subdued by the very small factor in the brackets. It is reasonable to conclude that the 2nd term in Eq (9) is dominant and therefore E_K ≈ E₂.

To estimate E₂ itself, an estimate for u^s(l) is needed. The inertial subrange contains only a small part of the total turbulence energy across the entire wavenumber range, 0 < k < ∞; a typical turbulence spectrum is the von Karman spectrum E_vk(k). The von Karman spectrum is essentially a fit between the low wavenumber energy spectrum (E(k) ∼ k⁴), and the high wavenumber inertial subrange spectrum (E(k) ∼ k^−5/3); the exact form is not important for our purposes; but see [16].

Let, the fraction of the turbulence energy in the inertial part of the spectrum be, I² = E_IN/E_TOT, where E_IN is the energy in the inertial range, and E_TOT is the total turbulence energy. Then if, and E_TOT ∼ U², then we obtain, , and u_ks ≈ IU. The choice of wavenumber which separates the inertial range from the large scales is somewhat arbitrary.

The wavenumbers that contribute to the sweeping of particle pairs at separation σ_l are in the range k₁ ≤ 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, (15) and using ε ∼ (u_ks)³/L ∼ (IU)³/L, this gives (16)

Using this in Eq (13), 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 (17) C_k is the constant of proportionality, which can depend upon p.

For p = 5/3, the residual error per unit timestep is, (18)

As p → 1, E_K → E_K(1) ≈ C_kU²/L ≈ Constant for σ_l/L < 1, 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, (19)

In this limit, nearly all the energy is contained in the largest scales and inertial subrange scaling is no longer applicable.

2.2 The integrated error

The incremental timestep error in the diffusivity E_K only tells us the error in KS due to the lack of dynamic sweeping per numerical time step. As the particle pair separations increase (on average), these errors will not only increase, but they will also accumulate over a period of time. The important question is, does this integrated error remain small?

We therefore need to calculate the integrated error, relative to the actual diffusivity, over the period while the pair separation remains within the inertial subrange. The limit that the pair separation approaches the upper end of the inertial subrange, σ_l ∼ L, implies that σ_l/η ≫ 1. In the limit of infinite inertial subrange this corresponds to σ_l/η → ∞.

Eq (17) provides a way of estimating an upper bound for the integrated residual error. Fig 2 shows the log-log plots of the factor for selected powers in the range 0 < p ≤ 3. The range over which is very short and close to x = 1; but over the rest of the range . Hence, to a good approximation, f(x) < x^(p−1)/2 for all x ≤ 1, and using this in Eq (17) with x = σ_l/L yields, (20)

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. The function , 0 < x ≤ 1, for selected powers 1 < p ≤ 3.

https://doi.org/10.1371/journal.pone.0189917.g002

The integrated residual error, , over a period of time is, (21)

Assuming the pair separation scaling , where t is the time and for some χ_p > 0, yields, (22)

If the pair diffusivity scales like, , for some γ_p > 0 then χ_p = 1/(1 − γ_p/2) is an exact relation.

The most important quantity is the relative integrated residual error with respect to the pair diffusivity, . Using the above expression for χ_p, and replacing the scaling with L by scaling with η, leads to (23)

For strict locality scaling γ_p = (1 + p)/2, and this becomes (24)

Since γ_p − 1 > 0, decreases with increasing pair separation for all p > 1. For p = 5/3, we have γ_p = 4/3 and we obtain, (25)

Thus, as σ_l/η → ∞ then decreases; and as σ_l/η → 0 then increases.

Even for non-local scaling, assuming that γ_p does not deviate too far from the local scaling, the above order of magnitude for is still approximately true. In fact, in KS we know, Section 3.3, that γ_p is slightly greater than locality so the errors will be slightly smaller than in Eq (25).

However, between these two asymptotic limits there must exist an intermediate range of separations, between 1 < σ_l/η < ∞, where the errors remain negligible. In reality, all turbulence spectra are finite in range. So what size of the inertial subrange do we need to approximate an ‘infinite’ subrange?

Even then, the crucial question is: is the intermediate range of scales where the KS sweeping errors are negligible wide enough for inertial range scaling to be actually observable in KS?

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.

The errors in Eqs (24) and (25) cannot be determined directly from simulations because only the KS particle trajectories are generated directly.

However, the effects of the error analysis can be observed in several important ways that will be described in the next section where we carry out KS calculations for the pair diffusivity.

3.1 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.

1. 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₁.
2. 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₀δ(k − k₀) at some low wavenumber k₀ < k₁, and such that E₀ is the energy in the von Karman spectrum in the range 0 < k < k₁.

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.

The analysis for Case 2 is similar to that which leads to Eq (24), except that U, L, and I are different. From I² ∼ E_IN/U², we obtain which remains relatively constant in the expressions (24) and (25). Thus, the differences in the errors between the two cases depends largely upon the different length scale L in the two cases, and assuming A_k = C_k UI ≈ constant we obtain, (26)

In Case 2, when large random scales are included in the simulations then L is about 10 times bigger than in Case 1. We may thereore expect the errors to be about 10 times bigger in Case2 than in Case 1.

The exact level of the error in Eq (26) depends upon the choice of cut-off scale between the sweeping and non-sweeping scales which is somewhat arbitray. It is best to view Eqs (24)–(26) as a general scaling law with parameters U, L, I, and the actual level of error can only be determined by simulations.

3.2 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). Although KS cannot capture the full dynamics of particle motion, it has the advantage that it can generate extremely large energy spectra that is far beyond current experimental or DNS capabilities, and therefore it can be used to examine some aspects of very high Reynolds number scaling laws. 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 [17], Turfus [18], Fung & Vassilicos [19], Malik & Vassilicos [20], Farhan et al. [21]. Meneguz & Reeks [22] carried out a DNS of inertial particle motion, and compared it to results from KS which they found to agree well with the DNS. Murray et al. [23] investigated inertial particle statistics using KS.

KS generates turbulent-like non-Markovian particle trajectories by releasing particles in flow fields that are incompressible by construction and which satisfy Eulerian statistics up to second order. A turbulent flow field realization is produced as a truncated Fourier series, (27) where N is a suitable number of representative wavemodes, typically hundreds for very long spectral ranges, k_η/k₁ ≫ 1. is a random unit vector ( and k_n = |k_n|). The coefficients A_n and B_n are chosen such that their orientations are randomly distributed and uncorrelated with any other Fourier coefficient or wavenumber, and their amplitudes are determined by , where E(k), k₁ ≤ k ≤ k_η, is the turbulent energy specturm. The angled brackets 〈⋅〉 denotes the ensemble average over many flow fields. This construction ensures incompressibility in each flow realization, ∇ · u = 0. The flow field ensemble generated in this manner is statistically homogeneous, isotropic, and stationary.

An important feature of KS is that unlike some other Lagrangian methods, by generating entire kinematic flow fields in which particles are tracked it 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 completely eliminated.

The energy spectrum E(k) can be chosen freely within a finite range of scales. In turbulent particle pair studies the interest is in Kolmogorov-like power law spectra, (28) C_E is a constant. The largest represented scale of turbulence is 2π/k₁, and the smallest is the Kolmogorov scale η = 2π/k_η. The constant is normalized such that the total energy contained in the range k₁ ≤ k ≤ k_η is 3(u′)²/2, where u′ is the rms turbulent velocity fluctuations in each direction. ε(p) is determined by integrating the spectrum, . (p = 1 is a singular limit which is not consider here.) v_η = (εη)^1/3 is the velocity micrcoscale, and τ_η = ε^−1/3η^2/3 is the Kolmogorov time micrcoscale.

The frequencies are chosen according to usual practice to be proportional to the eddy-turnover frequencies, i.e. . The choice of λ is somewhat arbitrary, but provided λ < 1 it does not affect the diffusion scaling itself—even frozen field with λ = 0 yields the same scaling [15]. λ = 0.5 is a common practice in KS which is also chosen here.

The distrbution of the wavemodes is geometric, k_n = k₁rⁿ⁻¹, with r = (k_η/k₁)^1/(N−1). The grid size in wavemode-space of the n^th wavemode is .

A particle trajectory is obtained by integrating the Lagrangian velocity W_L(t), (29)

The method of computing trajectories in KS is well established in the literature, see previous references above. Briefly, since an individual KS flow field realization is smooth and differentiable, Eq (29) can be readily integrated using any standard method such as Runga-Kutta or Predictor-corrector methods. The time step Δt should be smaller than any other time scale in the system—in this case the Kolmogorov time scale τ_k = 2π/k_η; thus we require Δt ≪ τ_k. In our simulations, particle trajectories are produced by integrating Eq (29) with a fixed time step of Δt ≈ 0.01τ_k in a fourth order Adams-Bashforth predictor-corrector method.

It is important to produce a very large ensemble of independent particle pair trajectories. Only eight independent particle pairs are released in any one given KS flow realization, each pair set part by more than an integral length scale. Each pair is initially released with a pair separation distance of l₀/η = 0.5. To obtain a true ensemble, this process is repeated in many thousands of KS flow realizations. (Releasing thousands of particle pairs in the same KS flow field will not produce the required independence which may produce biased statistics.) Pairs of trajectories are thus harvested over a large ensemble of flow realizations and pair statistics are then obtained from it for analysis.

The turbulent diffusivity 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²〉/dt, i.e. the derivative of the 〈l²〉, 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.

3.3 Results

KS simulations were performed with L = 1, k₁ = 1/L = 1, and k_η = 10⁶, C_E = 1.5 (Kolmogorov constant) and u′ = 1. There were 200 wavemodes per realization.

In Case 1 (swept frame of reference), E(k) = 0 for k < 1.

In Case 2, large scale random sweeping were added at the low wavenumber k₀ = 1/10, with E(k) = E₀δ(k − k₀). The energy in these sweeping scales, E₀, was equal to the energy contained in the von Karman turbulence spectrum for k < 1. k₀ 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.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. The turbulent diffusivity as log(K/(ηv_η)) against log(σ_l/η) obtained from KS.

From top to bottom, p = 1.01, 5/3, 3. Case 1 (red lines), Case 2 (green lines). No scale is show on the vertical axis because the energies in the two cases are different and Case 2 plots have been shifted vertically in order to compare the two cases directly. This does not affect the scalings (the slopes).

https://doi.org/10.1371/journal.pone.0189917.g003

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 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⁵, indicating that 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² is good.

For p = 3, the two cases overlap with a power-law scaling, γ₃ = 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.