Particle pair diffusion and persistent streamline topology in two-dimensional turbulence

Susumu Goto¹ and J C Vassilicos

Turbulence and Mixing Group, Department of Aeronautics, Imperial College, London SW7 2AZ, UK

¹Permanent address: Theory and Computer Simulation Center, National Institute for Fusion Science, 322-6 Oroshi-cho, Toki, 509-5292, Japan.

Email: susumu-goto@mail-box.jp and j.c.vassilicos@imperial.ac.uk

Received 23 March 2004
Published 25 June 2004

Contents

• 1. Introduction
• 2. Streamline topology: DNS observations
  • 2.1. Numerical scheme and parameters
  • 2.2. Flow field statistics
  • 2.3. Streamline topology
    • 2.3.1. In the local frame: zero-acceleration points and coherent vortices
    • 2.3.2. In the global frame: stagnation points and statistical persistence hypothesis
• 3. Physical picture of particle pair diffusion
  • 3.1. In the local frame
  • 3.2. In the global frame
  • 3.3. Geometrical (fractal-like) separation
• 4. Consistency of local and global viewpoints of particle pair diffusion
  • 4.1. Doubling time estimation from the local frame's viewpoint
  • 4.2. Doubling time estimation from the global frame's viewpoint
  • 4.3. Consistency
• 5. A simple model of particle pair diffusion
  • 5.1. Mean square separation
  • 5.2. PDF of separation
• 6. Numerical verification
  • 6.1. Mean square separation
  • 6.2. PDF of separation
  • 6.3. Doubling time
  • 6.4. Correlation functions
  • 6.5. Doubling time dependence on the stagnation point distribution
• 7. Particle pair diffusion in general cases
• 8. Concluding remarks
• Acknowledgments
• References

1. Introduction

Richardson's [1] study of particle^Note2  pair diffusion in homogeneous isotropic incompressible turbulence was motivated by measurements of turbulent diffusivities which showed them to vary over 11 or so orders of magnitude when the scale defining these diffusivities is varied from the smallest to the largest atmospheric eddy sizes. By fitting data available at the time, Richardson found that the effective diffusivity K is proportional to the power of the average size of the puff. The Reynolds-averaged advection diffusion equation for the average concentration of puffs would therefore depend on the size of the puff, an observation fundamentally at odds with the concept of Fickian diffusivity. Richardson recognized that a new approach was needed and introduced what, in modern terminology, is known as the probability density function (PDF) of particle pair separations Δ. He showed that the Reynolds-averaged advection diffusion equation with constant eddy diffusivity K₀ (independent of scale) implies that the equation for this PDF P(Δ, t) is

where d is the spatial dimension. He then recognized that, whereas the Fickian equation for the average concentration cannot be generalized, the PDF equation (1) can be generalized, and he postulated it to be

where F(Δ) is a function of Δ replacing 2K₀. This equation for P(Δ, t) is important because concentration fluctuations in turbulent flows can be calculated from integrations of particle pair PDFs [2, 3]. Richardson [1] chose F(Δ)~Δ^4/3 so as to recover K(Δ)≡(d/dt)Δ²~Δ^22/3 (averages being defined over many particle pairs) which seemed to, more or less, agree with the (mixed quality) observations at his disposal. This, perhaps fortuitous, choice of F(Δ) was later corroborated by Kolmogorov-type dimensional analysis which, based on the locality-in-scale hypothesis, leads to F(Δ)~ε^1/3Δ^4/3 and (d/dt)Δ²~ε^1/3Δ^22/3 where ε is the kinetic energy dissipation rate per unit mass [4, 5]. The locality-in-scale hypothesis states that the dominant contribution to the pair separation rate comes from eddies of size comparable to the pair separation. This hypothesis implies that the evolution of turbulent pair diffusion straddles the inertial range and is closely linked to the Kolmogorov-type power-law form of the energy spectrum.

New laboratory and computational techniques developed over the past 25 years (e.g. particle image velocimetry and direct numerical simulations (DNS)) are now reaching the point where they can be applied to the study of particle pair turbulent diffusion. Whereas these approaches still remain broadly inconclusive for three-dimensional turbulent pair diffusion, they do seem to offer tangible possibilities for the case of two-dimensional turbulence in the inverse energy cascade range. In this case, a turbulence with a broad inertial range can be simulated and can also be comprehensively visualized. Recent DNS by Boffetta and Sokolov [6] and laboratory measurements by Jullien et al  [7] of such flows have given rise to indications that they are characterized by a k ^- 5/3 form of the energy spectrum, by the Richardson Δ²~εt³ law (which can be obtained from (d/dt)Δ²~ε^1/3Δ^22/3), and perhaps also by a PDF P(Δ, t) that has the Richardson form arising from (2). Visualizations of particle pair trajectories by Jullien et al  [7] have revealed that pairs travel together for long times and then separate into sudden bursts. (This behaviour has also been observed in laboratory three-dimensional turbulence by Ott and Mann [8] where, however, the inertial range is very limited, if at all present.)

The question that is raised is this: is the particle pair separation mechanism by comparatively rare sudden bursts compatible with Richardson's PDF equation and with the locality-in-scale hypothesis which underpins Richardson's t³-law?

In an attempt to understand this separation mechanism by rare sudden bursts, one might draw attention to the streamline structure of the flow. In the extreme case where a velocity field is frozen in time, particle trajectories coincide with streamlines and pairs can be expected to travel together till they encounter a region of highly curved diverging streamlines around a straining (hyperbolic in two dimensions) stagnation point where they separate in a sudden burst.

In a turbulent flow, which is inherently unsteady, if the streamline structure is coherent enough in time it will then have a decisive imprint on trajectories. This point of view was first proposed qualitatively by Fung et al  [9] and Fung and Vassilicos [10]. Fung and Vassilicos in fact modified the locality-in-scale hypothesis as follows: the dominant contribution to the pair separation rate comes from straining stagnation points of size comparable to that of pair separation. They suggested a fractal cat's eyes within cat's eyes streamline structure of two-dimensional turbulence to illustrate how different hyperbolic stagnation points can be linked to streamlines of different sizes. This multiple-scale structure of streamlines was quantified by Dávila and Vassilicos [11].

However, streamlines and their persistence in time are not Galilean invariant. In this paper we introduce the idea that the frame where particle pair diffusion's dependence on streamline structure is studied must be the frame which, in some deterministic or statistical sense, maximizes the streamline persistence in time. As a consequence, on the basis of DNS of two-dimensional turbulence, we propose two frames, one local and one global, where we reformulate the locality-in-scale hypothesis in terms of the sudden burst mechanism. This allows us to derive the t³-law and a new PDF equation for particle pair separations from this mechanism. This PDF equation turns out to be a generalization of Richardson's PDF equation (2). We use our DNS to validate our approach.

The rest of this paper is organized as follows. In the following section, we summarize the numerical scheme and parameters of our DNS of inverse energy cascade two-dimensional turbulence (section 2.1), and we report the statistics of the velocity field where particles are tracked (section 2.2). Especially, we focus on the streamline structure and stagnation points of the field (section 2.3). In section 3 we propose two physical pictures of particle pair diffusion, one in a local frame (section 3.1), the other in a global frame (section 3.2), and we introduce our formulation of the locality-in-scale hypothesis (section 3.3). The consistency of the two pictures is established in section 4. In section 5 we derive the consequences of this hypothesis, specifically the t³-law of the mean square separation and an equation for the separation PDF, which we compare with Richardson's PDF equation. The predictions of our model are numerically verified in section 6. Finally, we show in section 7 how our model can be generalized to d-dimensional turbulence with an energy spectrum proportional to k ^- p, and we conclude in section 8.

