The classic description of the energy cascade in turbulence postulates that kinetic energy originates from forcing large-scale eddies, is subsequently transferred to smaller-scale eddies (forward cascade), and is eventually dissipated due to viscous e↵ects (Richardson 1922;Kolmogorov 1941). In a statistical sense, the sign and magnitude of third-order moments of velocity increments confirm this general direction of the energy cascade, as described by the 4/5 law governing the global average of the third-order longitudinal velocity increment (Kolmogorov 1941;Frisch 1995), h u L (`) 3 i ⌘ h([u(x + `) u(x)] • `/`) 3 i = 4 5 `h✏i, where h...i denotes global averaging, u L (`) is the longitudinal velocity increment and ✏ is the viscous dissipation rate, while the displacement `= |`| is assumed to be well inside the inertial range of turbulence. In this sense, the quantity 5 4 h u L (`) 3 i/ì s often interpreted as a measure of the energy flux going from scales larger than `to all smaller scales. Because turbulence is known to be highly intermittent in space and time (Kolmogorov 1962;Frisch 1995;Meneveau & Sreenivasan 1991) there has also been much interest in characterizing the local properties of the energy cascade, i.e. the fluctuations of the energy flux before averaging. However, without statistical averaging, the 4/5-law is less meaningful, e.g., the quantity 5 4 u 3 L /`cannot simply be interpreted as an energy flux locally in space and time. To enable such interpretation, it is necessary to consider explicit angular averaging over all possible directions of the vector `. Such formulations have been developed in prior works by Duchon & Robert (2000), Eyink (2002) and Hill (2001Hill ( , 2002a)). Duchon & Robert (2000) and Eyink (2002) use such equations to study the energy cascade and energy dissipation in the limit of zero viscosity. A review about extensions to the classic Kolmogorov equation is presented by Dubrulle (2019), specifically focusing on the Duchon & Robert (2000) local formulation. Hill (2001Hill ( , 2002a) ) developed a local version of the Kolmogorov equation in which the reference position x is symmetrically located halfway between the two points x + r/2 and x r/2 separated by r over which the velocity increment is computed. This equation, which we shall denote as the Kolmogorov-Hill (KH) equation (sometimes also called Karman-Howarth-Monin-Hill (Danaila et al. 2012;Yasuda & Vassilicos 2018) or Generalized Kolmogorov (Marati et al. 2004) equation), describes the evolution of the second-order (squared) velocity di↵erence, a measure of energy content of all scales smaller than |r| at a specific physical position x. As will be reviewed in §2, scale-space integration over r of the KH equation up to some scale `in the inertial range and without additional statistical averaging provides a localized description of the energy cascade process (Hill 2002b;Yasuda & Vassilicos 2018). The KH equation also includes e↵ects of viscous dissipation, viscous di↵usion, advection, and pressure. A number of prior works have studied various versions of the KH equation. For isotropic turbulence, Yasuda & Vassilicos (2018) quantified the variability of the energy flux that arises in this equation, while Carbone & Bragg (2020) considered a definition of mean energy flux approximated based on a solenoidal filtered velocity increments and examined its connections to average vortex and strain stretching rates. Besides applications to isotropic homogeneous flow, numerous studies have investigated the application of the statistically averaged KH equation to spatially non-homogeneous flows. For instance, in wall bounded flows, researchers have explored the energy cascade using a Reynolds decomposition to isolate e↵ects of mean shear and non-homogeneity (Antonia et al. 2000;Danaila et al. 2001Danaila et al. , 2004Danaila et al. , 2012;;Marati et al. 2004;Cimarelli et al. 2013). Investigations have also studied the energy cascade rates in boundary layer bypass transition (Yao et al. 2022), and flow separation (Mollicone et al. 2018). Furthermore, specific attention has been given to the study of inverse cascade in wake flows (Gomes-Fernandes et al. 2015;Portela et al. 2017) and at turbulent/non-turbulent interfaces (Zhou & Vassilicos 2020;Cimarelli et al. 2021;Yao & Papadakis 2023).
The notion of transfer, or flux, of kinetic energy across length-scales is of particular practical interest also in the context of large eddy simulation (LES). There the rate of energy cascade is commonly referred to as the subgrid or subfilter-scale (SGS, SFS) rate of dissipation. It is defined as the contraction between the subgrid stress tensor and the filtered strain-rate tensor and arises as a source term in the transport equation for subgrid/subfilter-scale kinetic energy (Piomelli et al. 1991;Meneveau & Katz 2000). This quantity characterizes the energy transfers between the resolved scale and the residual scale within the inertial range, which is also a local property (Eyink & Aluie 2009). The SGS dissipation is highly intermittent (Cerutti & Meneveau 1998), and can be both positive and negative locally, but on average, energy is known to be transferred from large scales to the residual scales (forward cascade). There is considerable literature on the subject starting from the seminal papers by Lilly (1967), Leonard (1975) and Piomelli et al. (1991). Some reviews include Meneveau & Katz (2000); Meneveau (2010); Moser et al. (2021).
Without averaging, it has been a common observation that the SGS/SFS dissipation can be negative which has often been interpreted as indicative of local inverse cascading of kinetic energy, i.e., energy transfer from small to large scales of motion ("backscatter" (Piomelli et al. 1991)). Borue & Orszag (1998) noted that the forward cascade occurs predominantly in regions characterized by strong straining, where the magnitude of negative skewness of the strain tensor and vortex stretching are large. Conversely, backscatter was observed in regions with strong rotation. The relationship between SGS dissipation and stress topology and stress strain alignment geometry was discussed and measured based on 3D PIV measurements by Tao et al. (2002). In a more recent study, Ballouz & Ouellette (2018) investigated the SGS tensor by considering the relative alignment of the filtered shear stress and strain tensors. They found that the energy cascade e ciency is quite low, a trend they attributed to energy being transferred largely between positions in physical space. Quantitatively, in expressing the subgridstress tensor as a superposition of all smaller scale Gaussian-filtered velocity gradients, Johnson (2020Johnson ( , 2021) ) was able to isolate the relative contributions of small-scale strain self-stretching and vortex stretching, finding both to be important.
It has been questioned whether it is the local quantity ⌧ ij Sij (where ⌧ ij and Sij are the subgrid-scale stress and resolved strain-rate tensors, respectively), or the work done by the SGS/SFS force, ũi @ j ⌧ ij (where ũi is the resolved velocity), that should be the genuine definition of local energy cascade rate. For instance Kerr et al. (1996) use the latter in their study of correlations of cascade rate and vorticity, and more recently Vela-Martín & Jiménez (2021) use both quantities in their analysis. Moreover the SGS force plays a central role for optimal LES modeling (Langford & Moser 1999). The SGS force is invariant to divergence-free tensor fields which therefore do not a↵ect the large-scale dynamics but addition of such a tensor field to ⌧ ij can certainly a↵ect the usual definition of subgrid-scale dissipation ⌧ ij Sij . By re-expressing the SGS stress and dissipation terms using an optimization procedure, Vela-Martín (2022) provided arguments that the often observed backscatter does not actually contribute to the energy cascade between scales but rather to the energy flux in the physical space, also suggesting that backscatter does not need to be explicitly modeled in LES.
As can be seen from this partial summary of the literature on backscatter and inverse cascade in the LES filtering approach, no consensus has been reached regarding the possible importance and physical interpretation of local backscatter using the definition based on the inner product of the subgrid stress and filtered strain-rate tensors. Also, the question of inverse cascade has not received much attention from the point of view of the local versions of the Kolmogorov equation in the structure function approach. Therefore, in this paper we first revisit the generalized local structure function formulation ( §2.1). We argue that in this formulation the term responsible for the energy cascade can be unambiguously interpreted as a flux of kinetic energy between scales since it appears inside a divergence in scale space. In this sense it di↵ers from the filtering formulation used in LES (reviewed in §2.2) in which typically a fixed filter scale is used and no change in scales is considered, thus making the concept of a "flux in scale space" less clearly defined and open to various interpretations.
With the definition of local cascade rate or energy flux clarified for the structure function approach, we perform a comparative study of both the structure function and the filtering approaches' energy flux terms in a relatively high Reynolds number DNS database of forced isotropic turbulence at a Taylor-scale Reynolds number of 1,250. The data analysis is greatly facilitated by the availability of these data in a new version of the Johns Hopkins Turbulence Database (JHTDB) System, in which python notebooks access the data directly (see appendix A). The comparisons involve various statistical properties of the energy flux. First, in §3 we provide comparisons of both quantities by means of simple statistical measures such as their mean values, joint probability density distributions and correlation coe cients, comparing both the two definitions of kinetic energy and kinetic energy cascade rate or flux. We then comparatively examine conditional averages based on the local molecular dissipation rate averaged over a ball of size `, specifically re-examining the Kolmogorov refined similarity hypothesis (KRSH) in §4. Then, in §5, we present comparative conditional averages of kinetic energy flux based on properties of the large-scale velocity gradient field such as the strain and rotation magnitudes, and the Q and R invariants. Particular attention is placed on events of local negative energy flux and whether or not such events can be considered to be of statistical significance. Overall conclusions are presented in §6.
In this section, both the structure function based (KH equation) and filtering (LES) energy equations are reviewed. We focus on the term representing energy cascade (energy flux) in each equation, and describe some of the prior e↵orts in the literature relating the structure function and filtering approaches.
The KH equation is a generalized Karman-Howarth equation that is directly derived from the incompressible Navier-Stokes equations without any modeling. Before averaging, the instantaneous KH equation with no mean flow and neglecting the forcing term reads (Hill 2001(Hill , 2002b)):
where u i = u i (x, r) = u + i u i is the velocity increment vector in the i th Cartesian direction over displacement vector r. The superscripts + and represent two points
x + r/2 and x r/2 in the physical domain that have a separation vector r i = x + i x i and middle point x i = (x + i + x i )/2 (see Fig. 1 (a)). The superscript ⇤ denotes the average value between two points, e.g., the two-point average dissipation is defined as ✏ ⇤ (x, r) = (✏ + + ✏ )/2, and ✏ ± here is the "pseudo-dissipation" defined at every point as ✏ = ⌫(@u i /@x j ) 2 (in Hill (2002b), an alternate expression involving the real dissipation was introduced (his Eq. 2.13) at the cost of including an additional pressure term). Note that throughout this paper, when referring to "dissipation" we will mean the pseudodissipation. Also, we will use r s = r/2 to denote the radial coordinate vector from the local "origin" x.
As remarked by Hill (2001Hill ( , 2002b) ) it is then instructive to apply integration over a sphere in r s -space up to a radius `/2, i.e. over a sphere of diameter `. The resulting equation is divided by the sphere volume V `= 4 3 ⇡(`/2) 3 and a factor 4, and Gauss' theorem is used for the r-divergence terms (recalling that @r = 2@r s ), yielding
(2.2) where S `represents the bounding sphere's surface of area S `= 4⇡(`/2) 2 and nj is the unit normal vector. Eq. 2.2 suggests defining a structure-function based kinetic energy at scale `according to
so that the first term in Eq. 2.2 corresponds to @ k sf,`/ @t. The 1/2 factor in front of the integral is justified since the volume integration over the entire sphere will double count the energy contained in u 2 i = (u + i u i ) 2 . Equation 2.2 thus describes the transport of two-point, structure function energy k sf,`, which represents energy within eddies with length scales up to ` (Davidson 2015) in both the length scale `and physical position x spaces. The last term in equation 2.2 represents r-averaged rate of dissipation with the radius vector r s = r/2 being integrated up to magnitude `/2,
(2.4)
As remarked by Hill (2001Hill ( , 2002b) ) this quantity corresponds directly to the spherical average of local dissipation at scale `and plays a central role in the celebrated Kolmogorov Refined Similarity Hypothesis, KRSH (Kolmogorov 1962).
The local energy cascade rate in the inertial range at position x and time t is defined as
where [..] S `indicates area averaging over the sphere of diameter `. We note that in this definition, `(x, t) represents the surface average of a flux that is defined positive if energy is flowing into the sphere in the r-scale space. The position is fixed at x and thus the quantity `(x, t) does not contain possible confounding spatial transport e↵ects. In terms of the overall average of Eq. 2.2, under the assumptions of homogeneous isotropic flow and statistical steady-state conditions, and for `in the inertial range of shows integration up to a ball of radius `in which pairs of points separated by distances r up to `are used as in the approach of Duchon & Robert (2000). For volume averaging, in (a) 3D integration over the vector rs is performed at fixed x while in (a) 3D integration over the vector r is performed at fixed x. For surface integrations, in (a) integration is done over the spherical surface of radius `/2 while in (b) is is done over a spherical surface of radius `.
turbulence, the unsteady, transport and viscous terms vanish. The pressure term is also zero due to isotropy and incompressibility. Therefore, Eq. 2.2 can be simplified and yields as expected
or equivalently [ u 2 i u j nj ] S `= 4/3 `h✏i, the 4/3-law Frisch (1995). In this paper the focus will be mainly on the flux term `with some attention also on the dissipation term ✏ `. Analysis of the time derivative, spatial advection terms and pressure terms is left for other ongoing studies. The viscous flux terms (in both spatial and scale spaces) are also not considered, since our present interest concerns the inertial range.
In this section, we review the transport equation of the subgrid-scale kinetic energy (Germano 1992) for k sgs,`⌘ 1 2 ⌧ ii , where ⌧ ij = g u i u j ũi ũj is the subgrid-scale stress tensor, where the tilde symbol (⇠) denotes spatial filtering of variables. The transport equation for k sgs,`r eads (Germano 1992)
(2.7) The last term is called the subgrid-scale rate of dissipation at position (x), and is often denoted as
(2.8)
For filtering, in the present work we consider a spherical-shaped sharp top-hat filter in physical space with a diameter equal to `. Therefore, for any field variable A(x), we define the filtered variable as Ã
r s . Note that each term in Equation 2.2 and in Equation 2.7 are thus evaluated at the same length scale. Terms in Equation 2.7 can be compared directly to terms in Equation 2.2, in particular the local dissipation terms are exactly the same, i.e,
(2.9)
Again, for homogeneous steady state turbulence in the inertial range (neglecting viscous di↵usion and resolved dissipation terms), upon averaging Eq. 2.7 simplifies to
which is similar to Equation 2.6 and thus on average certainly both definitions of energy cascade rate/flux agree with each other, i.e., h⇧
It is also of interest to compare the average value of the two definitions of kinetic energy used in both definitions of energy cascade rate/flux. In the inertial range of high Reynolds number turbulence, both hk sf,`i and hk sgs,`i can be evaluated based on the Kolmogorov r 2/3 law and k 5/3 spectrum, respectively. The result is (see Appendix B for details), hk sf,`i ⇡ 1.6 h✏i 2/3 `2/3 and hk sgs,`i ⇡ 1.2 h✏i 2/3 `2/3 . In other words, they are of similar order of magnitude but the SGS kinetic energy is slightly smaller.
In the present paper, we shall perform the data analysis and comparisons using the two approaches mentioned above (scale-integrated local KH and filtering formulations). However, it is useful at this stage to include some remarks regarding other structure function and energy definitions used in earlier works by Vreman et al. (1994), Constantin et al. (1994), Duchon & Robert (2000), Eyink (2002) and Dubrulle (2019). Those approaches typically focus on the structure function written at one of the endpoints instead of the midpoint. Duchon & Robert (2000) and Dubrulle (2019) focus on the two-point correlation quantity C(x, r) = u i (x)u i (x + r) (see Fig. 1 (b)). Local averaging over all values of r from r = 0 up to scale |r| = `at any given x then corresponds to the "mixed" energy quantity u i ũi /2 (denoted as E `in Dubrulle ( 2019)), and where the filtering is over a sphere of diameter 2`so as to combine two points with separation distances up to `. The quantity C(x, r) combines filtered and unfiltered velocities and hence it is more di cult to interpret for comparisons of structure function and LES filtering approaches. In its transport equation, Duchon & Robert (2000) show that a term similar to the third-order structure function term of Eq. 2.5 arises. However, in order for the structure function to correspond to scale `, one has to choose to integrate over a sphere of diameter 2`(the locally integrated dissipation rate would then be ✏ 2`) . In a spherical integration over r of powers of the velocity di↵erence [u i (x + r) u i (x)], only the first term is a↵ected by filtering or averaging over the spherical shell, while the center velocity u i (x) remains fully local. Note that in the scale-integrated local KH equation, the averaging a↵ects both end-point velocities in the same way, and both become averaged at scale `in a formally symmetric way.
An early connection between structure functions and filtering approaches was developed by Vreman et al. (1994). In the Vreman analysis, the structure function is defined based on the di↵erence of velocity u i (x + r) and the locally filtered velocity ũi centered at x. Spherical integration of (u i (x + r) ũi ) 2 over a sphere of radius `/2 then yields equivalence with the SGS kinetic energy at scale `. But (u i (x + r) ũi ) 2 does not equal the usual structure function definition, now due to a mixture of filtered and unfiltered quantities at two points even before local filtering.
Another interesting approach was presented in Constantin et al. (1994) and connected to the LES filtering approach by Eyink (1995), Eyink (2006) (equations 2.12-2.14). In fact as recounted in the review by Eyink & Sreenivasan (2006), early unpublished work by Onsager anticipated such expressions half a century prior. Written in terms of the sharp spherical filter we use here, the expression for the trace of the SGS stress reads
This equation represents an exact relationship between two-point structure functions and the subgrid-scale kinetic energy. But for the RHS to correspond to structure functions up to scale `, the integration must be done over a sphere of radius `and thus a filtering scale of 2`for the stress tensor in the filtering formulation. The suggested relationship then appears to be between subgrid-scale stress kinetic energy at scale 2`and structure functions up to two-point separations `but averaged over a local domain of size 2`, similarly as in the Duchon & Robert (2000) approach. Note that while each of the terms in Eq. 2.11 is also a mixture of filtered and unfiltered velocities, the subtraction cancels the local term and restores the fully filtered property inherent in the definition of ⌧ ii .
While not expecting qualitatively di↵erent results (except perhaps using the diameter instead of the radius as a name for "scale"), we here continue our focus on the more "symmetric" formulation by Hill, with fixed position x specified at the midpoint between two points separated by vector r whose magnitude then spans up to scale `(or integration radius r s up to radius `/2).
In this section, we provide comparisons of local kinetic energies in the structure function formalism, k sf,`, with that in the filtering formalism, k sgs,`. We also compare the local energy cascade rates `and ⇧ `. We consider data from a direct numerical simulation (DNS) of forced isotropic turbulence at R = 1,250 (the Taylor-scale Reynolds number) that used 8,192 3 grid points (Yeung et al. 2012) in a computational domain of size equal to (2⇡) 3 . The integral scale of the flow is L = 1.24, the velocity root-mean-square is u 0 = 1.58 and the mean dissipation is h✏i = 1.36. More details about the data and simulation parameters are available as supplementary material. The analysis is performed at four length-scales in the inertial range, `= {30, 45, 60, 75}⌘ where ⌘ = (⌫ 3 /h✏i) 1/4 is the Kolmogorov length scale, the value of which is 4.98 ⇥ 10 4 . Comparing to the transverse Taylor microscale (Pope 2000) g = u 0 (15⌫/2h✏i) 1/2 ⇡ 0.024, the four length scales are `= {0.62, 0.93, 1.24, 1.55} g , respectively.
To compute volume spherically filtered quantities such as k sgs,`a nd ⌧ ij (and filtered velocity gradient tensor to be discussed in §5), we fix the middle point coordinate x in the physical domain. Subsequently, we download data in a cubic domain using the JHTDB's cutout service in a cube of size equal to `3. The data are then multiplied by a spherical mask (filter) to evaluate local filtered quantities. Other quantities are obtained by utilizing pre-computed Getfunctions from the Johns Hopkins Turbulence database (JHTDB) including spatial interpolation and di↵erentiation, as explained in more detail in Appendix A. For surface averages such as `, we discretize the outer surface of diameter `into 500 points (for the largest `/⌘ = 75 case, 2000 points are used) that are approximately uniformly distributed on the sphere. The accuracy of this method of integration has been tested for the `/⌘ = 45 case by comparing the results from using 500 points to those using 2000 points, for a smaller testing subsample of 500 randomly chosen spheres. We verified that the di↵erence between the mean values of as well as the average of the absolute value of di↵erences were less that 1%. For volume averages such as ✏ `and k sf,`, we use 5 shells for `/⌘ = 30, 45, 60. The outermost shell comprises 500 uniformly distributed points, with a reduction in number of points towards the inner shells approximately maintaining the density. We tested 500 randomly chosen spheres to calculate ✏ `at `/⌘ = 45 using 5 shells and 10 shells. The di↵erence between the mean values of ✏ `as well as the average of the absolute value of di↵erences were less that 2%. For the larger length scale `/⌘ = 75, the number of shells was increased to 6; the accuracy is tested using the same method as employed for `/⌘ = 45. For all the calculations, data on the specified points are obtained from the database using 8th-order Lagrange spatial interpolation. We tested di↵erent spatial interpolation methods even without interpolation (using the closest grid point values), verifying that the averaged values of interest were essentially unchanged.
Overall mean values are obtained at the four scales and are plotted in Fig. 2 (a). The results for kinetic energy for the structure function approach are consistent with the analytical evaluation (see Appendix B). For the SGS kinetic energy, the numerical results fall below the theoretical inertial range prediction, due to the transfer function of top-hat filtering having a very di↵erent spectral signature compared to the structure function, and it emphasizes more the viscous range when integrating than the structure function operation, reducing the amount of SGS kinetic energy even at scales much larger than the Kolmogorov scale (see discussion in Appendix B).
Figure 3(a) shows the joint PDF of k sf,`a nd k sgs,`a t scale `= 45⌘. The correlation coe cient between both quantities is ⇢ kk = 0.97 (Fig. 2 (b)). The correlation coe cient is defined as ⇢ xy = h(x hxi)(y hyi)i/( x y ) where represents the variable's root mean square value. Similarly, Figure 3(b) shows the joint PDF of ⇧ `and `, also at scale `= 45⌘ for the same dataset. The correlation coe cient between both quantities is measured to be ⇢ ⇧ = 0.58 (Fig. 2 (b)), significantly lower than for the energies but still appreciable. It can be seen that negative values occur for both ⇧ `and `, although it appears that `has more variability and larger negative excursions than ⇧ `. As summarized in §1, the relevance of locally negative values of ⇧ `to the flow physics remains unclear, especially given the fact that upon averaging, the quantity becomes positive. Conversely, the quantity `has a clearer local interpretation, in the sense that locally negative values can clearly be interpreted as kinetic energy (local u 2 i /2)) showing a net flux out of a sphere of diameter `in scale space, i.e. becoming associated with energy at larger `, while its overall average is positive. An interesting question is whether negative values of ⇧ `or `survive under some type of statistical averaging. In the next sections we use conditional averaging to quantify the importance of negative values (inverse local cascade, or backscatter).
Motivated by KRSH and the fact that local viscous dissipation (small-scale) appears in both scale-integrated local KH equation and subgrid-scale kinetic energy equation, i.e., Eqn. 2.4 is identical to Eqn. 2.9, in this section, we compare conditionally averaged cascade rates/fluxes for both the structure function and filtering formulations, conditioned on ✏ `, i.e, h `|✏ `i and h⇧ `|✏ `i. According to KRSH (Kolmogorov 1962), the statistical properties of velocity increments depend on the local average dissipation within a sphere of scale `, rather than being determined by the globally averaged dissipation. Written in terms of the quantities of present interest, the KRSH would read h `|✏ `i = ✏ `since `is fully determined by the velocity increments envisioned in the KRSH. Loosely extending the KRSH arguments to the filtering formalism would suggest h⇧ `|✏ `i = ✏ `.
In order to assess this hypothesis, we evaluate the conditional averages based on the same dataset described before. Results for h `|✏ `i and h⇧ `|✏ `i are shown in Figure 4(a). Results for the four scales considered are included. As can be seen, the plot shows a close agreement between both h `|✏ `i and h⇧ `|✏ `i, with ✏ `. It is important to note that ànd ⇧ `are conditioned on the exact same values of ✏ `. The similarities and di↵erences observed in Figure 4 indicate that `and ⇧ `share many properties (same conditional averages) but they are not identical. For instance, it is clear from Fig. 3 (b) that the variance of `exceeds that of ⇧ `, even though their mean values are the same.
In general, the behaviors of both h `|✏ `i and h⇧ `|✏ `i confirm the validity of the KRSH in the present context. More detailed analysis of the KRSH for `and connections to Eq. 2.2 are reported in Yao et al. (2023). We also tested KRSH using the full viscous dissipation ⌫(@u i /@x j )(@u i /@x j +@u j /@x i ) instead of the pseudo-dissipation ⌫(@u i /@x j ) 2 when computing ✏ `. The largest di↵erence for h `|✏ `i is less than 1%. Additionally, the correlation coe cient between the two types of dissipation is 0.996.
Similarly, we evaluate the kinetic energies k sf and k sgs conditionally averaged on ✏ `. Results for k sf and k sgs are essentially indistinguishable, except for a constant o↵set consistent with the ratio of their mean values. In terms of their dependence on dissipation, we observe power-law scaling ⇠ ✏ `with ⇠ 0.79, slightly larger than the value 2/3 implied by standard Kolmogorov scaling. To verify the present data and analysis methods, we also evaluate the traditional longitudinal second order structure function conditioned on ✏ `, h u 2 L |✏ `i (where u L = u j nj ), for which the Kolmogorov scaling ⇠ ✏ 2/3 `according to the RKSH is well established (Stolovitzky et al. 1992). The result (shown as stars in Fig. 4(b)) indeed confirms the 2/3 scaling for this quantity. A more in-depth analysis and possible reasons for non-Kolmogorov scaling of k sf with ✏ `is left for future studies. At this stage, we simply note the similarity in scaling and overall behavior of k sf and k sgs .
In this section, motivated by large-scale properties of the flow that would be available in LES, we explore correlations between properties of the velocity gradient tensor filtered at scale `and the two definitions of energy cascade rate/flux. It is useful to cast the present comparative study of and ⇧ using analyses of the type that have been performed before in the context of LES. The velocity gradient tensor encapsulates information about fluid deformation and rotation and connections to the energy cascade have been studied extensively. Already Bardina et al. (1985) examined the impact of rotation on homogeneous isotropic turbulence (HIT) and observed that rotation decreases the dissipation (cascade) rate while increasing the length scales, suggestive of inverse energy cascade e↵ects. Goto (2008) investigated physical mechanisms underlying forward energy cascade and argued that forward cascade can be triggered in regions characterized by strong strain between two large-scale tubular vortices. The role of the filtered gradient tensor for energy cascade was first explored numerically in Borue & Orszag (1998) and experimentally in Van der Bos et al. (2002) building on the "Clark model" that approximates features of the subgrid-scale tensor using Taylor-series expansion. Recent studies by Johnson (2020Johnson ( , 2021) ) and Carbone & Bragg (2020) have significantly expanded on such analyses and examined the roles of self-strain amplification and vortex stretching driving the forward energy cascade process. For inverse cascade, a vortex thinning mechanism may be at play (Johnson 2021).
A first level of characterization of the properties of the velocity gradient tensor are its invariants. To characterize deformation and rotation, we evaluate the strain and rotation invariants from data, defined according to
where S ij and ⌦ ij are the symmetric and antisymmetric parts of the velocity gradient tensor A ij = @u i /@x j and the tilde denotes, as before, spherical tophat filtering over a ball of diameter `. For consistency with prior literature, these values will be normalized by the overall average hQ w i = 1 2 h⌦ 2 `i (equal to 1 2 hS 2 `i in homogeneous turbulence). A more detailed characterization of the statistics of velocity gradients involves the invariants Q and R (Vieillefosse 1982). It is well-known that the joint probability density function (JPDF) of Q-R exhibits a characteristic tear-drop shape (Chong et al. 1990;Meneveau 2011), from which flow topology information such as vortex stretching and compression can be inferred (Chong et al. 1990;Borue & Orszag 1998;Lüthi et al. 2009;Danish & Meneveau 2018). These two invariants (at scale `) are defined as usual according to
Under the assumption of restricted Euler dynamics (Meneveau 2011), the transport equation for the velocity gradient tensor leads to dQ `/dt = 3R `and dR
). The quantity R `can thus be considered as the (negative) rate of change of Q `and contains both vortex stretching and strain self-stretching mechanisms (Johnson 2021). In our comparative investigation of energy cascade rates, conditional averaging based on the four invariant quantities S 2 `, ⌦ 2 `, Q `and R `will be undertaken. We begin with qualitative visualizations of the fields in small subsets of the domains analyzed. Panel (a) and (b) of Figure 5 depict a sample instantaneous field of `and ⇧ `respectively, highlighting regions of strong local forward cascade (indicated by solid red circles) and strong inverse cascade (indicated by dashed circles). The correlation between these two variables is evident, the computed correlation coe cient between the snapshots is 0.64. On both panels (a) and (b) the fluxes are normalized by the global averaged dissipation h✏i. As already noted based on the joint PDFs, there are di↵erences between `and ⇧ `. The maximum magnitude of the positive cascade rate in `is about twice that of ⇧ `, while the magnitude of the negative cascade rate in `is about 3 to 4 times larger. Since h `i ⇠ h⇧ `i ⇠ h✏ `i ⇠ h✏i, the significant di↵erent maximum values indicates `is more variable and intermittent than ⇧ `. Also, `exhibits somewhat finer-scale spatial features. (1998) and recently Johnson (2021); Carbone & Bragg (2020)). We focus attention on the regions with negative energy cascade rates. Conditional averaging can elucidate the statistical significance of these regions. Specifically, we inquire whether there are largescale flow local features as characterized by the filtered velocity gradient invariants that are systematically accompanied by inverse cascade, i.e., negative `. Thus we perform conditional averaging of `based on the invariants S 2 `and ⌦ 2 `and repeat the analysis for the SGS energy flux quantity ⇧ `.
Figure 6 shows the joint conditionally-averaged `and ⇧ `based on S 2 `and ⌦ 2 `, denoted as h `|S 2 `, ⌦ 2 `i and h⇧ `|S 2 `, ⌦ 2 `i, respectively. The analysis is performed by computing averages over two million randomly distributed points x. In the presented results, ìs normalized by h✏i, while S 2 `and ⌦ 2 `are normalized by hQ w i = 1 2 h⌦ 2 `i. Panels (a), (b), (c) and (d) of figure 6 present the joint conditionally-averaged h `|S 2 `, ⌦ 2 `i at four di↵erent length scales, namely `= {30, 45, 60, 75}⌘, highlighting the dominance of the forward cascade by the extensive red region. This magnitude is many times larger than the maximum magnitude observed in the blue region, representing the inverse cascade. The red region covers a wide range of S 2 `and ⌦ 2 `values, consistent with the expectation that the global average would favor a forward cascade (h `> 0i). The highest positive values of h `|S 2 `, ⌦ 2 `i correspond to high strain rates and low rotation rates, and they decrease as the strain rate decreases. Interestingly, the inverse cascade appears explicitly in the lower-right corner of the plots, where the rotation rate is strong but the strain rate is weak. It is worth noting that the conditionally averaged values shown in Figure 6 reflect the combined outcome of the forward and inverse cascades. Consequently, in specific regions characterized by distinct strain and rotation rates, events with forward and inverse cascades can cancel each other out. Only in the lower right corner is there an indication of net inverse cascade when the cascade rate is defined using the structure function approach.
In panel (d), we superimpose dashed lines representing the isolines of Q `, with the Q `= 0 line indicating the condition of equal strain and rotation rates. The parallel dashed lines correspond to Q `= 10, 5, 0, 5, 10, and 15, respectively. The Q `= 15 line appears near the boundary separating the red and blue regions. However, the boundary of the blue region does not appear to align well with the Q `isoline. This observation suggests that Q `might be not enough to provide an adequate threshold for distinguishing the net forward and inverse cascade regions. Panels (d), (e), (f), and (g) in Figure 6 present results for the joint conditionallyaveraged ⇧ `based on S 2 `and ⌦ 2 `, corresponding to the same filter scales as panels (a), (b), (c), and (d). It is evident that trends for the positive cascade rate (red region) for ⇧ `closely resemble those of `, with the peak of the forward cascade occurring at a high strain rate and low rotation rate. The magnitude of the maximum forward cascade rate for ⇧ `is slightly weaker compared to that of `. The most significant di↵erence is that only a few instances of blue squares are observed in regions characterized by strong rotation and weak strain, indicating that the overall predominance of the forward cascade persists regardless of the local values of S 2 `and ⌦ 2 `. These results highlight some important statistical di↵erences between `and ⇧ `. We further evaluate the conditional averages of the energy fluxes `and ⇧ `, conditioned on either S 2 `and ⌦ 2 `individually.
This analysis is motivated by the work of Buaria & Pumir (2022), which highlighted di↵erent scalings of conditional averages with respect to strain rate compared to rotation rates even though one would expect similar results based on dimensional arguments. In figure 7, panel (a), we demonstrate that both `and ⇧ `exhibit a power-law relationship when conditioned on S 2 `, with a slope of 3/2. This scaling aligns with dimensional analysis and Kolmogorov scaling, h `|S 2 `i ⇠ [S 2 `]3/2 and h⇧ `|S 2 `i ⇠ [S 2 `]3/2 . In contrast, when conditioning `and ⇧ `on ⌦ 2 `(panel (b)), a di↵erent trend emerges, with much weaker dependence on rotation rate.
To develop a more detailed understanding of the inverse cascade region within h `|S 2 `, ⌦ 2 `i, we perform a further analysis by dividing the samples based on `> 0 and `< 0 for `= 45⌘. We then calculate the conditional average of these separated samples considering S 2 `and ⌦ 2 `, denoted as h `|S 2 `, ⌦ 2 `, `> 0i and h `|S 2 `, ⌦ 2 `, `< 0i. The results are presented in Figure 8. From panel (a), it can be observed that the forward cascade clearly increases with S 2 `, with the highest values of `concentrated in the upper-left corner. It increases also with ⌦ 2
`but less rapidly. Combined, the trend seems to be an increase roughly proportional to ⇠ S 2 `+ 0.75 ⌦ 2 `. Di↵erently, panel (b) illustrates that the inverse cascade is roughly proportional to ⇠ S 2 `+ 0.5 ⌦ 2 `, i.e. shallower isolines extending more in the horizontal direction than in the vertical compared to the forward cascade case shown in (a). This observation elucidates why the strongest red region in panel (b) of Figure 6 emerges at the largest S 2 `, while below this threshold, the forward cascade events progressively weaken and are gradually canceled out by the inverse cascade. Finally, in regions characterized by a weak strain rate and strong rotation rate, the inverse cascade becomes the dominant contribution.
Panel (c) and (d) display the distribution of the number of samples corresponding to positive and negative cascade rates in logarithmic scale (out of the 2 million samples (balls) considered). Our focus is specifically directed towards the bottom-right corner of the plots, which corresponds to the region where the inverse cascade is observed in panel (b) of Figure 6. Interestingly, we observe that at ⌦ 2 `/hQ w i ⇡ 40 and S 2 `/hQ w i < 10, the number of samples representing both inverse and forward cascade rates is roughly equivalent, falling within the range of 10 1 to 10 1.5 . This implies that within this region, the magnitude of the inverse cascade must be significant to achieve net negative values for the conditional average. Still, for the conditions with net inverse cascade, the number of occurrences for both forward and inverse cascade rates is quite small, on the order of only 1/10 5 of the total samples, indicating a very low frequency. We point out that in the extreme bins the conditionally averaged fluxes are unlikely to be fully converged. Therefore our observations are meant to be mostly qualitative in these regions. In particular, the inverse cascade regions depicted in Figure 6 are primarily attributed to rare but intense events. In the following subsection, we will show that inverse cascade can be better characterized by conditioning on Q `and R `invariants.
In the context of the ⌦ 2 `and S 2 `map shown in Figure 6, we observe the presence of a distinct inverse energy cascade in the region characterized by strong rotation but weak strain (corresponding to large Q values) for `. However, such observations did not hold for ⇧ `. But these results do not preclude the possibility that net forward and inverse cascade may be associated with other invariants of the filtered velocity gradient tensor Ãij .
Figure 9 shows the joint conditional averaged h `|Q `, R `i and h⇧ `|Q `, R `i at four di↵erent scales, namely, `= {30, 45, 60, 75}⌘. Across all figures, we can observe the distinctive teardrop shape pattern on the Q-R map, as reported in previous studies (Chong et al. 1990;Meneveau 2011). The black solid lines are the boundaries, separating the four quadrants based on the sign of Q and R. Notably, it becomes evident that both and ⇧ exhibit a strong and dominant inverse cascade in the quadrant characterized by Q > 0 and R > 0. It is useful to recall that the variable R `is associated to the rate of change of Q `(in fact assuming restricted Euler dynamics they are related by dQ `/dt = 3R ` (Cantwell 1992;Meneveau 2011). Therefore R `> 0 corresponds to a decreasing trend of Q `, i.e. vortex compression. We find that such events are associated to inverse energy cascade.. This observation is particularly interesting considering the absence of an observable inverse cascade for ⇧ `in the S 2 `and ⌦ 2 `map. Hence, these results indicate that the variables Q and R provide a more e↵ective characterization to identify inverse cascade using conditional averaging.
The region characterized by Q < 0 and R > 0, which corresponds to the straindominated region, exhibits the most pronounced forward cascade. This observation aligns with the findings of many prior analyses in the literature Borue & Orszag (1998) 2) . The data and editable analysis code that generated these joint PDFs (for the case at `= 45⌘) can be found at: https://cocalc.com/.../Figure9.
in Figure 9, further emphasizing that the strong local strain rate plays a crucial role in driving the forward energy cascade. We note that Borue & Orszag (1998); Van der Bos et al. (2002) display conditional averages weighted by the joint PDF of R and Q. In their results, there was hardly any indication of backscatter/inverse cascade in the Q > 0 and R > 0 "vortex compression" region, because the overall probability density of that region is smaller than the other regions. However, the unweighted conditional averaging represents the relevant values if the large-scale flow is in that particular state (Q > 0 and R > 0), and is therefore relevant to our analysis. We also notice a small blue region at Q `< 0, R `< 0 quadrant of h `|Q `, R `i (but not seen for h⇧ `|Q `, R `i). The occurrence of inverse cascade in strain-dominated, strain self-stretching regions appears intriguing. However the small number of samples ⇠ O(10) in the bin showing inverse cascade precludes us from ascribing much significance to this observation for now."
In a similar manner to Figure 8, we perform further conditional averaging also distinguishing positive and negative cascade rates. Figure 10 (a) and (b) present h `|Q `, R `, `> 0i and h `|Q `, R `, `< 0i at `= 45⌘. In the case of the inverse cascade, it is observed to occur in all four quadrants (panel (b)), with a more evenly distributed and symmetric presence in the upper two quadrants associated with Q > 0, i.e, the rotation-dominated regions. The characteristic teardrop shape is less prominent and exhibits a shorter tail compared to the forward cascade (panel (a)). Regarding the forward cascade, it is evident that it is most dominant in the Q < 0, R > 0 quadrant, consistent with Figure 9. However, in the Q > 0, R > 0 quadrant, the forward cascade is weaker and is overall canceled out by the stronger inverse cascade in that particular region.
Panels (c) and (d) display the distribution of number of samples of forward and inverse cascade rates, respectively. The shapes of the distributions align with Figure 10 (a) and (b), but a majority of the samples are concentrated at the center, corresponding to small values of Q and R. This observation confirms that the strong instances of inverse cascade and forward cascade observed in Figure 9 are primarily determined by infrequent but extreme events (intermittency). Note that the joint PDF in Fig. 10(d) is more left-right symmetric than that in Figure 10(c), suggesting a less non-Gaussian behavior of the flow in regions of inverse cascade than in the forward cascade regions.
Finally, to provide a visual impression of the spatial distribution of regions of negative `in the flow, in Fig. 11 we provide a 3D visualization of two instances in specific small subdomains of size 150 3 grid points (i.e. (225⌘) 3 out of the overall 8192 3 DNS domain. The selection of these subdomains was based on the condition that Q `/hQ w i > 15 and
> 15 at the center of each subdomain such that the center is at a strong vortex compression region. We then calculate the values of `, ⇧ `, Q `and R `at every second grid point. In panel (a) and (b) of Figure 11, the light blue regions correspond to isosurface of a large negative value of `/h✏i = 60, indicating the presence of an inverse cascade with significant magnitude. Clearly, we can see that the occurrence of strong inverse cascade is closely associated with the presence of the vortices. Panel (a) depicts that large negative `appears near the center and not at the core of the yellow tube, although one should recall that `(x, t) is defined locally as centered at x but represents the energy cascade into balls of diameter 45⌘, i.e. comparable to the diameter of the vortex (yellow region) shown. The blue region is also largely connected with the black isosurface, (R `= 20hQ w i 3/2 ) indicating a strong "vortex compression" region within the yellow tube. Interactive 3D versions of the figure that can be accessed following the links in the figure caption help elucidate the spatial structure. Panel (b) is an entirely di↵erent instance of similar conditions, showing a more broken up vortex and showing that `can also peak near the sides, and appear more scattered within the vortex. Coupled with the results shown in Fig. 9, the visualizations suggest that the strong inverse cascade occurs along the large scale vortices, in regions of these vortices in which R > 0, i.e. the vortex compression regions. We can also observe some yellow tubes, within which inverse cascade and compression are both absent. This is consistent to the statistics such that, when conditionally averaging in terms of Q `but irrespective of R, the inverse cascade becomes very week and almost non-existent. However, once one only considers R > 0 regions, inverse cascade can be clearly observed in high vortical regions.
We also show ⇧ `/h✏i = 20 (green-blueish isosurface) in the corresponding 3D subdomains shown as panels (c) and (d) of Figure 11. Clearly we can see that the greenblueish and black regions largely overlap within the yellow region in panel (c). In panel (d), the overlapping between green-blueish, yellow and black region occurs at the center and right-top region of the subdomain, indicating strong negative ⇧ `is also associated with strong vortex compression within high vortical region, consistent with Figure 9. However, the patterns of the green SGS flux regions are smoother, consistent to the 2D visualisation in Figure 5.
Caution must be expressed that visualizations only provide qualitative impressions and more quantitative analysis requires structure-based conditional averaging such as recently undertaken in Park & Lozano-Duran (2023). While such analysis is beyond the scope of the present paper, the conditional statistics presented in Fig. 9 already by themselves provide the strong statistically robust connection between cascade rate measures and features of the large scale velocity gradient tensor.
In this paper we explore, based on a DNS dataset of isotropic forced turbulence at a relatively high Reynolds number (R = 1250), local features of the energy cascade. We compare two common definitions of the spatially local rate of kinetic energy cascade at some scale `. The first is based on the cubic velocity di↵erence term appearing in the scale-integrated local KH equation, in the structure function approach. The second is the subfilter-scale energy flux term in the transport equation for subgrid-scale kinetic energy, i.e., as used in filtering approach often invoked in LES. Particular attention is placed on interpretation and statistical robustness of observations of local negative structurefunction energy flux and subfilter-scale energy flux. The notion and relevance of local inverse cascade or "backscatter" has been open to debates in the literature. We argue that the interpretation of `(x, t) as a spatially local energy flux appears unambiguous because it arises naturally from a divergence term in scale space. And, the symmetric formulation of Hill (2001Hill ( , 2002b) ) leads to the spherically averaged third-order structure function-based definition of a local cascade rate involving velocities at two points that are treated equally via angular averaging over the sphere.
The data confirm the presence of local instances where `(x, t) is negative, i.e. indicative of local inverse cascade events in 3D turbulence. Flow visualizations show that spatially the inverse cascade events are often located near the core of large-scale vortex structures. Comparable observations for the LES-based energy flux ⇧ `(x, t) (which also displays negative values at many locations in the flow as is well-known in the LES literature on "backscatter") show that ⇧ `(x, t) displays smoother and more blob-like features. Regarding the statistical significance of such observations, local observations from single realizations are extended using conditional averaging. Attention is placed first on relationships between the local cascade rate and the local filtered viscous dissipation rate ✏ `(x, t) that plays a central role in the classic KRSH Kolmogorov (1962). Results show that conditional averaging of both `(x, t) and ⇧ `(x, t) eliminates negative values and that the conditional averages in fact equate ✏ `(x, t) to very good approximation, entirely consistent with KRSH predictions.
The analysis then focuses on conditional averages of `and ⇧ `conditioned on properties of the filtered velocity gradient tensor properties, in particular four of its most invariants (strain and rotation rate square magnitude and the two Q R invariants). We find statistically robust evidence of inverse cascade as measured with `when both the large-scale rotation rate is strong and the large-scale strain rate is weak. When defined using ⇧ `, the conditional averaging based on large-scale strain and rotation rates does not lead to any significant average backscatter. When conditioning based on the R and Q invariants, significant net inverse cascading is observed for `in the "vortex compression" R > 0, Q > 0 quadrant. Qualitatively similar, but quantitatively much weaker trends are observed for the conditionally averaged subfilter scale energy flux ⇧ `. We recall that a multiscale decomposition of ⇧ `in terms of velocity gradients at multiple scales Johnson (2020Johnson ( , 2021) ) shows that ⇧ `< 0 appears associated with a vortex-thinning mechanism occurring at smaller scales interacting with large-scale strains.
In summary, present results show that locally negative values of kinetic energy fluxes at scale `are observed for both the structure function and filtering approaches, and at least for the structure function approach, the interpretation as a flux in scale space appears unambiguous. Regarding statistical robustness and the potential net impact of such local observations, conditional averaging reveals that both negative `and ⇧ `(representing inverse cascade) become statistically dominant mechanisms in regions where turbulent motions at scales larger than `exhibit a "vortex compression" behavior (R `> 0 and Q `> 0). However, the magnitude of inverse cascade in filtering approaches is weaker and negligible on the (S 2 `, ⌦ 2 `) map. For future work, it would be of interest to explore the sensitivity of results to Reynolds number, especially as it is expected that at higher Reynolds numbers the intermittency of the variables would increase. It would be also interesting to extend conditional averaging to more accurately reflect local energy distribution, entire flow structures and their possible connections to local inverse cascade mechanisms. Other pointwise quantities such as helicity can also be explored. It would also be instructive to connect present results with the multiscale decomposition of Johnson (2020Johnson ( , 2021) ) and thus be able to identify the small-scale mechanisms associated to local backscatter/inverse cascade events. And, further theoretically obtained exact relations between structure function and filtering approaches may yet be found. stored on volumes mounted to each user's SciServer container. The entirety of the isotropic8192 data set (8192 3 volume, 6 snapshots) in Zarr format is publicly available through Python Notebook on SciServer. Users can apply for a SciServer account freely and download the demo Notebook. In the Notebook, users can get access to pre-coded "Get" functions for arbitrary sets of points: GetPressure to retrieve and interpolate pressures, GetPressureGradient to retrieve and interpolate pressure gradient, and similarly GetPressureHessian, GetVelocity, GetVelocityGradient, GetVelocityHessian, GetVelocity-Laplacian, and GetCutout to read raw data for a user-specified box.
Demo codes for accessing data at user-specified arrays of points (in various sample geometrical configurations) are listed in the Notebook. The isotropic8192 datasets can be also accessed via the web-portal cutout service where the pyturb GetCutout function has replaced the legacy function for user queries (see https://turbulence.pha.jhu. edu/newcutout.aspx). JHTDB still provides and maintains other datasets (https: //turbulence.pha.jhu.edu/datasets.aspx) through legacy SQL systems with C, Fortran, Matlab, Python, and .Net interface. However, the aim is to transfer the existing datasets and any new coming datasets to the pyturb system in the future for faster and more stable services.
WolframAlpha online) and is given by
5/3 d = 544 ( 20/3), (B 6) with = k`and where (..) is the Gamma function. Evaluating and using C K ⇡ 1.6, the result is hk sgs,`i = 1 2 h⌧ ii i = 0.76 C K h✏i 2/3 `2/3 ⇡ 1.2 h✏i 2/3 `2/3 . (B 7)
As can be seen, the 3D integration needed to evaluate hk sf i involves the radius to the 8/3 power while that for hk sgs i involves the wavenumber to the -5/3 power. The former is thus much more strongly dominated by the large scale limit of integration (`/2) than the latter. As a result, the latter is more strongly a↵ected by the spectral behavior of turbulence at smaller scales, including the viscous range. This explains why the values of hk sgs i measured from DNS are significantly smaller than the prediction in Eq. B7, while the measurements of hk sf i agree well with the prediction in Eq. B3.