Multi-time structure functions and the Lagrangian scaling of turbulence

Structure functions, the statistical moments of velocity increments as a function of the spatial or time scale, are a fundamental tool to study fully developed turbulent flows and to assess how strong the intermittency is, i.e., how much the statistics deviate from the scaling predicted by Kolmogorov’s 1941 (K41) theory [1]. In the context of K41 theory, for very large Reynolds numbers and at scales (both spatial and temporal) that are small compared to the typical energy injection scales, the statistics of the flow should be self-similar and universal, and therefore independent of the forcing mechanisms. However, deviations from this prediction in the form of intermittency have been observed even at very high Reynolds numbers in homogeneous and isotropic turbulent flows, indicating that its occurrence is not a finite Reynolds effect [2]. Such deviations were also found in fractional, low order exponents of structure functions, which are less affected by extreme events and by departures from universality [3]. In order to quantify the deviations from self-similarity, the proper identification of the bounds of the inertial range in structure functions is key, and in the Eulerian frame the third-order moment of the longitudinal velocity increments has been recognized as a less ambiguous indicator of inertial-range scaling [4].

In the last decades the Lagrangian approach to turbulence, being the most natural for problems of mixing and transport, has received more attention, as both numerical and experimental methods became available to track individual particles at high rates and in large volumes [5, 6]. In spite of this, the Lagrangian characterization of turbulence and the understanding of intermittency in the Lagrangian frame remain challenging. Tracers’ velocity structure functions are known to display short ranges compatible with Kolmogorov scaling, with an amplitude and width that grows slowly with the Reynolds number [7, 8, 9, 10]. Furthermore, large-scale flows have been found to affect both small-scale Lagrangian and Eulerian turbulent statistics [11, 12, 13, 14, 15]. More recently, long Lagrangian trajectories were shown to retain information of the large-scale flow topology [16]. Evaluating how these effects impact on the observed scaling is important for the development of accurate models of Lagrangian turbulence [17, 18]. As an example, statistical theories derived directly from the Navier-Stokes equations without fitting parameters result in measurable relations for the exponents of Lagrangian velocity structure functions, which are expected to hold even in the absence of scaling [19]. It has been speculated that Lagrangian intermittency may result from the coupling of slow and fast time scales, arising from correlations between acceleration and velocity increments [20]. And it was also reported that the squared acceleration of tracers, coarse-grained over viscous time scales, recovers Gaussian statistics [21].

Many tools have been developed to improve the observed scaling and to reduce mean-flow and anisotropy effects induced by forcing mechanisms. In the Eulerian frame, extended self-similarity (ESS) [22] has been widely used to analyze scaling properties at finite Reynolds numbers. Angle averaging of structure functions has also been proposed as a way to recover isotropy in Eulerian statistics of arbitrarily forced flows [23, 24]. Indeed, performing an SO(3) decomposition to velocity correlation or structure functions has been used to disentangle anisotropic from isotropic contributions to the scaling in numerical simulations [24, 25, 26, 27]. In the Lagrangian frame, considering data conditioned to small spatial regions is also a usual practice to recover more homogeneous and isotropic statistics [13, 28], and also to reduce sweeping effects from the large-scale flow on the particles’ velocity statistics [15]. The Lagrangian acceleration spectrum has been used as a better proxy for inertial-range scaling, as it is less affected by the mean flow and finite Reynolds number effects [9]. Complementarily, the Hilbert-Huang transform methodology has been generalized and applied to single-particle velocity data to extract the hierarchy of the Lagrangian scaling exponents without resorting to structure functions or the ESS procedure [29].

Higher-degree velocity differences, for instance using wavelet analysis [30] and coarse-graining [31], have been proposed as unbiased quantities to study Eulerian and Lagrangian statistics, and as a way to detect high-order singularities in the velocity field. In the Eulerian frame, they provide a robust and straightforward technique to quantify intermittency in a signal with a steep power spectrum, such as those observed in surface wave turbulence [32]. A related tool, the so-called multi-point Eulerian structure functions, become useful as well to filter out the contribution of smooth and large-scale fields as in the case of studies of sub-ion-scale fluctuations in plasma turbulence, where the turbulent scaling is sub-leading and the power spectra are also steep [33, 34]. The importance of multi-point measurements has also been pointed out as a general way to distinguish between turbulent fluctuations and waves in space plasmas [35].