2. Streamline topology: DNS observations

One of the most important features of two-dimensional turbulence is the simultaneous cascading of energy (from smaller to larger scales) and enstrophy (from larger to smaller scales). Kraichnan [12] and Batchelor [13] predicted a k ^- 5/3 form for the energy spectrum in the inverse energy cascade range (hereafter we call this range the inertial range). Various authors have investigated the validity of this prediction in DNS of two-dimensional turbulence, e.g. [14]-[16], with conflicting conclusions. For example, Borue [14] did not find support for the k ^- 5/3 spectrum, whereas Sukoriansky et al  [15] did. Even though the energy flux is constant in the inertial ranges of both these DNS studies, the crucial difference is in their choices of large-scale drag. In this paper, we do not use the large-scale drag of Sukoriansky et al , but we follow Borue in modelling large-scale drag by a hyper-drag form. However, our hyper-drag is of much lower order than Borue's, which allows us to achieve good inertial range statistics. Specifically, our energy spectrum has a k ^- 5/3 scaling over a wide inertial range, and our energy flux is constant over the entirety of this range. The aim is to achieve a simulation where the energy dissipation rate can be calculated as accurately as possible, which is an important prerequisite for an accurate study of particle pair diffusion statistics.

2.1. Numerical scheme and parameters

We consider particle diffusion in an incompressible (∇u = 0) two-dimensional turbulent velocity field u(x, t), which is governed by the modified Navier-Stokes equation,

with periodic boundary conditions of period 2π in the two orthogonal directions. Here, p(x, t) is the pressure field, ρ( = 1) is a constant density, ν is a hyper-viscosity, α is a hyper-drag coefficient and f(x, t) is an external forcing vector field.

We solve numerically the vorticity equation, which is the curl of (3), by employing the Fourier spectral method (de-aliased by the -rule) for accurate spatial derivative estimations, and the 4th order Runge-Kutta method for time integration (the time increment is denoted by dt). We have used different numbers N² of numerical grid points ranging from 512² to 4096².

In the present DNS, we employ an effective forcing which acts only on the Fourier components in the wavenumber range between k_f and βk_f (β is a constant, slightly larger than 1, see table 1). We reset al l magnitudes of the Fourier components of vorticity in this wavenumber range at every time step to ensure that they are kept equal to the constant value (see table 1). However, the phases of these Fourier components are allowed to vary freely in time according to the vorticity equation. To obtain the k ^- 5/3 energy spectrum in the inverse energy cascade range, we set the forcing wavenumber k_f at a relatively large wavenumber, k_f = N/5. Note that the largest wavenumber of this system is N/3 because of the de-aliasing by the -rule.

Table 1. Numerical parameters: number of numerical grid points N2, forcing wavenumber kf, ratio of forcing wavenumber range β, magnitude of fixed Fourier component of vorticity , time increment of temporal integration dt, hyper-viscosity coefficient ν and its exponent m1, hyper-drag coefficient α and its exponent m2.

Run	N2	kf	β		dt	ν	m1	α	m2
A	5122	102	1.01	1	5×10 - 4	2.5×10 - 32	8	2.5	1
B	10242	205	1.005	0.79	2.5×10 - 4	7.5×10 - 37	8	2.5	1
C	20482	410	1.0025	0.63	1.25×10 - 4	2.0×10 - 41	8	2.5	1
D	40962	819	1.00125	0.5	6.25×10 - 5	6.0×10 - 46	8	2.5	1

Energy and enstrophy are fed into the flow at the forcing scale η_f = 2π/k_f, and they cascade towards larger and smaller scales, respectively. Hence, some dissipative terms at both the largest and smallest scales are needed to establish a statistically stationary state. For the small-scale enstrophy dissipation, we employ a high (specifically m₁ = 8)-order hyper-viscosity to ensure that the small-scale enstrophy dissipation does not contaminate inertial-range statistics.

For the large-scale energy dissipation, we employ first-order (m₂ = 1) hyper-drag. As shown in the next section, this hyper-drag produces the expected inertial range statistical signatures such as a constant energy flux and the k ^- 5/3 energy spectrum. Higher-order (m₂ > 1) hyper-drag leads to an accumulation of energy at the large scales which contaminates the inertial range and causes deviations from the k ^- 5/3 energy spectrum. This may be the reason why the energy spectrum observed in [14], where 8th-order hyper-drag was employed, is steeper than k ^- 5/3. The energy dissipation spectrum, being , becomes increasingly singular near k = 0 for increasingly large values of m₂. As a result, the numerical scheme fails to capture this dissipation accurately, and leads to a pileup of energy at large scales.

A second reason for choosing a small value of m₂ is to make sure that the integral length scale (defined in the usual way on the basis of the second-order longitudinal velocity correlation function) of the velocity field is significantly smaller than the system size L₀ = 2π. Otherwise, the periodic boundary condition may generate spurious effects such as large vortices being advected by their own periodic images.

Finally, we draw attention to the fact that when m₂ is equal or too close to zero, then the drag affects a larger portion of the would-be inertial range. In finite size simulations this can lead to an energy flux that is not constant over the entire range of relevant wavenumbers (see figure 2 of [16] for the resolution N² = O(1000²)). It is important that the energy flux should be constant over a substantial inertial range in order to accurately estimate the energy dissipation rate due to the large-scale drag.

Numerical parameters adopted in the present simulations taking into account the considerations just presented are listed in table 1.

We have run the smallest resolution simulations (run A) with several initial conditions, and invariably reached a statistically stationary state which appeared to be independent of the initial condition choice (however we have not done a full systematic study to confirm this statement). Run B was initialized by an interpolation of the stationary state of run A, run C by an interpolation of the stationary state of run B, and so forth for run D. In all our runs, the time required to reach a statistically stationary state is much longer (by a factor of order 100) than the eddy turnover time of the largest eddies.

In the statistically stationary state of the velocity field u(x, t), we track a large number, O(N²), of particles by integrating the advection equation

where x_p(t) stands for particle position. This numerical integration is carried out using the fourth-order Runge-Kutta scheme. The right-hand side of (4) is estimated by the 4²-point Lagrangian interpolation of the simulated velocity field. The time increment of the integration of (4) is twice that (dt) of the velocity field, since the fourth-order Runge-Kutta method requires the estimation of u(x, t) at times between increments.

2.2. Flow field statistics

Figures 1 and 2 respectively show temporal averages in the statistically stationary state of the energy spectrum,

and the energy flux function,

for runs A, B, C and D. Here, is the Fourier component of the velocity field, ∫dΩ denotes an integral over a thin shell k≤|k|≤k + dk, and T(k, t) is the energy transfer function defined by

where is the energy input rate by the external forcing. In these figures, all the quantities are non-dimensionalized by the r.m.s. u ' of a velocity component and the energy dissipation rate ε which we estimate from the constant value of the energy flux function in figure 2 (see table 2). Indeed, figure 2 demonstrates the constancy of the energy flux function in our simulations, and the negative value of this constant implies that the energy cascades inversely from small to large scales. As a consequence, the energy spectrum takes the Kraichnan-Batchelor form

