Extreme velocity gradients in turbulent flows

TL;DR

This study uses high-resolution simulations to analyze extreme velocity gradients in turbulent flows, revealing new scaling laws and challenging existing theories about the smallest flow scales.

Contribution

It provides novel empirical scaling laws for extreme velocity gradients and smallest scales in turbulence, based on unprecedented numerical resolution.

Findings

Velocity gradients scale as R_λ^0.775, larger than predicted.

Velocity increments at small scales can reach the r.m.s. of fluctuations.

Smallest length scale behaves as η R_λ^{-α}, contradicting existing theories.

Abstract

Fully turbulent flows are characterized by intermittent formation of very localized and intense velocity gradients. These gradients can be orders of magnitude larger than their typical value and lead to many unique properties of turbulence. Using direct numerical simulations of the Navier-Stokes equations with unprecedented small-scale resolution, we characterize such extreme events over a significant range of turbulence intensities, parameterized by the Taylor-scale Reynolds number ($R_\lambda$). Remarkably, we find the strongest velocity gradients to empirically scale as $\tau_K^{-1} R_\lambda^{\beta}$, with $\beta \approx 0.775 \pm 0.025$, where $\tau_K$ is the Kolmogorov time scale (with its inverse, $\tau_K^{-1}$, being the {r.m.s.} of velocity gradient fluctuations). Additionally, we observe velocity increments across very small distances $r \le \eta$, where $\eta$ is the Kolmogorov length scale, to be as large as the {r.m.s.} of the velocity fluctuations. Both observations suggest that the smallest length scale in the flow behaves as $\eta R_\lambda^{-\alpha}$, with $\alpha = \beta - \frac{1}{2}$, which is at odds with predictions from existing phenomenological theories. We find that extreme gradients are arranged in vortex tubes, such that strain conditioned on vorticity grows on average slower than vorticity, approximately as a power law with an exponent $\gamma < 1$, which weakly increases with $R_\lambda$. Using scaling arguments, we get $\beta=(2-\gamma)^{-1}$, which suggests that $\beta$ would also slowly increase with $R_\lambda$. We conjecture that approaching the limit of infinite $R_\lambda$, the flow is overall smooth, with intense velocity gradients over scale $ \eta R_\lambda^{-1/2}$, corresponding to $\beta = 1$.

Tables

Table 1. Table 1: Simulation parameters for the DNS runs used in the current work: the Taylor-scale Reynolds number ( R λ subscript 𝑅 𝜆 {R_{\lambda}} ), the number of grid points ( N 3 superscript 𝑁 3 N^{3} ), spatial resolution ( k m ​ a ​ x ​ η subscript 𝑘 𝑚 𝑎 𝑥 𝜂 k_{max}\eta ), ratio of large-eddy turnover time ( T E subscript 𝑇 𝐸 T_{E} ) to Kolmogorov time scale ( τ K subscript 𝜏 𝐾 \tau_{K} ), length of simulation ( T 𝑇 T ) in stationary state in terms of turnover time and the number of instantaneous snapshots ( N s subscript 𝑁 𝑠 N_{s} ) used for each run to obtain the statistics.

$R_{\lambda }$	$N^3$	$k_{max}\eta$	$T_E/{\tau }_K$	$T/T_E$	$N_s$
140	${1024}^3$	5.82	16.0	6.5	24
240	${2048}^3$	5.70	30.3	6.0	24
390	${4096}^3$	5.81	48.4	2.8	28
650	${8192}^3$	5.65	74.4	1.1	35

Equations

$${\tau }_{ext}={\tau }_K\times {R_{\lambda }}^{-\beta },\beta >0,$$

$$\partial \mathbf{u}/\partial t+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\nabla p/\rho +\nu {\nabla }^2\mathbf{u}+\mathbf{f},$$

$$f_X(x)=a\exp \left(-b{x}^c\right),$$

$$f_X(x_e)=a{R_{\lambda }}^{2\beta }\exp \left(-b{R_{\lambda }}^{2\beta c}{x_e}^c\right).$$

$$b{R_{\lambda }}^{2\beta c}=b^{\prime },a{R_{\lambda }}^{(2\beta +\delta )}=a^{\prime },$$

$$⟨\Sigma |\Omega ⟩{\tau }_K^2\sim (\Omega {\tau }_K^2{)}^{\gamma },\gamma <1.$$

$${\tau }_{ext}\sim {\eta }_{ext}/\delta u_{max}.$$

$$\delta u{\eta }_x/\nu \approx 1,$$

$${\eta }_{ext}\sim \eta {R_{\lambda }}^{-\alpha },\mathrm{where}\alpha =\frac{1-3{\mathrm{h}}_{\min }}{2(1+{\mathrm{h}}_{\min })}$$

$${\tau }_{ext}$$

$$\delta u_{max}$$

$${\eta }_{ext}\sim \eta {R_{\lambda }}^{-\alpha },\mathrm{with}\alpha =\beta -0.5,$$

$${\tau }_{ext}\sim {\tau }_K{R_{\lambda }}^{-1/2}.$$

$${\tau }_{ext}\sim \frac{\nu }{u^{\prime 2}}=\frac{\nu }{u_K^2}\frac{u_K^2}{u^{\prime 2}}\sim {\tau }_K{R_{\lambda }}^{-1},$$

$${\tau }_{ext}\sim {\tau }_K{R_{\lambda }}^{-\beta },\mathrm{with}\beta =\frac{1}{2-\gamma }.$$

$$\Delta x/\eta \approx {\left(R_{\lambda }/{R_{\lambda }}^{\ast }\right)}^{-\alpha },$$

Full text

I Introduction

Quantitative studies of turbulence in incompressible flows reveal that the averaged dissipation rate of turbulent kinetic energy, $\langle\epsilon\rangle$, is independent of kinematic viscosity, $\nu$, when $\nu\rightarrow 0$ or equivalently when the turbulence intensity, i.e., the Reynolds number, is very high FalkSreeni:06 ; KM.2013 . This empirical result, also known as the zeroth law of turbulence, implies that the amplitude of velocity gradients grows on average as $(\langle\epsilon\rangle/\nu)^{1/2}$. However, the fluctuations of velocity gradients are orders of magnitude larger than this average value, a phenomenon referred to as small-scale intermittency Frisch95 ; Sreeni97 . Such extreme events play a crucial role in numerous physical processes in both nature and engineering, e.g. turbulent dispersion falkovich01 , cloud physics shaw03 , turbulent combustion in jet engines Sreeni04 ; hamlington12 , and are also conjectured to be connected to regularity and smoothness of fluid equations leray ; Fefferman . Hence, understanding their formation and statistical properties is of central importance in developing a complete theory of turbulence Sreeni97 . The complexity of the problem is apparent in Fig. 1, which shows the structure the velocity gradients. The strongly intermittent nature of turbulence is clearly visible by the highly inhomogeneous distribution of the regions of very intense gradients (see also KM.2013 ). Fluid turbulence involves a wide range of spatial scales, from approximately the system size all the way down to the very finest scale, corresponding to the largest gradients. In this respect, it can be viewed as an emblematic example for other complex dynamical systems, where such extreme events are also observed Solli:2007 ; PhysRevX.8.011017 , including the climate system Rahmstorf:2011 , with its far-reaching implications.

Ever since Kolmogorov formulated and refined his seminal hypotheses K62 , intermittency in turbulence has been the subject of many studies Sreeni97 . In particular, detailed investigations demonstrate that the very large fluctuations in velocity gradients become more extreme with increasing Reynolds number MS91 ; Donzis:08 . While there have been theoretical proposals to describe quantitatively the Reynolds number dependence of velocity gradient fluctuations Paladin87 ; Sreeni88 ; Nelkin90 ; YS:05 , they have remained difficult to verify due to lack of reliable data. In fact, directly measuring the most intense fluctuations, experimentally or numerically, over a reasonable range of Reynolds number is a very challenging endeavor, as very high spatial and temporal resolution is required to accurately resolve such fluctuations. As recently pointed out PK+18 , this demand can be even stricter in numerical simulations than previously expected. Consequently, such high resolution investigations have been so far restricted to low Reynolds numbers Donzis:08 ; Schum+07 .

