Nonlinear amplification in hydrodynamic turbulence

Abstract Using direct numerical simulations performed on periodic cubes of various sizes, the largest being $8192^3$, we examine the nonlinear advection term in the Navier–Stokes equations generating fully developed turbulence. We find significant dissipation even in flow regions where nonlinearity is locally absent. With increasing Reynolds number, the Navier–Stokes dynamics amplifies the nonlinearity in a global sense. This nonlinear amplification with increasing Reynolds number renders the vortex stretching mechanism more intermittent, with the global suppression of nonlinearity, reported previously, restricted to low Reynolds numbers. In regions where vortex stretching is absent, the angle and the ratio between the convective vorticity and solenoidal advection in three-dimensional isotropic turbulence are statistically similar to those in the two-dimensional case, despite the fundamental differences between them.

dissipation (Tennekes & Lumley 1972;Chorin 1994). It has been widely reported that local depletion of nonlinearity occurs in developed turbulence (Pelz et al. 1985;Kit et al. 1987; Kraichnan & Panda 1988;Shtilman & Polifke 1989;Shtilman 1992;Tsinober, Ortenberg & Shtilman 1999;Bos & Rubinstein 2013). Based on this notion, it has been conjectured that turbulent flows spend a large proportion of time in neighbourhoods of Euler solutions, which are without dissipation and are weakly Beltrami in nature, i.e. velocity u and vorticity ω (≡ ∇ × u) are aligned with each other (Levich 1987;Moffatt & Tsinober 1992;Moffatt 2014). Such neighbourhoods are thought to be hospitable for flow structures to persist coherently in time (Hussain 1987). It has been demonstrated that there exists no direct correlation between regions of low dissipation and velocity-vorticity alignment (Kerr 1987;Rogers & Moin 1987;Speziale 1989;Wallace, Balint & Ong 1992), but it is unclear if the local 'Eulerization', i.e. the suppression of nonlinearity, occurs in the limit of vanishing viscosity. More importantly, how does the vortex stretching mechanism, which depends on the nonlinear character of the Navier-Stokes dynamics, change with decreasing viscosity or increasing Reynolds numbers? An exposition of these issues is the thrust of this paper.
Using a direct numerical simulation (DNS) database of isotropic turbulence that spans more than two orders of magnitude in box size, we particularly study regions with depleted nonlinearity for their energy cascade characteristics. The solenoidal projection of the nonlinear advection in the Navier-Stokes dynamics is used to quantify the strength of nonlinearity as a function of Reynolds number. Using the geometric alignment between the nonlinear advection and the vortex stretching, we explore how the fundamental energy transfer evolves with increasing Reynolds numbers. Our main findings are as follows: (i) we confirm previous low-Reynolds-number results that dissipation fluctuations are not influenced by local 'Beltramization' of the flow. Such regions have low nonlinearity and yet remain highly dissipative. (ii) We show that the depletion of nonlinearity is not an endemic feature of hydrodynamic turbulence in the vanishing viscosity limit, and that turbulence deviates further away from the Euler flow in this regard. (iii) Nonlinearity indeed grows with increasing Reynolds number, which causes the vortex stretching mechanism to become highly intermittent. (iv) Regions in three-dimensional (3-D) isotropic turbulence where vortex stretching is locally absent tend to resemble two-dimensional (2-D) flows in certain statistical aspects.
The outline of the rest of the paper is as follows. In § 2 we briefly describe the numerical method and provide some basic definitions. The results and discussion are contained in § 3. We summarize the main findings in § 4.

