Chaotic properties of a turbulent isotropic fluid

TL;DR

This paper investigates the chaotic behavior of turbulent isotropic fluids by measuring the Lyapunov exponent's dependence on Reynolds number through high-resolution numerical simulations, revealing scale-independent divergence growth.

Contribution

It establishes a quantitative relationship between the Lyapunov exponent and Reynolds number in turbulent flows using direct numerical simulations.

Findings

Lyapunov exponent scales as Re^{0.53}

Trajectory divergence grows uniformly across scales after transient

Linear divergence rate depends on energy forcing rate

Abstract

By tracking the divergence of two initially close trajectories in phase space in an Eulerian approach to forced turbulence, the relation between the maximal Lyapunov exponent $\lambda$, and the Reynolds number $Re$ is measured using direct numerical simulations, performed on up to $2048^3$ collocation points. The Lyapunov exponent is found to solely depend on the Reynolds number with $\lambda \propto Re^{0.53}$ and that after a transient period the divergence of trajectories grows at the same rate at all scales. Finally a linear divergence is seen that is dependent on the energy forcing rate. Links are made with other chaotic systems.

Equations

$$\lambda \sim \frac{1}{\tau }\sim \frac{1}{T_0}R{e}^{\alpha },\alpha =\frac{1-h}{1+h}.$$

$${\partial }_t\mathbf{u}=-\nabla P-\mathbf{u}\cdot \nabla \mathbf{u}+\nu \Delta \mathbf{u}+\mathbf{f},\nabla \cdot \mathbf{u}=0,$$

$$\bm{\hat{f}}(\bm{k},t)=\begin{array}[]{c}(\epsilon/2E_{f})\bm{u}(\bm{k},t)\ \text{if}\ 0<|\bm{k}|<k_{f};\\ 0\ \text{otherwise}\end{array}\ ,$$

$$E(k,t)=\frac{1}{2}{\int }_{|\mathbf{k}|=k}d\mathbf{k}|\hat{\mathbf{u}}(\mathbf{k},t){|}^2,$$

$$E_d(k,t)=\frac{1}{2}{\int }_{|\mathbf{k}|=k}d\mathbf{k}|{\hat{\mathbf{u}}}_1(\mathbf{k},t)-{\hat{\mathbf{u}}}_2(\mathbf{k},t){|}^2,$$

Full text

In Press Physical Review Letters 2018

††preprint: APS/123-QED

Turbulence displays chaotic dynamics Bohr et al. (2005) and ideas from chaos theory find many different applications in turbulence including the dispersion of pairs of particles Taylor (1921); Richardson (1926); Salazar and Collins (2009); Biferale et al. (2005), the presence of Lagrangian coherent structures Haller (2015), turbulent mixing Ottino (1990), turbulent transitions Eckhardt et al. (2007) and predictability Lorenz (1963); *Lorenz1969; Aurell et al. (1997); Leith (1971); Leith and Kraichnan (1972); Boffetta et al. (2002). Chaos has been seen and applied in systems as diverse as quantum entanglement, where the classical dynamical properties are linked to the quantum counterparts Jalabert and Pastawski (2001); Jacqoud and Petitjean (2009), planetary dynamics Laskar (1990), and biological systems May (1974).

Using the Eulerian approach, we track the divergence of fluid field trajectories, which initially differ by a small perturbation. We do a model independent analysis, evolving the Navier-Stokes equations for three dimensional homogeneous isotropic turbulence (HIT) using direct numerical simulation (DNS). The Eulerian approach to the study of the chaotic properties of turbulence has received only limited numerical tests prior to this Letter. Amongst approximate models, there have been EDQNM closure approximations Metais and Lesieur (1986) and shell model studies Crisanti et al. (1993a); Aurell et al. (1996); Yamada and Saiki (2007). Amongst exact DNS studies, there have been some in two dimensions Kida et al. (1990); Boffetta et al. (1997); Boffetta and Musacchio (2001) and single runs in three dimensions at comparatively small box sizes Deissler (1986); Kida and Ohkitani (1992), all more than a decade and a half ago. This Letter tests the theory of Ruelle Ruelle (1979) relating the maximal Lyapunov exponent $\lambda$ and $Re$ in DNS of HIT in a Eulerian sense. The paper also examines the time history of the divergence and finds a uniform exponential growth rate across all scales at an intermediate time and to show a linear growth for late time in three dimensional HIT. The simulations are also the largest yet for measuring the Eulerian aspects of chaos in HIT for DNS, performed on up to $2048^{3}$ collocation points and reach an integral scale Reynolds number of $6200$. This allows a more accurate measurement of the $Re$ dependence of $\lambda$.