In this work, we characterize the dependence of the extreme velocity gradients on the Reynolds number, and illuminate the underlying physical processes. To this end, we use high resolution direct numerical simulations (DNS) of isotropic turbulence, based on highly accurate Fourier pseudo-spectral methods. To accurately resolve the extreme gradients, all our simulations were carried out with a small-scale resolution at least 3-4 times higher than typical turbulence simulations, along with appropriate temporal resolution PK+18 . Going up to grids of $8192^{3}$ points, we have obtained results at Taylor-scale Reynolds number (${R_{\lambda}}$) ranging from 140 to 650.

In order to characterize the gradients, we first consider the probability density functions (PDFs) of square of the norm of strain and vorticity, which represent the symmetric and skew-symmetric components of the velocity gradient tensor respectively. They are analogous to dissipation rate and enstrophy and have the same mean values (within a prefactor) in isotropic turbulence, given by $1/\tau_{K}^{2}$, where $\tau_{K}=(\nu/\langle\epsilon\rangle)^{1/2}$ is the Kolmogorov time scale Donzis:08 . Consistent with previous works MS91 ; Yeung12 , we observe that the PDFs of these quantities, when normalized by $\tau_{K}$ exhibit tails that become broader with increasing ${R_{\lambda}}$. By further characterizing these PDFs, we demonstrate that their tails can be collapsed very well over the range of ${R_{\lambda}}$ covered here, when instead normalized by the time scale:

[TABLE]

which implies that the strongest gradients in the flow correspond to a time scale $\tau_{ext}$ which increasingly decreases with respect to $\tau_{K}$ as ${R_{\lambda}}$ increases (and hence the strongest gradients in the flow grow as $\tau_{K}^{-1}{R_{\lambda}}^{\beta}$). Numerically, we find that $\beta\approx 0.775\pm 0.025$. The tails of the PDFs of velocity increments $\delta u_{r}$, normalized by the Kolmogorov velocity scale $u_{K}(=(\nu\langle\epsilon\rangle)^{1/4})$, over distances $r\leq\eta$ (where $\eta=(\nu^{3}/\langle\epsilon\rangle)^{1/4}$ is the Kolmogorov length scale), also become broader when ${R_{\lambda}}$ increases. On the contrary, when normalized by the r.m.s. of velocity fluctuations $u^{\prime}$, the tails grow very slowly. With the understanding that the most intense gradients in the flow occur with velocity increments of order $u^{\prime}$ over a scale $\eta_{ext}$ we conclude that $\eta_{ext}\sim\eta{R_{\lambda}}^{-\alpha}$, with $\alpha=\beta-\nicefrac{{1}}{{2}}$, represents the smallest scale in the flow. The collapse of the tails of PDFs of $\delta u_{r}$, when normalized by either $u_{K}{R_{\lambda}}^{\beta}$ or $u^{\prime}{R_{\lambda}}^{\alpha}$, supports these findings. Comparisons with existing theoretical predictions Paladin87 ; YS:05 , point to difficulties in explaining our data. However, these theories utilize the phenomenological definition that the smallest scales in the flow correspond to a local Reynolds number of unity Frisch95 , which, contrary to the numerical results of Jimenez93 ; Jimenez98 and also our own, does not appear to be satisfied at the location of intense gradients, as utilized in current work to characterize the smallest scales.

Consistent with earlier works Siggia:81 ; Jimenez93 ; Ishihara09 , we find that the structures corresponding to the largest velocity gradients to be vortex tubes. We do not find extreme events in strain and vorticity to be colocated Donzis:08 ; Yeung15 . Conditional averaging shows that intense strain is always likely to be accompanied by equally intense vorticity. However, intense vorticity is found to be accompanied by relatively less intense strain, with an approximate power law dependence corresponding to exponent $\gamma<1$, which very slowly increases with ${R_{\lambda}}$. With the interpretation that $\eta_{ext}$ is the radius of most intense vortex tubes, we use simple scaling arguments to relate it to the conditional strain, and thereby relate $\gamma$ to $\beta$. This suggests that $\beta$ would also slowly increase with ${R_{\lambda}}$. We conjecture that the limit $\beta=1$ (and $\alpha=0.5$), as predicted by some intermittency theories, would only be realized for ${R_{\lambda}}\rightarrow\infty$.

The rest of the manuscript is organized as follows. In Section II, we describe our numerical methods. Our numerical results concerning the scaling of extreme velocity gradients are presented in Section III. The structure of regions of very intense velocity gradients is investigated in Section IV. Section V contains a discussion, comparing our results with existing theories, and then providing an an alternative description connected to flow structure examined in Section IV. We briefly discuss the implications of our results on future DNS and experiments in Section VI. Finally, we present our conclusions in Section VII.

II Numerical approach and database

The present work is based on DNS of the incompressible Navier-Stokes equations

[TABLE]

where $\mathbf{u}$ is the velocity field (satisfying $\nabla\cdot\mathbf{u}=0$), $p$ is pressure, and $\mathbf{f}$ is the forcing term used to maintain a stationary state EP88 ; DY2010 . The equations are solved utilizing a massively parallel implementation of Rogallo’s pseudo-spectral algorithm Rogallo , whereby the aliasing errors are controlled by a combination of truncation and phase-shifting PattOrs71 . We use explicit second-order Runge-Kutta scheme for time integration, with the time step $\Delta t$ subject to a constraint for numerical stability expressed in terms of the Courant number, $C=\frac{\Delta t}{\Delta x}\left(||\mathbf{u}||_{1}\right)_{max}$, where $||\cdot||_{1}$ represents the $L^{1}$-norm and the maximum is taken over all ($N^{3}$) grid points. The flow simulated is homogeneous and isotropic with periodic boundary conditions, on a cubic domain of $(2\pi)^{3}$ for all cases.

As stressed earlier, appropriate numerical resolution of the small scales is crucial to our study of extreme velocity gradients. Spatial resolution in pseudo-spectral DNS is typically measured by the parameter $k_{max}\eta$, where $k_{max}=\sqrt{2}N/3$ is the largest wavenumber resolved and $\eta$ is the Kolmogorov length scale. Equivalently, one can use the ratio $\Delta x/\eta$ ($\approx 2.96/k_{max}\eta$), where $\Delta x=2\pi/N$ is the grid spacing. Most turbulence simulations, aimed at reaching high Reynolds number, are in the range $1\leq k_{max}\eta\leq 2$ Yeung15 ; Ishihara16 . However, resolution studies have shown that such a resolution is inadequate for studying extreme events in velocity gradients Schum+07 ; Donzis:08 ; PK+18 . Hence, we have consistently used $k_{max}\eta\approx 6$ in all the runs shown here. Additionally, we have also used a Courant number of $0.3$, instead of 0.6 in previous studies e.g. Donzis:08 ; Yeung12 ; Yeung15 , as it was recently found that the latter led to spurious over-prediction of the gradients PK+18 . Resolution studies presented in PK+18 and our own tests confirm that the resolution used here is adequate to address the questions asked in this work.

The database used here and the corresponding simulation parameters are listed in Table 1. The Taylor-scale Reynolds numbers (${R_{\lambda}}$) considered here are similar to those in some previous works Donzis:08 ; Yeung12 , but with a much higher small-scale resolution as emphasized earlier.

These high resolution simulations were recently used in PK+18 . In the present work, we simply restarted these runs (which were already in a stationary state) and extended them to substantially longer times to greatly improve statistical convergence. We list the length of the current simulation $T$ in terms of the large-eddy turnover time $T_{E}$. The statistical results shown here were obtained by analyzing $N_{s}$ instantaneous snapshots for each run. Whereas a long simulation is typically desirable for sampling accuracy, finite resources have limited the value of $T$ for the highest resolution runs. Nevertheless, since our focus is on highly intermittent velocity gradients, one can improve sampling by simply analyzing more snapshots for a given simulation length as the Reynolds number increases. This is justified, both by the increase of the ratio of time scales $T_{E}/\tau_{K}$ with ${R_{\lambda}}$, and also by the increasingly smaller time scales associated with the extreme events, as discussed in the manuscript.