In this work we construct and characterize multi-time differences (or high-degree differences) of the Lagrangian velocity and study their moments, employing particle data from two turbulent swirling flows: Taylor-Green direct numerical simulations, and a von Kármán laboratory experiment. In both systems sweeping is observed to affect the particles’ statistics, as the tracers are affected by the large-scale counter-rotating eddies that make up the mean, macroscopic flow. At this point, it is worth noting that these flows display a mean flow (or circulation) with a complex geometry, together with large- and small-scale fluctuations. This contrasts with the paradigmatic case of homogeneous and isotropic turbulence (HIT). Distinguishing between mean circulation and large-scale fluctuations is not trivial, specially in the Lagrangian framework and for flows in which the mean flows are not stationary and average to zero in volume [36, 37]. How to disentangle these contributions in realistic flows is a long-standing problem in turbulence (see, e.g., a method proposed by Batchelor [38] and a recent application [39]). Our main motivation is then to reduce the contribution resulting from large-scale flow components, which in our flows entangle mean circulation and random sweeping (the case of HIT, without a mean flow, will be also considered), and to evaluate Lagrangian intermittency once this effect has been alleviated. When compared to their two-times counterparts, we observe that multi-time Lagrangian structure functions display a significant reduction of the large-scale flow contribution on the tracers’ statistics. This allows us to identify a time scale at which contamination from the mean flow becomes dominant, and which separates two different Lagrangian scaling ranges. The results from the multi-time analysis also suggest that Lagrangian intermittency is not a result of large-scale flow effects. The method presented here is suitable for other types of flow geometries, as a priori knowledge of forcing mechanisms or boundary conditions is not required for its application.

In this Section we introduce the multi-time technique, following the presentation for the multi-point Eulerian structure functions given in Cho [40]. The goal is to remove slow (or large-scale) variations in a signal, in order to isolate the turbulent-like fluctuations.

For the sake of consistency with our Lagrangian case (and without loss of generality) we take that signal to be a time series of the particle velocity, denoted herein by $v(t)$. This signal is considered to be composed of short-scale (i.e., fast) fluctuations $v_{S}(t)$ with correlation time $\tau_{S}$, and large scale (i.e., slow) fluctuations $v_{L}(t)$, so that $v(t)=v_{L}(t)+v_{S}(t)$. Firstly, the fast fluctuations are assumed to have turbulent-like statistics, in the sense that its power spectral density follows some power law scaling, with the largest power present at time scale $\tau_{S}$, and with less power present as the time scale decreases. Secondly, the mean of $v_{S}(t)$ is taken to be equal to $0$ when averaging over time intervals larger that $\tau_{S}$, i.e., $\langle v_{S}(t)\rangle_{\Delta t>\tau_{S}}=0$. Central differences using multiple times (or “points” in the time series) are then defined as

and so on. Given a number of times $n$, the corresponding structure function of order $p$ is then

where $\langle\cdot\rangle$ indicates an average over the time $t$, and $\Lambda^{(2\text{pt})}=1$, $\Lambda^{(3\text{pt})}=1/3$, and $\Lambda^{(4\text{pt})}=1/10$ (see Appendix A of Ref. [40] for a general expression of $\Delta v_{\tau}^{(n\text{pt})}$ and $S_{{p}}^{n\text{pt}}$). Strictly speaking, the velocity differences defined above are high-degree differences, which are always centered around the same time instant $t$. Since we are interested in studying high-order moments of said differences, and for consistency with the name used in plasma physics in the Eulerian frame of “multi-point structure” functions, we refer to the structure functions defined in Eq. (4) as multi-time structure functions.