over the range where Π(k, t) is constant. The constant C_K is found to be approximately 7 irrespective of the resolution. This value compares well with previously obtained DNS values (e.g. C_K = 6.0±0.4 [16]) and the experimental value 5.5-7 obtained by Paret and Tabeling [17].

Table 2. Statistics of the simulated turbulence: energy dissipation rate ε, r.m.s. of a velocity component u ', integral scale , and its ratio to the forcing scale ηf = 2π/kf.

Run	ε	u '
A	3.0	3.2	0.38	6.3
B	1.4	2.5	0.32	10
C	0.60	1.9	0.29	19
D	0.20	1.3	0.23	30

As mentioned above, an important advantage in using hyper-drag in equation (3) is the constancy of the energy flux in the inertial range. This implies that we can accurately estimate the dissipation ε caused by the large-scale drag, which is a great advantage as the estimation of ε is in general highly non-trivial both in simulations and in the laboratory. The accuracy of our estimation is confirmed in figure 3 by the -law of the third-order longitudinal structure function,

in the inertial range (see derivation in [18]). Although this relation is clearly observed in our simulations, its scaling range is somewhat narrower than the inertial range observed in figures 1 and 2.

Figure 3. (a) Third-order longitudinal structure function for runs A, B, C and D. The dotted line indicates the -law. (b) Same but compensated.

2.3. Streamline topology

One might expect turbulent particle pair diffusion to occur as a result of encounters with high straining regions where large velocity gradients separate the pairs. However, in the case of stationary velocity fields, pairs can also separate when they encounter regions where strain is low but streamlines are highly curved and divergent. Of course, turbulence is not stationary but it is not delta-correlated in time either, and therefore effects of persistent streamline topology on particle trajectories can be important. In general, if the velocity field is stationary within some timescale then the topology of particle trajectories are solely determined by the topology of the streamlines. However, the persistence of streamlines is a frame-dependent concept. Consider a stationary velocity field: when subjected to a Galilean transformation it becomes non-stationary, and its streamlines lose their persistence. Inversely, given a time-dependent velocity field, it should be possible to choose a frame which maximizes the persistence of streamlines. In the following subsections, we introduce two such frames, one local (section 2.3.1) in the sense that it is concerned with specific particles, the other global (section 2.3.2) in the sense that it is concerned with statistical ensembles of particles.

Along with streamline persistence, the other important property of streamlines which is fundamental to our study of turbulent diffusion is the multiple-scale nature of streamline topology. A convenient tool for the study of scale-dependent flow structure is coarse-graining. In this paper, we coarse-grain by sharply low-pass filtering the energy spectrum at the cut-off wavenumber k_c. We plot in figure 4 coarse-grained vorticity fields for the cut-off wavenumbers , which are all in the inertial range. Because the enstrophy spectrum is proportional to k^1/3 in this range, and because k_c is smaller than k_f, the peak of the coarse-grained enstrophy spectrum is at k_c. This qualitatively explains the size of the vortices observed in figure 4.

Figure 4. Coarse-grained vorticity field and elliptic (yellow dots) and hyperbolic (black dots) zero-acceleration points. The cut-off wavenumber kc is in (a)-(d), respectively. Shade of red/blue indicates magnitude of positive/negative vorticity. The size of the area plotted is L0/4 = π/2. Run C.

In the following subsection, we define frames which maximize persistence of streamlines at different coarse-graining scales, η_c = 2π/k_c.

2.3.1. In the local frame: zero-acceleration points and coherent vortices Typical temporal evolutions of coarse-grained vorticity fields for and in the statistically stationary state are shown in figure 5 (see also movie 1a and movie 1b). The time lag between successive panels in figure 5 is equal to 0.075 which is about half the eddy turnover time of the largest vortices. The lifetimes of the observed vortices, in movie 1a, movie 1b and figure 5, appear to be longer than their eddy turnover times. Indeed the r.m.s. vorticity of these coarse-grained fields is estimated to be about 33 for and 20 for , which correspond to eddy turnover times 0.03 and 0.05, respectively.

Figure 5. Temporal evolution (0≤t≤0.375) of vorticity field coarse-grained after time integration. The cut-off wavenumber is (left row, movie 1a) and (right row, movie 1b). Shade of red/blue indicates magnitude of positive/negative vorticity. Yellow/black spheres are elliptic/hyperbolic zero-acceleration points. Run C. The size of the area plotted is L0/4 = π/2. The eddy turnover time is .

The observed long lifetime of these vortices suggests that in the frame moving with one of these vortices, streamlines around its centre are persistent. In fact, this local frame moving with the vortex is the frame where persistence of streamlines is maximized. To unambiguously determine this frame (how to unambiguously define the velocity of a vortex?) we need to introduce zero-acceleration points, x_a, defined by a(x_a, t) = 0, as we now explain.

In figures 4 and 5, the zero-acceleration points are plotted together with the magnitude of vorticity. We observe that some zero-acceleration points coincide with centres of vortices, and some others are located within zero-vorticity lines (white zones in figures 4 and 5). To classify these two types of zero-acceleration points, we condition them according to the sign of

which is reminiscent of the Weiss criterion [19]. Zero-acceleration points corresponding to negative/positive values of ∇a are coloured yellow/black in figures 4 and 5. As seen in these figures, many, though not all, zero-acceleration points with negative or positive acceleration divergence correspond to vortex centres or straining strips of zero vorticity, respectively. We call the former elliptic (vortical) zero-acceleration points, and the latter hyperbolic (straining) zero-acceleration points.

We first examine elliptic zero-acceleration points. It is observed (see movie 1a, movie 1b and figure 5) that elliptic zero-acceleration points at centres of vortices move with these vortices for times comparable with the vortex lifetimes.^Note3  This observation allows us to capture the motion of vortices by following the motion of elliptic zero-acceleration points at their centres, and thereby unambiguously determine the frame moving with the vortex in which streamline structure is approximately persistent. More precisely, the frame can be determined as the one which moves with the velocity u(x_a, t) at the elliptic zero-acceleration point, since in two-dimensional flow vortices can be expected to move approximately with the fluid velocity at their centre, i.e. u(x_a, t). It might be worth mentioning, in passing, that the determination of the frame with maximized streamline persistence on the basis of zero-acceleration points is consistent with the fact that at such points ∂u/∂t vanishes in the frame moving with the velocity u(x_a, t). In this frame, such points are stagnation points, and ∂u/∂t = 0 follows from a = ∂u/∂t + u∇u. Note that the acceleration field is Galilean invariant, and therefore the zero-acceleration points are also Galilean invariant.

We now examine hyperbolic zero-acceleration points. It is observed that some, though not all, hyperbolic zero-acceleration points residing in straining strips of zero vorticity have a long lifetime (see movie 1a, movie 1b and figure 5). This suggests that, in the frame moving with such a hyperbolic zero-acceleration point, the streamline structure closely surrounding it is persistent. Hyperbolic zero-acceleration points are located between vortices, therefore their relatively long lifetimes may correspond to finite lifetimes of coherent vortex clusters (as indeed observed in some cases in movie 1a, movie 1b).

