Self-organization and transition to turbulence in isotropic fluid motion driven by negative damping at low wavenumbers

We observe a symmetry-breaking transition from a turbulent to a self-organized state in direct numerical simulation of the Navier-Stokes equation at very low Reynolds number. In this self-organized state the kinetic energy is contained only in modes at the lowest resolved wavenumber, the skewness vanishes, and visualization of the flows shows a lack of small-scale structure, with the vorticity and velocity vectors becoming aligned (a Beltrami flow).

In the course of studying the direct numerical simulation (DNS) of forced isotropic turbulence, we have found that self-organized states form at low Taylor-Reynolds numbers, in simulations which were invariably turbulent at large Taylor-Reynolds numbers, thus hinting at the possibility of a transition to turbulence. We observed depression of the nonlinear term in the Navier-Stokes equation (NSE) during the formation of the self-organized state. Visualization of the flow showed that in the self-organized state the velocity field u and vorticity field u  ω = × are aligned in real space, forcing the nonlinear term u ω × in the NSE to vanish. This is known as a Beltrami field [1,2]. As it is also a condition for maximum helicity, and the initial state has zero helicity, it follows that the transition to the self-ordered state is symmetry-breaking. Moreover, the flow shows only large-scale structure.
We begin by discussing the details of our 'numerical experiment'. The incompressible forced NSE was solved numerically, using the standard fully de-aliased pseudospectral method on a 3D periodic domain of length L 2 box π = , resulting in the lowest resolved wavenumber k L 2 1 min box π = = . We note that our simulation method is both conventional and widely used, so we give only some necessary details here. Full details of our numerical technique and code validation may be found in [3].
We found that a choice of 32 3 collocation points proved to be sufficient in order to resolve the Kolmogorov dissipation scale, as all our simulations in this particular investigation were at very low Reynolds numbers. We did however verify that the formation of the self-organized state was also observed at a higher resolution (64 3 ). All simulations were well resolved, satisfying k 2.16 max η ⩾ [4], where η denotes the Kolmogorov dissipation scale. The maximum time the simulations were evolved for was t 1000 s = , or, in terms of initial largeeddy turnover times t , where U denotes the initial rms velocity and L the initial integral length scale. Initial Taylor-Reynolds numbers R λ ranged from R 2.61 = λ to R 3.78 = λ , where λ denotes the Taylor microscale. A summary of simulation details is given in table 1.
The initial conditions were random (Gaussian) velocity fields with prescribed energy spectra of the form is the Fourier transform of the vorticity field x ( ) ω , was negligible for all simulations. The system was forced by negative damping, with the Fourier transform of the force f given by is the instantaneous velocity field (in wavenumber space). The highest forced wavenumber, k f , was chosen to be k f = 2.5. As E f was the total energy contained in the forcing band, this ensured that the energy injection rate was constant W ε = . Negative damping was introduced to turbulence theory by Herring [5] in 1965 and to DNS by Machiels [6] in 1997. Since then it has been used in many different investigations to simulate homogeneous isotropic turbulence (for example, [4,[7][8][9][10][11][12]) and particularly in the benchmark simulations of Kaneda and co-workers [13], with Taylor-Reynolds numbers up to R 1200 = λ . It has also been studied theoretically by Doering and Petrov [14]. It should be emphasized that the forcing in the DNS does not relate the force given in (2) to the velocity field in the NSE at the same instant in time. The force added at time t is calculated from the velocity field at the previous time step t t d − , where in practice an intermediate predictor-corrector step is taken. This has the important consequence that, because of nonlinear phase mixing during this additional step, the force and the velocity field in the NSE at a given instant in time are not in phase, as long as the nonlinear term is active. Although this forcing procedure does not have an explicit stochastic element, it is nevertheless generally regarded as adding a random force to the fluid. As the velocity field at t = 0 consists of a (pseudo-)randomly generated set of numbers, the forcing rescales this set of numbers and feeds it back into the system at the first time step; and so on.
In order to investigate aspects of turbulence at high Reynolds number, which are beyond currently available computing power, our invariable practice is to keep W ε constant, and reduce ν. This corresponds to taking the limit of infinite Reynolds number [15]. We have carried out investigations up to R 435 = λ with regards to the behaviour of the dimensionless dissipation rate [12] and the behavior of the second-order structure function [11]. During these investigations, we observed the peculiar asymptotic behavior of the fluid at low Reynolds numbers as shown in figure 1, which was confirmed using the publicly available code hit3d [16,17].
At such low Reynolds numbers, we found that initially the simulations developed in the usual way, with clear evidence of nonlinear mixing. But, after long running times, the total energy ceased to fluctuate, and instead remained constant, while the skewness dropped to zero. Once this had happened, the kinetic energy was confined to the k = 1 mode 4 and the nonlinear transfer had become zero. The cascade process was thus absent and no small-scale structures were being formed. Hence the system had self-organized into a large-scale state.
Similar results were found both with our code and with hit3d as shown in figure 1(a). Figure 1(b) shows the evolution of the velocity derivative skewness. Note that it drops to zero at long times, indicating that the observed state is Gaussian and hence is not turbulent. As the velocity derivative skewness is not recorded by hit3d, a comparison using this quantity was not possible.
Since the energy spectrum in the self-organized state is unimodal at k = 1, we can predict the asymptotic value E t E ( ) = ∞ in this state from the energy input rate W ε and the viscosity ν using the spectral energy balance equation for forced isotropic turbulence where E k t ( , ) and T k t ( , ) are the energy and transfer spectra, respectively, and is the work spectrum of the stirring force. By invoking stationarity ( W ε ε = , where ε denotes the dissipation rate), we obtain for the total energy after self-organization All runs for all values of the initial Taylor-Reynolds number stabilized at this value of E(t).
Note that the same value is obtained both with our code and with hit3d as shown in figure 1(a). The realizations obtained from hit3d and from our code only differ in the random initial condition while all other parameters such as the viscosity and forcing rate are identical. As can also be seen in figure 1(a), although both realizations attain the same asymptotic value of E(t), they do not do so at the same rate. The development of the large-scale state was accompanied by a depression of nonlinearity, as shown in figure 2(a) where the transfer spectrum T k t ( , ) is plotted at different 4 That is, for the six wavevectors ( 1, 0, 0), (0, 1, 0) ± ± and (0, 0, 1) ± .
times. Again, it should be emphasized that at early times T k t ( , ) has the behaviour typically found in simulations of isotropic turbulence (see e.g. figure 9.3, p 293 of [18]), but at later times tends to zero. At the same time, the energy spectrum E k t ( , ), tends to a unimodal spectrum at k = 1 as can be seen in figure 2(b). Figure 3 shows snapshots of the system before and after self-organization. The arrows indicate the velocity field u x ( ) and vorticity field x ( ) ω at different points in real space. We , where U is the initial rms velocity and L the initial integral scale. the Taylor-Reynolds number in the self-organized state, ν the kinematic viscosity, R λ the initial Taylor-Reynolds number and t max the time the simulations were evolved for.   for a coefficient α (which must have dimensions of inverse length). Vector fields satisfying (6) are eigenfunctions of the curl operator and therefore helical. It should be noted that the nonlinear term in the NSE vanishes identically if the velocity and vorticity fields are aligned. From the image on the right hand side of figure 3, it can clearly be seen that the flow is in an ordered state, as velocity and vorticity vectors are everywhere parallel to each other, giving visual evidence of the Beltrami property (6) of the large-scale field. Furthermore, we observe a lack of small-scale structure of the flow, which reflects the measured unimodal energy spectrum in the self-organized state. The Beltrami condition has also been verified by analyzing helicity spectra, showing nonzero helicity at k = 1 and zero helicity at all higher wavenumbers. We found that the final states always satisfied the Beltrami condition. However, final states occur where vorticity and velocity are aligned (positive helicity) or antialigned (negative helicity). Which of the two possible helicity states is chosen in the selforganized state depends on the random initial conditions. If a particular final state occurs, then its mirror image will also be a possible asymptotic state that may occur, as there is no systematic preference for a particular configuration.
A short film showing the transition to the self-organized state can be found in the supplemental material (available online at stacks.iop.org/jpa/48/25FT01/mmedia). Two interesting points may be observed from the film. First, the system shows behaviour similar to 'critical slowing down' (see e.g. [19]). Secondly, we observe a transition from a nonhelical state to a helical state. That is, the formation of the helical self-organized state is symmetrybreaking. This latter aspect can also be seen in figure 3.
Our results may be peculiar to our particular choice of forcing. We note that as the nonlinear term tends to zero, the phase-mixing also vanishes, and the negative damping becomes in phase with the velocity field. Accordingly there is then no mechanism to restore randomness, and in this simple picture, the self-organized state must be stable. A further implication of this, is that negative damping would be rather difficult to achieve in a real physical system. Accordingly we are currently exploring the effect of other types of forcing and this work will be reported in due course.
However, it is arguable that the combination of the NSE and negative damping is an interesting dynamical system in its own right. For initial values of the Taylor-Reynolds number greater than about 25 (corresponding to a steady-state value of about 40), it reproduces isotropic turbulence; while for initial R λ less than about 5, it makes a spontaneous, symmetry-breaking transition to a Beltrami state. In the past, there have been mathematical conjectures about the possibility of Beltrami fields [20,21] (and references therein) in turbulence; and numerical studies have suggested that such fields can occur in a localized way [22][23][24]. To the best of our knowledge, this present investigation reports the first instance of a spontaneous transition to a stable Beltrami state for an entire flow field.
Furthermore, the emphasis in DNS has always been on achieving ever-larger Reynolds numbers, in order to study asymptotic scaling behaviour. Yet in the case of the dimensionless dissipation, for example, the emphasis is on low Reynolds numbers, and indeed on looking for universal behaviour under these circumstances. Accordingly, the present work also draws attention to the surprising fact that stirred fluid motion at low R λ may not actually be turbulent.
One of the few practical applications of isotropic turbulence is in geophysical flows. Our results could possibly be relevant to atmospheric phenomena such as 'blocking' [25], for instance. Recent numerical studies in magnetohydrodynamics show the development of similar organized states [26]. This suggests a wider relevance of this type of result, not least because it was achieved in a different system using a different forcing scheme from our own.
Investigation into the nature of this transition continues. In view of recent successes of a dynamical systems approach to the transition to turbulence in wall-bounded shear flows, work is currently being carried out taking a similar approach, and this will be submitted for publication in due course.