III Scaling of extreme velocity gradients

III.1 PDFs of vorticity and strain

In order to study small-scale intermittency, in this sub-section, we characterize the velocity gradient tensor by its two quadratic invariants Sreeni97 , namely $\Omega=\omega_{i}\omega_{i}$, where $\mathbf{\omega}=\nabla\times\mathbf{u}$ is the vorticity, and $\Sigma=2s_{ij}s_{ij}$, where $s_{ij}$ is the strain rate tensor defined as $s_{ij}=(\partial u_{i}/\partial x_{j}+\partial u_{j}/\partial x_{i})/2$. The former is the enstrophy and the latter is simply the dissipation divided by viscosity i.e. $\Sigma=\epsilon/\nu$. In isotropic turbulence, as considered here, $\langle\Omega\rangle=\langle\Sigma\rangle=1/\tau_{K}^{2}$, where $\tau_{K}$ is the Kolmogorov time scale.

Investigating the extreme events in $\Omega$ or $\Sigma$ amounts to focusing on the outmost parts of the (wide) tails of their PDFs. Issues of statistical convergence makes a precise determination of these quantities extremely difficult. With the available data, we estimated the statistical error in each bin, by looking at the sample to sample fluctuations across various snapshots used to determine the statistics. We kept only bins with an error less than $20\%$ compared to the mean. These PDFs, converged with respect to both the small-scale resolution PK+18 and statistical sampling, allowed us to determine the properties of the extreme events, as presented next.

Fig. 2 shows the PDF of (a) $\Omega$ and (b) $\Sigma$, normalized by their mean value, $1/\tau_{K}^{2}$, at various Reynolds numbers. We primarily observe that as the Reynolds number increases, the tails of these PDFs get wider and extend to much higher values. Equivalently, the likelihood of finding a value of $\Omega\tau_{K}^{2}$ or $\Sigma\tau_{K}^{2}$ larger than a given large value increases with the Reynolds number. This is expected and consistent with previous studies. One notices, however, that the part of the PDFs, corresponding to events smaller than about 10 times the mean, appear to approximately collapse for different Reynolds numbers. This can be seen in the insets of Fig. 2, which show a zoomed in version.

The existence of increasingly large fluctuations, as shown in Fig. 2, leads us to ask how large are the extreme gradients and how quickly do they grow with increasing Reynolds number. We propose to answer this by rescaling the PDFs. Namely, we use a different time scale, $\tau_{ext}=\tau_{K}\times{R_{\lambda}}^{-\beta}$ (see Eq. (1)), to rescale the extreme values. Denoting $f_{\Omega}(\Omega_{e})$ and $f_{\Sigma}(\Sigma_{e})$ as PDFs of $\Omega_{e}=\Omega\tau_{ext}^{2}$ and $\Sigma_{e}=\Sigma\tau_{ext}^{2}$ respectively, Fig. 3 shows ${R_{\lambda}}^{\delta}f_{\Omega}(\Omega_{e})$, and ${R_{\lambda}}^{\delta}f_{\Sigma}(\Sigma_{e})$. The factor ${R_{\lambda}}^{\delta}$ provides a measure of how rare the largest fluctuations of $\Omega_{e}$ or $\Sigma_{e}$ are, when ${R_{\lambda}}$ increases. As shown in Fig. 3, using $\beta\approx 0.775$ and $\delta\approx 4.0$, the wide tails of rescaled PDFs are almost perfectly collapsed. This indicates that while the average events in $\Omega$ and $\Sigma$ scale as $\tau_{K}^{-2}$, the most extreme events behave like $\tau_{K}^{-2}{R_{\lambda}}^{2\beta}$.

The exponents $\beta$ and $\delta$, used in Fig. 3 to collapse the large tails of the PDFs, can be also empirically determined by utilizing a functional form of the tails of PDFs of $\Omega$ and $\Sigma$. While theories have proposed several functional forms for the entire range of PDFs benzi91 ; Sreeni97 ; yakhot06 ; wilczek09 , the stretched exponential function is known to empirically fit the tails of the PDFs very accurately MS91 ; KSS92 ; zeff:2003 ; Donzis:08 ; PK+18 . Since the tails of the PDFs in Fig. 3 collapse, we use the following stretched exponential functional form to mathematically verify the value of $\beta$:

[TABLE]

where $x=\Omega\tau_{K}^{2}$ or $\Sigma\tau_{K}^{2}$ (or alternatively $\Omega/\langle\Omega\rangle$ and $\epsilon/\langle\epsilon\rangle$ respectively in the notation of Donzis:08 ; PK+18 ) and $a$, $b$, $c$ are the fitting parameters. Applying a change of variable $x_{e}=x\times(\tau_{K}/\tau_{ext})^{2}$, where $x_{e}=\Omega_{e}$ or $\Sigma_{e}$, the PDF of $x_{e}$ becomes

[TABLE]

The collapse shown in Fig. 3 implies that

[TABLE]

such that the constants $b^{\prime}$ and $a^{\prime}$ are independent of ${R_{\lambda}}$. Thus, the dependence of $b^{1/c}$ as a function of ${R_{\lambda}}$, provides a direct access to $\beta$.

To determine the coefficients, we simply fit the logarithm of the PDF to the functional form $(\log{a}-bx^{c})$ of Eq. (3). We choose the fitting window to be $x\geq 50$, which sufficiently excludes the region around the mean value, where the PDFs appear to collapse for various ${R_{\lambda}}$ (as shown in insets of Fig. 2). We also explicitly checked by extending the fitting range to smaller values, but found that the results remained virtually unchanged. The determination of the three parameters $a$, $b$ and $c$ then leads to a non-linear regression. However, since non-linear regression can be very sensitive to the initial guess values for the fitting parameters – especially the value of the exponent $c$ in this particular case – determining the parameters directly in such a manner can result in significant error Fletcher87 . On the other hand, if the value of $c$ is known beforehand, then a very robust fit can be obtained, since the fitting procedure reduces to a linear regression to determine only $a$ and $b$. The exact values of $a$ and $b$ would obviously substantially differ for different values of $c$, but this would not matter if they all provide the same value of $\beta$ (which as shown next, is the case).

The values of the exponent $c$ in previous numerical studies Donzis:08 ; PK+18 were found to be close to the range $0.23-0.25$, with a possible scatter within $0.19-0.29$ and no clear dependence on Reynolds number (e.g. see Table 4 of Donzis:08 ). Keeping this in mind, we therefore fit the PDF by assuming fixed values of $c$, ranging from $0.19$ to $0.29$ in increments of 0.02, and determine the parameters $b$ and $a$ for PDFs $\Omega$ and $\Sigma$ for all available ${R_{\lambda}}$. Note that a wider range for $c$ may be considered, but this chosen range falls within the error obtained from a naive non-linear regression and hence for values outside the chosen range, the quality of fit starts deteriorating. To provide a measure, for the chosen values of $c$, the coefficient of determination ($R^{2}$) was greater than $0.995$ for every fit. Additionally for each $c$, the resulting values of $a$ and $b$ are always obtained with greater than 95% confidence, resulting in negligible error bars. In fact, these values are even found to be quite insensitive to minor variations in the fitting window, e.g., our fits compare extremely well with those of Donzis:08 , who considered a fitting window of $5\leq x\leq 100$ for ${R_{\lambda}}\leq 240$. In this regard, we make a note that the results of Donzis:08 can only be trusted for ${R_{\lambda}}\leq 240$, since the higher ${R_{\lambda}}$ runs were affected by resolution issues, as reported in PK+18 . Nevertheless, the excellent quality of fit is evident in Fig. 3 (also see Fig.9 of Donzis:08 which is also in near perfect agreement with our fits).