Streamlines around an elliptic and a hyperbolic zero-acceleration point in the frame moving with each one of them are plotted in figures 6(c) and (d) respectively. Streamlines in the original global frame are shown in figure 6(b), which is a magnification of figure 6(a). It is indeed observed in these local frames (not in a global frame) that these points are respectively elliptic and hyperbolic stagnation points with corresponding streamline structures surrounding them. From our argument and evidence presented above, we expect these streamline structures to be significantly persistent in these frames. We also observe that the streamline structures in figures 6(c) and (d) to be similar, thus suggesting that the elliptic and hyperbolic zero-acceleration points marked in these plots move approximately together.

Figure 6. Streamlines with elliptic (circles) and hyperbolic (squares) zero-acceleration points of the coarse-grained velocity field. . Run C. The corresponding vorticity field is shown in figure 4(d). (a) Streamlines in the global frame introduced in section 2.3.2. (b) Magnification of the square sub-region in (a). (c) Streamlines in the local frame moving with the red marked elliptic zero-acceleration point. (d) Streamlines in the local frame moving with the red marked hyperbolic zero-acceleration point.

In summary, many (though not all) elliptic zero-acceleration points of the coarse-grained field are located at the centres of coherent vortices, while many (though not all) hyperbolic ones reside along boundaries between the vortices. Coherent vortices seem to survive with sufficiently longer lifetimes than their own eddy turnover times, and so also elliptic zero-acceleration points located at their centre. On the other hand, some hyperbolic zero-acceleration points are created by a combination of a few vortices, and therefore their lifetime is determined by that of the coherent cluster of vortices, which seems, from casual observation of our movies, to be comparable to the eddy turnover time of the vortices. In the local frames moving with a zero-acceleration point, this point is a stagnation point and streamlines around it are approximately persistent. Incidentally, elliptic or hyperbolic zero-acceleration points with very short lifetimes, i.e. blinking points in movie 1, cannot play any role in the particle diffusion because streamlines around them are not persistent even in the frame moving with them.

2.3.2. In the global frame: stagnation points and statistical persistence hypothesis Since the evolution of particle pair separations is frame-independent, it is possible to calculate statistics of such separations by following every single pair in a different local frame. The local frame picture developed in the preceding subsection is therefore appropriate for the calculation of particle pair diffusion. However, as pointed out in [11] and as examined in detail in the present paper, statistics of particle pair diffusion should depend on the global frame streamline topology. For example, the scaling of mean square separation of particle pairs is related to the fractal dimension of the spatial distribution of stagnation points (see (50) in section 7). However, streamlines and stagnation points are not Galilean invariant, and it is crucial to identify the appropriate global frame where statistics of stagnation points and particle pair diffusion can be calculated together. The important assumption made here is that such a global frame exists, that it is the one which maximizes streamline persistence in some statistical sense, and that this persistence is long enough for global streamline topology to strongly influence turbulent particle diffusion. We call this assumption the statistical persistence hypothesis.

As described in the preceding section, maximization of streamline persistence is achieved by choosing a local frame which moves with a vortex, and the velocity of this local frame can be determined by the fluid velocity, u(x_a, t), at a zero-acceleration point located at the centre of the vortex. As illustrated in figure 6, in a global frame zero-acceleration points do not coincide with stagnation points. This implies that we lose some persistence of streamlines when we move from a local to a global frame. Nevertheless, our assumption implies that it is possible to identify a global frame where the loss of streamline persistence is statistically minimized. We expect the velocity of this global frame to be equal to the average velocity of local frames, u(x_a, t), over all zero-acceleration points with finite lifetimes. In isotropic turbulence, this global frame is the same as the one where the mean fluid velocity vanishes, and is independent of the cut-off scale, η_c = 2π/k_c. This is the global frame adopted here.

The statistical persistence hypothesis is strong, and its justification is left for future studies. Nevertheless, it implies that pair separation events caused by highly curved streamlines around hyperbolic stagnation points dominate pair diffusion statistics. Hence, the statistics may be captured by the spatial distribution of stagnation points. We count the number of stagnation points of coarse-grained velocity fields for different cut-off scales and different runs A, B, C and D, and plot the results in figure 8. This figure supports, in the case of and d = 2, the formula

with

where and η denote respectively the outer and inner length scales of the scaling range, L₀ = 2π is the system size, and p is the exponent in E(k)~k ^- p(1 < p < 3). This general formula has previously been shown to be in agreement with kinematic simulations, three-dimensional DNS turbulence and wind-tunnel turbulence data [11]. In the case of figure 8, the outer length scale is the integral scale of the second-order longitudinal velocity correlation function, and the cut-off length scale η_c serves as the inner length scale η. Equation (11) is the quantitative expression of the cat's eyes within cat's eyes picture of streamline topology in [10] which we confirm here in figure 7. It is worth mentioning that the number of elliptic stagnation points must be equal to the number of hyperbolic points plus one in two-dimensional flows according to Euler's relation [21].

Figure 7. Streamlines and stagnation points of the coarse-grained velocity field in the global frame for different cut-off wavenumbers: (a) , (b) , (c) and (d) . (a), (b) and (c) are respectively magnifications of the square sub-regions in (b), (c) and (d). Square sub-regions are chosen around a cluster of stagnation points. Run C.

Figure 8. Number of stagnation points in the entire computational domain (i.e. in an area equal to L02) of the coarse-grained velocity field as a function of cut-off length scale ηc = 2π/kc. Solid circles, run A; open circles, B; solid squares, C; open squares, D. The straight line is . is the integral scale.

This global picture is quite powerful, and predictions based on it are shown in section 4 to be consistent with those based on the local frame picture.

Before closing this subsection, we compare the number of stagnation points to that of zero-acceleration points. We count the number N_a of zero-acceleration points in the coarse-grained field by the use of a Newton-Raphson method. The result is plotted in figure 9, and is well fitted by the scaling,

where η_c = 2π/k_c is the cut-off length scale and L₀ = 2π is the system size. C_a is a constant which we estimate numerically to be C_a≈0.7 (the coefficient is slightly larger for hyperbolic zero-acceleration points than for elliptic ones). The scaling (13) is significantly steeper than (11) for d = 2 and , and implies that the mean distance, , between zero-acceleration points is nearly equal to η_c irrespective of the cut-off length scale. For comparison, the mean distance, L₀/N_s^(1/d), between stagnation points is proportional to (for ) which is the Taylor scale when η is the Kolmogorov scale.

3. Physical picture of particle pair diffusion

3.1. In the local frame

We now explain how particle motion in a local frame, where streamlines are approximately persistent, can be described in terms of streamline topology and two elementary types of motion. Recall that the local frame depends on the coarse-graining scale; and that it moves with an elliptic zero-acceleration point located at the centre of a vortex at that coarse-graining scale. In this frame, persistent streamlines, which approximately coincide with particle trajectories, are elliptic around the zero-acceleration point (as observed in figure 6(c)), and particles swirl around it. This trapping by a vortex is the first elementary type of particle motion, and we call this vortex a patron of the particle. (Other works, see e.g. [22], have also related one particle diffusion to vortex trapping.) We expect the majority of particles to belong to a patron at each scale and at any one time.

However, at each coarse-graining scale, there exist hyperbolic zero-acceleration points between patrons (as observed in figure 4), and the local persistent streamlines are hyperbolic around such points (see figure 6(d)). Hence, when a particle belonging to a patron encounters a hyperbolic point of the coarse-graining scale, it may leave the patron for another adjoining patron of the same scale. This patron change is the second elementary motion of particles.

