Structure function tensor equations in inhomogeneous turbulence

Exact budget equations for the second-order structure function tensor $\langle \delta u_i \delta u_j \rangle$ are used to study the two-point statistics of velocity fluctuations in inhomogeneous turbulence. The Anisotropic Generalized Kolmogorov Equations (AGKE) describe the production, transport, redistribution and dissipation of every Reynolds stress component occurring simultaneously among different scales and in space, i.e. along directions of statistical inhomogeneity. The AGKE are effective to study the inter-component and multi-scale processes of turbulence. In contrast to more classic approaches, such as those based on the spectral decomposition of the velocity field, the AGKE provide a natural definition of scales in the inhomogeneous directions, and describe fluxes across such scales too. Compared to the Generalized Kolmogorov Equation, which is recovered as their half trace, the AGKE can describe inter-component energy transfers occurring via the pressure-strain term and contain also budget equations for the off-diagonal components of $\langle \delta u_i \delta u_j \rangle$. The non-trivial physical interpretation of the AGKE terms is demonstrated with three examples. First, the near-wall cycle of a turbulent channel flow at $Re_\tau=200$ is considered. The off-diagonal component $\langle -\delta u \delta v \rangle$, which can not be interpreted in terms of scale energy, is discussed in detail. Wall-normal scales in the outer turbulence cycle are then discussed by applying the AGKE to channel flows at $Re_\tau=500$ and $1000$. In a third example, the AGKE are computed for a separating and reattaching flow. The process of spanwise-vortex formation in the reverse boundary layer within the separation bubble is discussed for the first time.


Introduction
Since the early days of fluid mechanics, understanding turbulence fascinates scholars, enticed by the goal of identifying the key mechanisms governing turbulent fluctuations and eventually determining the mean flow. This is essential for developing and improving RANS and LES turbulence models, useful in engineering practice. Most turbulent flows of applicative interest, in particular, are challenging because of their anisotropic and inhomogeneous nature.
Among the several approaches pursued so far to address the physics of inhomogeneous and anisotropic turbulence, the two most common ones observe the flow either in the space of scales, or in the physical space. In the scale-space approach, the characteristic  Kim et al. (1987). Left, adapted from Kim et al. (1987): one-dimensional energy spectra versus streamwise wavenumber κx, at two wall distances. Continuous, dashed and dotted lines refer to streamwise, spanwise and wall-normal velocity fluctuations. Right, adapted from Mansour et al. (1988): terms in the budget equation for u ′ 1 u ′ 1 , with notation as in the original paper. P11: production; ǫ11: dissipation; Π11: velocity pressure-gradient term; T11: turbulent transport; D11: viscous diffusion.
shape and size of the statistically most significant structures of turbulence are deduced from two-point second-order statistics. A spectral decomposition of the velocity field can be employed to describe the scale distribution of energy, while spatial correlation functions are used to characterise the shape of the so-called coherent structures (Robinson 1991;Jiménez 2018). Since a turbulent flow contains eddies of different scales, the power spectral density of turbulent fluctuations is a gauge to the actual eddy population, and provides useful information to develop kinematic models of turbulence capable to explain some of its features. One such model rests on the attached-eddy hypothesis by Townsend (1976), and predicts self-similar features of turbulent spectra in wall-bounded flows (Perry & Chong 1982). Two-points correlations of velocity fluctuations are the inverse Fourier transform of power spectra. They emphasise the spatial coherence of the largest and strongest turbulent fluctuations, and have been, for instance, employed to describe the streaky structure of near-wall turbulence (Kline et al. 1967), to identify largescale structures in high-Reynolds number flows (Smits et al. 2011;Sillero et al. 2014) or to describe the structural properties of highly-inhomogeneous separating and reattaching turbulent flows (Mollicone et al. 2018;Cimarelli et al. 2018).
In the physical-space approach, it is possible to characterise the spatial organisation of production, transfer and dissipation of the turbulent kinetic energy associated with the temporal fluctuations of the three velocity components. The tools of choice are the exact single-point budget equations for the components of the Reynolds stress tensor and of its half-trace, the turbulent kinetic energy k. This approach has been successfully applied to canonical wall-bounded flows and, more recently, to more complex turbulent flows. For the former, the main focus has been the inhomogeneity and anisotropy induced by the wall (Mansour et al. 1988) and the effect of the Reynolds number (Hoyas & Jiménez 2008) on the Reynolds stress budgets. For the latter, the Reynolds stress production and transport phenomena have been studied in free shear layers and recirculation bubbles (Mollicone et al. 2017;Cimarelli et al. 2018;Cimarelli et al. 2019b), where local nonequilibrium results in significantly different physics.