The dependence of the coefficient $b^{1/c}$ on ${R_{\lambda}}$ at various values of $c$ is shown in Fig. 4. The data points correspond to curve fits for both $\Omega$ and $\Sigma$ (thus giving two sets of points for each $c$). For the sake of clarity, the values of $b^{1/c}$ are divided by their corresponding values at ${R_{\lambda}}=650$, which imposes that all the curves shown in Fig. 4 pass through $1$ at ${R_{\lambda}}=650$ (since higher ${R_{\lambda}}$ provides a larger fitting range and hence can be expected to be most robust). This also allows us to directly compare the data points for every $c$ value considered. We find that all sets of points superpose reasonably well, and remarkably point to a similar power law in ${R_{\lambda}}$. This collapse demonstrates that the determination of the exponent $\beta$ is not very sensitive to the precise value of $c$, at least with the available data. However, weak deviations from scaling cannot be ruled out, especially if an even larger range of ${R_{\lambda}}$ is considered in future. In fact, we will later (in Section V) present arguments supporting a very weak growth of $\beta$ with ${R_{\lambda}}$.

We would like to further clarify that ideally the best curve fits to Eq. (3) may lead to a dependence of $c$ on the Reynolds number ${R_{\lambda}}$. However, choosing a fixed $c$ substantially improves the quality of curve fit and also minimizes the sensitivity to the fitting range. In fact, a fixed value of $c$ also helps in determining the scaling without ambiguity, since if $c$ is a function of ${R_{\lambda}}$, the constant $b^{\prime}$ in Eq. (5) will also become a function of ${R_{\lambda}}$, which would recursively require additional non-linear fits to obtain $\beta$. The minor deviations for different $c$ values at various ${R_{\lambda}}$, can also be possibly explained by this. Nevertheless, the good collapse seen for such a wide range of $c$ values provides clear evidence that the scaling proposed provides a compelling description of our data, at least over the available range of ${R_{\lambda}}$. We will see that this is also further supported by results shown in section III.2. Finally, by fitting a power law through the obtained data points (marked by dashed line), we obtain $\beta=0.775\pm 0.025$ (where the error bar takes into account the variation across different $c$ values), which was used in scaling the PDFs in Fig. 3. The same procedure for the parameter $a$ (not shown) gives $\delta\approx 4.0$, with deviations of approximately 5-10%. Thus, systematically characterizing the PDFs of $\Omega$ and $\Sigma$, we are able to mathematically determine that the strongest gradients in the flow grow as $\tau_{K}^{-1}{R_{\lambda}}^{\beta}$, with $\beta=0.775\pm 0.025$. We again emphasize that obtaining such a result required statistically well converged PDFs (in turn requiring adequate spatial and temporal resolutions) over a wide enough range of Reynolds numbers.

III.2 PDFs of velocity increments

In order to further validate the scaling obtained in Section III.1, we next investigate the PDFs of velocity increments. In simplified notation, velocity increments are given as $\delta u_{r}=u(x+r)-u(x)$, where the separation distance $r$ can be either in the direction of $u$ (longitudinal increments) or perpendicular to $u$ (transverse increments). Over very small distances, the velocity differences essentially reduces to the velocity gradients (within a constant factor), aside from systematic but small errors introduced by use of finite differencing. Hence, we can expect velocity increments over small distances to show the same scaling as derived earlier. However, to extract information about the gradients, one still needs to ensure that $r$ is sufficiently small. We note in this respect that the high resolution of our runs, $k_{max}\eta\approx 6$, effectively allows us to calculate increments over a very small distance ($r\approx\eta/2$). A benefit of using velocity increments is that their 1D surrogates can be also obtained and verified using experiments EB2014 .

Fig. 5a shows the PDF of $\delta u_{r}$, normalized by the Kolmogorov velocity scale $u_{K}$, for $r/\eta\approx 0.5$, at various Reynolds numbers. Both the longitudinal and transverse components are shown (in dashed and solid lines respectively). Since we are interested only in the magnitude of the increments, we take the absolute value of $\delta u_{r}$. Consistent with the results shown in Fig. 2, the initial part of the PDFs, corresponding to moderate events, superpose very well (see inset of Fig. 5a) and as Reynolds number increases, the tails start growing. Both the longitudinal and transverse increments show the same behavior, with the transverse component being expectedly larger Frisch95 .

To further characterize the velocity increments, we next consider the PDFs of $\delta u_{r}/u^{\prime}$, where $u^{\prime}$ is the r.m.s. of velocity fluctuations. Fig. 5b shows the PDFs of $\delta u_{r}/u^{\prime}$ corresponding to $r=\eta/2$, at various Reynolds numbers and for both longitudinal and transverse components. The corresponding PDFs for $r=\eta$ are shown in the inset. Even at such a small separations, we observe that the velocity differences can be as high as $u^{\prime}$, for the entire range of Reynolds number considered. This appears to be consistent with the observation of Jimenez93 ; Jimenez98 at ${R_{\lambda}}\lesssim 170$. Additionally, it appears that while the probability density for a given $\delta u_{r}/u^{\prime}$ decreases with increasing ${R_{\lambda}}$, the extent of $\delta u_{r}/u^{\prime}$ itself slowly increases.

Following similar ideas as in Fig. 3b, we determine the PDFs of $\delta u_{r}/u_{K}$, rescaled by ${R_{\lambda}}^{-\beta}$, shifted by a factor ${R_{\lambda}}^{\delta}$. The result is shown in Fig. 5c corresponding to $r/\eta\approx 0.5$ for both longitudinal and transverse components. Note for the PDFs of $\delta u_{r}/u^{\prime}$ this corresponds to rescaling by ${R_{\lambda}}^{-\alpha}$, with $\alpha=\beta-0.5$, since $u^{\prime}/u_{K}\sim{R_{\lambda}}^{1/2}$ Frisch95 (the importance of $\alpha$ is discussed later in Section V). We find that the rescaled PDFs collapse very well for different Reynolds numbers. The scatter towards the very end of the tails (especially for the transverse component) can be attributed to lack of statistical convergence for the endmost bins. In the inset of Fig. 5c, we repeat the exercise, but now for PDFs corresponding to $r/\eta\approx 1$. The superposition, although comparatively worse with respect to $r/\eta=0.5$, still remains very good.

An alternative approach to investigate velocity increments is to consider the quantity $\delta u_{r}/r$, which for sufficiently small $r$, is a proxy for the velocity gradient and thus also independent of $r$. Hence, it is tempting to use the same scaling based on $\tau_{ext}$ as before, to collapse the tails of various PDFs of $\delta u_{r}/r$ (as done in Fig. 5c). Fig. 6 shows such rescaled PDFs for several values of $r/\eta$, at ${R_{\lambda}}=140$ and $650$. Once again, we find that the tails of the curves for $r/\eta\approx 0.5$ collapse very well for both ${R_{\lambda}}$. However at $r/\eta\approx 1$, the curve for ${R_{\lambda}}=650$ starts to deviate from this collapse. At $r/\eta\approx 2$, the curve for ${R_{\lambda}}=650$ significantly deviates from that for ${R_{\lambda}}=140$. This provides yet further evidence that the resolution of $\Delta x/\eta\approx 0.5$ is adequate for both ${R_{\lambda}}$. On the other hand, $\Delta x/\eta\approx 1$, while adequate for ${R_{\lambda}}=140$, is insufficient for ${R_{\lambda}}=650$. Additionally, we also see that as $r/\eta$ grows, the deviations of the curves at ${R_{\lambda}}=650$ increase faster than those at ${R_{\lambda}}=140$. This also provides a hint that the smallest length scale in the flow is actually smaller than $\eta$ and additionally decreasing with increasing ${R_{\lambda}}$. This observation will be further analyzed and discussed in Section V.

IV Structure of regions of intense vorticity and strain

In order to gain some understanding on the structure of regions of intense gradients, we first use flow visualization. Vorticity arranged in tube-like structures has been repeatedly seen in DNS over a very large range of Reynolds numbers – from very low at ${R_{\lambda}}\approx 45$ Siggia:81 all the way to ${R_{\lambda}}\approx 1100$ Ishihara09 (although the small-scale resolution in these simulations was limited). Whether such tubes carry the most intense regions of velocity gradients in the flow, however, has been questioned by a recent study Yeung15 , which suggested that the largest values of $\Omega$ and $\Sigma$ appear colocated and without any coherent structure. It is important to note that these observations may have been affected by the numerical artefacts documented in PK+18 . One of the motivations of the present work is to revisit the issue. We again stress that the present visualizations are based on DNS at much higher spatial resolution than previously available.