Now let us describe particle pair diffusion in terms of the two elementary motions that we have just identified. Suppose that the pair separation is Δ. This means that both particles belong to a same patron of scale Δ. As described in section 2.3.1, this patron is defined by an elliptic zero-acceleration point of the coarse-grained field with cut-off scale η_c = Δ. As long as they belong to the same patron, their distance remains O(Δ). However, when one of the particles encounters a hyperbolic zero-acceleration point and changes its patron, their distance suddenly becomes larger than Δ. However, they now belong to a common patron at a larger scale, ξΔ, where ξ > 1. The subsequent pair evolution is similar but now determined by zero-acceleration points at the larger scale ξΔ.

This particle-with-patron picture is observed in our simulation shown in figure 10 and movie 2, where we follow the two blue particles. First (from t = 0 to 0.225) they are trapped by a common patron. Notice the yellow elliptic zero-acceleration point at the patron's centre and the long time persistence of this point and patron. Notice also how the pair is trapped by this patron until one of the particles enters the influence zone of an adjoining hyperbolic zero-acceleration point at time t = 0.225. The change of patron can also be seen with more detail in figure 11 where we also plot the local frame streamlines at t = 0.3. Since streamlines are hyperbolic around the red square in figure 11, one of the particles leaves the patron to eventually be trapped by another one, and for the pair to remain trapped by a larger scale patron.

Figure 10. Temporal evolution of a particle pair (large blue balls; the small balls map the trajectory of each particle between t - 0.075 and t) and of the coarse-grained vorticity field (shade of red/blue indicates magnitude of positive/negative vorticity). Yellow/black dots are elliptic/hyperbolic zero-acceleration points. Run C. The size of the area plotted is L0/4 = π/2. The eddy turnover time is . See also movie 2 which is without zero-acceleration points.

This scale-by-scale separation based on the particle-with-patron picture provides a concrete physical mechanism of particle pair diffusion in the inertial range of two-dimensional turbulence. Note that this picture is described based on the persistent streamlines rather than the high straining regions, although strain is likely to take large values around some hyperbolic zero-acceleration points.

3.2. In the global frame

The statistical persistence hypothesis corresponds to the picture where the velocity field is frozen during a typical time scale of particle pair separation and the particle trajectories approximately coincide with streamlines of the frozen field. In other words, a particle pair separates in high-curvature regions of streamlines, which is likely to be observed around hyperbolic stagnation points of this frozen field. Here, we assume that most intense particle separation events occur when pairs encounter hyperbolic stagnation points. Note that hyperbolic stagnation points are associated with a cut-off scale η_c; they are the stagnation points which disappear when the cut-off scale increases above η_c. Pairs of separation Δ are most likely to have their separation increased to ξΔ(ξ > 1) by stagnation points of scale η_c = Δ because, on account of equation (11), they are more numerous than larger scale stagnation points. Stagnation points of smaller scales cannot increase the separation of the pair.

The pair distance remains ξΔ until the pair encounters a hyperbolic stagnation point of length scale η_c = ξΔ.

3.3. Geometrical (fractal-like) separation

Both pictures described in sections 3.1 and 3.2 require that the particle pair separation Δ evolves in a series of bursts as follows:

with a constant ξ ( > 1). This fractal-like separation is quite reasonable when the turbulent field is statistically self-similar. The assumption of a multiplicative process with constant ξ reflects the self-similar multiple-scale structure of the turbulence and of turbulent diffusion, and is expected to hold in the inertial range. This assumption represents our reformulation of the locality-in-scale hypothesis. We denote by T_ξ(ξⁿΔ₀) the characteristic time for pair separations to increase from ξⁿΔ₀ to ξ^n + 1Δ₀. T_ξ can be interpreted to be the doubling time introduced by Artale et al  [23] and Boffetta et al  [24]. So hereafter, we refer to T_ξ as the doubling time although this nomenclature is strictly only appropriate for ξ = 2. An important difference with the original use of doubling time concept is that ξ is not a free parameter here, but a constant charactering the multiplicative process of particle pair turbulent diffusion. We determine this constant using DNS in section 6.

4. Consistency of local and global viewpoints of particle pair diffusion

In the previous section, we have proposed two physical pictures of turbulent particle diffusion, one in the local (section 3.1) and one in the global (section 3.2) frame. The purpose of this section is to consider the consistency of the two viewpoints. For this purpose, it is helpful to address a generalized d-dimensional (d = 2 or 3) turbulent field, where the energy spectrum E(k) is proportional to k ^- p and not restrict ourselves to the special case of d = 2 and . The typical (largest) velocity scale being denoted by u ' and the inner/outer length scales being denoted by η and , the energy spectrum in the range can in general be expressed as

with a constant C_K. Note that this function reduces to (8) when on account of Taylor's relation . Here, the exponent is assumed to be

The exponent p must be larger than 1 to ensure finite energy in the limit of , and smaller than 2 for the enstrophy spectrum to be an increasing function of the wavenumber. To check the consistency of the two pictures, we use them both to estimate the doubling time T_ξ, and show that they lead to the same result irrespective of p and d.

4.1. Doubling time estimation from the local frame's viewpoint

As mentioned in section 3.1, the doubling time T_ξ(Δ) must scale with the eddy turnover time of vortices of size Δ. This eddy turnover time of vortices of scale 2π/k can be estimated as (k³E(k)) ^- 1/2 as long as the exponent p is between 1 and 2, see (16). Thus the doubling time T_ξ(Δ) may be expressed by

because we assume the energy spectrum to be given by (15).

4.2. Doubling time estimation from the global frame's viewpoint

For the physical picture in the global frame (section 3.2), the doubling time T_ξ(Δ) is estimated as the average time required for pairs of separation Δ to encounter a hyperbolic stagnation point of scale larger than Δ.

We therefore need to calculate the mean distance between such hyperbolic stagnation points, which is

where we have used (11). In the global frame, where the mean velocity vanishes, the r.m.s. particle velocity is u '. Hence, the doubling time T_ξ(Δ) may be estimated by

which implies

Here, C_ξ is a constant proportional to C_s ^- 1/d.

4.3. Consistency

The local frame observation that, on average, particles remain in patrons of size Δ for a time T_ξ(Δ) given by (17) implies that, in the global frame, particle pairs of separation Δ should travel together for the same time without encountering hyperbolic stagnation points. Having shown, in the global frame, that pairs of separation Δ travel together for a time T_ξ(Δ) given by (20), we are led to conclude that the two time scales, (17) and (20), must be equal.

It is interesting to remark that the equivalence between (17) and (20) implies (12), i.e. p + (2/d)D_s = 3. This is remarkable because this relation between p and D_s has been established mathematically, and it has been validated in three-dimensional DNS and wind tunnel turbulence [11] as well as here in two-dimensional DNS turbulence. The mathematical arguments that have been used to obtain p + (2/d)D_s = 3 rely on Gaussianity of velocity increments (see [11, 25] and section 16.2 of [26]), but it is by no means clear that Gaussianity is a necessary condition for this relation, as the turbulence validations indeed indicate.

5. A simple model of particle pair diffusion