For a chaotic system, an initially small perturbation $|\delta\bm{u}_{0}|$ should grow according to $|\delta\bm{u}(t)|\simeq|\delta\bm{u}_{0}|e^{\lambda t}$ where $t$ is time. It is theoretically predicted that the Lyapunov exponent should depend on the Reynolds number according to the rule Ruelle (1979); Crisanti et al. (1993b)

[TABLE]

The Holder exponent, $h$, is given by $|\bm{u}(\bm{x}+\bm{r})-\bm{u}(\bm{x})|\sim Vl^{h}$, where $V$ is the rms velocity, $l$ the size of the eddy, $Re=VL/\nu$ the integral scale Reynolds number, $L=(3\pi/4E)\int(E(k)/k)dk$ the integral length scale, $E$ the energy, $\nu$ the viscosity, $T_{0}=L/V$ the large eddy turnover time, $\tau=(\nu/\epsilon)^{1/2}$ the Kolmogorov time scale, and $\epsilon$ the dissipation rate. In the Kolmogorov theory, $h$ is predicted to be $1/3$ and so $\alpha$ is predicted to be $1/2$ Ruelle (1979); Crisanti et al. (1993b); Kolmogorov (1941).

Some of the new results found in this Letter from the Eulerian approach are inaccessible to the Lagrangian approach, such as the linear growth rate of the divergence at late times which has no direct Lagrangian counterpart. The paper also highlights different results from the two approaches. For instance, within the Lagrangian approach, the relation $\lambda\tau\approx const$ has been found before in tracer particles Biferale et al. (2005); Bec et al. (2006) and for infinitesimal volume deformation Girimaji and Pope (1990). Furthermore, these results suggest that $\lambda\tau$ decreases slightly with Reynolds number, and that due to intermittency corrections this implies $\alpha<0.5$ Crisanti et al. (1993a); Bec et al. (2006). As will be shown, we find that in the Eulerian approach $\lambda\tau$ increases slightly with Reynolds number, which is consistent with our result that $\alpha>0.5$. There is nothing which says the Lyapunov exponent in the Eulerian and Lagrangian frames should be the same. An example is ABC flow in which the Lyapunov exponent in the Lagrangian frame is positive but in the Eulerian frame is non-positive Dombre et al. (1986). The prediction of Ruelle for turbulence does not distinguish between Eulerian and Lagrangian frames Ruelle (1979).

We perform DNS of forced HIT on the incompressible Navier-Stokes equations using a fully de-aliased pseudo-spectral code in a periodic cube of length $2\pi$

[TABLE]

where $\bm{u}$ is the velocity field, $P$ the pressure, $\nu$ the viscosity and $\bm{f}$ the external forcing. The density was set to unity The . The primary forcing used was a negative damping scheme which only forced the low wavenumbers (large scales), $k_{f}=2.5$, according to the rule

[TABLE]