Fig. 7 shows a collection of instantaneous snapshots from the ${R_{\lambda}}=650$ run, focusing on the region of most intense gradients. Since vorticity is in general larger than strain, the domains are chosen such their centers correspond to the maximum value of $\Omega$. The various panels show different contour thresholds (indicated as C in sub-figure caption) in cyan for $\Omega\tau_{K}^{2}$ and in red for $\Sigma\tau_{K}^{2}$. In the first panel (Fig. 7a), a domain of $301^{3}$ grid points or $(150\eta)^{3}$ is shown with the contour threshold of 50 for both vorticity and strain. Structures consisting of clusters of vortex tubes, qualitatively similar to e.g. Ishihara09 , are readily seen. The spacing between neighboring tubes widely varies: some vortices are relatively isolated, others seem to be more strongly interacting with their surrounding. Large values of strain are mostly located around large vortices, a phenomenon noticed many times (see Yeung15 and references therein).

Panel (b) zooms into the region of most intense gradients, showing a domain of $(50\eta)^{3}$ with a contour threshold of 300. The structure is composed of two closely interacting vortex tubes, wrapped around by intense strain. In panels (c) and (d), contour levels are successively increased to 500 and 1500 respectively and we also further zoom in to show a domain of $(25\eta)^{3}$ in (d). The vortex tube structure becomes very distinct, whereas the region occupied by strain reduces substantially in (c) and completely disappears in (d). Note, the largest value of $\Omega\tau_{K}^{2}$ is equal to about 3000 (at the center of the domain in each panel). In comparison, the largest value of $\Sigma\tau_{K}^{2}$ is about 1800, located in top left corner of domain shown in (c) – and hence no coherent strain region is visible in (d). Although not explicitly shown here, we also confirmed that the velocity increments around the center of each panel Fig. 7 correspond to far tails of PDF of $\delta u_{r}$ as shown in Fig. 5, i.e., $\delta u_{r}\simeq u^{\prime}$.

We analyzed many such flow fields corresponding to different snapshots and virtually all of them show a qualitatively similar behavior, i.e., as the contour thresholds are increased tube-like vorticity structures become prominent and high-strain regions shrink and disappear at lower values than the high-vorticity regions. While not directly evident in Fig. 7, we also find that the locations of maximum values of vorticity and strain are typically separated by at least $10-20\eta$ and never coincident, e.g. in Fig. 7c, These observations confirm that the largest values of $\Sigma$ are much smaller than the large values of $\Omega$ and the regions for large values of $\Sigma$ and $\Omega$ are not co-located. Hence, we conclude that visualizations in Donzis:08 ; Yeung15 were also affected by resolution issues reported in PK+18 . We remark in this respect that new independent tests at ${R_{\lambda}}=1300$ and $k_{max}\eta=3$, with a time step twice smaller than in Yeung15 - although not shown here - confirm that the qualitative aspect of the regions of extreme vorticity/strain are similar to that shown in Fig. 7.

In order to quantify the relation between strain and vorticity, we next consider their conditional expectations with respect to each other – shown in Fig. 8 for various ${R_{\lambda}}$. For low values of $\Omega$ or $\Sigma$, the conditional dependencies are very weak, i.e., strain and vorticity appear to be decorrelated. However, for conditional values greater than unity, i.e., the mean value, the conditional expectations clearly increase, seemingly showing a power law. Comparison with a dashed line of slope 1 (on log-log coordinates), suggests that $\langle\Omega|\Sigma\rangle\sim\Sigma^{1}$. In contrast

[TABLE]

This implies that intense events in strain are always likely to be accompanied by equally strong events in vorticity, whereas the strain is comparatively weaker in very intense vortices. This appears to be consistent with the earlier observations of vorticity being more intermittent than strain Siggia:81 ; Chen97 ; Donzis:08 and ultimately concerns with the inter-relationship of vorticity and strain, which is still an open question in turbulence. Note that Donzis:08 shows a similar plot as Fig. 8, however their curve for $\langle\Sigma|\Omega\rangle$ spuriously approaches a slope of 1 (for large $\Omega$) because of resolution issues PK+18 .

Interestingly, Fig. 8b also suggests that the exponent $\gamma$ slowly increases with ${R_{\lambda}}$. By fitting approximate power laws, we find that $\gamma$ varies from $0.60-0.72$ over the range of ${R_{\lambda}}$ considered here (see inset of Fig. 8b), although the variation appears to get weaker as ${R_{\lambda}}$ increases. This naturally leads to the question of what the limit of ${R_{\lambda}}\to\infty$ entails, which our data is unable to answer conclusively. Theoretical considerations suggest that $\gamma=1$ in the large ${R_{\lambda}}$ limit Sreeni97 ; He98 ; Nelkin99 . Given the very slow increasing trend of $\gamma$, it is evident that extremely high ${R_{\lambda}}$ would be necessary to realize $\gamma=1$, if at all possible (a simple sigmoidal or power law extrapolation suggests $\gamma=0.99$ would be realized for ${R_{\lambda}}\gtrsim 20000$). Therefore, the differences between strain and vorticity are expected to persist, even at the highest turbulence levels on earth. A fundamental understanding of Eq. (6) from first principles, i.e. a determination of the strain acting on a given vortex, resulting from the tangle of vortices as shown in Fig. 7 is still an open question in turbulence. As we will suggest in Section V.2, the power law dependence on the strain conditioned on vorticity in fact provides a way to understand the scaling exponent $\beta$.

V Theoretical considerations

In this section, we discuss the observation that the extreme velocity gradient fluctuations scale as $\tau_{ext}^{-1}=\tau_{K}^{-1}\times{R_{\lambda}}^{\beta}$. We first compare our result with existing theories, especially the multifractal model, and thereafter provide a new description for the observed scaling which relates to the structure of the flow discussed earlier in Section IV.

In simplest terms, extreme velocity gradients result from large velocity differences over a very small length scale; the largest velocity gradient in the flow can be written as proportional to $\delta u_{max}/\eta_{ext}$, where $\delta u_{max}$ is largest velocity difference over the smallest length scale $\eta_{ext}$ Nelkin90 , which given the flow structure, can be physically interpreted as the radius of the smallest vortex tube Jimenez93 . Our notation that the largest velocity gradients scale as $\tau_{ext}^{-1}$ therefore implies

[TABLE]

Thus, the question how large the gradients can grow entails answering how large $\delta u$ can become over the smallest length scale $\eta_{ext}$. Based on earlier resolution studies Donzis:08 ; PK+18 and also on the results presented in Section III.2, the smallest scale $\eta_{ext}$ can be defined by the resolution at which the PDFs of the gradients have converged – which for the present range of ${R_{\lambda}}$, gives $\eta/2\leq\eta_{ext}\lesssim\eta$. Notice that one could formally define scales smaller than $\eta_{ext}$, but given a smooth velocity field, the velocity increments at such scales will simply decrease linearly with the scale size, with respect to those at $\eta_{ext}$ (as also demonstrated by the velocity increment PDFs for $r/\eta=0.5$ and $1$ at ${R_{\lambda}}=140$ in Fig. 6). Thus, any length scale smaller than $\eta_{ext}$ would show the same scaling as $\eta_{ext}$ itself and would be immaterial for the purpose of present study.

V.1 Comparisons with existing theories

It is natural to interpret our results using existing theories. To this end, we begin by reviewing the multifractal model, which provides explicit predictions concerning the smallest scales in the flow. In the multifractal model, as well as in some other phenomenological approaches, a recurring concept is that of the fluctuating local viscous cutoff scale, say $\eta_{x}$, defined such that the velocity increment over this scale has a local Reynolds number of unity Paladin87 ; YS:05 :

[TABLE]

which essentially results from equating the viscous time scale $\eta_{x}^{2}/\nu$ to the convective time scale $\eta_{x}/\delta u$. In the multifractal framework, the velocity increment over a distance $r$ is given as $\delta u_{r}/u^{\prime}\sim(r/L)^{h}$, where $L$ is the energy-injection scale, and $h$ is the local Hölder exponent within an interval $\left[h_{min},h_{max}\right]$ such that a fractal set $D(h)$ can be determined for every $h$. Thereafter, following the derivation of Paladin87 , the smallest scale in the flow can readily be obtained corresponding to the minimum Hölder exponent