Let’s now consider a random signal with a given and fixed variance, and with a slow trend. If for $\tau\gtrsim\tau_{S}$ the difference associated to the slow (large-scale) variations is smaller than the one associated to the fast (small-scale) variations, i.e., $|\Delta v_{\tau,L}^{(n\text{pt})}|<|\Delta v_{\tau,S}^{(n\text{pt})}|$, then the structure function using $n$ times is able to substantially remove the contributions corresponding to the slow modulation. In particular, it can be shown that a structure function using $n$ times (or points) can completely remove slow polynomial contributions of order $n-2$ [40]. Focusing on the case $p=2$, this implies that while structure functions using a small number of times will still show contamination from the slow variation, for sufficiently large $n$ and for $\tau\gtrsim\tau_{S}$ a plateau with amplitude proportional to the squared standard deviation of the small scale fluctuations $\sigma_{S}^{2}$ should become visible for our random signal. To illustrate this behavior we consider a synthetically generated random signal superimposed on a slow variation (in this example, a 1st degree polynomial), as shown in Fig. 1(a). The inset shows the power spectral density of the fast fluctuations $v_{S}(t)$ in log-log scale, which displays a power-law decay for $f\tau_{S}\geq 1$. In Fig. 1(b), the second-order structure function calculated using 2 and 3 times are shown, compensated by $2\sigma_{S}^{2}$. Note that the two-times second-order structure function, $S_{{2}}^{2\text{pt}}$, presents a plateau only at a very reduced range of time lags. In contrast, the three-times counterpart, $S_{{2}}^{3\text{pt}}$, successfully removes the effect of the slow variations, displaying a significantly broader plateau which extends beyond a decade and whose amplitude is proportional to $\sigma_{S}^{2}$.

While multi-point structure functions (or related techniques using, e.g., wavelets) were used with Eulerian data in plasma physics and wave turbulence [32, 40], with a few exceptions they have not been used to analyze Lagrangian statistics [30], and higher-order moments of the statistics were not considered. We thus consider two different datasets (albeit with similar large-scale flows in spite of differences in the boundary conditions and the forcing mechanism [37, 15]) to study whether multi-time methods can reduce the effect of the large-scale flow in Lagrangian statistics, and recover previously observed scaling laws. The two datasets stem from two sources: a von Kármán (VK) laboratory experiment, and Taylor-Green (TG) direct numerical simulations (DNSs).

The VK experiment consists of two facing disks of diameter $D=19$ cm, each fitted with $8$ straight blades, separated by a distance of $20$ cm and contained in a square cross-section cell with side of $20$ cm. As the propellers counter rotate, they agitate the fluid, which in this case is water. This generates two large counter-rotating circulation cells, which produce (on average) a strong shear layer in the mid-plane of the propellers, as well as a secondary circulation in the axial direction. In this way we obtain a three-dimensional (3D) fully developed turbulent flow in a confined region of space of dimensions $(20\times 20\times 20)~{}\text{cm}^{3}$. The flow has a mean, macroscopic structure which is anisotropic: the large scale structures in the horizontal direction (i.e., parallel to disks’ plane) are larger that the structures in the axial direction. For this study we consider two values for the rotation frequency of the disks $f_{0}^{\text{VK}}$, namely $0.83$ Hz and $1.25~{}$Hz (or equivalently, $50$ rpm and $75$ rpm, respectively). These choices were made to have Reynolds numbers in the experiment of the same order as those obtained in the DNSs described below. The flow is seeded with tracer particles, which are neutrally buoyant (density $1$ g cm^-3) polyethylene microspheres of diameter $d=250-300~{}\mu$m (Cospheric). The measurement of their dynamics is done using Particle Tracking Velocimetry (PTV), employing two high-speed cameras at sampling frequency $f_{S}^{\text{VK}}\approx 2/\tau_{\eta}$, with $\tau_{\eta}=(\nu/\varepsilon)^{1/2}$ being the Kolmogorov time scale. The 3D individual evolution of each particle is tracked in a region of size $(16\times 16\times 16)~{}\text{cm}^{3}$ around the geometrical center of the cell, whose size is comparable to the entire experimental volume. It comprises both the shear-dominated region of the flow and the areas closer to the propellers, where energy is injected. After several realizations of the experiment there are $\mathcal{O}(10^{4})$ trajectories with a mean duration of $0.34/f_{0}^{\text{VK}}$. From the individual trajectories, the instantaneous velocity is then computed after applying a Gaussian filter. For further details on the experimental setup and the measurement technique see Ref. [15] (see also [41, 12, 28] for studies on other VK experiments).