where $E_{f}$ is the energy in the forcing band and $\bm{u}(\bm{k},t)$ is the Fourier coefficient of field $\bm{u}$. This well tested forcing function Linkmann and Morozov (2015); Kaneda and Ishihara (2006) allows the dissipation rate, $\epsilon$, to be known a priori. We set $\epsilon$ to 0.1 for all runs unless otherwise stated. A full description of the code, including the forcing, can be found in Yoffe (2012). The Reynolds number quoted throughout this Letter is the integral scale Reynolds number, $Re$, which was changed by varying $\nu$. The simulations were well resolved, with $k_{max}\eta>1$ for all simulations, where $k_{max}$ is the largest wavenumber in the simulation and $\eta$ the Kolmogorov length. $T_{0}$ and $L$ vary between simulations. Over resolved simulations, with $k_{max}\eta\gg 1$, were performed to test if the box size had a statistically significant effect on the results and this was not the case. All simulations parameters are given in the Supplementary Material.

To implement the perturbation, a copy of the evolved field $\bm{u}_{1}$ was made and perturbed slightly to create field $\bm{u}_{2}$. This perturbation was achieved by not calling the forcing function at one particular timestep. This meant that the perturbation would be in the band of wavenumbers $0<|\bm{k}|<k_{f}$ and would depend non-trivially on the field itself by Eq. (3). The difference field $\delta\bm{u}=\bm{u}_{1}-\bm{u}_{2}$ was then calculated. Fields $\bm{u}_{1}$ and $\bm{u}_{2}$ were then evolved independently and the statistics of $\delta\bm{u}$ were tracked. The same realisation of the external forcing is used on both fields. The key statistic measured was the energy spectrum of the field, $E(k,t)$, which in Fourier space is defined by

[TABLE]

with total energy, $E(t)=\int_{0}^{\infty}dkE(k,t)$. Analogously, we define the energy of the difference spectrum, $E_{d}(k,t)$ as

[TABLE]

which is useful in assessing the degree of divergence of two fields at a particular scale. We then similarly define $E_{d}(t)=\int_{0}^{\infty}dkE_{d}(k,t)$ as the total energy in the difference spectrum. By inspection we can see that $|\delta\bm{u}(t)|=(2E_{d}(t))^{1/2}$.

After a statistically steady state of turbulence was reached, perturbations were made for a range of Reynolds numbers from $Re\approx 10$ to $Re\approx 6200$ at box sizes from $64^{3}$ to $2048^{3}$. We found that the growth of $|\delta\bm{u}|$ best fit an exponential $\exp(\lambda t)$. We multiply $\lambda$ by $T_{0}$ to non-dimensionalize the simulation time. A plot of $Re$ vs. $\lambda T_{0}$ is shown in Fig. 1. From the data we find a good fit to the functional form $\lambda T_{0}=CRe^{\alpha}$ with $\alpha=0.53\pm 0.03$ and constant $C=0.066\pm 0.008$, in reasonable agreement with the theory value prediction Ruelle (1979). Previous results from a shell model analysis relying on a phenomenological multifractal model to extract a fit gave $\alpha=0.459$ Crisanti et al. (1993a), whilst other Lagrangian results have suggested $\alpha<0.5$ Bec et al. (2006). We cross-checked the $Re$ dependence using an alternative DNS implementation of HIT described in Chumakov (2008), which gave a result within one standard error of ours (see Supplementary Material). In a Lagrangian study Bec et al. (2006) a decrease in $\lambda\tau$ was associated with $\alpha<0.5$. As is shown in the inset in Fig. 1 our data shows an increase in $\lambda\tau$ with $Re$, which agrees with $\alpha>0.5$ found here. This shows at least one difference between the Eulerian and Lagrangian approaches, which may have some significant underlying reason worth exploring in future work.