We now construct a simple model of turbulent particle pair diffusion based on the physical picture developed in the previous sections. The key ingredient of this model is that particle pair separates in stepwise bursts as described in (14). The central quantity that determines the statistics of particle pair separations is the doubling time T_ξ(Δ). We have shown in the previous section that the two viewpoints in sections 4.1 and 4.2 based on different (local and global) frames lead to the same estimation of T_ξ(Δ) given by (17) or equivalently (20). Now we return to the case of two-dimensional (d = 2) turbulence with . In this case,

where C_ξ is a constant proportional to C_s ^- 1/2. The only unknown physical parameters in this model are the coefficient C_ξ and the multiplicative constant ξ characterizing, respectively, the time and length scales of the multiple-scale stepwise evolution (14).

5.1. Mean square separation

Neglecting fluctuations around expression (21) for the doubling time T_ξ, we can estimate mean square separation at time t as follows. Suppose that at time t the pair separation is Δ = ξⁿΔ₀, where Δ₀ is the initial separation. Then,

which can be easily rewritten, in terms of the well-known Richardson constant G_Δ, as

with

The averaging in (23) comes from the fact that T_ξ in (22) are average time scales. Note that G_Δ~C_ξ ^- 3~C_s^3/2 and that expression (23) tends to the t³-law,

in the limit of

It is worthwhile, in passing, to observe that the doubling time T_ξ(Δ) may be regarded as the separation's correlation time. In the large-t limit, (22) leads to

which is a constant irrespective of time t. This observation explains why, as is often stated [7], the particle pair diffusion has a long memory of the order of t (see section 6.4).

5.2. PDF of separation

In the previous subsection we have derived the Richardson t³-law for mean square separations. In this section, we extend our approach and derive the full PDF of pair separations.

We define the probability Q_n(t) for a particle pair to be separated by a distance between Δ_n = ξⁿΔ₀ and Δ_n + 1. By definition, ∑_n = - ∞^∞Q_n = 1. Based on the physical picture described in section 3, we construct the equation for Q_n(t) as

Here the first (second) term on the right-hand side is proportional to the probability per unit time for the distance between a particle pair to change from Δ_n - 1 to Δ_n (from Δ_n to Δ_n + 1). The coefficients b_n must be proportional to the inverse doubling times, T_ξ(Δ_n) ^- 1. Therefore, using (21), b_n is expressed as

where B≡C_BC_ξ ^- 1 ε^1/3 is a constant (independent of Δ_n).

The above equation (28) leads to a new equation for the PDF P(Δ, t) of separations between particle pairs. First, note the relationship between Q_n(t) and P(Δ_n, t);

where α is defined by

Then, assuming the change of Q_n with n to be small, a Taylor expansion up to second order of the right-hand side of (28) leads to

for P(Δ, t). Omitting the second order of the Taylor expansion leads to (32) without the second term on the right-hand side. The solution of this first-order approximation leads to Richardson's t³-law (23) of Δ², but if P(Δ, 0) is a delta function then it remains so for ever. The second-order term allows P(Δ, t) to evolve away from a delta function as seen in (34). Keeping higher orders leads to small corrections which we hope do not accumulate upon time integration. Note that the new equation (32) has two physical parameters: the coefficient B of the doubling time, and the length ratio ξ (or α = log ξ). Our equation (32) is a generalization of the Richardson diffusion equation (2), which in two-dimension (d = 2) is

where C_R is a constant proportional to ε^1/3 and related to G_Δ. This generalization is expressed by the presence of parameter α which quantifies the multiplicative particle pair diffusion process assumed in the model. Our equation becomes identical to that of Richardson's only for a specific value of .

The solution of (32) with the initial condition P(Δ, 0) = 2δ(Δ) is given by

where A is a time-independent constant obtained from

with

It is interesting that the t³-law for the mean square separation is consistent with (32) irrespective of the parameter α. Indeed, the solution (34) leads to

with

Note that the Richardson constant G_Δ, as defined in (23) and (25), is expressed by

which is proportional to C_s^3/2 because C_ξ~C_s ^- 1/2.

Using (37), it is easy to show that the PDF of separation Δ normalized by its r.m.s. is time-independent, i.e.

has the form

with a prefactor

Note that the similarity form (41) depends only on α, and not on C_ξ and C_B.

6. Numerical verification

6.1. Mean square separation

Temporal evolutions of mean square separations Δ² are plotted in figure 12(a) for different initial separations Δ₀ = 2^jη_f (j = - 5,  - 4, ..., - 1 plotted with blue curves; j = 0, 1, 2 plotted with black curves). As reported in [6], a strong dependence on initial separation is observed. Indeed, when Δ₀ < η_f (blue curves in figure 12), the effects of initial separation contaminates the entire range of times and scales. It is possible to find a particular value of Δ₀ < η_f (e.g. Δ₀ = η_f/8 in figure 12(a)) for which the data might be fitted by a t³-law, but such a fit is spurious because it is dependent on Δ₀ being carefully chosen. Different power laws fit best for curves with smaller or larger initial separations, e.g. Δ₀ = η_f/4 and η_f/16 in figure 12(a). A similar observation has been made with kinematic simulations [27]. This observation has some bearing on reported claims of a t³-law obtained for one specific value of Δ₀( < η_f) [6, 7]. This same observation may also cast some doubt on the prediction [2] that the temporal evolution of the mean square separation during times when Δ(t) < η_f is exponential. We cannot confirm this exponential behaviour with the present DNS; the semi-logarithmic plots of Δ² versus time do not reveal clear straight lines.

Figure 12. (a) Temporal evolution of mean square separation. Thick black curves, Δ0≥ηf; thin blue curves, Δ0 < ηf. The dashed line indicates t3. Run C. The eddy turnover time is . (b) Normalized by the initial separation Δ0 and the energy dissipation rate ε. Red curve, prediction (43) with GΔ = 6.9.

The model proposed in the preceding section, as well as the Richardson diffusion equation and t³-law, are valid in the inertial range, i.e. for . Our DNS results plotted in figure 12(b) indicate that (23), and the implied collapse of the data where Δ² is normalized by Δ₀² and t is normalized by Δ₀^2/3ε ^- 1/3, cannot be observed for Δ₀ < η_f but can be observed for Δ₀≥η_f as long as Δ² is smaller than . In figure 12(b) we also compare the DNS results with our model's prediction (23) rewritten with the above normalization as

We plot the function (43) for the case where G_Δ = 6.9 with a red curve in figure 12(b), and found it to be consistent with our numerical results except in the vicinity of t = 0; although the discrepancy in the small t region may be emphasized in this logarithmic plot, it is fairly small. Note that dΔ²/dt = 0 in DNS, while (43) has a finite derivative at t = 0. The estimated value G_Δ≈6.9 is about twice that obtained in [6] (G_Δ≈3.8) for a t³ fit of the temporal evolution of Δ² for a particular (unjustified) choice of Δ₀( < η_f). A comparable value (≈3) can be obtained if we fit the third blue curve from the bottom in figure 12(a) with G_Δεt³. From laboratory measurements with ratio smaller than ours and that used in [6], the value 0.5 has been suggested [7] from another arbitrary value Δ₀(  η_f).