Structure function tensor equations in inhomogeneous turbulence
3 Typical results ensuing from the two approaches above are exemplified in figure 1, where key plots from Kim et al. (1987) and from Mansour et al. (1988) are reproduced. Both diagrams stem from the analysis of the same DNS database for a low-Re turbulent channel flow. The two leftmost plots are one-dimensional turbulent energy spectra as function of the streamwise wavenumber, each computed at a specific distance from the wall. The right plot shows the wall-normal behaviour of the terms appearing in the budget of the 1,1 component of the Reynolds stress tensor.
Despite their fundamental importance, both approaches suffer of some limitations. Indeed, it is well known since Richardson (1922) that turbulence is a truly multi-scale phenomenon, where fluctuations of different spatial extent non-linearly interact through energy-cascading mechanisms. Even more so, in inhomogeneous flows these interactions vary in space significantly, leading to a transfer of momentum between different spatial locations. The single-point budget equations for the Reynolds stresses do not contain information about the scales involved in such energy fluxes, and therefore miss the multiscale nature of turbulence. The spectral decomposition and two-point spatial correlations do discern the different scales, but fail to provide direct information on their role in the processes of production, transfer and dissipation of k, and therefore lack a dynamical description of turbulent interactions.
These limitations are overcome when space and scale properties of turbulence are considered jointly. For example, to recover the scale information Lumley (1964), Domaradzki et al. (1994) and more recently Mizuno (2016) and Lee & Moser (2019) analysed spectrally decomposed budget equations for the Reynolds stresses. They observed inverse energy transfers from small to large scales, supporting substantial modifications of the Richardson scenario in wall-bounded flows. Unfortunately, however, spectral analysis does not allow a definition of scales in statistically inhomogeneous directions, such as the wall-normal one in wall-bounded flows. Hill (2001), Danaila et al. (2001), Hill (2002) and Dubrulle (2019) proposed a complementary approach, free from this restriction, and generalized the Kolmogorov (1941) description of the energy transfer among scales from isotropic flows to inhomogeneous flows. The Generalized Kolmogorov Equation or GKE (see for example Danaila et al. 2004;Marati et al. 2004;Rincon 2006;Cimarelli et al. 2013Cimarelli et al. , 2015Cimarelli et al. , 2016Portela et al. 2017) is an exact budget equation for the trace of the so-called second-order structure function tensor, i.e. the sum of the squared increments in all three velocity components between two points in space. This quantity is interpreted as scale energy, and provides scale and space information in every spatial direction, regardless of its statistical homogeneity. The present work discusses the Anisotropic Generalized Kolmogorov Equations (AGKE), which extend the scale and space description of the GKE, limited to scale energy. The goal is to describe each component of the structure function tensor separately, thus capturing the anisotropy of the Reynolds stress tensor and of the underlying budget equations. This provides a complete description of energy redistribution among the various Reynolds stresses. The AGKE identify scales and regions of the flow involved in the production, transfer and dissipation of turbulent stresses, thus integrating the dynamical picture provided by single-point Reynolds stress budgets with the scale information provided by the spectral decomposition. The relationship between the second-order velocity increments and the two-point spatial correlation functions can be exploited to identify the topological features of the structures involved in creation, transport and destruction of turbulent stresses. This endows the kinematic information provided by the spatial correlation functions with additional dynamical information from exact budget equations.
The present work aims at introducing the reader to the AGKE and to their use via example applications to inhomogeneous turbulent flows. The paper is structured as fol- Figure 2. Sketch of the quantities involved in the definition of the second-order structure function. x = X − r/2 and x ′ = X + r/2 are the two points across which the velocity increment δu is computed.
lows. First, in §2 the budget equations for the structure function tensor are presented and provided with a physical interpretation, and the numerical datasets used in the example flows are described in §2.2. Then AGKE are applied to canonical turbulent channel flows. In particular, §3 focuses on the near-wall turbulence cycle of a low-Re channel flow. The energy exchange among the diagonal terms of the structure function tensor via the pressure-strain term is discussed, and the complete AGKE budget of the off-diagonal component is described for the first time. Then, §4 demonstrates the capability of the AGKE to disentangle the dynamics of flows with a broader range of scales by considering the outer cycle of wall-turbulence in channel flows at higher Reynolds numbers. Finally, §5 considers the separating and reattaching flow over a finite rectangular cylinder, and shows how the AGKE do in such highly inhomogeneous flows. The paper is closed by a brief discussion in §6. Additional material is reported in three appendices. The complete derivation of the AGKE and their complete form, both in tensorial and component-wise notation, are detailed for reference in Appendix A. Appendix B lists the symmetries of the AGKE terms in the specialised form valid for the indefinite plane channel. Appendix C describes the computation of the velocity field induced by the ensemble-averaged quasistreamwise vortex, employed in §3.

Anisotropic Generalized Kolmogorov Equations (AGKE)
Let us consider an incompressible turbulent flow, described via its mean and fluctuating velocity fields, U i and u i respectively, defined after Reynolds decomposition. The Anisotropic Generalized Kolmogorov Equations or AGKE are exact budget equations for the second-order structure function tensor δu i δu j , derived from the Navier-Stokes equations. The operator · denotes ensemble averaging, as well as averaging along homogeneous directions, if available, and over time if the flow is statistically stationary. The structure function tensor features the velocity increment δu i of the i-th velocity component between two points x and x ′ identified by their midpoint X = (x + x ′ ) /2 and separation r = x ′ − x, i.e. δu i = u i (X + r/2, t) − u i (X − r/2, t). (In the following, unless index notation is used, vectors are indicated in bold.) In the general case, δu i δu j depends upon seven independent variables, i.e. the six coordinates of the vectors X and r and time t, as schematically shown in figure 2, and is related (Davidson et al. 2006;Agostini & Leschziner 2017) to the variance of the velocity fluctuations (i.e. the Reynolds stresses) and the spatial cross-correlation function Structure function tensor equations in inhomogeneous turbulence 5 as follows: is the sum of the single-point Reynolds stresses evaluated at the two points X + r/2 and X − r/2 at time t, and is the two-point spatial cross-correlation function. The AGKE contains the structural information of R ij ; however, for large enough |r| the correlation vanishes, and δu i δu j reduces to V ij , whereas the AGKE become the sum of the single-point Reynolds stress budgets at X ± r/2.

