Experimental and numerical study of the Lagrangian dynamics of high Reynolds turbulence

We introduce an original acoustic method to track the motion of tracer particles in high Reynolds number turbulent flows. We present in detail the experimental technique and show that it yields a measurement of the Lagrangian velocity variations of single particles, resolved across the inertial range turbulence. Second-order quantities such as the velocity autocorrelation function and time spectrum are in agreement with Kolmogorov 1941 phenomenology. Higher-order quantities reveal a very strong intermittency in the Lagrangian dynamics. Using both the results of the measurements and of direct numerical simulations, we show that the origin of intermittency can be traced back to the existence of long-time correlations in the dynamics of the Lagrangian acceleration. Finally, we discuss the role played by vortices in the Lagrangian dynamics.

have also been developed to model the dispersion of particles in turbulent flows. In contrast, most experimental studies are still about Eulerian quantities, because of the technical difficulties associated with the implementation of particle measurements. One difficulty is to single out a fluid particle. Two options have been investigated: giving a real fluid particle a different physical parameter from the surrounding fluid (temperature, colour etc) or seeding the flow with small neutrally buoyant particles. The first attempts trace back to the 1920s with Taylor's work to relate Lagrangian velocity statistics to heat dispersion-the width of the thermal wake behind a heated wire in a wind tunnel is related to the autocorrelation function of the Lagrangian velocity [4]- [6]. This relation is restricted to the autocorrelation and is perturbed by diffusion effects and by the non-stationarity of the turbulence in the wind tunnel. Most of the recent experimental work in Lagrangian dynamics uses solid particles that are tracked with optical techniques. One is then confronted with the second technical difficulty in Lagrangian studies: obtaining a measurement of the motion resolved in time and space. Film cameras were used in the 1970s [7,8], but they were restricted to only a few images per second and hence to low Reynolds numbers. At the end of the 1980s, the development of computers and digital cameras gave access to moderate Reynolds numbers and to more complex statistical quantities [9]. The 3D particle tracking velocimetry technique [10] is especially suited for Lagrangian studies but so far has been used only with standard camera rates (i.e. 25 frames s −1 ) which limit the measurements to moderate Reynolds numbers [11]. The main results of all of these experiments so far are that the Lagrangian velocity distribution is Gaussian and that the velocity autocorrelation coefficient is likely to decrease exponentially with the Reynolds number. More recent experiments use ultra-fast detectors to access particle accelerations with a very accurate time resolution [12]- [14], but the measurements are limited to very short time scales (a few Kolmogorov times τ η ). They have observed a highly non-Gaussian distribution of acceleration components [15].
In parallel to these few experiments, the development of direct numerical simulations (DNS) has produced Lagrangian data from simulated turbulent flows. Although still restricted to moderate Reynolds numbers, they yield a lot of information as all dynamical variables are accessible [16]- [18]. These numerical experiments have revealed that (i) the Lagrangian velocity autocorrelation function has an exponential decrease with a time T L characteristic of the large-scale motion; (ii) the direction of the Lagrangian acceleration decorrelates within a few Kolmogorov times while the magnitude remains correlated for times of order T L and (iii) the time increments of Lagrangian velocity have probability density functions (PDFs) which evolve from a Gaussian distribution at integral scale to one with wide stretched exponential tails near τ η .
It is the aim of the method presented here to investigate experimentally the above findings, using a resolved measurement that covers the entire inertial range of time scales of Lagrangian motion, in high Reynolds number flows. To this end, we have developed a new experimental technique based on the use of ultrasound [19]. The Lagrangian velocity is extracted from the Doppler frequency shift of ultrasound scattered by the solid particles-this step required the development of a high-resolution spectral analysis algorithm [20]. The main advantage of ultrasound is that, using only one emitter and one receiver, it is possible to have a large measurement volume without degrading the velocity resolution. The velocity is directly estimated from the Doppler shift, and not after differentiation of the particle position along its trajectory. It means that the velocity resolution and the time resolution are independent of the size of the measurement volume, which is not the case with digital cameras (for a given pixel resolution, increasing the measurement volume degrades the spatial resolution and thus the The von Kármán swirling flow. The flow is generated by two counterrotating discs of radius 9.5 cm, with a 18 cm separation. The fluid is water, enclosed in a cylindrical tank of radius 10 cm so that the total volume is 9 l. The fluid temperature is regulated to ambient temperature by an external recirculation loop. The discs are driven by two 1 kW dc motors whose rotation frequency is kept constant. A measured trajectory of a 250 µm particle is displayed. velocity resolution). As this method is quite new, we shall compare our essential results with DNS data of Lagrangian motion.
The experimental technique (including acquisition and signal processing) and the numerical simulation are described in section 1. In section 2, we compare the results of our measurements to the dimensional predictions of the Kolmogorov 1941 theory. The third section is devoted to the analysis of the Lagrangian acceleration. Section 4 reports the analysis of the intermittency of the increments in time of the Lagrangian velocity. Our findings suggest a dynamical origin which is analysed further in section 5 with the help of the DNS data.