We find that an initial perturbation must adopt a particular spectrum, described below for $E_{d}(k)$, before $E_{d}(k)$ grows uniformly at all scales and maintains this profile during exponential growth. This particular spectrum is shown in Fig. 2 for a run with $Re\approx 2500$. The spectrum of $E_{d}(k)$ has three main characteristics; at low $k$ there is an approximately $k^{3}$ power law dependence, at intermediate $k$ $E_{d}(k)$ has a peak between the peaks of $E(k)k^{2}$ and $E(k)k^{3}$, and for high $k$ there is an exponential dependence on wavenumber, which we approximate as $E_{d}(k)\sim\exp(-Sk)$. Our DNS show that this exponential slope becomes flatter with increasing $Re$ according to a power law, this dependence is very strong and is shown in Fig. 3 which plots the relationship between $Re$ and the magnitude of the exponential slope, $S$. Thus, as $Re$ becomes large, $E_{d}(k)$ becomes flat for wavenumbers higher than the peak. The difference spectrum at low $k$ for an EDQNM approximation was found to be $k^{4}$ Metais and Lesieur (1986), whilst in a single run of DNS it was $k^{2}$ with large error Kida and Ohkitani (1992). Similar difference spectra as ours at all scales have been seen in atmospheric models Vannitsem (2017).

To understand the origin of the peak in $E_{d}(k)$, it is useful to look at the theory of Ruelle (1979), where it is assumed that the maximal Lyapunov exponent is inversely proportional to the smallest characteristic eddy time, which is the Kolmogorov time $\tau$. Naively we might expect that the peak of $E_{d}(k)$ to be $k_{\eta}$, the wavenumber corresponding to $\eta$, which is the Kolmogorov length scale with $\eta=(\nu^{3}/\epsilon)^{0.25}$. This is not observed. Instead, we can define a frequency for eddies at wavenumber $k$ of $f(k)k$ where $f(k)=\sqrt{E(k)k}$ Hinze (1975). This would make the divergence dominated by the eddies of the size of the peak of $E(k)k^{3}$, which is close to the observed peak of $E_{d}(k)$.

It is also interesting to plot the growth of $E_{d}(k)/\langle E(k)\rangle$ for selected wavenumbers as is done in Fig. 4, for the run with $Re\approx 2500$ on box size $1024^{3}$, with angled brackets representing a steady state average. The perturbation was performed at the forcing wavenumbers, $k<k_{f}$. There are three stages of growth. The first stage is a transient stage during which the characteristic $E_{d}(k)$ spectra is adopted. For the low wavenumber perturbation, the large scales remain close for at least one $T_{0}$, waiting until the small scale divergence has reached a certain size, as seen before in one dimensional atmospheric models Lorenz (2006). This is the cause for the different behaviour of $k=1,2$ in Fig. 4 compared to the other wavenumbers. In our simulations $E_{d}(k)\sim t^{2}$ for the small scales when the perturbation was made at low wavenumber. If the perturbation is made at high wavenumber, the large scales do not remain close and there is an initial convergence of the fields, as seen in 2D turbulence, suggesting a common behaviour Boffetta et al. (1997). If the perturbation is made at low wavenumber then there is no initial convergence.

Note that, although the plot in Fig. 4 is of one particular initial state and initial perturbation vector, we find that the presence of these three stages appears to be independent of the form of the perturbation made and initial state. Only the initial transient stage depends on the form of the perturbation. Perturbations made at high wavenumber exhibited the same form in the latter two stages as those made at low wavenumber. This suggests it is a characteristic feature of the difference field evolution.

The second stage is the exponential growth stage, where it is notable that all scales grow at the same exponential rate and this exponent is the same as the maximal Lyapunov exponent. In test simulations, forcing was performed at intermediate wavenumbers so that wavenumbers lower than the inertial range could be simulated. These simulations also showed the same exponential growth rate at every scale, including those larger than the forcing scale. This suggests it is not a feature of the well known forward cascade of energy in turbulence. This scale independent growth has also been seen in quasi-geostrophic turbulence in a channel McWilliams and Chow (1981), atmospheric models Vannitsem and Nicolis (1997); Vannitsem (2017), and other systems of non-linear equations Gao et al. (2006); Bohr and Christensen (1989). We now also measure it in a large turbulent simulation. In Fig. 4 this stage is relatively short but can be extended arbitrarily by having a smaller perturbation, these checks also showed our perturbation could be considered infinitesimal.