Budget equations
The budget equations for δu i δu j describe production, transport and dissipation of the turbulent stresses in the compound space of scales and positions, and fully account for the anisotropy of turbulence. For a statistically unsteady turbulent flow, these equations link the variation in time of δu i δu j at a given scale and position, to the instantaneous unbalance among production, inter-component transfer, transport and dissipation. The full derivation starting from the Navier-Stokes equations is detailed in Appendix A, and Appendix B mentions the symmetries that apply in the plane channel case. The AGKE can be cast in the following compact form (repeated indices imply summation): For each (i, j) pair, φ k,ij and ψ k,ij are the components in the space of scales r k and in the physical space X k of a six-dimensional vector field of fluxes Φ ij , and are given by: . (2.7) Here δ ij is the Kronecker delta, ν is the kinematic viscosity, the asterisk superscript f * denotes the average of the generic quantity f between positions X ± r/2, and ǫ ij is the pseudo-dissipation tensor, whose trace is the pseudo-dissipation ǫ. The sum of the equations for the three diagonal components of δu i δu j reduces to the Generalized Kolmogorov Equation (Hill 2001). Each term contributing to the fluxes in equations (2.5) and (2.6) can be readily interpreted in analogy with the single-point budget equation for the Reynolds stresses (see e.g. Pope 2000) as the mean and turbulent transport, pressure transport and viscous diffusion. φ ij describes the flux of δu i δu j among scales, and turbulent transport is the sole nonlinear term. ψ ij describes the flux of δu i δu j in physical space, and all its terms but the viscous one are nonlinear. The source term ξ ij describes the net production of δu i δu j in space and among scales; it is similar to the one appearing in the GKE, but additionally features a pressure-strain term, involved in the energy redistribution process between different components of turbulent stresses. Each term in equation (2.4) informs on the spatial position X, scale r and time t at which production, transport and dissipation of Reynolds stresses are statistically important.
The diagonal components of δu i δu j are positive by definition, and their budget equations inherit the interpretation proposed by Marati et al. (2004) and Cimarelli et al. (2013) for the GKE: they are analogous to scale energy, and the AGKE enables their discrimination into the separate diagonal components of the Reynolds stress tensor. The non-diagonal components, however, can in general assume positive or negative values, also when the sign of u i u j can be predicted on physical grounds. For these components, ξ ij has the generic meaning of a source term, which can be viewed as production or dissipation only upon considering the actual sign of δu i δu j at the particular values of (X, r). In analogy with the concept of energy cascade, paths of δu i δu j in the (X, r) space represent fluxes of Reynolds stresses through space (X) and scales (r) at time t. The shape of the paths is determined by ψ ij (space fluxes) and φ ij (scale fluxes).

Simulations and databases
As anticipated in §1, the AGKE analysis below stems from the post-processing of velocity and pressure fields obtained via Direct Numerical Simulations (DNS) of two flows. The former is the turbulent plane channel flow, whose inner and outer turbulent cycles will be discussed in §3 and §4 respectively. The latter is the separating and reattaching flow around a finite rectangular cylinder, discussed in §5.
The turbulent channel flow simulations have been carried out for the present work via the DNS code introduced by Luchini & Quadrio (2006). The incompressible Navier-Stokes equations are projected in the divergence-free space of the wall-normal components Reτ  of the velocity and vorticity vectors and solved by means of a pseudo-spectral method, as in Kim et al. (1987). Three database are used, with friction Reynolds number Re τ = u τ h/ν of Re τ = 200, 500 and 1000. Here h is the channel half-height, and u τ = τ w /ρ is the friction velocity expressed in terms of the average wall shear stress τ w and the density ρ. The size of the computational domain is L x = 4πh and L z = 2πh in the streamwise and spanwise directions, discretised by N x = N z = 256, 512 and 1024 Fourier modes (further increased by a factor 3/2 for de-aliasing). In the wall-normal direction the differential operators are discretised via fourth-order compact finite differences using respectively N y = 256, 250 and 500 points collocated on a non-uniform grid. Further details are provided in table 1. In this table and throughout the whole paper, quantities denoted with the superscript + are given in viscous units, i.e. normalised with u τ and ν.
The database for the flow around around a finite rectangular cylinder is taken from the DNS study by Cimarelli et al. (2018), where the information on the numerical setup can be found. A rectangular cylinder of length 5h, thickness h and indefinite span is immersed in a uniform flow with free-stream velocity U ∞ aligned with the x direction. The Reynolds number is Re = U ∞ h/ν = 3000. The streamwise, wall-normal and spanwise size of the computation domain is (L x , L y , L z ) = (112h, 50h, 5h). The leading edge of the cylinder is located 35h past the inlet of the computational box. The fluid domain is discretised through a Cartesian grid consisting of 1.5 · 10 7 hexahedral cells. The average resolution in the three spatial direction is (∆x + , ∆y + , ∆z + ) = (6.1, 0.31, 5.41).
The AGKE terms are computed with an efficient code specifically developed for the present work, which extends a recently written code for the computation of the GKE equation . The symmetries described in Appendix B are exploited to minimise the amount of memory required during the calculations. Each term of equations (2.5), (2.6) and (2.7) is decomposed into simpler correlation terms, which are then computed as products in Fourier space along the homogeneous directions, with huge savings in computing time. For maximum accuracy, derivatives in the homogeneous directions are computed in the Fourier space, otherwise a finite-differences scheme with a five-points computational stencil is used. Finally, a parallel strategy is implemented (see Gatti et al. 2019, for details). The calculation receives in input the fluctuating velocity field for each snapshot of the databases. It outputs δu i δu j , the flux vectors ψ ij and φ ij , and the various contributions to the source term ξ ij as in equation (2.7) for each of the six different second-order structure functions, and in the whole physical and scale space.
The statistical convergence of the data is verified by ensuring that the residual of equation (2.4) is negligible compared to the dissipation, production and pressure-strain terms.  Table 2. Maximum values for diagonal terms of δuiδuj + , its source ξ + ij , absolute pressure strain |Π + ij | and production P + ij and positions in the (r + y , r + z , Y + )-space.