[TABLE]

where $\eta$ is the Kolmogorov length scale. It also follows

[TABLE]

which implies $\beta=2\alpha$.

Earlier works have suggested that $h_{min}=0$ Paladin87 ; Nelkin90 , which gives $\alpha=0.5$, $\beta=2\alpha=1$ and hence $\delta u_{max}\sim u^{\prime}$. The value of $\beta\approx 0.775$ derived earlier can be obtained by using $h_{min}\approx 0.06$, and would additionally imply $\alpha=\beta/2\approx 0.39$. However, $h_{min}\approx 0.06$, or rather any non-zero positive value of $h_{min}$, implies $\alpha<0.5$, and hence suggests that $\delta u_{max}/u^{\prime}$ would decrease with increasing ${R_{\lambda}}$. In fact, $h_{min}>0$ also suggests that the range of $\delta u_{r}/u^{\prime}$ decreases with ${R_{\lambda}}$ for a fixed $r/\eta$. Our observation in Fig. 5b, which demonstrates that the range of $\delta u_{r}/u^{\prime}$ does not show any sign of decreasing with ${R_{\lambda}}$ at $r/\eta\lesssim 1$ – and rather appears to be slowly increasing, does not unambiguously support the decay of $\delta u_{max}/u^{\prime}$ implied by the theory. At the same time, since $\eta/2<\eta_{ext}\lesssim\eta$ over the present range of ${R_{\lambda}}$, the extent of PDFs in Fig. 5b implies $\delta u_{max}\gtrsim u^{\prime}$ for the strongest gradients. We notice that the probability density for $\delta u\gtrsim u^{\prime}$ appears to decrease slowly with ${R_{\lambda}}$, however, the excellent superposition of the PDFs in Fig. 5c indicates that the decay is at best algebraic (decreasing as ${R_{\lambda}}^{-\delta}$), and therefore, the probability should remain finite even as ${R_{\lambda}}\rightarrow\infty$. The above suggests $h_{min}=0$ to allow $\delta u_{max}\sim u^{\prime}$, leading to $\beta=1$ (for $\alpha=0.5)$ as suggested by Paladin87 ; Nelkin90 , but at odds with the observed value of $\beta=0.775$ (and the corresponding $h_{min}\approx 0.06$).

In our view, the inconsistency above is a result of the assumptions built into the extension of the multifractal theory, originally developed to describe the inertial scales, to far dissipative scales. In particular, the definition in Eq. (8) obtained by equating the convective and dissipative time scales, while reasonable for inertial range (where the rate of energy transfer across scales can be assumed to be constant), does not appear justified at smallest scales, where dissipation dominates. This is readily observed in Fig. 5a, where the local Reynolds number from the tails (corresponding to strongest gradients residing in vortex tubes) increases steadily with the ${R_{\lambda}}$ (much strongly than $\eta_{ext}$ decreases with ${R_{\lambda}}$). In fact, earlier works based on DNS at ${R_{\lambda}}\lesssim 170$, have already have shown that the local Reynolds numbers corresponding to the vortex tubes where the intense gradients are localized are much larger than unity, and appear to scale differently from the multifractal prediction Jimenez93 ; Jimenez98 . Our numerical results at significantly higher ${R_{\lambda}}$ (and also higher small-scale resolution) further strengthen this conclusion and puts into question the (phenomenological) criterion that the smallest scale in the flow can be determined by a local Reynolds number of order unity and hence, also the relation $\alpha=\beta/2$. However, an additional remark is necessary in this context. While we associate the smallest scales of motion with the extreme gradients (which appear to reside in vortex tubes), the theoretical constructs resulting from multifractal considerations are based on the local scaling, where the smallest scales simply correspond to the minimum Hölder exponent $h_{min}$, without any explicit connections to the flow structure. As a result, there is no guarantee that the scales resulting from $h_{min}$ actually correspond to the structures observed in Fig. 7, calling for some caution when comparing our results with the multifractal theory.

An alternative description, which also utilizes Eq. (8), is that of Yakhot and Sreenivasan YS:05 . In their approach, the even moments of the velocity increment ($\delta u)^{2n}$ (or the structure functions), each correspond to a unique dissipative length scale $\eta_{n}$, such that the smallest possible scale in the flow corresponds to $n\to\infty$. Thereafter, again utilizing Eq. (8), $\eta_{n}$ can be related to the anomalous inertial range scaling exponents of structure functions and $n\to\infty$ results in a similar prediction as that of multifractal theory with $h_{min}=0$, i.e., $\alpha=0.5$. However, once again, such an approach does not appear as justified given the lack of evidence for Eq. (8) to define the smallest scales. In fact, numerical results of Donzis:08 ; Schum+07 ; schum07sub all suggest that the smallest scales in fact grow weaker than the prediction of Yakhot-Sreenivasan theory (and hence also that of multifractal theory for $h_{min}=0$). This observation is once again reinforced by the results presented in current work.

V.2 Alternative description in light of strain-vorticity dynamics

In view of the apparent shortcomings of intermittency theories reviewed in the previous subsection, we propose, in order to reconcile our observation concerning the very large velocity differences and the exponent $\beta<1$, a different description that directly relates to the flow structure explored in Section IV.

Based on Fig. 5b, we propose that the strongest gradients correspond to $\delta u_{max}\sim u^{\prime}$ over the smallest scale $\eta_{ext}$. This is in line with observations of Jimenez93 ; Jimenez98 and also $h_{min}=0$ as postulated by Paladin87 . Note, this assumption also essentially implies that the PDFs of velocity (and hence velocity increments) are bounded leray . Thereafter, substituting $\delta u_{max}\sim u^{\prime}$ and $\tau_{ext}=\tau_{K}{R_{\lambda}}^{-\beta}$ into Eq. (7) gives

[TABLE]