Numerical simulations and definitions
We consider a constant-density turbulent flow of a viscous fluid in a tri-periodic cubic domain with the velocity fluctuation u governed by the Navier-Stokes equations, ∂u/∂t + u · ∇u = −∇p/ρ + ν∇ 2 u + f , (2.1) and the continuity equation, ∇ · u = 0. Here, p is the pressure, ρ is the density and ν is the kinematic viscosity of the fluid. The forcing acceleration f , based on the Ornstein-Uhlenbeck random process, is used to obtain a statistically steady, homogeneous and isotropic state (Eswaran & Pope 1988). Data from DNS spanning a wide viscosity range have been used for this analysis (see table 1).  Table 1. Salient DNS parameters: the grid resolution N on a side of the box; the Taylor-microscale Reynolds number R λ = u λ/ν, where u is the root-mean-square velocity, the microscale λ is given by λ 2 = 15νu 2 / and = (ν/2) i,j (∂u i /∂x j + ∂u j /∂x i ) 2 is the energy dissipation rate; the number of independent snapshots N r used for data averaging; and the spatial resolution parameter k max η = √ 2Nη/3, where η = (ν 3 / ) 1/4 is the Kolmogorov length scale. The box length L 0 = 2π, the integral scale L ≈ L 0 /5 and the large-eddy time scale T E = L/u = 0.8. Also given are the second moments of / and relative helicity density (see (2.4)). Error bars indicate 95 % confidence intervals based on Student's t distribution. The error bars for cos 2 θ are insignificant and hence not reported. method, which by now is standard, and the DNS database, see Donzis & Yeung (2010) and Yeung, Zhai & Sreenivasan (2015). We point out that the large range of box sizes allows us to span a large enough range in the Reynolds number, so as to be able to ascertain Reynolds-number variations; the highest microscale Reynolds number of 1300 is high enough to indicate conditions of the asymptotic state.
We now provide a few definitions before proceeding further. We define the Lamb vector as (Wu, Zhou & Fan 1999) Λ ≡ ω × u (2.2) and the helicity density as h ≡ ω · u. (2. 3) The Lamb vector and helicity density are connected by |Λ| 2 + h 2 = |u| 2 |ω| 2 . We define the relative helicity as the cosine of the angle θ between the fluctuating velocity and vorticity vectors as As for any other vector, one can decompose the Lamb vector into irrotational and solenoidal parts (Helmholtz decomposition) and write where α is a scalar field and β is a solenoidal vector field (∇ · β = 0), with no loss of generality. Recasting the advection term in (2.1) as u · ∇u = Λ + ∇(|u| 2 /2), we note that ∇ × β is the solenoidal projection of the advection term in (2.1) (Moffatt 1990;Shtilman 1992;Tsinober, Yeung & Vedula 2001). By taking the curl of the Lamb vector, using (2.5) and utilizing the fact that β is solenoidal (Speziale 1989), we get where ω · ∇u is the vorticity stretching term that gives rise to the energy cascade from the large scales to the small, and u · ∇ω is the convective vorticity term. Equation (2.6) is also illustrative of how the nonlinear advection term effects the vortex stretching mechanism through its solenoidal component.