The third stage is the late time saturation stage, the details of which depend on the size of the inertial range. At late times, the growth of $E_{d}(t)$ enters a linear stage before saturation, which is entered as soon as $dE_{d}(t)/dt\approx\epsilon$. This implies that the threshold energy is $E_{d}\approx\epsilon/2\lambda$. If this energy is greater than the saturation of the difference, then the growth of the difference is exponential until it saturates. $E_{d}$ for runs at $Re\approx 130$ and $Re\approx 800$ are shown in the inset of Fig. 5, where late time starts at $t\approx 20$ for $Re\approx 130$ and $t\approx 7$ for $Re\approx 800$.

By varying the rate of dissipation we can see the dependence of this linear growth rate on $\epsilon$, which is the energy input rate for a statistically steady state system. A plot of $\epsilon$ against $dE_{d}(t)/dt$ for late times is shown in Fig. 5. The values here are not normalized and we find $dE_{d}(t)/dt=1.12\epsilon$. $dE_{d}(t)/dt$ is really a quantification of the rate of separation of trajectories in phase space, which is related to information creation, i.e. Kolmogorov-Sinai (KS) entropy. If it is possible to interpret $dE_{d}(t)/dt$ as the KS entropy, we can relate our results with corollary (2.2) of Ruelle (1982) which shows that the upper bound of the KS entropy in an isothermal fluid in equation (2.9) of Ruelle (1982) is related to the dissipation.

The findings of linear growth in $E_{d}$ at late time in a 2D DNS of turbulence were justified on the basis that there is a characteristic timescale for the eddies $\tau(k)\sim k^{-2/3}$ Boffetta and Musacchio (2001), which is in agreement with the definition of our frequency $f(k)k$. However, in our data we find instead that $\tau(k)\sim k^{-1/3}$. This linear growth at late times does not have a clear Lagrangian counterpart. For high $Re$ the exponential growth phase may be very brief and so the majority of the divergence will be dominated by the linear growth, which only depends on the dissipation. In this way the divergence of two velocity field trajectories may be universal in the Kolmogorov sense at high $Re$.

We have found that, if one scale $E_{d}(k)$ diverges exponentially, then all scales do so. This could indicate the presence of a turbulent regime. If there is no turbulent regime, then there are no scales which diverge exponentially in the Eulerian framework. This is different to the Lagrangian case. Instead of associating the inverse Lyapunov exponent with Kolmogorov time $\tau$, a slight reinterpretation of Ruelle’s theory is to associate the characteristic time with $l_{T}/V$ where $l_{T}$ is the Taylor microscale, which only exists if an inertial range exists (see Supplementary Material for data). This would also give $\alpha$ close to 0.5. This quantity uses the largest velocity and smallest length scale exclusive to turbulence to achieve the smallest time scale.

In summary, we have shown that the degree of chaos for forced HIT appears to be uniquely dependent on the large scale Reynolds number according to the law $\lambda T_{0}\sim Re^{0.53}$. Divergence does not occur at all scales until the velocity field difference spectrum adopts a characteristic form. After this spectrum is adopted, the normalized energy difference spectrum $E_{d}(k)/\langle E(k)\rangle$ grows similarly for all wavenumbers at intermediate times. Due to the shape of the spectrum, the smallest length scales will become decorrelated long before the largest length scales, as has been predicted before Lorenz (1969). At the large scales, predictability for a fixed tolerance should be possible for much longer than at the smallest scales. The late time growth of $E_{d}(t)$ was found to be linear and approximately equal to the energy input rate.

This Letter has made thorough numerical demonstrations of the links between chaos and turbulence in a Eulerian context, and so by extension relates turbulence to other chaotic processes and might provide a different perspective for their study. In chaos containing multiple length and time scales, applying ideas from turbulence may be especially fruitful because we have seen similar features here in turbulence to those found in chaotic systems which are not considered turbulent Kandrup and Sideris (2003); Vannitsem (2017); Bohr and Christensen (1989); Gao et al. (2006). There are interesting similarities between the linear growth behavior found in this paper and others Ruelle (1982); Bianchi et al. (2017), which we will examine in the future.