For the TG simulations, we performed DNSs of the incompressible Navier-Stokes equations, where the turbulence is sustained by an external volumetric forcing $\mathbf{F}$. The equations are solved in a 3D cubic, $2\pi$-periodic domain using a parallel pseudospectral method with the GHOST code [42, 43], with spatial resolutions of $512^{3}$, $768^{3}$, and $1024^{3}$ grid points. The mechanical forcing $\mathbf{F}$ is based on the TG flow [44], with Cartesian components

and with forcing wave number $k_{F}=1$. This forcing consists of a periodic array of counter-rotating large-scale vortices. Note that it is anisotropic: it is similar in the horizontal ($x$, $y$) direction but it injects no energy into the $z$ component of the velocity. As a result of the symmetries that the flow presents in a statistical sense, while still being anisotropic, the full domain can be split into eight cells of size $\pi^{3}$. In each of these cells, the resulting flow strongly resembles the large-scale geometry of the VK flow: two counter-rotating large-scale vortices separated by a shear layer in the midplane; these similarities with the experiment hold for both the Eulerian and Lagrangian descriptions, and have been exploited in several studies (see [37, 15] and references therein). In the turbulent steady state of this flow, $10^{6}$ Lagrangian particles are evolved in time along with the flow according to

where $\mathbf{x}_{P}$ is the position of the Lagrangian tracer at time $t$, and $\mathbf{u}(\mathbf{x}_{P},t)$ is the velocity of the fluid element at position $\mathbf{x}_{P}(t)$. The instantaneous position, velocity and acceleration of each particle are computed as it evolves. Note these particles are ideal tracers, as they do not interact with each other nor do they affect the flow’s dynamics.

In order to characterize both flows, and considering their anisotropy, we use the one-component r.m.s. tracers’ velocity $U$, associated to motions in the horizontal plane. That is, we define $U$ as

where $v_{i}$ is the tracers’ velocity component along the $i$ direction, and the brackets $\langle\cdot\rangle$ denote averages over time and trajectories (in previous studies we observed that this way of estimating the r.m.s. velocity in the TG and VK flow results in better comparisons between non-dimensional parameters in DNSs and experiments [15]; taking into account the three components of the velocity does not change significantly the estimation of $U$, resulting in values $4$ to $8\%$ smaller). We also define the integral length scale $L_{0}$ based on the scale where energy is injected into the systems. In the experiments, $L_{0}^{\text{VK}}=D=0.19$ m, the diameter of the disks. In the TG DNSs, the squared norm of the perpendicular components of the forcing wave vector (with respect to $\hat{z}$) satisfy $\mathbf{k_{F\perp}}\cdot\mathbf{k_{F\perp}}=2$, so that $L_{0}^{\text{TG}}=2\pi/|\mathbf{k_{F\perp}}|=2\pi/\sqrt{2}$ (we use $\mathbf{k_{F\perp}}$ for consistency with the experiment, where we only consider the diameter of the disks as the energy injection scale, ignoring the vertical flow dependence). In this way, an integral time can then be defined as

The energy injection rate $\varepsilon$ is directly available in the DNSs. In the VK flow we estimate it as