Flow characteristics
Our measurements are performed in a confined von Kármán swirling flow, used as a model flow for homogeneous and isotropic turbulence studies [12], [21]- [23]. The fluid is water, set into motion by two counter-rotating discs (figure 1) with radius 9.5 cm, set 18 cm apart. The discs may be smooth or with eight blades of height 5 mm. The flow is enclosed in a cylindrical shell of inner radius 10 cm so that the total volume is 9 l. The temperature is regulated to ambient by an external cooling recirculation loop. The discs are driven by two 1 kW dc motors whose rotation frequency is kept constant. The values of the parameters are displayed in table 1. In a typical experiment (such as #3), the rms amplitude σ of velocity fluctuations is 1 m s −1 . The integral length scale L is close to 4.5 cm, the Taylor length scale λ is around 600 µm and the Kolmogorov dissipation length scale η is 20 µm. The Lagrangian integral time scale T L is 20 ms and the Kolmogorov time scale τ η is 0.2 ms. The average dissipation heat is measured by Table 1. Parameters of the different experiments. is the frequency of rotation of the discs, σ the standard deviation of the Lagrangian velocity, L the integral length scale, the flow power consumption per unit mass, η the Kolmogorov dissipation length, τ η the Kolmogorov dissipative time scale, λ the Taylor microscale, R λ = σL/ν the turbulent Reynolds number, S the number of Lagrangian trajectories analysed, corresponding to a number N 0 of data points, T L the Lagrangian integral time scale, T 2 the outer time scale in Sawford stochastic equation, σ a the standard deviation of the Lagrangian acceleration and C 0 the Kolmogorov constant in the Lagrangian velocity autocorrelation function. Error bars in the parameters can be estimated from the number of significant digits in the reported values. estimating the thermal power removed by the cooling system. The external cooling recirculation flow rate is tuned to stabilize the internal fluid temperature (in the middle of the flow) to the external ambient temperature so that there is no thermal flow between inside and outside the tank. The heat that is removed by the cooling system comes only from the turbulent dissipation.
In this way, one obtains an estimation of the global dissipation in the experiment. The problem with this method is that the dissipation is not homogeneous in the tank. For example, there is a strong shear between the back of the disc and the wall. In case of smooth discs, this shear may be the dominant source of dissipation as the entrainment by the disc is done through the boundary layer, which is not a very efficient process. This will lead to an overestimate of the dissipation in the measurement volume. Even between the discs the flow is not fully homogeneous. Thus, the estimate of the dissipation obtained from the global heating may differ strongly from the actual dissipation rate in the measurement volume. Another way to estimate the dissipation is to use the relationship where u is the standard deviation of the Eulerian flow velocity. The latter quantity is estimated from the Lagrangian velocity measurement, as the moments of the velocity are equal when estimated from Eulerian or Lagrangian measurements [24], so that u = σ. The integral length scale L can be estimated from the Eulerian autocorrelation function of the velocity. We have measured it in the centre of the flow, using air and a geometry very close to that used in water. The integral length scale L is found to be close to 4.5 cm and independent of the Reynolds number [23], in agreement with [14]. This allows the estimation of the dissipation rate as reported in table 1. One can see that, in most cases, the dissipation heat estimated from the heating is too large. In this paper, we will use the value of the dissipation rate computed from equation (1) as the actual value for in the measurement volume (this is the reason why numerical values may differ slightly from previously published results on the same experiment [25]).

The experimental technique
2.2.1. Principle of measurement. Fluid particles are singled out by the use of a small number of solid polystyrene spheres that are advected by the flow. Their density is 1.06, so that the advection by the turbulent flow greatly dominates any density effect (the buoyant forces are less than a thousand times smaller than the inertial forces). In all experiments (except #5), the particle diameter is 250 µm. This is smaller than the Taylor length scale, so we expect the particle to follow the fluid motion in the inertial range. We will return to this point in the next section and discuss it in detail. The measurement technique is based on the Doppler effect of the ultrasound scattered by the solid particle.Acoustic waves are emitted and detected by transducers located 18 cm off the axis of the cylindrical vessel as shown in figure 2(a). The emitter produces continuous monochromatic ultrasound at frequency 2.5 or 2.8 MHz (figure 3); the corresponding wavelength is around 0.6 mm in water. With square transducers of 2 mm sides, diffraction produces an emission beam with an aperture angle close to 35 • . In this way, a large volume of the flow is insonified and probed (the emission beam has a diameter close to 10 cm at the centre of the flow). When a polystyrene sphere enters the emission beam, it scatters the ultrasound. The scattered sound frequency is shifted by Doppler effect so that the instantaneous pulsation is where ω 0 is the emission angular frequency, q = k scat − k inc is the scattering wavevector and v(t) the instantaneous velocity of the particle. The scattered sound is recorded by an array of transducers. The array is used to improve the signal-to-noise ratio (SNR) of the measurement: the transducer detects independent realizations of the noise but identical Doppler signals. The receiver array is made of 9 square transducers, 2 mm each side, lying on a cross and separated by 0.5 mm ( figure 2(b)). The measurement domain, being the intersection of the two emitter and receiver beams, has a typical size of 10 cm. The measurement volume is located around the centre of the tank so that the average velocity is zero in the volume. Its size is larger than the integral length scale (L = 4.5 cm) so that one can expect a significant number of particles (a) Transverse cut of the experimental vessel at equal distance from the discs. The dark grey region around the vessel corresponds to a sound absorbing resin used to damp the reflections. The measurement region is the intersection of the emission and the reception beams. (b) Geometry of the transducer array: each transducer is a square with a 2 mm side and the spacing between transducers is 100 µm.
remaining in the ultrasonic beam for a long enough time to study the Lagrangian dynamics of the flow up to its large time scales. The results presented below show that the Lagrangian integral time scale is indeed smaller than the typical sweeping time of the particles in the measurement volume. We observed that the average duration of the recorded tracks is typically more than twice as long as the estimated Lagrangian integral time scale. An emitter-receiver pair yields the measurement of one component of the Lagrangian velocity. Two such pairs (figure 2(a)) are used in order to record two components of the velocity of the tracer particle. Both components lie in the transverse plane of the tank. It is important to note that when measuring the transverse components, we do not insonify any moving wall such as the rotating discs that generate the flow. This means that the echoes off the wall are not frequency-shifted (they are damped by over 30 dB by an adequate coating of the inner walls of the cylinder). Even if the echoes are over Principle of measurement. A monochromatic ultrasound wave is emitted by one transducer. Its size (2 mm) is comparable with the sound wavelength (3 MHz corresponds to a 0.5 mm wavelength), so that diffraction ensures that a large volume of the flow is insonified. As a solid particle enters the emission beam it scatters ultrasound with a Doppler shift that mirrors the particle velocity. The scattered sound is recorded by one receiving transducer or by an array of transducers.
three orders of magnitude larger than the particle contribution, they can be removed in the data processing using a digital stop-band filter at the emission frequency. In this way, one detects the signal scattered by the particle as long as the particle velocity is not strictly zero.

Acquisition.
Two velocity components are measured by recording the signal received by two arrays of nine transducers each. In order to avoid cross-talk problems, two distinct emission frequencies are used, equal to 2.5 and 2.8 MHz respectively. Because of the echoes on the wall, the signal has a large dynamic range. The largest part of the signal comes from the echoes and covers roughly 60 dB in magnitude. The scattering contribution from the particle has an amplitude that varies with the particle position in the ultrasonic beam (the amplitude of the beam decreases with distance to the transducer and also laterally due to the diffraction pattern) and its detection requires an extra 60 dB of dynamic range in the recorded transducer signal. The requirements for the acquisition device are as follows: 120 dB of dynamic range (to resolve the scattering signal) and a digitizing rate greater than 6 MHz (to sample the ultrasound accurately), for each of the 18 receiving transducers. As the Doppler shift induced by the particle motion is small compared with the frequency of the incoming wave ( ω/ω ≈ σ/c 0 , the flow Mach number), frequency multiplexing is an alternative to the use of 18 independent high-speed, high-resolution digitizers [26]. The output from the transducers in each array are thus combined into a single signal. The principle of operation of the frequency multiplexer is shown in figure 4: the central frequency of the signal received by the nine transducers is shifted into nine distinct spectral bands. To this end, nine very stable clock signals are synthesized at nine different frequencies from a master The output is summed with the output of the other channels and send to the acquisition device. High-level performance is achieved using ultra-fast and ultra-low noise components.  One can easily see the nine channels. Each channel is made of a sharp peak corresponding to the sound emission frequency. At the bottom of this peak is the part corresponding to the particle signal. It spans about 15 kHz (in this example the spectrum is computed over long time interval, so that the particle signal appears as a wide band and not as a single Doppler line), with an amplitude of 15 dB in this case. clock signal at 40 MHz. These clock signals are multiplied by the analogue acoustic signals coming from the nine ultrasonic receivers. The nine outputs are then summed and digitized by a Agilent e1430A 10 MHz, 23 bits digitizer. Using ultra-fast and ultra-low noise components, the multiplexing device satisfies both the requirements of rapidity and large dynamic range. The digitizer computes a numeric heterodyne detection at the centre frequency of the nine frequency bands. A typical power spectrum of the output of the multiplexing device is shown in figure 5.
One can see the nine spectral bands, each corresponding to one transducer. Every band is made of one large and narrow peak, due to the coupling emitter-receiver and echoes, and a wide, lowamplitude band due to the scattering by the particles. The individual spectra look alike because they measure the same scattered acoustic signal produced by a single particle moving in the measurement zone. After digital sampling, the first signal processing operation is to separate the nine signals, using numerical filtering and heterodyne detection. In its final form, the signal from each element of the transducers array correspond is complex and displays a spectrum centred about 0 Hz. The resulting data set is then analysed in two ways: (i) From each transducer signal, the Doppler shift from the sound scattered by a single particle moving in the measurement zone is extracted. This measures one component of the particle velocity. The advantage of the array detection is to obtain nine different realizations of the noise and, thus, to reduce the effective noise level by a factor of 3 compared with the use of a single transducer. (ii) The phase shift between the different transducers can be used to extract the direction of emission of the scattered sound. One line of transducers acts like a linear antenna so that by using the two arrays, one obtains the position of the particle.
We note here that the particle position and the particle velocity are obtained independently, a feature that provides a calibration and consistency check of the signal.

Signal processing
2.3.1. Principle. The velocity information is coded into the frequency modulation of the scattered ultrasound signal. This frequency modulation is highly non-stationary as the velocity may fluctuate on time scales as small as the Kolmogorov time scale (which is less than 1 ms), or more probably the eddy turnover time corresponding to the size of the tracer particle. At an emission frequency equal to 2.5 MHz, the typical order of magnitude of the Doppler frequency shift for a velocity fluctuation equal to 1 m s −1 is 3000 Hz. To estimate this frequency one needs a time window over which the signal may be considered to be stationary, i.e. spanning less than 1 ms. If one requires a 10% resolution on the frequency estimation (300 Hz) then the product 1 ms × 300 Hz = 0.3 is less than 1. It means that a detection based on Fourier analysis would not be efficient enough to satisfy this resolution requirement because of the uncertainty principle. The problem is made even more complex by the low SNR. Better performance can be obtained using parametric estimation. Indeed, in order to improve the resolution of the frequency estimation, one needs to make some hypotheses about the structure of the signal. After pre-processing filtering operations that remove the contributions from the echoes, the signal contains only the sound scattered by the particle, which can be written as The signal is the sum of M complex exponential functions that are generated by every particle present in the measurement volume. Each one of these exponentials has a complex amplitude a m which fluctuates slowly due to the position of the particle in the measurement region (and due to the filters used to remove the low-frequency component due to the echoes), and a frequency modulation f m (t) that is proportional to the velocity component. The signal also contains an amount of noise n(t) that comes mostly from the electronics and which will be assumed to be white, Gaussian and i.i.d. (independent, identically distributed). A very small number of particles is introduced into the flow, so that M < 2 for most events. The case M 2 is not conceptually different but is practically much more difficult because two (or more) particles may have the same velocity component at the same time so that several spectral components may have the same frequency. It raises problems such as changing the order of the estimation (M), and dealing with ambiguities to attribute a spectral component to each particle. As we are interested in singleparticle measurement, we lower the concentration of particle to have at most one particle in the measurement region for most events. The estimation of the different parameters (M, a m , f m and the noise variance) is done using an approximated maximum likelihood (AML) algorithm. This algorithm has been described in detail in a previous paper [20] and the reader is referred to it for technical details. The next subsection gives an example of extraction of the particle position and velocity from the scattered signal.

Particle velocity.
The acoustic signal (heterodyned at null frequency) is made of alternating periods of time with (at least) one particle in the measurement region (we call this period an event) and periods without any particles. The first step of the processing is to select the events. This is done by thresholding the sound intensity. An example of a typical experimental event is presented in figure 6(a). The real part of the signal is displayed. At small times, the particle has not yet entered the measurement region, so that the signal is essentially noise. At t close to 10 ms, a particle enters the measurement volume and one observes the slow evolution of the amplitude corresponding to the motion of the particle inside the measurement volume. The frequency modulation due to changes in its velocity are also easily observed. Figure 6(b) displays a time-frequency picture (spectrogram) of the same signal. One column of the picture corresponds to a small time window of the signal centred around the corresponding time. The grey levels correspond to the amplitude of the power spectrum in this window. One can see a line of high energy that corresponds to the time evolution of the Doppler frequency shift due to the particle motion. It is obvious that the resolution of this technique is not high enough to extract the fast oscillations that are observed about t = 80 ms. The solid line is the output of the AML algorithm, using a time window of 15 consecutive samples. It can be seen that the algorithm extracts the frequency modulation even when it displays fast oscillations. The time response is restricted in part by the length of the time window and by the amplitude of the model noise [20,26]. The resolution of the AML algorithm in estimating the Doppler shift is expected to be of the order of a few per cent of the velocity variance for the results presented below.

Particle position.
The same algorithm can be used to extract the position of the particle. Indeed, let us consider N transducers in a line array, separated by a distance d. A monochromatic source at frequency ω, set at infinity in a direction that makes an angle θ with the line of transducers, creates an amplitude vector on the transducer array given by where c is the speed of sound. One sees that such a source induces a spatial frequency on the line of transducers of Its measurement is equivalent to the estimation of the Doppler shift discussed above, and the AML algorithm can be used to extract the angle θ directly. Our set-up in which each transducer array is made of two perpendicular lines of transducers allows the measurement of two such angles with each cross array and, thus, defines a direction in space. Using two cross arrays, one obtains the position of the particle as the intersection of the two lines given by each of the arrays.
An example of such a trajectory is given in figure 7. The resolution of the position measurement has not been studied in detail. The resolution of the x and y components is expected to be a fraction of a millimetre and the z component resolution is about 1 mm due to the use of only four transducers in one of the arrays. A good position resolution would require the use of more transducers as it depends mainly on the aperture of the array.

Numerical simulation
Several of our results will be compared with the analysis of DNS data which will also be used to link the properties of Lagrangian motion to the Eulerian characteristics of the flow field. This section describes the numerical method. Numerical simulations refer to the integration of the incompressible Navier-Stokes equations in a cubic domain with 2π-periodic boundary conditions, using a parallel distributedmemory pseudo-spectral algorithm. The time evolution of Fourier modes is performed by a (second-order) leap-frog scheme. The large-scale kinetic-energy forcing is adjusted at each time step by scaling the amplitudes of modes 1.5 k < 2.5 uniformly (phases are left to fluctuate freely), such as to compensate exactly for losses due to the viscous dissipation. This forcing scheme permits a rapid relaxation to stationarity [27]. Some Eulerian statistics of Run2 are reported in [28]; in particular, the Kolmogorov 4/5 law (including the viscous contribution) for isotropic and homogeneous turbulence has been checked carefully.
Fluid-particle trajectories are resolved in time by a second-order Runge-Kutta scheme [16]. Cubic spline interpolation is used to evaluate Lagrangian quantities from Eulerian fields. For Run1, 5000 particles, uniformly distributed at the initial time, are tracked for about 50T E (T E is the Eulerian integral-scale eddy turnover time). For Run2, 10 000 particles are tracked for about 15T E . Along the fluid-particle trajectories, the velocity v, the gradient of pressure ∇p, the dissipation ν v and all derivatives of the velocity ∂ j v i are recorded at each time step. The parameters of the numerical simulations are given in tables 2 and 3.

Kolmogorov's 1941 theory
The Kolmogorov 1941 theory (K41) addresses mainly the Eulerian velocity field and relies on the von Kármán-Howarth relationship [29]. In the case of one-particle Lagrangian velocity, there is no equivalent theoretical relation but one can derive dimensional predictions in the same spirit as the K41 theory. For example, one expects the Lagrangian second-order structure function to behave linearly in both the time scale τ and the mean dissipation rate in the inertial range: Parameters and Eulerian characteristic scales of simulated flows: u rms is the root-mean-squared velocity, ν the kinematic viscosity, ε the mean dissipation, η ≡ (ν 3 /ε) 1/4 the Kolmogorov dissipative scale, λ ≡ 15νu 2 rms /ε the Taylor micro-scale and R λ ≡ u rms λ/ν is the Taylor-scale Reynolds number.
Run1 128 Table 3. Eulerian and Lagrangian characteristic times of simulated flows. T E ≡ L/u rms is the integral-scale eddy turn-over time, T L the Lagrangian velocity correlation time and τ η ≡ (ν/ε) −1/2 the Kolmogorov dissipative time. δt is the time step of the simulation. It follows from this relation that the power spectrum of the Lagrangian velocity should scale as This section is devoted to a comparison of our measurements with these dimensional K41 predictions. Figure 8 displays the PDF of one component of the Lagrangian velocity at R λ = 810. It is compared with a Gaussian function of the same variance. The Gaussian shape for the Lagrangian velocity PDF is expected from the Eulerian measurements. Indeed, for homogenous and isotropic turbulence, the theory predicts that any single point statistical quantities should be identical for the Lagrangian or Eulerian fields [24]. There is a general agreement that the Eulerian velocity has a Gaussian distribution for homogeneous and isotropic turbulence [29]. Thus, the same distribution is expected in the Lagrangian case. One can see that, the agreement is very good, except for strongly negative values. The departure from the Gaussian shape observed for negative values may be due to the non-homogeneity of the flow as there is a mean recirculation of the fluid (generated by the differential rotation of the discs and by the centrifugal pumping). The decrease in probability for velocities close to 0 m s −1 is due to the pre-processing filtering operations that remove the contributions from the echoes at zero frequency shift, but which also remove the signal part. This cannot be avoided as the echoes create a large spectral component around zero frequency whose magnitude is typically three orders of magnitude larger than the particle signal. It makes the extraction of the particle Doppler impossible for zero velocities. For the data set presented in figure 8, the non-accessible velocities are less than 0.1 m s −1 . For other Reynolds number, the observations are the same once rescaled by the velocity variance. Note that the probability of observing a zero velocity is not strictly zero, as the spectral estimation algorithm is coupled with a Kalman filter that partially reconstructs the zero-crossings of the frequency modulation. We do not expect this gap to have a significant effect on the analysis of the Lagrangian dynamics. Indeed, we show in section 5.1 that the PDF of velocity increments for time lags larger than the integral scale do not display this type of gap and are almost perfectly Gaussian. In addition, the observed intermittency properties, as shown below, are very close to those derived from DNS data.

Lagrangian velocity autocorrelation function
We consider two estimators of the autocorrelation coefficient R L (τ) of one component v of the Lagrangian velocity. The first stems from its definition, where is the average over all recorded events. σ(τ) 2 is the velocity variance of recorded tracks whose duration is longer than τ. These tracks are the only ones to contribute to the estimation of R L (τ) as obviously v(t) and v(t + τ) cannot be measured for tracks shorter than τ. Because of the finite size of the measurement volume, the length of the recorded tracks displays a non-trivial distribution and we observe the following measurement bias: the particles that remain the longest in the measurement volume are the slowest. This means that as the value of τ is increased, one selects particles that are slower and slower. The evolution of σ 2 (τ) is displayed in figure 9. The decrease in the variance with increasing width τ of the time interval is small but visible. We expect the normalization of equation (8) by σ(τ) 2 instead of the total velocity variance to help correct for the bias. . Velocity variance σ 2 (τ) of the recorded tracks whose duration is longer than τ. Experiment #3, at R λ = 810.
A second way to estimate the autocorrelation coefficient is to compute the second-order structure function D L It is related to the autocorrelation coefficient, for a stationary process, by The results of the two estimators are displayed in figure 10(a). The two curves display small but significant differences. The direct estimate decreases faster and becomes negative at τ = 0.05 s, whereas the structure function estimate remains positive for times up to integral time scale and larger. The estimate based on the second order structure function displays an almost perfect exponential decrease within the estimation fluctuations. The difference of shape comes from the measurement bias discussed above. The normalization by σ 2 cannot correct the bias in the first estimate as it becomes negative (and will remain negative when divided by a positive value). This is not the case for the second estimate as the second-order structure function takes only positive values and is indeed unbiased by the normalization; in the following, we use the second estimator in the computation of R L (τ). One can see here the importance of the size of the measurement volume compared with the integral length scale. In our case, the linear size of the measurement volume is about twice the integral length scale L, and there is still a small effect on the measurement. This effect is expected to be much more significant for smaller volumes. The negative values of the first estimator can be interpreted in the following way: particles that circle inside the measurement volume have a tendency to remain a long time in it. This creates a bias towards an anti-correlation of the Lagrangian velocities. The bias does not have a large effect on our results especially for those based on the structure functions as they can be corrected for these biases.
After correction for the bias, one can conclude that the autocorrelation function for homogeneous turbulence displays an exponential decrease for large time scales. This is in agreement with what has been suggested by previous experiments at lower Reynolds numbers [8]- [11], [30]. One defines the integral time scale T L as the characteristic time of the exponential decrease of the Lagrangian velocity autocorrelation. In the experiment at R λ = 810, one obtains  (8), (--) estimation using the second-order structure function; (· · · · · ·) exponential decrease fitted to the solid line. Inset: same plot in semi-logarithmic scale. (b) The dashed and dash-dotted lines correspond to the autocorrelation of two orthogonal components of the velocity vector. The solid line corresponds to the cross-correlation between these two components.
T L = 22.4 ms, which is of the same order of the frequency of the blades (17.4 ms: disc rotation frequency/number of blades). In the case of smooth discs (experiment #1), one obtains T L = 42.5 ms compared with the period of the disc rotation, T = 57.1 ms. It can be seen that this Lagrangian integral scale is of the same order of magnitude as the characteristic time of flow forcing (energy injection). We have also observed that when the geometry of the flow is fixed (experiments #2-4) the product of the velocity variance by the Lagrangian integral time scale defines an integral length scale that is constant (σT L ∼ 22.5 ± 0.6 mm for all three experiments). Figure 10(b) displays the autocorrelation of two orthogonal components of the velocity together with their cross-correlation. The cross-correlation remains null at all times (within the estimation uncertainty). We note that these correlations are determined less accurately than in figure 10(a) because the necessity of matching the events from the two acquisition channels decreases the amount of available data (the number of samples is decreased by a factor 3 compared with figure 10(a)). Figure 11 shows the autocorrelation coefficient for the smallest values of the time lag τ. One can observe a clear departure from the exponential behaviour (which appears linear at these time scales). The experimental curve is rounded and displays a horizontal asymptote at τ = 0. This behaviour is expected from the well-known physical argument that the curvature of the velocity autocovariance is equal to the acceleration variance. As this acceleration variance is finite for finite Reynolds numbers, the velocity autocorrelation has to display a finite curvature at the origin. This introduces another time scale [31]. To estimate it, we fit the experimental curve by the expression proposed by Sawford in his second-order stochastic model [2], For the measurement at R λ = 810, the fit gives T 2 = 0.4 ms. This value is 60 times smaller than the integral time scale T L and twice the estimated Kolmogorov time scale. In [14], a value of T 2 ≈ 0.7 τ η is reported. One interpretation is that T 2 is a characteristic time scale related to the dissipation. In our case, our interpretation is that due to its size, the particle does not respond perfectly to the flow solicitations and acts as a low-pass filter. Thus, the value of T 2 estimated in our measurements is equal to 3τ η and should be considered a time characteristic of the response of the particle. Nevertheless, our observations show than only two time scales are involved in the velocity autocorrelation: the integral scale T L related to the energy input, and a small scale linked with the energy dissipation (as T 2 is clearly in the dissipation range).
In the limit of large Reynolds numbers, in the framework of the Kolmogorov theory, one expects the second-order structure function to behave as and C 0 to be a universal constant. A first way to estimate the numerical value of C 0 is to find the maximum C 0 of the function which should tend to C 0 as the Reynolds number tends to infinity. This function is displayed in figure 12 for R λ = 810 and the values of the maximum for the various experiments are given in table 1. The different curves in figure 12 correspond to two different sizes of the particles, the Sawford model [2] (with a 0 = 6 from [14] and with C 0 = 4.1 to reproduce the observed velocity autocovariance), and the exponential fit of the velocity autocovariance function. The Sawford model is expected to be close to the real curve. Even for the smaller particle, our estimation is significantly lower than the Sawford model curve at the smallest times which leads to the estimate C 0 = 3.4, whereas the Sawford model gives C 0 = 3.8. This is due to the time response of the particle which is not fast enough to follow the flow at its smallest time scales. The maximum is seen to be even lower for the larger particle. It means that in our experiments, C 0 depends more on the particle cut-off frequency than on the Reynolds number and makes any comparison between experiments difficult. One way to circumvent this problem is to try to extrapolate to infinite Reynolds number. As the smallest time scale T 2 is related to the Kolmogorov time scale, it will vanish at infinite Reynolds number so that one expects the second-order structure function to be In that case, the C 0 parameter becomes simply This method to estimate C 0 from finite Reynolds number measurements is also found in [32]. The estimated values of C 0 are displayed in table 1. Except for the experiment with smooth disks, the value of C 0 is close to 4, which is consistent both with the estimation from atmospheric measurements [33] and with the values reported in [32]. For a given flow geometry, there is no dependence on the Reynolds number on this estimate of C 0 , within the precision of the experiment. The higher value obtained for the case of smooth disks may be due to the fact that the flow geometry is slightly different. Note that in equation (14), C 0 is defined only in terms of large-scale quantities so that it may be dependent on the structure of the flow.

Power spectrum
The estimation of the power spectrum of the Lagrangian velocity is difficult because of the structure of the signal. The raw data consist of a set of events of varying length (imposed by the residence time of the particle in the measurement volume). The common estimator of the power spectrum is based on the Fourier transform. One computes the spectrum over a window, of, say, 1024 samples, and then averages this spectrum by sliding the window over the signal. In our case, the estimation is reduced to the events that are larger than the length of the window. In particular, if one wishes to resolve the smallest frequencies, one has to increase the window size. It means that fewer and fewer events are used for the estimate. Another estimate of the power spectrum is obtained from the Fourier transform of the (biased) estimate of the autocovariance function. Even if the two estimates are supposed to have the same result as the number of samples tends to infinity, they yield different results in practice. The first estimator is well adapted at high frequencies and the second one gives better results at low frequencies (because the autocovariance function is noise-sensitive at high frequencies). Combining the two estimators in their respective domain of best confidence, one obtains the curve displayed in figure 13, which is compared with the Fourier transform of the fit given by equation (10): One observes that the agreement between the two curves is excellent. Four regions can be identified in the experimental curve. At low frequency (f < 1/2πT L ) the spectrum is flat (the Lagrangian velocity becomes uncorrelated at times larger than T L ). In the second region (1/2πT L < f < 1/2πT 2 ), the spectrum has a nearly algebraic decay with an exponent close to −2, consistent with (15) and with the Kolmogorov 1941 prediction (equation (7)). The third regime (f > 1/2πT 2 ) has a faster fall-off than the inertial range power law until the signal reaches the noise level-for f 1 kHz in figure 13. At the highest frequencies (in figure 13 for f > 1.5 kHz), one can observe oscillations that trace back to the demodulation algorithm (whose cut-off frequency is about 1.1 kHz with the settings of figure 13). An inertial range is thus developing in between the frequencies 1/T L and 1/T 2 ≈ 1/τ η . Its width is thus proportional to the turbulent Reynolds number R λ of the flow, since T L /τ η ∼ L √ ε/σ √ ν = R λ . This, again, is observed in figure 13. The extension of the inertial range is narrower than in the Eulerian case,  (15)). The dashed straight line is an guide to the eye for an ω −2 , K41 scaling. and it would require much higher Reynolds numbers to observe a clear power-law scaling in the spectrum.
These first measurements (PDF and autocorrelation coefficient) are in full agreement with the few existing theoretical predictions that come from symmetry arguments and from Kolmogorov's 1941 theory. We can estimate the Kolmogorov constant C 0 used in models of Lagrangian dispersion. From an experimental point of view, the results presented in this section show that the acoustic technique developed here yields a measurement of the Lagrangian motion of single particles over more than two decades of frequencies (from 3 Hz to 1 kHz) with a dynamic range in excess of 60 dB (when the signal reaches the noise level). This level of performance is comparable with most hot-wire measurements used in Eulerian studies.

The Lagrangian acceleration
Even if the particle is not a perfect Lagrangian tracer, its response time is close to the Kolmogorov time scale. Thus, we expect its acceleration to display most of the properties of the actual Lagrangian acceleration. We first describe the experimental observations of the particle acceleration and then the results of numerical simulations which give the actual Lagrangian acceleration, albeit at a lower Reynolds number.

Experimental measurements
4.1.1. Particle size effect. As the particle is larger than the Kolmogorov scale, the smallest scales of the inertial range will be partially filtered out by the particle response to fluid solicitations. This effect has been studied in detail in a previous paper [19] for the case of a particle settling in a fluid at rest. It was shown there that the particle response has a characteristic time where r p is the particle radius, d the density of the particle relative to that of the fluid and a p a typical value for the acceleration of the particle. Let us assume that this behaviour is valid for a small particle in a turbulent flow. On the other hand, using a Kolmogorov scaling, the characteristic time at a length scale r should be proportional to τ(r) = −1/3 r 2/3 (17) (up to a multiplicative constant of order unity). Equating these two expressions with r = r p , one can extract an expectation for the acceleration variance. The result is given in table 1 as σ a−exp together with the measured acceleration variance. The two results are in good agreement considering the crudeness of the argument. The mismatch increases with the Reynolds number or the size because the implicit hypothesis that the particle is small compared with turbulence length scales becomes less and less accurate. Nevertheless, we can conclude so far that the filtering observed at small time scales comes is due to the finite size of the particle. The particle does not react perfectly to solicitations of the flow that occur at scales smaller than its radius. The experiment at Reynolds number 1000 can be compared with the data by Voth et al [14] at the same Reynolds number and a particle size of 450 µm, which has roughly the same ratio to the Kolmogorov length scale (about 25). They report a value of a 0 lying between 1.5 and 2, whereas in our experiment #4 we get only 0.5. The difference between the two experiments cannot be explained by the uncertainty on the estimation of the dissipation rate alone. Figure 14 shows the PDF of the acceleration at R λ = 810. It is strongly non-Gaussian and has a flatness factor equal to 19. This curve must be compared with the Lagrangian acceleration measured from the motion of very small tracer particles [15]; in this work, flatness factors up to 55 have been reported. In our case, the distribution is not as wide because the particle does not respond to the smallest scales of the flow. Numerically, we observe at R λ = 810 a particle acceleration variance equal to 280 m s −2 , and very intense acceleration events can be detected, as intense as 4000 m s −2 (i.e. 14 times the acceleration variance or 400 times larger than gravity!). Example of such intense acceleration events are displayed in figure 15. We have observed that these very intense events look either like impulsions ( figure 15(a)) or oscillations (figures 15(b) and (c)). Figure 16(a) shows the autocorrelation coefficient for one component of the particle acceleration vector. The correlation decreases very rapidly and becomes negative for times larger than 9τ η . This zero crossing is expected since the integral of the autocorrelation must be zero. Thus, the correlation takes small negative values and tends to zero at infinite time lags. On the same graph, the cross-correlation of two components of acceleration is displayed. It can be seen that they are uncorrelated at all times. The zero crossing of the autocorrelation occurs for time lags of about eight Kolmogorov times. This is four times longer results than reported from other experiments or numerical simulations [12,14,16,34] due to the time response of the particle. In [12,14], the particle tracer is much smaller (about one Kolmogorov length scale) and thus displays a much shorter time response. We cannot study directly the amplitude of the acceleration vector since we have access to only two of its components, but we can make a rough estimate of its dynamics by studying the absolute value of the acceleration components. The autocorrelation and cross-correlation coefficients are displayed in figure 16(b). Their behaviour is quite different from that of the signed components. First, the decrease is much slower. The absolute value remains correlated for times up to the integral scale. Secondly, the cross-correlation is no longer zero and remains very close to the autocorrelation.

Particle acceleration autocorrelation.
The above observations are consistent with Pope's proposition [35] to write the Lagrangian acceleration as where A(t) is the amplitude and e(t) the direction vector, these two quantities being independent. His numerical observation is that e has a typical decorrelation time that scales as τ η . This is confirmed by experiments [14] that show a decrease of the (signed) acceleration autocovariance over time close to τ η . The DNS ofYeung and Pope [16] andYeung [34] shows that the acceleration amplitude remains correlated over a time that increases with R λ . Our observation is that the decrease in the correlation of the acceleration amplitude covers the entire inertial range, ending at the integral scale. The link between the amplitude of one acceleration component |a i (t)| and the acceleration modulus A(t) is not straightforward except under the hypothesis that both e i and |e i | have a decorrelation time of the order of τ η . In this case, and assuming the independence of e and A, the autocorrelation of one acceleration component maybe rewritten as Then for times larger than a few τ η one obtains Under the same set of assumptions, the cross-correlation would also be proportional to R L A in the inertial range. The hypothesis of the short time correlation of e i has been verified by Yeung and Pope [16] in a low Reynolds number DNS, but the same hypothesis for |e i | is ad hoc and can only be checked from a full three-component measurement of the Lagrangian acceleration or from DNS data. We will return to this point in the next section.
We conclude from the experimental results that the Lagrangian acceleration presents complex dynamics that involves both the dissipation time scale for the evolution of its direction and the integral time scale for the dynamics of its amplitude. This behaviour violates the K41 postulate that the acceleration characteristics should be independent of the integral time scale.

Numerical simulations
Even though our simulations are restricted to relatively low Reynolds numbers (R λ = 75 or 140), they have the advantage of giving access to full 3D quantities, in particular the acceleration amplitude. In addition, one gets the true fluid particle acceleration in this case. We will show that the features described in the previous section concerning the particle acceleration are also observed in the numerical fluid acceleration. We also check the hypothesis developed in the previous section. Figure 17 shows the autocorrelation coefficients of the (signed) acceleration components. The curves are qualitatively similar to the experimental autocorrelation functions for the particle acceleration except that they reach zero in a time of 2.05τ η at R λ = 75 and 1.95τ η at R λ = 140. These values are similar to the one reported by Yeung and Pope [16] and Yeung [34]. The components of the Lagrangian acceleration are uncorrelated.

Amplitude and direction of the Lagrangian acceleration.
We now turn to the behaviour of the acceleration modulus. Figure 18 shows the autocorrelation coefficient of the modulus of the acceleration vector, to which we compare the autocorrelation of the absolute value of one acceleration component and the cross-correlation of the absolute value of two components. For times larger than 1.5τ η , the three curves can be superimposed after multiplication by a numerical factor. This in agreement with equation (20). To check the corresponding hypothesis, we study the behaviour of the unit vector e giving the direction of the Lagrangian acceleration. In figure 19(a) it is observed that the (signed) components have a short correlation time (the zero crossing occurring at 3.7τ η ) and they are not cross-correlated, nor are they correlated with the amplitude of acceleration. When the same functions are computed for the absolute value of the components of e, one observes in figure 19(a) that the correlation time is close to 1.4τ η (the autocorrelation is always positive in that case). In addition, the absolute value of the components of e is not correlated with the amplitude a. These observations justify the hypothesis needed to obtain equation (20). To measure the autocorrelation of the magnitude of the Lagrangian acceleration in the inertial range, one only needs a 1D measurement of the acceleration. If one has a closer look at figure 18, one can see that the curves corresponding to the two different Reynolds numbers are not superimposed when the time is scaled either by T L or by τ η . One can only partially rescale the curves: for long times (τ > 2T L ) the curves seem to scale with T L so that the full damping of the correlation occurs only at times of the order of the large time scale. For very short times (τ < 2τ η ), the scaling seems to be with τ η . In the inertial range, the two curves are not found to scale either with the dissipation nor the integral time scales. The question of the dependence of this new time scale T as a function of the Reynolds number cannot be addressed here except by saying that T/T L seems to decrease with R λ while T/τ η seems to increase with R λ .

Characteristic times.
It is instructive to compare the characteristic times for the velocity, acceleration and acceleration modulus. This is done in figure 20 where the autocorrelation coefficients for these quantities are plotted in the same graph. One observes there that the acceleration modulus shares the long time correlation established for the velocity. This again emphasizes that a scale separation between velocity and acceleration is not observed, in contradiction with the Kolmogorov 41 hypothesis. The Lagrangian acceleration displays features that are definitely linked with the integral scale.

Lagrangian intermittency
The Lagrangian acceleration has a strongly non-Gaussian behaviour, as observed in our experiments and reported in [13]. On the other hand, it is also well established that the Lagrangian velocity has Gaussian statistics. Thus, one expects to observe intermittency, i.e. a progressive evolution of the Lagrangian velocity time increments from Gaussian statistics at large time scales to non-Gaussian distributions at dissipative time scales. This section is devoted to the analysis of this behaviour.

PDF of velocity increments
To show and characterize this Lagrangian intermittency, we study the evolution of the statistics of the Lagrangian velocity time increments as the time scale τ varies. The experimental PDFs for τ v are presented in figure 21(a) for τ varying from 0.75τ η to 1.75T L . At scales larger than the integral scale, the measured PDF is Gaussian. As the time scale decreases, the PDFs develop stretched exponential tails. For the smallest time scale, one recovers the shape of the acceleration PDF. The tails do not become as wide as that of the acceleration PDF [13] (dashed line) because of the particle size effect. Nevertheless, the consistency of our data with this measurement is quite obvious. Even though the DNS is done at moderate Reynolds numbers, intermittency is well observed there, as shown in figure 21(b), where the PDFs of the increments at R λ = 140 are displayed.

The flatness
A first way to quantify the observed intermittency, as the change of shape of the velocity increment PDF, is to study the evolution of the flatness factor, i.e. the normalized fourth moment: This is shown in figure 22 for the DNS data and the experimental measurement. At large scales the flatness factor is close to 3 as expected for a Gaussian distribution. The flatness then increases as the time scale decreases. The DNS data show that it increases for time scales down to 0.1τ η . For example, at R λ = 140, the flatness is only 14 at τ = τ η and saturates at a value equal to 22 at the smallest scales. This increase in the flatness for time scales below the Kolmogorov dissipative time is noteworthy. For instance, it means that the flatness factor computed from the optical measurements in [13] is probably slightly underestimated, as the effective cut-off corresponds to a scale close to 0.2τ η . Our measurements have the same behaviour, except that the increase in the flatness with decreasing scale shows an inflection at τ = 3τ η so that at τ = τ η it has almost saturated to a value close to 18. However, because of the finite time response of the tracer particles used here, this saturation value is much lower than what is expected at this Reynolds number (≈55; from [15]).

Structure functions
Intermittency is traditionally quantified by the evolution of the successive moments of the velocity increments with scale, i.e. using the structure functions: Note that the PDFs of τ v are symmetric due to the isotropy of the flow: a particle trajectory and a trajectory symmetric under any reflection are equiprobable. This is in contrast with the Eulerian case where the scale under consideration is a length scale. In this case, the von Kármán-Howarth-Monin relation states that the third moment of the Eulerian velocity increments is non-zero and related to the energy transfer rate. The experimental structure functions are displayed in figure 23(a) for values of p from 1 to 6. We do not have a sufficient amount of data to estimate higher orders. For odd values of p, the structure functions are estimated using the absolute values: because the signed moments are zero due to the symmetry of the distributions. Figure 24(a) shows the same quantities computed from the DNS data, for p up to 10. Such a high value is allowed by the small Reynolds number (R λ = 75 is displayed) and the high statistics gathered in the numerical simulation. In both the DNS and the experiment, clear power laws D L p (τ) ∝ τ ζ p are not observed in natural unit plots. In the DNS case, this could be justified by the low Reynolds number, but even for the experiment at R λ = 1000 no true scaling can be found. Note that this could be expected from the exponential behaviour of the second-order structure function, which is not rigourously compatible with a linear scaling in the inertial range. In such situations, it has now became common to study scaling properties using the hypothesis of extended self-similarity (ESS) [36], i.e. looking for power laws in the evolution of one structure function relative to another: 1.92 ± 0.14 9 1.92 ± 0. In the Eulerian case, the reference structure function is of third order, as it is assumed to scale linearly with the increment. In the Lagrangian case, the second-order structure function plays a similar role and is used as a reference. Figures 23(b) and 24(b) show the relative behaviour of the structure functions versus D L 2 . One observes a clear power-law scaling; the corresponding exponents ξ L p are given in table 4 (for the DNS data, relative scaling can be observed up to very high orders as shown in figure 24(b)). The values of measured exponents are given in table 4 and plotted in figure 25; as expected one observes a clear departure from the dimensional behaviour ξ L p = p/2. It is found that the exponents from the DNS data are equal to the ones computed from the measurements (within estimation errors), despite the very large difference in Reynolds numbers. As observed in Eulerian measurements [37], the Lagrangian relative scaling exponents do not seem to depend on Reynolds number. Another aspect of intermittency which has been the object of much debate over the years [29] is the function form of ξ p as a function of p. Given the limited data set available in our study, the matter cannot be resolved here. The experiments allow an estimate of the first six orders, but for p 6, there is no significant difference between the experimental and numerical exponents. Within experimental error, the ξ p (p) curve, for p < 6, is well fit by a quadratic law which should be considered a next order improvement over the linear dependance assumed in a Kolmogorov K41 argument. For higher orders, the above equation predicts that the exponent ξ p should decrease with increasing order. This is not consistent with the fact that the velocity should remain bounded. Indeed, it is not observed in our data, particularly in the numerics where orders up to 10 have been computed. Our observation is that the measured values of the exponents saturate for p > 7.

Multifractal random walk
As an alternative, a dynamical approach to intermittency has been proposed in [38]. In that framework, the Lagrangian velocity fluctuations are modelled by a stochastic equation of the Langevin type in which the force term has long-range correlations. Specifically, the stochastic process is a multifractal random walk (MRW) [39]; it is built such that the autocovariance of the logarithm of the velocity increments at the smallest time scale decreases logarithmically with increasing increment interval which reaches zero at a large time scale T . This ingredient is shown to lead to a log-normal multifractal process in which the correlation coefficient λ 2 1 is exactly equal to the intermittency parameters λ 2 in equation (26). The main advantage of the model is that the statistical description of intermittency in linked to specific features in the Lagrangian dynamics of fluid particles.
In the DNS, the smallest time scales are resolved (τ min < τ η ), and we can study here in more detail the functional form of the autocovariance of the logarithm of the acceleration amplitude ( figure 26(a)). The continuous lines show the actual acceleration modulus, whereas the dashed lines show the absolute value of one component of the acceleration. One observes that the acceleration displays almost one decade of logarithmic behaviour in the inertial range at R λ = 140. In the dissipation range, the decay is slower. The decay to zero is not as abrupt as in the MRW model but is a smooth transition. The coefficient of the logarithmic decay is computed to be λ 2 1 ≈ 0.14, a value that is in good agreement with the intermittency parameter computed from the low-order scaling exponents, λ 2 = 0.115 ± 0.01. For the experimental measurement, instead of the acceleration, which we cannot estimate, we use velocity increments for the smallest time scale we actually resolve. We also cannot use the amplitude of the vector so we have to use the absolute value of the velocity increment. Its behaviour is close to that of the acceleration in the case of the DNS ( figure 26(a)). The main effect of this replacement is to reduce the interval in which the logarithmic behaviour is observed. The experimental curves are plotted in figure 26(b) together with the logarithmic fits for two experiments: experiment #1 in which the discs have no blades (R λ = 510), and experiment #3 with blades (R λ = 810). In both cases, more than one decade of logarithmic decrease is observed. The decay rate λ 2 1 is close to the intermittency exponent. One important difference is the value of T . In the case of discs with blades, T is three times the integral time T L . This is different from the DNS where T is twice T L for both cases. However, in the case of the flow driven with smooth discs, T is 20 times T L which is a very long time (note that our experimental technique allows us in that case to follow a particle during such long time intervals). We thus conclude that the Lagrangian acceleration dynamics can involve time scales that can be one order of magnitude larger than the integral time scale which is characteristic of the forcing.
Such long time correlations have yet to be explained, on theoretical or phenomenological grounds. One picture comes easily to mind. If a particle is caught in a vortex, then its motion consists of a centripetal acceleration. One component of the acceleration will display oscillations with a very short period. The acceleration amplitude, in contrast, will remain high as long as the particle is trapped inside the vortex, which can be as long as lifetime of the vortex. This (not unique) picture is supported by some experimental observations: the acceleration component is indeed short-time-correlated and its amplitude is long-time-correlated. In that case, the correlation time of the acceleration magnitude reflects in part the lifetime of the vortices present in the flow. The acceleration component has a time scale of the order of the Kolmogorov time, so this suggests that the vortices would have a length scale close to the dissipative scale, which is indeed true for the case of vorticity filaments. On the other hand, their lifetime would then be distributed up to the large time scales [21,40]. Note that, in the above picture, filamentary vortices would contribute strongly to the statistics particle acceleration, but not necessarily to the statistics of the velocity itself which may be dominated by regions of weak acceleration. Indeed, the filaments are not associated with high-velocity events (in terms of velocity standard deviation).

Acceleration, dissipation and vorticity
As reported in previous work, information about the dynamics of Lagrangian motion can be gained by studying how the flow enstrophy and dissipation vary along particles paths. These quantities are extremely difficult to measure experimentally; one measurement has been reported very recently by Zeff et al [41] in a small volume fixed in space, but their variations along Lagrangian trajectories are only accessible through direct numerical simulations [16]. Using the DNS data, this section is devoted to the study of the particle acceleration conditioned on the Eulerian flow characteristics encountered during its motion: enstrophy, dissipation and vorticity.

Enstrophy and dissipation.
We first consider the joint evolutions of the acceleration modulus (A), enstrophy (ω 2 ) and dissipation ( ), and begin with an analysis of their time correlations along Lagrangian trajectories. Figure 27 displays the correlation coefficients: and so on for the other pairs; the brackets stand for time average. These curves are very similar to those published byYeung and Pope [16]. The cross-correlations (acceleration-vorticity and acceleration-enstrophy) have significant values, peaking at 0.5 for a zero time delay and decreasing only very slowly. It takes a time of the order of the Lagrangian integral time scale for the correlation coefficient to drop to half of its maximum value. These correlation functions are also observed to be asymmetric. The acceleration-vorticity cross-correlation peaks for negative values of the time lag and the decrease of the correlation is slower for negative values of τ. This means that the high-acceleration events tend to precede the high-vorticity events. This is consistent with the filament picture. If a particle is attracted to a vortex then the acceleration is higher on the edge of the vortex core, whereas the vorticity is greatest in the core. The vorticity-dissipation correlation peaks for positive values of τ and in the same way the decrease in the correlation is slower for τ > 0. This is consistent with observations by Gotoh and Rogallo [42] and with recent measurements by Zeff et al [41], who report the vorticity lagging behind the dissipation in Eulerian measurements. Again this can be explained in the vorticity filament picture, the filaments being surrounded by a region of high dissipation that the particle must cross before reaching the high vorticity core. The acceleration-dissipation curve is harder to interpret. It is maximum for τ > 0 and the correlation decreases faster for positive values of τ than for negative values. Note that long-time dynamics is observed in all these cross-correlation functions. One thus observes that the long correlations that constitute an essential feature of the Lagrangian dynamics are also found in the enstrophy and dissipation fields.
Large enstrophy or dissipation regions are probable candidates for the large changes observed in the momentum of the Lagrangian tracers, leading to the highly non-Gaussian Lagrangian acceleration statistics observed in experiments. It is thus of interest to investigate the relationship between the acceleration of the particle and the fluctuations of enstrophy or strain encountered during Lagrangian motion. One way to do this is to compute joint probabilities. For two independent variables, one should observe, for example, where P(A, ω 2 ) is the joint PDF of A and ω 2 and P(A) is the PDF of the acceleration A. An estimation of the possible dependency is given by the study of the ratio of the left to the right components of the equation above, i.e. in the above example If the two quantities are independent then this ratio should be equal to unity. Figures 28(a), (b) and (c) show this ratio for the pairs of variables (A, ω 2 ), (A, ) and (ω 2 , ), respectively, at R λ = 75 (plotted only for value of the joint probabilities larger than 10 −6 ). It is easily seen that the three quantities are not independent. For the three pairs of variables, the probability of observing simultaneous high values of the two quantities is much higher than would be the case if they were independent. The highest ratio is observed for the pair (A, ω 2 ) for which the ratio reaches values over 50. The pair (A, ) shows values as high as 30, whereas the pair enstrophy-dissipation reaches values over 40. In all three cases the ratio increases along the diagonal of the picture. In contrast, the probability ratio to observe a small value of one of the variables while the other takes a large value is much less than 1. In agreement with previous numerical simulations [16] and recent measurement at a lower Reynolds number [41], this also indicates that high values of enstrophy and dissipation happen jointly. Another way to show the non-independence of the three quantities is to compute the conditional averages of one variable given another. We call, for example, a| the average of the acceleration magnitude a conditioned on a given value of the dissipation . The different combinations are displayed in figure 28(d). As expected from the previous observations, the conditional averages are not constant. All of them are increasing functions of the conditioning variable confirming that the occurrence of large values of the three variables are highly correlated.

Vorticity filaments.
The above observations support the hypothesis of a strong influence of vorticity on the acceleration dynamics. This suggests that the acceleration of the Lagrangian particle has its strongest variations when it enters regions of high vorticity. We call these regions 'filaments' (with no attempt to define them more precisely) and only associate their vorticity modulus and the local direction of their axis with these vortical structures. Doing so, it can be shown that the vortices that have the largest influence on the Lagrangian motion are those with axes perpendicular to the trajectory of the particle. To wit, we first compute the acceleration components parallel and perpendicular to the local vorticity field along the trajectory: a ⊥ = a − a || (33) and the acceleration induced by the filaments: (Note that 1 2 ω is the local rotation rate of the fluid particle; in the simple case of a solid rotation in the core of a single vortex, this quantity is the acceleration of the particle.) An example of time variations of these quantities along a Lagrangian trajectory is shown in figure 29. It can be seen that a ω reproduces the occurrences of the large acceleration events (the variance of acceleration for the DNS represented in this picture is 0.09). When the acceleration values are comparable with the variance, a ω does not follow the same evolution as the particle acceleration, as these quantities are expected to be comparable only for high vorticity events. The acceleration perpendicular to the vorticity vector is seen to be almost equal to the full acceleration, especially for large accelerations. The acceleration parallel to vorticity displays much smaller fluctuations and is comparable with the full acceleration only when the latter is of low amplitude (of the order of the variance). The features observed in single trajectories are robust at a statistical level. This is first shown in figure 30, which shows the PDFs of the four quantities. It is seen that the PDF of a ω captures the wide tails observed in the PDF of the full acceleration component, confirming the importance of the vorticity filaments. The acceleration parallel to the vorticity vector displays a much narrower  distribution than the components of acceleration lying in the plane perpendicular to the vorticity vector, which show the exact same tails as the full acceleration. As a second test, we show, in figure 31, the time-correlation functions of (a, a ⊥ , a || ) together with the cross-covariance of a ω and a. The cross-covariance of a ω and the full acceleration peaks at half the acceleration variance showing that a ω reproduces a significant part of the acceleration. This cross-correlation is only about 50% a ω and reproduces the behaviour of the acceleration only for large acceleration events (note that since the autocovariance is a second-order moment, it is not very sensitive to the tails of the distributions).

Conclusion
We have described an experimental technique for the measurement of the Lagrangian velocity of a fluid particle. The velocity is estimated from the Doppler frequency shift induced on sound scattered by a small solid tracer advected by the flow. The use of high-performance acquisition electronics and of a high-resolution spectral analysis algorithm yields a measurement with a fast time resolution (less than a ms) and a large dynamic range of velocity fluctuations (about 60 dB). When compared with optical tracking techniques, one important advantage of this acoustic technique is that it samples a very large measurement volume in both time and space. Particles are tracked for durations that exceed several integral time scales of the flow. Together with the time resolution of the measurement it yields records of Lagrangian velocity variations at all time scales in the inertial range and up to the largest time scales of the flow. The observations deduced from the measurement have also been compared with direct numerical simulations of Lagrangian trajectories. The agreement is excellent and thus validates the experimental measurement technique with which high Reynolds number flows can be studied. In addition, the numerical simulations complement the measurements because they compute some quantities that are not accessible to the experiment, such as the vorticity and dissipation rate.
We have observed that the velocity field is Gaussian with an exponential correlation in time. Second-order quantities such as the Lagrangian velocity spectrum are observed to follow dimensional predictions à la Kolmogorov (1941). The dynamics of the Lagrangian acceleration is shown to have a mixed behaviour in terms of the time scales involved. The direction of acceleration evolves over a time scale close to the Kolmogorov time scale, i.e. the dissipative scales. In contrast, the acceleration magnitude displays correlations over a much longer time scale. In particular, the correlation of the logarithm of acceleration decays to zero in a time of the order of the integral scales. This time scale depends on the flow under consideration: for the DNS it is 2T L , for the flow entrained by discs with blades it is 3T L and for smooth discs it is much larger than the integral scale (close to 20T L ). We have also observed that the Lagrangian dynamics displays a very strong intermittency, defined as the change in the functional form of the probability density of velocity increments. Structure function exponents relative to the second-order moment are found to be independent of the flow and Reynolds number. We believe that the long time correlations observed in the Lagrangian dynamics are the origin of the strong intermittency observed in the Lagrangian velocity time series. Dynamically, we observe that the Lagrangian dynamics is strongly influenced when the particles enters regions of high local vorticity (filaments). The analysis of 1 2 ω × v definitely shows that these intense events are responsible for the wide tails of the distribution of the acceleration components. performed using IBM SP computers at the Centre National Informatique de l'Enseignement Supérieur (CINES).