where we have used $u^{\prime}/u_{K}\sim{R_{\lambda}}^{1/2}$ from classical scaling estimate Frisch95 . The value $\beta=0.775\pm 0.025$ found numerically leads to $\alpha=0.275\pm 0.025$, which is also the value used in Fig. 5c to collapse the PDFs of $\delta u_{r}/u^{\prime}$ for $r=\eta/2<\eta_{ext}$. Note, the convergence of $\delta u_{r}/r$ to the velocity gradient ensures that the PDFs have a well-defined limit for $r\leq\eta_{ext}$ (and hence the PDFs at $r<\eta_{ext}$ can be simply obtained by linearly rescaling the PDF at $r=\eta_{ext}$ by a factor $r/\eta_{ext}$, which essentially is the same as ${R_{\lambda}}^{\alpha}$ for $r=\eta/2$.

This value of $\alpha\approx 0.275$ can be further verified by considering the PDFs shown in Fig. 6. While for $r\leq\eta_{ext}$, $\delta u_{r}/r$ converges to the velocity gradient, systematic deviations arise for $r>\eta_{ext}$. These deviations from the gradient can be accordingly quantified by analyzing the higher order terms in a Taylor series expansion of $\delta u_{r}/r$ and can be shown to be approximately proportional to $r/\eta_{ext}$ Donzis:08 ; YS:05 . Since $\eta_{ext}/\eta$ decreases when ${R_{\lambda}}$ increases, at a given value of $r/\eta$, the deviations from the PDFs, especially in the tail, from their limiting form at $r/\eta_{ext}\ll 1$ also increases, as clearly seen in Fig. 6. In addition, we find that the deviations of the PDF tails at fixed values of $r/\eta_{ext}$ to be independent of ${R_{\lambda}}$. For $\alpha\approx 0.275$, $\eta_{ext}/\eta$ decreases by approximately $1.53$ between ${R_{\lambda}}=140$ and $650$. In contrast, taking the value of $\alpha$ predicted by the multifractal theory, $\alpha=\beta/2\approx 0.39$, leads to a variation of $1.82$ in the ratio $\eta_{ext}/\eta$. The $r/\eta$ values at ${R_{\lambda}}=140$ shown in Fig. 6 are chosen, as close as possible, within these factors (of $1.53$ and $1.82)$, compared to the $r/\eta$ values for ${R_{\lambda}}=650$ (the curves corresponding to 1.53 are shown in dashed-dotted lines, whereas curves corresponding to 1.82 are shown in dashed lines). As visible, the PDFs corresponding to the factor of $1.53$ between the two ${R_{\lambda}}$ cases collapse remarkably well (especially as $r/\eta$ increases), hence providing an alternative means to verify $\alpha\approx 0.275$.

Interestingly, alongside some scaling arguments to evaluate $\alpha$, Eq. (12) leads to two different limits (which incidentally also correspond to previously reported cases in literature). The first limit corresponds to simply assuming that the smallest length scale in the flow is the Kolmogorov length scale, i.e., $\eta_{ext}=\eta$. Using this, we get $\alpha=0$ and $\beta=0.5$, and hence

[TABLE]

This result was derived in Jimenez93 ; Jimenez98 , based on the analysis of DNS data at relatively low Reynolds numbers (${R_{\lambda}}\lesssim 170$), which the present work greatly improves upon. The second limit consists in taking into account the extreme fluctuations of the velocity gradients. Physically, the smallest scale in a flow can be thought to result from a balance between viscosity $\nu$ and strain $\Sigma$ ($=2s_{ij}s_{ij}$, as defined earlier), which leads to the expression of the length scale: $\eta_{ext}\simeq(\nu^{2}/\Sigma)^{1/4}$, familiar in a number of contexts Burgers48 . Assuming $\eta_{ext}=\eta$, which leads to Eq. (13), amounts to a mean field approximation, consisting in replacing the strain by its averaged value (as $\eta$ is calculated from the mean dissipation). Taking into account the large fluctuations of $\Sigma$ results in $\eta_{ext}$ being smaller than $\eta$ Sreeni88 . In this regard, the second limit can be simply derived by evaluating $\eta_{ext}$ using the maximum value of strain ($\Sigma_{max}$), i.e., $\eta_{ext}=(\nu^{2}/\Sigma_{max})^{1/4}$. Thereafter, using $\Sigma_{max}\sim\tau_{ext}^{-2}$ and $\delta u_{max}\sim u^{\prime}$ based on earlier results, it follows from Eq. (7)

[TABLE]

where we have used $\nu/u_{K}^{2}=\tau_{K}$ and $u^{\prime}/u_{K}\sim{R_{\lambda}}^{1/2}$. This implies $\beta=1$ and $\alpha=0.5$ from Eq. (12), as also predicted by intermittency models discussed earlier Paladin87 ; Nelkin90 ; YS:05 . However, this is not completely surprising, as defining $\eta_{ext}$ based on $\Sigma_{max}$ with $\Sigma_{max}\sim\nu\tau_{ext}^{-2}$ also leads to $\beta=2\alpha$, which is essentially the multifractal prediction, and in conjunction with $\beta=\alpha+0.5$ derived in Eq. (12) gives $\beta=1$ (and $\alpha=0.5$). Additionally for this scenario, the local Reynolds number, which can be written as $\eta_{ext}^{2}\nu^{-1}\tau_{ext}^{-1}$ using Eq. (7), comes out to be constant as inherently assumed in intermittency theories discussed earlier.

The numerically observed value of $\beta\approx 0.775$ lies between $\beta=0.5$ and $1$, which suggests that $\eta_{ext}$ results from a strain, intermediate between the two limits considered before. In fact, this is precisely what we observed in Section IV. As noted earlier, Eq. (6) (and Fig. 8b) suggests that the strain acting on a very intense vortex tube is significantly weaker than naively expected by postulating $\Sigma\propto\Omega$. A simplified estimate consists in substituting $\Sigma_{max}$ in the argument leading to Eq. (14) by $\tau_{K}^{-2}(\tau_{K}^{2}\Omega_{max})^{\gamma}$, as suggested by Eq. (6). Thereafter, we get

[TABLE]

The limits of $\beta=0.5$ and $1$ correspond to $\gamma=0$ and 1 respectively. In view of the weak dependence of $\gamma$ shown in the inset of Fig. 8b, Eq. (15) suggests a dependence of $\beta$ on ${R_{\lambda}}$. The values of $\gamma$ observed over the range of ${R_{\lambda}}$ studied here, $0.60\lesssim\gamma\lesssim 0.72$, implies a variation of $\beta$ in the range: $0.72\lesssim\beta\lesssim 0.78$, which is quantitatively consistent with $\beta\approx 0.775$ determined empirically in Section III. The weak variation of $\beta$ implied by Eq. (15) may also explain the slight deviations from scaling seen for ${R_{\lambda}}=140$ in Fig. 4. In fact in Fig. 4, considering only data points at ${R_{\lambda}}=140$ and 240, the slope corresponds to $\beta\approx 0.73$, which appears to be remarkably consistent with that obtained from $\gamma$ for these ${R_{\lambda}}$. On the other hand, the interesting possibility that $\gamma\rightarrow 1$ when ${R_{\lambda}}\rightarrow\infty$ would then suggest, in view of Eq. (15), that $\beta\rightarrow 1$, as originally expected by some theories (albeit corresponding to a constant local Reynolds number much larger than unity). However, the very slow variation of $\gamma$ shown the inset of Fig. 8b would indicate that $\beta=1$ would be attained at extremely large values of ${R_{\lambda}}$, likely larger than practically relevant. In this regard, the scope of existing predictions in understanding finite Reynolds number scaling appears to be severely limited.

Whereas the prediction of the exponent $\gamma$ and its dependence on ${R_{\lambda}}$ is a very challenging task, we briefly note that the cascade model of She and Leveque SL94 presents a similar idea, though with shortcomings. The model postulates that the locally averaged dissipation field $\epsilon_{r}$ at a scale $r$ and the corresponding moment ratios: $\epsilon_{r}^{(p)}=\langle\epsilon_{r}^{p+1}\rangle/\langle\epsilon_{r}^{p}\rangle$, are related to the hierarchy of complex structures in the flow. The most singular structures correspond to $\epsilon_{r}^{(\infty)}$, which in the phenomenology of SL94 obeys the power law dependence: $\epsilon_{r}^{(\infty)}\simeq\langle\epsilon\rangle\left(L/r\right)^{\mu}$, with $\mu=2/3$. While their original arguments were postulated for inertial scales, if one were to extend the cascading process down to smallest scale, i.e., $r=\eta_{ext}=\left(\nu^{3}/\epsilon_{r}^{(\infty)}\right)^{1/4}$, it leads to $\alpha=(3\mu)/{(8-2\mu)}$. Using $\mu=2/3$ as proposed by She-Leveque then gives $\alpha=0.3$, which is close to our current prediction of $\alpha=0.275\pm 0.025$. However, they also suggest $h_{min}=1/9$ within the multifractal formalism, which using Eqs. (9)–(11), gives $\alpha=0.3$, $\beta=0.6$ and $\delta u_{max}\sim u^{\prime}{R_{\lambda}}^{-0.2}$, which are clearly inconsistent with our data. Ultimately, phenomenological descriptions (as those of SL94 ) are at best weakly connected to flow structures, and typically assume a constant value of the exponents such as $\beta$ and $\alpha$, thus ignoring any possible dependence on ${R_{\lambda}}$, as suggested from Eq. (15) and Fig. 8. Hence, it appears that the closeness of $\alpha$ between the She-Leveque model and our current result is only fortuitous.

In conclusion, our analysis of the most intense vortex structures observed in the flow relates the exponent $\beta$ with the properties of the strain acting on vortices, and in particular with the exponent $\gamma$ defined by Eq. (6). The weak variation of $\gamma$ with ${R_{\lambda}}$, see Fig. 8, implies that $\beta$ should increase with ${R_{\lambda}}$. A very natural conjecture is that the symmetry between strain and vorticity, clearly broken at finite ${R_{\lambda}}$, will be restored as ${R_{\lambda}}\rightarrow\infty$, and that the exponents $\gamma$ and $\beta$ both tend to $1$, corresponding to earlier predictions Paladin87 ; YS:05 . Understanding the ${R_{\lambda}}$-dependence of the strain acting on intense vortex tubes appears as an essential question in this regard, that deserves renewed theoretical attention.

VI Implications for simulations and experiments

The identification of the smallest scale $\eta_{ext}$, characteristic of the largest velocity gradients in the flow, which decreases faster than $\eta$ when ${R_{\lambda}}$ increases, has some obvious consequences for the resolution constraints required in both DNS and experiments. In DNS, it is typical for most studies based to be performed with a $k_{max}\eta$ or $\Delta x/\eta$ held constant across the range of ${R_{\lambda}}$ simulated (e.g. see Ishihara16 ; BSY.2015 ). On the other hand, in experiments, the resolution, determined by the probe size or the data acquisition frequency, often gets worse as ${R_{\lambda}}$ increases EB2014 . The present results, however, show that in studies focused on intermittency, one must continuously improve $\Delta x/\eta$ as ${R_{\lambda}}$ is increased to adequately resolve the smallest scales, i.e., $\Delta x/\eta_{ext}$ should be held constant across various simulations. In fact, this suggestion was also put forward by YS:05 , though their criterion was stricter than the present numerical results suggest.

The simulations presented here suggest, based on PDFs of various components of the velocity gradient tensor, that a resolution of $\Delta x/\eta\approx 1$ is sufficient for ${R_{\lambda}}=240$, but barely insufficient for ${R_{\lambda}}=390$ (this is also evident from Fig. 6). An earlier resolution study at ${R_{\lambda}}\leq 240$ Donzis:08 , also supports this. In fact, in Donzis:08 , the authors also explored a practical approach to determine the necessary resolution based on calculating the error between $p$-th order structure function and its analytic behavior for small distances (obtained from Taylor series expansion), such that the result is a function of $p$ and ${R_{\lambda}}$. However, such expressions are limited to small values of $p$, since the derivation retained only a small number of terms in the Taylor expansion and hence also cannot be generalized to estimate $\eta_{ext}$.

Nevertheless, based on previous and current results, it follows empirically that $\Delta x/\eta$ to accurately resolve the velocity gradients should be

[TABLE]

where ${R_{\lambda}}^{\ast}\approx 300$ is the reference Taylor-scale Reynolds number, at which $\eta_{ext}\approx\eta$. The above relation provides a practical resolution criteria for future simulations at even larger problem sizes than considered here. While $\alpha\approx 0.275$ for the current range of ${R_{\lambda}}$, we anticipate newer simulations at higher ${R_{\lambda}}$ would progressively update $\alpha$ and also quantify its dependence on ${R_{\lambda}}$ (though given the slow growth of $\alpha$ with ${R_{\lambda}}$, a very substantial range of ${R_{\lambda}}$ might be required). The constraint provided by Eq. (16) should also apply to experimental investigations, which are currently capable of providing data at much higher ${R_{\lambda}}$ compared to DNS EB2014 ; SHREK . A simple estimate suggests that $\eta/\eta_{ext}\gtrsim 3$ corresponding to these laboratory experiments at ${R_{\lambda}}\approx 6000-10000$. While this correction is unlikely to affect the dynamics in wind tunnel experiments EB2014 , it might enhance the quantum effects in liquid-He experiments at SHREK . However, more quantitative studies, even at relatively lower ${R_{\lambda}}$, would be useful, since currently resolving even $\eta$ in such high ${R_{\lambda}}$ experiments is an outstanding technical challenge.

VII Conclusions

Using very well-resolved DNS of isotropic turbulence, both in space and time, at Taylor-scale Reynolds number ${R_{\lambda}}$ ranging from 140 to 650, we have characterized the extreme fluctuations of the velocity gradients. In particular, we focused on the square of vorticity, $\Omega=\omega_{i}\omega_{i}$, and strain, $\Sigma=2s_{ij}s_{ij}$ (synonymous with enstrophy and dissipation), which have the same mean value (equal to $1/\tau_{K}^{2}$). Whereas the PDFs of $\Omega\tau_{K}^{2}$ and $\Sigma\tau_{K}^{2}$ superpose well around their mean values, the extents of their tails strongly grows as ${R_{\lambda}}$ increases. We find that these tails can be empirically collapsed by using a smaller time scale $\tau_{ext}$, defined as $\tau_{ext}=\tau_{K}\times{R_{\lambda}}^{-\beta}$, implying that the extreme velocity gradients in the flow grow as $\tau_{ext}^{-1}$. The numerical results indicate that $\beta\approx 0.775\pm 0.0025$.

The above result is further validated by analyzing the PDFs of velocity increments $\delta u_{r}$ at distances $r$ equal to or less than the Kolmogorov length scale $\eta$. Our results show that $\delta u_{r}$ can be as large as the velocity r.m.s. $u^{\prime}$, and slowly increases with ${R_{\lambda}}$. The excellent superposition of the rescaled PDFs of ${R_{\lambda}}^{-\alpha}\delta u_{r}/u^{\prime}$, with $\alpha=\beta-\nicefrac{{1}}{{2}}$ over the range of ${R_{\lambda}}$ covered in this study, suggests that the largest velocity gradients consist of velocity increments of approximately $u^{\prime}$ over a size $\eta_{ext}\sim\eta\ {R_{\lambda}}^{-\alpha}$. The existence of scales smaller than $\eta$ is consistent with previous ideas, although our results do not quantitatively support the existing phenomenological theories. In particular, the assumption that extreme events correspond to a local Reynolds number of unity does not appear to be justified in vortex tubes, where extreme gradients are found to reside – as revealed by flow visualizations in Figs. 1 and 7 and also consistent with previous studies Jimenez93 ; Ishihara09 .

Further analysis of flow structures around intense gradients in Fig. 7 reveals that vorticity and strain are generally not spatially colocated, as suggested in some earlier studies Donzis:08 ; Yeung15 . Conditional averaging shows that strain acting on intense vorticity, is on average weaker than the vorticity, and shows an approximate power law behavior given by Eq. (6). It is important to note that the structures shown in Fig. 7 taken in isolation do not lead to a strong vorticity amplification. This suggests that the stretching, necessary to create the very intense velocity gradients, in a representation of the velocity field in terms of a Biot-Savart equation PumSig87 , could originate from a non-local mechanism. The observation that the strain acting on intense vortices is significantly weaker than corresponding vorticity (as reflected in exponent $\gamma<1$ in Eq. (6)) can be viewed as a consequence of this non-locality. Using scaling analysis, we are able to quantitatively relate $\beta$ with the exponent $\gamma$. The weak increase in $\gamma$ with ${R_{\lambda}}$ suggests the same for $\beta$. This leaves open the possibility that $\beta$ could asymptote to $1$ (and $\alpha$ to $0.5$), in the limit of ${R_{\lambda}}\rightarrow\infty$ – a simple extrapolation of our data suggests that $\beta\geq 0.99$ would require ${R_{\lambda}}\gtrsim 20,000$. However, such high Reynolds numbers might not be feasible experimentally or numerically. To conclude, explaining the scalings discussed here, especially in the light of $\gamma$, remains an outstanding theoretical challenge. Much remains to be learned by analyzing well-resolved data from even higher Reynolds numbers than considered here, from both DNS and experiments.

Acknowledgments

We gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) for funding this project by providing computing time on the GCS supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC), where the simulations reported in this paper were primarily performed. This work was also partly supported by the advanced supercomputing resources provided to author P.K.Y. by the National Center for Supercomputing Applications (NCSA) at the University of Illinois at Urbana-Champaign and the Texas Advanced Computation Center (TACC) at the University of Texas at Austin which provided access to supercomputers Blue Waters and Stampede2 respectively. D.B., A.P. and E.B. acknowledge support from EuHIT–European High-performance Infrastructure in Turbulence, which is funded by the European Commission Framework Program 7 (Grant No. 312778). P.K.Y. was supported by National Science Foundation (NSF) Grants ACI-1036170 and 1640771 under the Petascale Resource Allocations program. We thank K.R. Sreenivasan for helpful comments on an earlier draft of the manuscript.