where $\alpha$ is a non-dimensional proportionality constant. We use the TG DNSs to estimate $\alpha$ by computing $\varepsilon/(U^{3}/L_{0})$, and observe that it varies by less than $5\%$ in the range of Reynolds numbers explored. Thus, we use $\alpha=1.69$, which corresponds to the TG DNS value using $1024^{3}$ grid points, i.e., for the simulation with the largest Reynolds number available.

We can then estimate two Reynolds numbers based on Lagrangian quantities, in order to compare experiments and DNSs on an equal footing. An integral Reynolds number based on the tracers’ velocity $U$ is defined as

where $\nu$ is the kinematic viscosity. This number should not be confused with a Reynolds number of the flow at the particle size, as we only consider tracers in this study. Also, using the definition of $\varepsilon$ in Eq. (9), the Taylor-based Reynolds number can be estimated as

Finally, besides the integral time $T_{0}$ we can also estimate a Lagrangian correlation time directly from the tracers’ dynamics. The particles’ velocity correlation function $R_{L}^{(i)}(\tau)=C_{v}(\tau)/C_{v}(0)=\langle v_{i}(t)~{}v_{i}(t+\tau)\rangle/\langle v_{i}^{2}(t)\rangle$ is computed for each of the Cartesian components of the tracers’ velocities. The one-dimensional Lagrangian correlation time $\tau_{L}^{(i)}$ is then estimated from its corresponding correlation function as the time of the first zero-crossing. Then, for each dataset, the correlation time is defined as $\tau_{L}=(\tau_{L}^{(x)}+\tau_{L}^{(y)})/2$. Note that we base this estimation on the correlation function of the horizontal components of the velocity, as the largest structures in the flow lie in the $x-y$ plane, and as a result the correlation time of the tracers is dominated by these large-scale eddies. In the experiments, the correlation function is also affected by the different lengths of the trajectories: the shorter the curves, the faster the decay observed in $R_{L}(\tau)$. In particular, in the experiments where particles can get in or out of the measurement volume, the trajectories do not have the same extension in time, and the distribution of lengths of trajectories has an exponential decay with the trajectory length. To take into account the effect of this on the estimation of the correlation time, we sub-sampled the numerical data to obtain a similar distribution of trajectories lengths. To this end we randomly picked a subset of trajectories from the $10^{6}$ trajectories available in each run, and each one was cropped from an initial time chosen randomly, up to a final time, so that the lengths of all cropped trajectories followed the same distribution of lengths observed in the experiments. From this synthetically sub-sampled data we computed the velocity correlation function and the time lag of the first zero crossing, and compared its value with $\tau_{L}$, the zero crossing of $R_{L}(\tau)$ using the complete dataset with the long trajectories. This allowed us to estimate the finite-time (or finite sampling) correction to the time $\tau_{L}$ measured in the experiments.

All the resulting parameters for the two experimental and the three numerical datasets are listed in table 1. Unless stated otherwise, the results presented in the following sections correspond to the largest Reynolds numbers available in each setup: the VK experiment at $75$ rpm, and the TG DNS with $1024^{3}$ grid points.

In this work we studied multi-time structure functions (or high-order differences) of the Lagrangian velocity, and their moments, considering particle data from two turbulent flows: a von Kármán experiment, and Taylor-Green direct numerical simulations. In previous studies it was observed that the poor agreement between theory and data in the scaling of the Lagrangian inertial range [7, 25, 9, 10] may be the result of large-scale flow effects [11, 13, 14, 20]. In a recent study, we pointed out how the effect of sweeping by the macroscopic swirling flow (common to both flows considered in the present study) might be responsible for the anomalous scaling observed in the Lagrangian second-order structure function [15]. In particular, we found a scaling proportional to $\tau^{2/3}$ for time lags $\tau\lesssim\tau_{L}$, consistent with sweeping of fluid elements by the large-scale eddies, and indicative of the fact that sweeping plays a relevant role in the particles’ evolution even in the Lagrangian frame. Motivated by these results, in this study we considered the multi-time Lagrangian structure functions as a way to quantify the effect of sweeping on Lagrangian statistics, and to disentangle its contribution from inertial-range scaling, further extending the analysis to also consider Lagrangian intermittency. In the Eulerian case, multi-point structure functions were already shown to reduce the effect of sweeping and to provide better results for the determination of steep scaling laws (see, e.g., [32, 30, 40]).