The temporal evolutions of Δ² in the cases where the value of Δ₀ is larger than η_f reveal clear universality consistent with the model prediction (23), and equivalently (43), but we cannot observe the t³-law clearly. From figure 12(b) we infer that is not enough for the t³-law to be unambiguously observed. This implies that in the present DNS might be too small to capture the Richardson law. Yet, as suggested by some authors [8], the time-shifted plot may be helpful. Equation (43) suggests a time shift equal to G_Δ ^- 1/3. The results, where we adopt the value G_Δ = 6.9, are plotted in figure 13. In this figure, we plot four curves all corresponding to Δ₀ = η_f but for different resolutions. As the resolution increases (i.e. as the ratio increases) the t³-law becomes clearer.

6.2. PDF of separation

The main result of our model is the functional form of the PDF of separations. This PDF has the similarity form (41) when the separation is normalized by its r.m.s. We plot (41) in figure 14 together with DNS results of PDFs at times when Δ(t) is smaller than the integral scale. All the DNS PDFs collapse well, and they also are in good agreement with (41) with α = 1.3 which fits the data best. This fit is poor at the PDF tails, but the deviation becomes smaller as increases, and may therefore be due to finite-size effects. Incidentally, when plotted linearly (figure 15) the deviations of the PDF tails are negligible.

Figure 14. PDF of separation Δ normalized by its r.m.s. The initial separation is Δ0 = ηf. Solid circles, t = 0.05; open circles, 0.1; solid triangles, 0.15; open triangles, 0.2; solid squares, 0.25; open squares, 0.3; crosses, 0.35; Only the data for times when are shown. The solid curve is the model prediction (41) with α = 1.3; the dotted curve is the Richardson prediction corresponding to α = 1.5. (a) Run A. The eddy turnover time is . (b) Run B. . (c) Run C. . (d) Run D. .

The value α = 1.3 is surprisingly close to the value corresponding to the Richardson diffusion equation in two-dimensional turbulence. Note that α = 1.3 implies ξ≈3.6 because α = logξ.

6.3. Doubling time

The central concept underpinning our model is the doubling time, which is estimated in section 4 to be given by (21). In this section, we verify this formula by DNS. We estimated the average doubling time T_ξ(Δ) for ξ = 3.6. Recall that ξ is not a free parameter in our model, but a physical constant determined by the turbulence. In the preceding subsection, ξ≈3.6 is determined by the PDF of separations for our two-dimensional DNS turbulence. In figure 16, we plot T_ξ = 3.6(Δ) for three different initial separations Δ₀(≥η_f), and for two different values of (runs C and D). Note that T_ξ = 3.6(Δ) is plotted for Δ = 1.2ⁿΔ₀ (n = 0, 1, 2, ...) instead of ξⁿΔ₀ to make the plots denser. This figure confirms (21) as long as Δ is in the inertial range. Deviations can be observed but only at the very early stages of the temporal evolutions (Δ(t)  4Δ₀), beyond which the average doubling time becomes independent of Δ₀ and the resolution if it is normalized by ε ^- 1/3. Our data allow us to determine the constant C_ξ in (21) as C_ξ≈1.1.

Figure 16. Average doubling time. ξ = 3.6. Solid circles, Δ0 = ηf; open circles, 2ηf; solid squares, 4ηf. The solid line is Tξ = 1.1ε - 1/3Δ2/3. (a) Run C. (b) Run D.

From the prediction (39) of our model, C_B in b_n = C_B(T_ξ(Δ_n)) ^- 1 (see (29)) can be determined in the case where ξ≈3.6. We have determined G_Δ≈6.9 for these simulations (see section 6.1) and therefore use (39) with ξ≈3.6 (α≈1.3) to determine C_B≈1.3 because C_ξ≈1.1.

6.4. Correlation functions

In this subsection, we examine the correlation function

of the Lagrangian velocity difference vector, δv(t) = u(x⁽¹⁾(t), t) - u(x⁽²⁾(t), t), of a particle pair. This correlation function is important because it determines the temporal evolution of mean square separation through the exact relation

The correlation function H(t, τ) estimated by DNS (run C) is plotted as a function of the time increment τ in figure 17(a) for different times t and different initial separations Δ₀. It is seen in this figure that the correlation time is longer for larger separations.

Figure 17. (a) Correlation function of velocity differences for different times and initial separations. Red symbols, Δ0 = ηf; blue, 2ηf; black, 4ηf. Solid circles, t = 0.05; open circles, 0.1; solid triangles, 0.15; open triangles, 0.2; solid squares, 0.25. Run C. The eddy turnover time is . (b) Time increment is normalized by the auto-correlation time τc(t).

As is mentioned in section 5, the doubling time can be regarded as the auto-correlation time of the separation process. Hence, the correlation time τ_c(t) of the particle pair at the distance Δ may be scaled as

using (23). According to this argument, the correlation function H(t, τ) should have a similarity form

We plot the correlation function estimated by DNS in figure 17(b), where the time increment τ is normalized by the predicted correlation time τ_c(t) with G_Δ = 6.9. This figure strongly supports the similarity (47); H(τ/τ_c) is identical irrespective of the time t or the initial separation Δ₀. Incidentally, the scaling τ/τ_c of the correlation time tends to τ/t in the large t limit, as has been observed in [7, 27].

6.5. Doubling time dependence on the stagnation point distribution

The global frame picture in section 3.2 suggests that the statistics of particle pair diffusion are related to the spatial distribution of stagnation points. To verify this relationship, we follow [11] and consider particle motions governed by

instead of (4). Here, V is an unphysical velocity which acts only on particles but not on the flow. In other words, the particles integrated by (48) are not fluid elements. Note that this addition of a uniform velocity is not a Galilean transformation and does not affect the persistence of the streamlines even though it modifies them. Hence, the global frame where the statistical persistence of streamlines is maximized is the same as when V = 0. It is the frame where the spatial average of u(x, t) vanishes. The statistical persistence hypothesis in this frame implies that particle trajectories are now influenced by streamlines of u + V. The number N_s of the stagnation points of u + V in this frame is plotted in figure 18(a) which shows that N_s scales as (11) with and d = 2, but the coefficient C_s decreases as V = |V| increases. The average doubling time T_ξ(Δ) is plotted in figure 18(b). The results in this figure support the scaling (21) of T_ξ(Δ) with an increasing prefactor C_ξ for increasing V. We note that this scaling extends over a smaller range of Δ for larger value of V. At the time of writing, we did not know how to explain the shrinking of this range from the global frame's viewpoint. However, from the local frame's point of view, it may be due to the disruption of the particle-with-patron picture at scales where the swirling velocity (E(Δ ^- 1)Δ ^- 1)^1/2 around the patron is exceeded by V.

Figure 18. (a) Number of stagnation points of the velocity field u + V. Solid squares, V = 0; open circles, V = 0.5u '; solid circles, V = u '. The solid lines are with Cs = 0.6, 0.5, 0.3. Run C. (b) Average doubling time for different V. Symbols are the same as in (a). ξ = 3.6. The solid lines are CξΔn2/3 with Cξ = 1.1, 1.3 and 1.7.

The results in figure 18 imply a negative correlation between the number of stagnation points and the doubling times of particle pair diffusion. This correlation is qualitatively consistent with the global frame picture of particle pair diffusion, where particle pairs separate only when they encounter a hyperbolic stagnation point. For a quantification of this effect, we plot C_ξ as a function of C_s in figure 19. This plot is consistent, though in a narrow scaling range, with the scaling C_ξ~C_s ^- 1/2 suggested in section 4.2.

7. Particle pair diffusion in general cases