Results and discussion
3.1. Relative helicity density Table 1 lists the second moment of the dissipation and the mean square of the relative helicity density, see (2.4), for R λ = 38 to 1300. With increasing Reynolds number, R λ , the dissipation fluctuations become more intense, thereby resulting in the larger second moment. As R λ increases, the scale separation -between the integral scale L and the Kolmogorov scale η -increases, resulting in a decrease in the mean-square relative helicity, cos 2 θ . As the isotropic velocity and vorticity vectors become more uncorrelated with increasing scale separation, the probability density function (p.d.f.) of cos θ becomes essentially flat at higher R λ (all values of cos θ are equally likely) as seen in figure 1, and so cos 2 θ tends towards 1/3. This geometric decorrelation between vorticity and velocity does not preclude local regions of preferential alignment, corresponding to high helicity. It has been demonstrated previously that weakly dissipative regions are not necessarily highly helical (Pelz et al. 1985;Kerr 1987;Wallace et al. 1992), since a flow can be locally Euler yet non-Beltrami (Tsinober 1990b). However, the converse scenario of high helicity implying low dissipation has been conjectured to be true (Moffatt 1985(Moffatt , 1986Moffatt & Tsinober 1992). The rationale is that highly helical regions are characterized by strong velocity-vorticity alignment (|cos θ | ≈ 1) which, in turn, implies small Lamb vector magnitudes (|Λ| ≈ 0) and hence a depletion of nonlinearity, manifesting in decreased energy transfer and dissipation of energy.
In order to examine whether high helicity implies low dissipation, we condition, at R λ = 1300, the first two moments of dissipation on the relative helicity density, and show the result in figure 2. The conditional moments of dissipation are close to the respective unconditional means, with the maximum difference (<5 %) occurring at cos θ ≈ ±1. Other small-scale quantities such as the enstrophy density (vorticity-squared), although not displayed here, also show no dependence on the velocity-vorticity alignment, at the various R λ examined (see table 1). This result confirms the theoretical arguments of Speziale (1989) that strongly helical regions where |cos θ| ≈ 1 can be associated with small Lamb vector magnitudes |Λ| but with a large solenoidal projection of the advection term, i.e.large |∇ × β| (see (2.5)), thereby leading to increased downscale energy transfer and increased dissipation due to (2.6). The conclusion is that the solenoidal projection of the advection term (∇ × β) is important to the energy transfer mechanism, rather than the velocity-vorticity alignment (Speziale 1989).

Nonlinear amplification
Taking the scalar product of Λ with itself, and noting that (∇α) · (∇ × β) = 0, (2.5) yields |Λ| 2 = |∇α| 2 + |∇ × β| 2 , (3.1) assuming homogeneity. Now consider the ratio of the solenoidal to irrotational parts of (3.1) (Shtilman 1992): The ratio R is a measure of the nonlinear reduction, as it quantifies the relative departure of the Lamb vector from its potential part. We note that, in a steady-state Euler flow, Λ = ∇(−p − 1 2 u 2 ) and the solenoidal part is zero, thereby giving R = 0. A Gaussian solenoidal vector field also has a substantial potential part, |∇α| 2 / |Λ| 2 ≈ 0.7 or, equivalently, R ≈ 0.43 (Tsinober 1990a). Figure 3 shows R for various Reynolds numbers. At R λ = 38, R = 0.27 ± 0.01, which is consistent with results in Shtilman (1992). This means that, at that low Reynolds number, almost three-quarters of the mean square |Λ| 2 gets its contributions from the potential part, with the solenoidal part making up for only a quarter. Because of this imbalance, it was thought that the former tends to depress its nonlinear strength (Kraichnan & Panda 1988;Shtilman 1992). However, as R λ increases, we find that R increases, although only logarithmically, as shown in the inset of figure 3. (Such logarithmic changes are common for all intermittent processes in turbulence (see Sreenivasan & Bershadskii 2006).) At R λ = 1300 the potential and solenoidal parts contribute almost equally. This suggests that the Navier-Stokes dynamics strengthens the solenoidal part of the advection term at higher R λ , in comparison to its behaviour at low Reynolds numbers, where the tendency is to suppress nonlinear advection (Kraichnan & Panda 1988).

Vortex stretching mechanism
The quantity ω · ∇u represents the enhancement of vorticity by stretching and tilting of vortex lines and is responsible for the forward cascade dynamics in 3-D turbulence   (Tsinober 1990a). The inset shows that the DNS data collapse onto the curve R = 0.18 ln R λ − 0.38, obtained using a least-squares fit. Error bars, most often subsumed by the symbol thickness, indicate 95 % confidence intervals.   against the logarithm of the normalized wavenumber kη. Symbols ( , blue), (•, red), ( , maroon) and ( , purple) correspond to R λ = 140, 400, 650 and 1300, respectively. The dashed lines show exponent ϕ(k) ∼ k υ for k 1/η, using a least-squares fit. The inset shows the wavenumber k at which ϕ(k) is maximum, as a function of R λ . The drop-off at the highest wavenumber k max is due to finite resolution effects. (b) Scaling exponent υ of the vortex stretching spectrum, as a function of R λ . The arrow at υ = 5/3 shows the self-similar estimate. Error bars indicate 95 % confidence intervals. (Chorin 1994). Figure 4(a) shows the vortex stretching spectrum ϕ(k), which is related to its variance as ϕ(k) dk = |ω · ∇u| 2 . The spectrum exhibits a power-law behaviour ϕ(k) ∼ k υ up to kη ≈ 0.1, followed by a peak in the dissipation range. As seen in the inset of figure 4(a), the peak wavenumber shifts further into the dissipation range with increasing Reynolds number, indicating that larger and larger wavenumbers are excited within the vortex stretching mechanism as the viscosity decreases. The scaling exponent υ, for kη 1, obtained from a least-squares fit, is shown to gradually decrease as a function of R λ in figure 4(b). A consequence of the self-similar theory of Kolmogorov (1941), which assumes that the local velocity over a distance r scales as r 1/3 , is that the vorticity stretching spectrum scales as ϕ(k) ∼ k 5/3 (Ishihara et al. 2003). The DNS exponent at R λ = 1300 differs from the self-similar estimate of υ = 5/3 by about 8 %, which is presumably related to intermittency effects in the vorticity stretching spectrum.
The anomaly displayed by the vortex stretching spectrum is reflected in the deviation from Gaussianity of the p.d.f. of the vortex stretching magnitude. Figure 5 shows the p.d.f. of the vortex stretching magnitude at different R λ . The p.d.f.s exhibit increasingly wider tails with increasing R λ , reflecting the fact that intense fluctuations of ω · ∇u occur with greater probability at higher R λ . In contrast, the likelihood of fluctuations around the mean decreases with Reynolds number, as can be seen in the inset of figure 5. The inference is that vortex stretching becomes increasingly intermittent, with low to moderate vortex stretching events interspersed by more and more extreme events, occurring with greater likelihood with increasing R λ , similar to other small-scale quantities such as the energy dissipation (Sreenivasan & Meneveau 1988;Sreenivasan & Antonia 1997).
In 3-D turbulence, the three terms in (2.6) form a degenerate triangle with 0 ≤ γ ≤ π and Ψ ≥ 0. Figure 6(a) shows that the tendency for the vectors u · ∇ω and ∇ 2 β to be  Figure 7. The conditional probability density function of the ratio Ψ of the magnitudes of ∇ 2 β and u · ∇ω ((2.6) and (3.4)), conditioned on γ = π at R λ = 140 ( , blue), 650 ( , maroon) and 1300 ( , purple). The inset shows the expanded view around Ψ = 1, which corresponds to 2-D flows, to reveal the increasing R λ trend for the peak values of P.
antiparallel increases with R λ . When they are antiparallel, the vortex stretching magnitude fluctuations attain a minimum, as seen in figure 6(b). In order to examine the ratio of the magnitudes of ∇ 2 β and u · ∇ω when they are antiparallel, we evaluate in figure 7 the conditional p.d.f. of Ψ , conditioned on the event γ = π. The ratio peaks at unity, with the peak increasing with R λ , as shown in the inset. As expected, this conditional ratio approaches unity with increasing R λ , indicating that the conditional p.d.f. tends essentially to a large peak at unity for very large R λ . The likelihood that γ = π (corresponding to Ψ = 1, which occurs when γ = π) increases as the logarithm of the Reynolds number (see figure 8 corresponding to γ = π from figure 6(a) as a function of the Reynolds number). The likelihood that γ = π is seen to follow a logarithmic trend in R λ (Sreenivasan & Bershadskii 2006). The conditional p.d.f. peaks of figure 7 also follow a similar R λ trend (although not shown here). This suggests that, in the asymptotic R λ → ∞ limit, γ = π and Ψ = 1 in vast regions of 3-D isotropic turbulence, resembling 2-D turbulence (Boffetta & Ecke 2012).

Conclusions
Tsinober (1990a) has shown that there is a depletion of nonlinearity purely from kinematic constraints. Based on the computations of Kraichnan & Panda (1988) and Shtilman (1992) at R λ ∼ 30, it was thought that the Navier-Stokes dynamics has a tendency to suppress this suppression further. It was argued that a turbulent flow was more likely to be in regions of phase space that contain fixed-point Euler solutions. High-helicity regions with favourable velocity-vorticity alignment were put forth as a plausible reason for the nonlinear suppression (Levich 1987;Levich & Shtilman 1988). Since then, results at higher R λ , such as Speziale (1989) and Wallace et al. (1992), have found no significant anticorrelation between helicity and dissipation in isotropic turbulence. Since the Navier-Stokes dynamics depends on Reynolds number, it seemed reasonable to examine the nonlinear depression at higher Reynolds numbers, covering a range of Reynolds numbers.
Using a DNS database with a Taylor-scale Reynolds number between 38 and 1300, we calculate the ratio of moments of the potential and irrotational parts of the Lamb vector. This ratio is shown to increase as the logarithm of the Reynolds number. For R λ = 1300, the moments of the potential part are almost equal to those of the irrotational part. This result indicates that nonlinear suppression is not a characteristic feature of high-Reynolds-number Navier-Stokes turbulence as previously thought, and that nonlinear amplification continues to be strong at higher R λ . A plausible implication of this strong presence is that turbulent flows might spend less time around fixed-point Euler solutions than previously perceived.
The growth of the solenoidal part of the advection term has a major impact on the vortex stretching mechanism. The Laplacian of β, where ∇ × β is the solenoidal part of the advection term in the Navier-Stokes equations, is antiparallel and nearly equal in magnitude to the vortex convection u · ∇ω, almost everywhere. A consequence of this geometric alignment is that some vortex stretching manifests itself nominally for most of the time and that intense stretching events occur only intermittently. A related observation here is that, in vast regions of a fully isotropic 3-D turbulent flow, the angle and the ratio between the convective vorticity and solenoidal advection vectors is the same as that in 2-D flows, and that this resemblance becomes stronger with increasing Reynolds number. A possible implication is that more profound connections (and differences) exist between the statistical properties of 3-D and 2-D turbulence than what is already appreciated (Biferale, Buzzicotti & Linkmann 2017;Goldenfeld & Shih 2017;Xia & Francois 2017). Quite apart from the presence of intermittent vortex stretching, which controls the energy cascade dynamics in three dimensions, it is worth remarking that other differences exist between turbulence in three and two dimensions (Bos 2021). An exploration of these connections in the context of small-scale turbulence theory will be undertaken as future work.