To this end, we first considered the multi-time structure functions of order $2$, which we computed for 3 and 4 times, and compared to the standard 2-times $S_{{2}}^{2\text{pt}}$. Increasing the number of times leads to broader and better-delimited scaling laws, from which we found a scaling proportional to $\tau$ for time-lags smaller or of the order of the Lagrangian correlation time. It is worth noting that this scaling, compatible with the Lagrangian K41 prediction, was not present in the (standard) 2-times second-order structure function, except for a very narrow range at significantly smaller time scales. Our results show that the range of time scales for which the multi-time structure functions display the most significant differences with $S_{{2}}^{2\text{pt}}$ corresponds to the time-scale at which the effects of sweeping by the large-scale flow become dominant over the turbulent inertial-range contributions. This allowed us to derive a phenomenological expression for the scaling of $S_{{2}}^{2\text{pt}}(\tau)$, and to identify a time scale at which the two effects become comparable.

The scaling of higher-order Lagrangian structure functions $S_{p}(\tau)$ (with $p>2$), pose a considerable challenge in experiments and simulations due to their heightened sensitivity to finite-Reynolds-number effects and finite width of their inertial range. Multi-time structure functions of order $p$ ranging from 1 to 5, using $2,3$, and $4$ times, showed an increasing width of the range of time scales compatible with power-law scaling, and a reduced contamination from the mean flow. The values obtained from the scaling exponents are, within uncertainties, in good agreement with those reported in the literature for different flows [41, 49, 50, 51], which are often obtained using extended self similarity as a way to compensate for the narrow scaling range and the poor statistics. The agreement of the scaling exponents independently of the number of times considered indicates that the large-scale flow is not the direct cause behind departures from self-similarity in Lagrangian scaling. Intermittency is an inherent property of Lagrangian turbulence, which may still result from some other subtle coupling between slow and fast time scales, or may be associated, e.g., with the time irreversibility reported in tracers’ trajectories in turbulent flows [52]. Besides this, the multi-time structure functions provide for a direct and flow-agnostic method to improve scaling without the need of extended self similarity, and results in an effective broadening of the inertial range, leading to a more precise determination of the associated exponents.

A study of the behavior of multi-time statistics on the Reynolds number, using all our datasets spanning microscale Reynolds number in the range $\textrm{Re}_{\lambda}\in[215,530]$, confirmed that the sweeping range is independent of Re, a result in agreement with a phenomenological argument and with the fact that the intensity of the large-scale flow is controlled mostly by the boundary conditions and the energy injection mechanism (even though the mean flow geometry is expected to depend weakly on Re). For structure functions using 3 and 4 times, an increase in the Reynolds number results in a widening of the inertial range, as expected. For all values of Re considered, the multi-time structure functions also considerably decrease the contamination of the inertial scaling.

Can we expect scaling to improve indefinitely as more times are used in computation of the structure functions? Our results indicate that if a given structure function displays inertial-range scaling, increasing the number of times $n$ improves the scaling. But although it may seem advantageous to then use a larger number of times, in practice datasets have finite length in time, both in simulations and in experiments. As $n$ is gradually increased, the range of available time lags shrinks drastically. This can also result in noisier structure functions and scaling exponents due to the lack of statistical convergence associated to the reduced number of data employed in their determination. As a consequence, there is an optimal value for the number of times that should be considered, which depends on the amount of data available.

As a tool, multi-time Lagrangian structure functions succeed in reducing the sweeping effect of the large-scale flow on the Lagrangian statistics of tracers. The technique is also suitable for other types of flow geometries, as a priori knowledge of the forcing mechanisms or boundary conditions is not required for its application. These structure functions also have a relation with the energy spectrum, providing an interesting link with other tools used in the characterization of turbulence.