Example: the near-wall turbulence cycle
A turbulent channel flow at Re τ = 200 is considered in the following. The mean velocity vector is U (y) = {U (y), 0, 0}, directed along the streamwise direction x = x 1 and varying only with the wall-normal coordinate y = x 2 ; z = x 3 is the spanwise direction, and u = u 1 , v = u 2 and w = u 3 indicate the three fluctuating velocity components. Since y is the only direction of statistical inhomogeneity, δu i δu j (Y, r) and all AGKE terms are function of the physical space only through the spatial coordinate Y = (y + y ′ ) /2, while still depending upon the whole scale vector r. Similarly, spatial transport of δu i δu j occurs along Y through the only nonzero component of the spatial flux ψ ij = ψ Y,ij .
The GKE for the scale energy δu 2 i has been thoroughly discussed in literature, (see e.g. Marati et al. 2004;Cimarelli et al. 2013Cimarelli et al. , 2015Cimarelli et al. , 2016, and different interpretations and visualisation techniques have been suggested. For this reason, in the following we only address the new information offered by the AGKE. This includes the analysis of the anisotropic scale-energy redistribution operated by the pressure-strain terms, and that of the budget equation for −δuδv . The analysis is also restricted to the subspace r x = 0: this is motivated by the turbulent vortical structures in channel flow being predominantly aligned in the streamwise direction. Such structures typically induce the largest negative correlation of velocity components for r x = 0 and characteristic values of r z . A classic example are the so-called near-wall streaks, for which r + z ≈ 60. As a consequence of (2.1), the local maxima of, for instance, δuδu and terms appearing in its budget equation also occur for r x ≈ 0. Note that in the r x = 0 space the terms of the AGKE are not defined below the Y = r y /2 plane, owing to the finite size of the channel in the wall-normal direction.

Scale-energy redistribution by pressure strain
The pressure-strain term Π ij redistributes energy among the diagonal components of δu i δu j . Hence, at different scales and positions this term can be a source or a sink depending on its sign. To better understand its behaviour and link it to physical processes, it is instructive to briefly analyse the scales and position at which δuδu , δvδv , δwδw and their sources ξ ij are important.
The position and the intensity of the maxima, hereinafter denoted with the subscript m, of the diagonal components of δu i δu j and of the associated ξ ij and Π ij are reported in table 2. δuδu , δvδv and δwδw peak at small scales within the buffer layer, similarly to δu 2 i (Cimarelli et al. 2016), with δvδv m located further from the wall. The anisotropy of the flow is denoted, for instance, by δwδw m being much lower than δuδu m and occurring at r z = 0 and small r y , whereas the other maxima occur at r z = 0 and r y = 0. This difference is explained by the quasi-streamwise vortices populating the near-wall cycle (Schoppa & Hussain 2002): they induce negatively correlated regions of spanwise fluctuations at r y = 0 and of streamwise and wall-normal fluctuations at r z = 0. The region of negative source terms partially coincides with the one of the source term in the GKE (see e.g. Cimarelli et al. 2016). As in the GKE, negative sources are observed at the lower boundary Y = r y /2, and in the whole channel height at r y , r z → 0: viscous dissipation dominates near the wall and at the smallest scales. However, the regions of large positive sources vary significantly among the three diagonal components (see table 2). This is due to the different nature of the positive source of the three diagonal components of δu i δu j . Indeed, in a turbulent channel flow the streamwise fluctuations are fed by the energy draining from the mean flow (i.e. by the production term P 11 ), whereas the cross-stream fluctuations are produced by the redistribution processes (i.e. the pressure-strain term Π 22 and Π 33 ). This explains also the larger order of magnitude of ξ 11,m . Unlike the GKE, the scale and space properties of this energy redistribution can be extracted from the AGKE (see equation 2.7). Figure 3 plots the pressure-strain term for the diagonal components, with values and positions of their maxima as reported in table 2. The figure shows the location of the pressure-strain maximum in absolute value together with the maximum production. Large values of P 11 occur near the plane Y + = r + y /2+14, except for the smallest scales in the region r + y < 30 and r + z < 20. On the other hand, Π 11 is negative almost everywhere, showing that the streamwise fluctuations lose energy at all scales to feed the other components. In particular, large negative values of Π 11 , albeit much smaller than P 11 , are seen near the plane Y + = r + y /2 + 24, except for the region r + y , r + z < 20. This brings to light the dominant scales and wall distances involved in the process of redistribution of δuδu towards the other components, and discriminates them from those involved in its production. On the contrary, at the smallest scales where viscous dissipation is dominant production and redistribution are not observed.
The pressure-strain terms of the cross-stream components, Π 22 and Π 33 , are positive almost everywhere; they show a positive peak near the wall and remain larger than dissipation in different regions of the r x = 0 space. Their maxima are located in the vicinity of the plane Y + = r + y /2 + 40 for Π 22 and Y + = r + y /2 + 14 for Π 33 , where Π 11 is negative. Hence, at these scales and wall-normal distances δuδu loses energy towards δvδv and δwδw . Moreover, Π 22 is negative in the very near-wall region, Y + < r + y /2+5, owing to the non-penetration wall boundary condition which converts δvδv into δuδu  Table 3. Maximum value for −δuδv + , maximum and minimum for the source ξ + 12 , minimum for the pressure strain Π + 12 and maximum of the production P + 12 and positions in the (r + y , r + z , Y + )-space.
and δwδw . Indeed, here Π 11 and Π 33 are positive. This phenomenon is known as the splatting effect (Mansour et al. 1988), and shows no scale dependency. Different values of Π 22 and Π 33 imply an anisotropic redistribution of the streamwise fluctuations to the other components. Owing to the incompressibility constraint, the following relationship holds: Hence, Π 22 /Π 11 = Π 33 /Π 11 = −0.5 corresponds to isotropic transfer of energy from the streamwise fluctuations towards the other components. In figure 3a the isosurface Π 22 /Π 11 = −0.5 is shown. The inner side at small scales of this surface is characterised by Π 22 /Π 11 < −0.5, and thus by Π 22 > Π 33 (as long as Π 11 < 0). Hence, at small scales the pressure strain preferentially redistributes streamwise energy to the vertical fluctuations.
On the contrary, on the outer side of the surface Π 33 > Π 22 holds, implying that at larger scales the streamwise energy is preferentially redistributed towards spanwise fluctuations.

Scale-by-scale budget of the off-diagonal term −δuδv
The only off-diagonal term associated with a nonzero component of the Reynolds stress tensor is −δuδv which, unlike the diagonal terms, is not definite in sign. Therefore, −δuδv and its fluxes cannot be interpreted in terms of energy and energy transfer. −δuδv describes the statistical dependence or, more precisely, the correlation between δu and δv and, for large r, the mean momentum transfer. Concepts as production and dissipation only apply to the source term ξ 12 after the sign of −δuδv is taken into account.

Intensity, production and redistribution
The off-diagonal term −δuδv and its budget are plotted in figure 4, and corresponding quantitative information is reported in table 3. As shown by figure 4a, −δuδv is positive almost throughout the entire physical/scale space except at very small separations (r + z = 0, r y ≤ 10) for Y + < 50. The largest positive values of −δuδv are in the buffer layer at 15 ≤ Y + ≤ 60, at spanwise scales 40 ≤ r + z ≤ 80 and vanishing r y . A second, less prominent local maximum of −δuδv is located near the r z = 0 plane.
The source term ξ 12 , plotted in figure 4b, is dominated by the (positive) production term P 12 and the (negative) pressure-strain term Π 12 (see equation (A 17) for their definitions). Indeed, the viscous pseudo-dissipation D 12 plays a minor role, as in the single-point budget for −uv (see e.g. Mansour et al. 1988). Large positive and negative values of ξ 12 define two distinct regions in the buffer layer (figure 4d). The positive peak corresponds to spanwise scales 10 ≤ r + z ≤ 50, while the negative one to small scales (r + z ≈ 0). Moreover, ξ 12 is negative in a portion of the Y + = r + y /2 plane, implying that turbulent structures extending down to the wall are inactive in the production of −δuδv .
It is worth noting that ξ 12 strongly varies with spanwise separation, as seen in the r y = 0 plane (figure 4c; see also figure 5). In comparison to the global picture obtained from single-point analysis of −uv in the buffer layer (here recovered in the limit r z → L z /2) where the source term is slightly negative, one can additionally appreciate the existence of a large positive peak of ξ 12 at r + z = 20 and a negative one at r + z = 70 (figure 5b). Indeed, P 12 and Π 12 are of the same order of magnitude throughout the r y = 0 plane, but reach their extreme values at different spanwise scales, see figure 4c. In particular large values of P 12 are found at (r + z , Y + ) ≈ (30, 17), whereas large negative values of Π 12 are found at (r + z , Y + ) ≈ (60, 16). The structural interpretation of these findings is discussed below in §3.2.3.

Fluxes
The transfer of −δuδv in space and among scales is determined by the flux vector (φ y , φ z , ψ), and is visualised via its field lines. These field lines can be grouped in two families. The lines of the first family enter the domain from the channel centerline, Y = h, and descend towards the wall; they can be further grouped in sets I, II and III as shown in figure 4b. The second family only contains set IV, and is visible in the zoomed figure 4d; its field lines are confined to the near-wall region, and connect the positive and negative peaks of ξ 12 .
Various quantities can be tracked along representative field lines, as done in figure 6. The position along a field line of length ℓ in the (r z , r y , Y ) space is described by the normalised curvilinear coordinate s = 1 ℓ ℓ 0 ds with ds = dr 2 z + dr 2 y + dY 2 . (3. 2) The values of r y , r z and Y (see figure 2) are plotted in the left column of the figure; the

Structure function tensor equations in inhomogeneous turbulence
13 central column plots the evolution of −δuδv , ξ 12 , P 12 , Π 12 and the pseudo-dissipation D 12 along the line; the right column plots the evolution of the correlation coefficient ρ ij defined by where repeated indices do not imply summation. R ij is linked to δu i δu j by equation (2.1).
The top and central panels of figure 6 illustrate the evolution of various quantities along representative lines of set I and II. Both lines are qualitatively similar: they highlight a transfer of −δuδv from the centerline to the near-wall region, through first decreasing and then increasing wall-normal scales. At the centreline they are parallel to the Y axis, consistently with the AGKE symmetries (see appendix B). However, lines of set I are attracted by the negative peak of ξ 12 towards smaller r z , while those of set II are repulsed from the positive source peak towards larger r z . Lines of set III are not shown for the sake of brevity, since they pass through regions of large separations and are characterised by almost zero correlation, see equation (3.3). On the other hand, lines of sets I and II exist at smaller r y and r z and, as shown in the upper-right and central-right panels of figure  6, are characterised by finite levels of correlation. Along lines of set I and II, −δuδv increases from zero at the centerline (due to the AGKE symmetries) to reach a positive peak in the near-wall region. Similarly, ξ 12 shows a negative/positive peak when the lines of set I/II approach the near-wall region as the pressure-strain/production overcomes the production/pressure-strain.
The evolution of the correlation coefficients −ρ 12 and −ρ 21 (recall that ρ ij = ρ ji for i = j, see equation (3.3)) is used to extract information about the turbulent structures involved in production, transfer and dissipation processes highlighted along the lines. As shown in the left-top and left-central panels of figure 6, at values of the curvilinear coordinate s > 0.75 corresponding to Y + < 60, lines of set I intersect positive −ρ 12 and −ρ 21 for small r z and r y , while those of set II intersect negative correlations at larger r + z ∼ 50 and smaller r y . For both sets, this is consistent with the flow field induced by near-wall quasi-streamwise vortices, creating positive and negative cross-correlation at values of separation in agreement with the present analysis; positive −ρ 12 is associated to u and v fluctuations at the same-side of the vortices (i.e. small r z ), whereas negative −ρ 12 is associated to opposite-side fluctuations. Hence, we relate the peaks of P 12 and Π 12 (and consequently of ξ 12 ) along the lines of set I and II to such structures.
The lines of set IV, shown in figures 4b and 4d and in the bottom panels of figure  6, behave differently. The field lines originate in the lower boundary of the domain at (r + y , r + z , Y + ) = (6, 15, 3). Along their path they first intercept the positive peak of ξ 12 at small r y where −δuδv is maximum. Then, they pass through the negative peak of ξ 12 , located at smaller r z and larger r y , where −δuδv is smaller. Eventually, they again vanish in the lower boundary of the domain.
Focusing on the correlation coefficient −ρ 12 , lines of set IV intersect a positive value along their complete extension. In detail, the lines first intersect small values of −ρ 12 for r + z ≈ 20 and Y + < 5 and then larger −ρ 12 for smaller r + z and larger Y + . Hence, this set of lines highlights a transfer of −δuδv between the small uv-structures created in the viscous sublayer by the wall boundary condition (Sillero et al. 2014) and the turbulent structures of the near-wall cycle.

Structural properties of wall turbulence
To connect the main statistical features of −δuδv in the buffer layer to the turbulent structures that populate it, we compute the −δuδv AGKE budget from the velocity field induced by the ensemble-averaged quasi-streamwise vortex. Such vortex, visualised in figure 7a, represents the characteristic near-wall coherent structure in the average sense. The procedure which extracts the ensemble-average vortical structure from the DNS database is very similar to the one presented by Jeong et al. (1997), which is slightly modified here to focus on the structures in the buffer layer only. Details of the procedure are provided in Appendix C.
The ensemble-averaged velocity field is shown in figure 7b in a z + − y + plane passing through the vortex centre. The corresponding −δuδv , normalised by its maximum in the r x = 0 space, is shown in figure 7c in the r x = r y = 0 plane. −δuδv computed for the average structure shows a remarkable agreement with the same quantity computed for the turbulent channel flow. In particular, its maximum occurs at r + y , r + z , Y + = (0, 52, 25), i.e. nearly the same location r + y , r + z , Y + = (0, 53, 30) observed for the full velocity field (see table 3). Figure 7d shows the production P 12 and the pressure-strain Π 12 normalised with the maximum production in the r x = 0 space. Again, the average quasi-streamwise vortex represents well the typical r z scales of production and pressure-strain of −δuδv . The peak of P 12 occurs at (r + z , Y + ) = (39.2, 20.0) while the minimum of Π 12 is located at (r + z , Y + ) = (52.3, 19.0), i.e. at a larger spanwise scale, similar to what figure 4c shows for the full velocity field.

Example: the outer turbulence cycle
Thanks to its ability to account for scales also in directions of statistical inhomogeneity, the AGKE becomes increasingly informative as the range of turbulent scales widens. For the turbulent channel flow, Re τ is the ratio between the outer geometrical lengthscale h and the inner viscous lengthscale ν/u τ . Hence, for increasing values of Re τ , the typical scales of the autonomous near-wall cycle discussed in §3 are constant in viscous units but shrink compared to h. Meanwhile, a whole new hierarchy of motions starts to appear: they include larger scales in the logarithmic region and form the so-called outer cycle (see, for instance, Cossu & Hwang 2017). The wall-normal extent of such motions is typically not accounted for by other frameworks for the analysis of scale transfers, but can be easily studied by the AGKE.
A comparative AGKE analysis for a channel flow at the three different values of Re τ = 200, 500 and 1000 is presented below. The main features of the DNS databases have been already introduced in §2.2. The profiles of mean velocity and variance of velocity fluctuations at all values of Re considered in the following are reported in figure 8, which confirms the full agreement of such statistics with the database available from Lee & Moser (2015). Figure 9a shows the contour ξ 11 = 0 in the (r z , Y ) plane at r x = r y = 0. Taking r y = 0 is equivalent to the classic approach, where only wall-parallel separations (or wavenumbers in the spectral analysis) are considered. Three different regions of net energy source ξ 11 > 0, enclosed by the isoline ξ 11 = 0, can be distinguished. The first region, which collapses for all values of Re with viscous scaling, corresponds to the net production of δuδu within the near-wall cycle, already described in §3, and takes place at all spanwise separations. The second region of ξ 11 > 0 is found for approximately r + z ≤ h + and Y + ≤ 0.6h + . Here the left boundary of the contour ξ 11 = 0 represents the cross-over value of r z , for a given Y , separating larger production scales from smaller inertial scales. The cross-over scale increases approximately linearly with the wall distance, in agreement with the overlap layer predictions of the attached-eddy model (Townsend 1976). Cimarelli et al. (2015) carry out a detailed analysis of the scaling properties of this second source region, albeit in terms of δu i δu i , while Marusic & Monty (2019) discuss the attached-eddy model and its implications. This second region of ξ 11 > 0 is observed also with the analysis based on one-dimensional premultiplied spectral budgets (see, for instance, figure 5  one of the near-wall cycle. It is also interesting to note that this region, albeit weak and confined to a tiny range of spanwise scales and wall-normal positions, is already apparent at Re τ = 200, something that can not be observed as easily from one-dimensional spectra. Only for the largest value Re τ = 1000 considered here, a third region of ξ 11 > 0 appears, with spanwise scales 2h + ≤ r + z ≤ 3h + and values of Y + pertaining to the logarithmic layer. This third region is related to the production by additional largescale turbulent features, whose statistical footprint cannot be predicted by using the attached-eddy hypothesis (Marusic & Monty 2019). These motions have been named superstructures (Hutchins & Marusic 2007) when found in boundary layers and Large Scale Motions (LSM) or Very Large Scale Motions (VLSM) (Guala et al. 2006) when observed in turbulent channels, pipes and plane Couette flows. Henceforth we will adopt the acronym LSM, disregarding the slight differences in the definition of the three terms given in literature. LSM are important for two main reasons. First, their relative contribution to the total turbulent kinetic energy and Reynolds shear stress rapidly increases with Re τ (Ganapathisubramani et al. 2003), making LSM one of the main players in the outer cycle and thus an obvious target for flow control. Second, LSM modulate the inner cycle (Mathis et al. 2009) and superpose to the near-wall turbulence (Hoyas & Jimenez 2006), thus causing the failure of exact viscous scaling for several statistical quantities, such as for example the wall-normal profiles of the streamwise and spanwise velocity fluctuations. Figure 9b focuses on the Re τ = 1000 case, and illustrates how the AGKE can naturally consider scales in the wall-normal inhomogeneous direction, something particularly useful to describe the volume-filling LSM. Contours of ξ 11 at Re τ = 1000 are plotted in the (r y , Y ) plane for r x = 0 and r + z = 2300, i.e. the spanwise scale at which LSM have been observed in Figure 9a. The results reveal the wall-normal distribution of the net positive source, i.e. net production of δuδu , occurring at the scales of the LSM throughout the channel. Positive ξ 11 is observed for 150 ≤ Y + ≤ 0.5h + at wall-normal scales in the range 0 ≤ r + y ≤ 400, while the bottom part of the contours runs parallel to the line Y + = r + y /2 + 150, indicating that the wall-normal scales related to LSM are selfsimilar, contrary to the spanwise ones. The wall-normal location and scale at which ξ 11 is active agrees remarkably well with the wall-normal extent of LSM measured by Madhusudanan et al. (2019) utilising high-Re DNS data and linearised Navier-Stokes equations subject to stochastic forcing. Interestingly, positive ξ 11 at the LSM spanwise scale occurs also for r + y ≈ 1.7h + and Y + ≈ h + (see figure 9b), indicating that δuδu is also produced at very large wall-normal scales at the centerline and thus that large-scale negative correlation of the streamwise velocity fluctuations is produced across the two channel halves.

Example: separating and reattaching flows
The separating and reattaching flow over a rectangular cylinder with length-to-height ratio of 5 is a popular benchmark for bluff-body aerodynamics (Bruno et al. 2014), known as BARC. It is considered here as an example of complex flow with two inhomogeneous directions and multiple separations and reattachments. Various flow structures are known to exist in different parts of the main recirculating bubble, and recently it has been suggested ) that streamwise-and spanwise-oriented vortices populate the attached and detached portion respectively of the reverse boundary layer.
The snapshots used below for the AGKE analysis of the BARC flow are taken from the DNS study by Cimarelli et al. (2018). Figure 10 visualises the mean and instantaneous velocity fields. Three recirculation zones are present: a large-scale primary bubble originating from the leading-edge separation, a separation in the wake and a smaller secondary recirculation within the primary bubble. Separating and reattaching flows often feature the simultaneous presence of small scales, related to turbulent motions, and large scales, related to shedding of large-scale vortices. A full understanding of their interaction would be of paramount importance for the correct prediction and control of the flow (Kiya & Sasaki 1983;Cherry et al. 1984;Kiya & Sasaki 1985;Nakamura et al. 1991;Tafti & Vanka 1991). In particular, transition in the leading-edge shear layer is strongly affected by such multi-scale interactions: a region with negative turbulence production has been identified (Cimarelli et al. 2019a), which leads to overwhelming difficulties with turbulence closures (Bruno et al. 2014). A key role is played by the turbulent structures advected within the main recirculating bubble, which trigger the transition of the leadingedge shear layer that in turn creates them, thus effectively belonging to a self-sustaining cycle. Remarkably, these structures appear to be quasi-streamwise vortices at the beginning of the reverse boundary layer and, while working their way upstream, become spanwise vortices. However, this process is far from being fully understood, and the AGKE will be used to clarify it. Note that, since statistical homogeneity only applies to the spanwise direction and time, all two-point statistics involved in the AGKE are now function of the separation vector r, and the two spatial coordinates X = (x + x ′ ) /2 and Y = (y + y ′ ) /2. In the figures that follow, lengths and velocities are made dimensionless with the free-stream velocity and the cylinder height. We start with the component δvδv , since it is the most obvious proxy for the local alignment of turbulent structures; in fact a streamwise structure would be revealed by a local maximum of δvδv at r x = 0 and a finite r z , whereas a spanwise structure implies a local maximum at finite r x and r z = 0. In figure 11 the pressure-strain term Π 22 is shown in the (X, Y, r z ) space that embraces the whole primary bubble for r x = r y = 0. Π 22 is first observed to mark clearly the outer edge of the bubble. Within the bubble, Π 22 is highly scale-and position-dependent, and it differs from channel flow as discussed in §3. For instance, along the reverse attached boundary layer, i.e. for −0.8 ≤ X ≤ 1 and Y ≤ 0.75, Π 22 shows an evident positive peak at small spanwise scales (r z < 0.1) even very near the wall, whereas in the channel flow the splatting effect leads to negative Π 22 (see figure 3 in §3). Therefore, in this region Π 22 feeds clearly identified spanwise scales which are compatible with streamwise-aligned vortices. However, closer to the detachment of the reverse boundary layer (i.e. −0.8 ≤ X ≤ −1.2), an abrupt change takes place: Π 22 becomes positive at every spanwise separation, suggesting that once detached the reverse boundary layer is no longer populated by streamwise vortices.
Further insight on the local structure of turbulence in the detachment zone is obtained by looking at δvδv and δwδw in the (X, r x , r z ) space, shown in figure 12 for (Y, r y ) = (0.64, 0). Identifying spanwise-oriented structures requires considering scales r x along the inhomogeneous streamwise direction. Indeed δvδv locally peaks at (X, r x , r z ) = (−0.95, 0.3, 0), i.e. exactly at the X position where the boundary layer detaches and for a specific streamwise scale. This confirms the suggestion by Cimarelli et al. (2018) that  spanwise-oriented structures are indeed present. δwδw too exhibits a local maximum for finite r x , precisely at (X, r x , r z ) = (−1.13, 0.65, 0). However, the streamwise extent of this peak is larger than that for δvδv . Moreover, δwδw increases within the secondary recirculation bubble, where it features a non-monotonic behaviour in r z , while δvδv does not. Hence, the detached reverse boundary layer and, in particular, the secondary recirculation bubble appear to be populated by a broader range of structures than just spanwise-oriented vortices. The process behind the formation of spanwise-aligned structures is addressed in figure 13, which shows the production terms P 11 and P 22 in the same (X, r x , r z ) space of figure 12. P 11 has a local maximum at (X, r x , r z ) = (−1.05, 0.35, 0). At these scales, the streamwise fluctuations drain energy from the mean shear and feed δvδv , which has been connected at such scales to spanwise structures. The process is described by the pressure-strain terms: at these scales indeed it is found (not shown) that Π 11 < 0 and Π 22 > 0. Similarly, P 22 is negative everywhere, with a relative minimum in same range of scales where P 11 is maximum. Thus, P 22 reconverts the energy δvδv received via pressure strain back to the mean flow, thereby sustaining the detachment of the reverse boundary layer.
Hence, within the limits of this necessarily brief example, the AGKE successfully confirm the literature suggestion that spanwise-oriented structures exist at the detachment of the reverse boundary layer. Moreover, they reveal that these structures do not simply derive from the upstream streamwise-oriented ones simply via a gradual reorientation. Instead, their appearance is rather abrupt, mediated by pressure-strain redistribution but mainly driven by local positive and negative production.

Concluding discussion
Exact budget equations for the components of the second-order structure function tensor δu i δu j (X, r) have been considered. Because of its close relationship with twopoint velocity correlations and spectra, δu i δu j is interpreted as scale Reynolds stress. In this spirit, the budget equations, that we name Anisotropic Generalized Kolmogorov Equations (AGKE), describe production, transport and dissipation of the scale Reynolds stresses in the combined physical and scale space.
Compared to the Generalized Kolmogorov Equation (GKE), which is half the trace of the AGKE and thus describes scale energy only, the AGKE fully account for the anisotropy of the structure function tensor, and allow the description of purely redistributive processes like pressure-strain. They are a powerful tool to complement energy spectra of turbulent fluctuations and spectral Reynolds stress budgets (see, for instance, Mizuno 2016; Lee & Moser 2019), to which they add two major features: i) scales are defined along directions of statistical inhomogeneity; and ii) fluxes are defined in the space of scales. Thanks to the former feature, scale properties of turbulence can be assessed also along the wall-normal direction of wall-bounded turbulent flows and, in general, in complex turbulent flows. Thanks to the second feature, fluxes of δu i δu j across all scales and in physical space can be clearly recognised. Thus, beside the identification of scales acting as donors or receivers of scale Reynolds stresses, already possible within the framework of spectral Reynolds stress budgets, the AGKE allow to quantify the local direction of the fluxes of δu i δu j throughout the whole (X, r) space, informing on the different physical processes underlying the transfer of scale Reynolds stress in space or through scales at different spatial positions in the flow.
The AGKE have been demonstrated via three examples. With a low-Re turbulent plane channel flow, the near-wall turbulence cycle has been observed and described in terms of the AGKE, thanks to its multi-dimensional and multi-component information. The pressure-strain term of the diagonal components of δu i δu j is analysed to identify scales and positions involved in the inter-component energy redistribution processes. Moreover, the budget equation for the off-diagonal component −δuδv , the other important element that the AGKE adds to the GKE, is presented and discussed. In contrast to the energetic interpretation of the diagonal components, the scale Reynolds shear stress −δuδv is not positive definite, and is rather interpreted as statistical proxy for coherent structures and related to the production of δuδu . The main transport mechanisms are identified via the combined analysis of the AGKE terms and of the correlation levels along typical transport patterns in the physical and scale space.
Channel flows at higher Re (up to Re τ = 1000) are also considered in order to demon-strate the AGKE on flows characterised by a broader range of scales with particular focus on the outer cycle of wall-turbulence. The range of scales and positions responsible for the net production of streamwise turbulent fluctuations in the outer layer are identified.
In particular, the presence of two well-separated self-regenerating cycles belonging to scales attached to the wall and to very-large scale motions are unequivocally detected in a quantitative way. Finally, the separating and reattaching flow over a finite rectangular cylinder is considered as a test case with two inhomogeneous directions. The AGKE describe how streamwise-oriented structures in the reverse boundary layer within the main recirculation bubble become spanwise-oriented structures in the detachment region. The pressurestrain and production terms show that the spanwise structures form abruptly near the detachment, rather than being gradually reoriented.
The AGKE are a tool with several potential applications. Thanks to the relationship between δu i δu j and the unresolved stresses , the AGKE can be useful to develop large-eddy turbulence models. Indeed, Cimarelli & De Angelis (2014) already used the GKE a posteriori to improve modeling, and the AGKE could further this approach, by fully accounting for anisotropy, an essential property of wall-bounded turbulent flows. For canonical turbulent flows at large values of Re, the AGKE seem apt to comprehensively describe the large-scale structures involved in the outer regeneration cycle (Hwang & Cossu 2010) and their modulating effect (Mathis et al. 2009) onto near-wall turbulence. Such structures, characterised by a large wall-normal extent (Hutchins & Marusic 2007), may be involved in a non-negligible transfer of δu i δu j across wall-normal scales, which is captured by the AGKE but escapes either the spectral Reynolds stress budgets and the analysis based upon structure function alone (Agostini & Leschziner 2017). Similarly, in plane Couette flow the AGKE could be used to study the transfer from small to large scales, resulting from the interaction of small near-wall structures with large scales further from the wall, which has been experimentally observed by Kawata & Alfredsson (2018) only for the Reynolds shear stress but not for the normal components. The AGKE can also be used to study how turbulent wall-bounded flows are modified by drag reduction .
Beside their application to canonical flows, the present paper demonstrates that AGKE can provide significant contributions in the study of all those complex flows, such as a backward-facing step, a three-dimensional turbulent boundary layer, flows over complex surfaces, with shear layers and with separation, where anisotropy and inhomogeneity are important.
where ρ is the fluid density, and p the pressure. The two points x and x ′ are independent: hence v i and p only depend on x, while v ′ i and p ′ only depend on x ′ , and The Reynolds decomposition of the velocity field is now introduced: v i = U i + u i where U i = v i denotes the mean velocity and u i the fluctuations. The two equations become: By subtracting equation (A 5) from (A 6) and using the following relations, derived from the independence of x and x ′ , an equation for the velocity increment δu i = u ′ i − u i is obtained: (A 7) By adding and subtracting u k (∂δu i /∂x ′ k )+u k (∂δU i /∂x ′ k )+U k (∂δu i /∂x ′ k )+U k (∂δU i /∂x ′ k ) to the left-hand side and observing that

equation (A 7) becomes
(A 8) Equation (A 8) multiplied by δu j is now summed to the same equation, with the i-index switched to j-index and after multiplication by δu i . We then use incompressibility and