The global frame picture (section 3.2) of particle pair diffusion is based on the statistical persistence hypothesis of streamline topology. This hypothesis allows us to predict the doubling time T_ξ(Δ), which is the central Lagrangian quantity of our reformulation of the locality-in-scale hypothesis in section 3.3 and of the model described in section 5, as a function of the number of stagnation points in the global frame. Our approach is fundamentally based on the scaling law (11), which leads to (20).^Note4  Hence, the statistical persistence hypothesis allows us to generalize our approach to any value of d and D_s. The strength of this approach is that Eulerian statistics of stagnation points may be more accessible than Lagrangian statistics.

With the form (20) of the doubling time, the argument of section 5.1 can be repeated to give

with

where .

Similarly, the argument in section 5.2 can also be generalized and leads to,

with a coefficient (C_B is a dimensionless constant) and an exponent

which characterizes the spatial density of stagnation points in the d-dimensional space. The solution of this equation for initial condition P(Δ, 0) = 2δ(Δ) is

with a normalization factor

Here, . With solution (53) we recover (49) with Δ₀ = 0, and a generalized Richardson constant

where the constant G₀ is defined by

We obtain, from (53) and (49) with Δ₀ = 0, a similarity form of the PDF, which is independent of time t and the Richardson constant G_Δ, specifically

with

All the above equations (49), (51), (53), (54), (55), (56), (57) and (58) respectively reduce to (23), (32), (34), (35), (39), (38), (41) and (42) when and d = 2 (or D_s = 2 and d = 3) and therefore .

An important point of this generalization is (50) which relates the Lagrangian quantity γ to the fractal dimension D_s of spatial distribution of stagnation. This relationship between γ and D_s has been confirmed by a series of kinematic simulations by systematically changing the exponent p of the energy spectrum so as to change D_s [11]. Note that the relation (12) is obeyed by these kinematic simulations and that, with (50), it implies

which is known from Kolmogorov-type dimensional analysis based on the locality-in-scale hypothesis [10].

It is also important to clarify the relationship between our PDF equation and Richardson's diffusion equation. The PDF equation (51) is the same as Richardson's diffusion equation (2), only if

D_s = 2 and are indeed observed in three- and two-dimensional DNS of homogeneous isotropic turbulence, respectively.

8. Concluding remarks

This paper puts forward a new approach to turbulent particle pair diffusion based on the notion of persistent streamline topology. We have identified two frames (one local and one global) where streamlines are approximately or statistically persistent (see section 3). The local viewpoint is based on the idea that the velocity field in the frame moving with a coherent vortex is approximately stationary. We defined such local frames unambiguously in terms of zero-acceleration points. The global viewpoint, on the other hand, is based on the statistical persistence hypothesis, which states that the persistence of streamlines is maximized in a statistical sense in the global frame moving with the average velocity of all local frames.

These local and global viewpoints lead to the same characteristic time and length scales for particle pair diffusion (see section 4). The simple model, developed in section 5, based on these time and length scales, describes well the statistics of particle pair diffusion in the inverse energy cascade range of two-dimensional turbulence (see section 6). The two viewpoints are complementary. The local frame viewpoint provides a concrete picture (the particle-with-patron picture) of particle pair diffusion (illustrated in figure 10), but cannot be used to predict overall statistics of particle diffusion in terms of an easily measurable Eulerian quantity. The global frame viewpoint is based on the statistical persistence hypothesis, which we are currently unable to prove or disprove, but leads to predictions such as (50) which relates the exponent of the mean square separation of particles (a Lagrangian quantity) to the fractal dimension of the spatial distribution of stagnation points (an Eulerian quantity).

Acknowledgments

The authors acknowledge support from the Japanese Ministry of Education, Culture, Sports, Science and Technology, and from the Royal Society of London.

References

[1]

Richardson L F 1926 Proc. R. Soc. A 110 709 
CrossRef

[2]

Batchelor G K 1952 Proc. Camb. Phil. Soc. 48 345 
CrossRef

[3]

Durbin P A 1980 J. Fluid Mech. 100 279 
CrossRef

[4]

Obukhov A 1941 Izv. Akad. SSSR, Ser. Geogr. Geofiz. 5 453 

[5]

Batchelor G K 1950 Q. J. R. Meteorol. Soc. 76 133 
CrossRef

[6]

Boffetta G and Sokolov I M 2002 Phys. Fluids 14 3224 
CrossRef

[7]

Jullien M-C, Paret J and Tabeling P 1999 Phys. Rev. Lett. 82 2872 
CrossRef

[8]

Ott S and Mann J 2000 J. Fluid Mech. 422 207 
CrossRef

[9]

Fung J C H, Hunt J C R, Malik N A and Perkins R J 1992 J. Fluid Mech. 236 281 
CrossRef

[10]

Fung J C H and Vassilicos J C 1998 Phys. Rev. E 57 1677 
CrossRef

[11]

Dávila J and Vassilicos J C 2003 Phys. Rev. Lett. 91 144501 
CrossRefPubMed

[12]

Kraichnan R H 1967 Phys. Fluids 10 1417 
CrossRef

[13]

Batchelor G K 1967 Phys. Fluids 12 233 (Suppl. II) 

[14]

Borue V 1994 Phys. Rev. Lett. 72 1475 
CrossRefPubMed

[15]

Sukoriansky S, Galperin B and Chekhlov A 1999 Phys. Fluids 11 3043 
CrossRef

[16]

Boffetta G, Celani A and Vergassola M 2000 Phys. Rev. E 61 R29 
CrossRef

[17]

Paret J and Tabeling P 1997 Phys. Rev. Lett. 79 4162 
CrossRef

[18]

Bernard D 1999 Phys. Rev. E 69 6184 
CrossRef

[19]

Weiss J 1991 Physica D 48 273 
CrossRef

[20]

Miura H and Kida S 1997 J. Phys. Soc. Japan 66 1331 
CrossRef

[21]

Flegg H G 1974 From Geometry to Topology  (London: The English Universities Press) 

[22]

Provenzale A 1999 Ann. Rev. Fluid Mech. 31 55 
CrossRef

[23]

Artale V, Boffetta G, Celani A, Cencini M and Vulpiani A 1997 Phys. Fluids 9 3162 
CrossRef

[24]

Boffetta G, Celani A, Crisanti A and Vulpiani A 1999 Europhys. Lett. 46 177 
IOPscience

[25]

Orey S 1970 Z. Wahrscheinlichkeitstheorie 15 249 
CrossRef

[26]

Falconer K 1990 Fractal Geometry, Mathematical Foundations and Applications  (New York: Wiley) 

[27]

Nicolleau F and Vassilicos J C 2003 Phys. Rev. Lett. 90 024503 
CrossRefPubMed

Notes

Note2  Except where otherwise stated, by particles we mean `fluid elements' in this paper.

Note3  This could have been expected from the fact that the pressure field has a local minimum at the centre of the vortex (defined as a region with concentrated vorticity and swirl), because of the balance between the centrifugal force and the pressure gradient, and from a≈ - (1/ρ)∇p if the viscous terms and external forcing are negligible. Such considerations might be generalizable to three-dimensional turbulence (see [20]).

Note4  As shown in section 4, (20) is identical to the eddy turnover time of vortices of size Δ if the formula (12) between D_s, the fractal dimension of the spatial distribution of stagnation points, and p, the exponent of energy spectrum, is valid.