Invariants of the velocity-gradient tensor in a spatially developing inhomogeneous turbulent flow

Tomographic particle image velocimetry experiments were performed in the near field of the turbulent flow past a square cylinder. A classical Reynolds decomposition was performed on the resulting velocity fields into a time invariant mean flow and a fluctuating velocity field. This fluctuating velocity field was then further decomposed into coherent and residual/stochastic fluctuations. The statistical distributions of the second and third invariants of the velocity-gradient tensor were then computed at various streamwise locations, along the centreline of the flow and within the shear layers. These invariants were calculated from both the Reynolds-decomposed fluctuating velocity fields and the coherent and stochastic fluctuating velocity fields. The range of spatial locations probed incorporates regions of contrasting flow physics, including a mean recirculation region and separated shear layers, both upstream and downstream of the location of peak turbulence intensity along the centreline. These different flow physics are also reflected in the velocity gradients themselves with different topologies, as characterised by the statistical distributions of the constituent enstrophy and strain-rate invariants, for the three different fluctuating velocity fields. Despite these differing flow physics the ubiquitous self-similar ‘tear drop’-shaped joint probability density function between the second and third invariants of the velocity-gradient tensor is observed along the centreline and shear layer when calculated from both the Reynolds decomposed and the stochastic velocity fluctuations. These ‘tear drop’-shaped joint probability density functions are not, however, observed when calculated from the coherent velocity fluctuations. This ‘tear drop’ shape is classically associated with the statistical distribution of the velocity-gradient tensor invariants in fully developed turbulent flows in which there is no coherent dynamics present, and hence spectral peaks at low wavenumbers. The results presented in this manuscript, however, show that such ‘tear drops’ also exist in spatially developing inhomogeneous turbulent flows. This suggests that the ‘tear drop’ shape may not just be a universal feature of fully developed turbulence but of turbulent flows in general.

Tomographic particle image velocimetry experiments were performed in the near field of the turbulent flow past a square cylinder. A classical Reynolds decomposition was performed on the resulting velocity fields into a time invariant mean flow and a fluctuating velocity field. This fluctuating velocity field was then further decomposed into coherent and residual/stochastic fluctuations. The statistical distributions of the second and third invariants of the velocity-gradient tensor were then computed at various streamwise locations, along the centreline of the flow and within the shear layers. These invariants were calculated from both the Reynolds-decomposed fluctuating velocity fields and the coherent and stochastic fluctuating velocity fields. The range of spatial locations probed incorporates regions of contrasting flow physics, including a mean recirculation region and separated shear layers, both upstream and downstream of the location of peak turbulence intensity along the centreline. These different flow physics are also reflected in the velocity gradients themselves with different topologies, as characterised by the statistical distributions of the constituent enstrophy and strain-rate invariants, for the three different fluctuating velocity fields. Despite these differing flow physics the ubiquitous self-similar 'tear drop'-shaped joint probability density function between the second and third invariants of the velocity-gradient tensor is observed along the centreline and shear layer when calculated from both the Reynolds decomposed and the stochastic velocity fluctuations. These 'tear drop'-shaped joint probability density functions are not, however, observed when calculated from the coherent velocity fluctuations. This 'tear drop' shape is classically associated with the statistical distribution of the velocity-gradient tensor invariants in fully developed turbulent flows in which there is no coherent dynamics present, and hence spectral peaks at low wavenumbers. The results presented in this manuscript, however, show that such 'tear drops' also exist in spatially developing inhomogeneous turbulent flows. This suggests that the 'tear drop' shape may not just be a universal feature of fully developed turbulence but of turbulent flows in general. The smallest scales present in a turbulent flow have long been thought to be well approximated as statistically homogeneous and isotropic (Kolmogorov 1941;Batchelor 1953). The general topology of these fine scales may be shown to depend on the invariants of the velocity-gradient tensor (VGT) which may be split up into a symmetric and skew-symmetric component, respectively the strain-rate and rotation tensors, in which u i denotes the fluctuating components of velocity from a classical Reynolds decomposition. These invariants are defined as the coefficients in the characteristic equation for the VGT of the form ξ 3 + Pξ 2 + Qξ + R = 0. (1. 2) The first invariant, P, is the negative trace of the VGT (P = −a ii ) and is thus identically zero for an incompressible flow. Hence, the generalised topology of the flow may be described by the invariants Q and R, defined as in which ω i = ijk ω jk are the components of the vorticity vector. Q may thus be considered the local excess of swirling over strain rate, with its constituent invariants Q ω being simply half the magnitude of the enstrophy whilst Q s = − /4ν, where is the rate of dissipation of turbulent kinetic energy. R may be considered the local excess of strain-rate production (self-amplification) over enstrophy production (vorticity stretching due to the interaction between strain rate and rotation).
The joint probability density function (p.d.f.) between Q and R, f QR (Q, R), is observed to make self-similar 'tear drop' shapes in a variety of fully developed turbulent flows including homogeneous isotropic turbulence, mixing layers and wall-bounded flows (e.g. Soria et al. 1994;Blackburn, Mansour & Cantwell 1996;Tsinober 2009;Buxton & Ganapathisubramani 2010). The ubiquity of this 'tear drop'-shaped joint p.d.f. has led to it being described as a universal feature of fine-scale turbulent motions (Chacin & Cantwell 2000;Elsinga & Marusic 2010).
The state of the VGT itself may be broadly divided up into four sectors of the Q-R space according to the signs of R and the discriminant of (1.2), ∆. For an incompressible flow (P = 0) the line ∆ = Q 3 + (27/4)R 2 separates purely real solutions to (1.2) (∆ < 0), physically interpreted as the flow being locally straining only, from one real and a complex conjugate pair of roots, for which the flow is locally swirling (∆ > 0). The four sectors may thus be defined as I : ∆ > 0, R < 0 is referred to as stable foci which is the sector primarily responsible for enstrophy amplification, II : ∆ > 0, R > 0 is referred to as unstable foci and primarily responsible for enstrophy attenuation, III : ∆ < 0, R > 0 is referred to as unstable nodes and IV : ∆ < 0, R < 0 is referred to as stable nodes (Soria et al. 1994). These sectors are illustrated in figure 6(a). Cantwell (1992) made the assumption that the cross-derivatives of the pressure field and viscous diffusion were negligible in the evolution equation for the VGT and was hence able to model some of the phenomenology of fine-scale turbulence. However, Cantwell (1993) showed that whilst this restricted Euler analysis Invariants of the velocity-gradient tensor in a spatially developing flow 3 may account for certain features of the 'tear drop'-shaped joint p.d.f., it was unable to account for the particular form of this statistical distribution.
The vast majority of studies examining the statistical distribution of the invariants Q and R have done so in fully developed turbulent flows, in which the spectrum of turbulent kinetic energy is akin to the model spectrum proposed by Pope (2000). More recently, however, Gomes-Fernandes, Ganapathisubramani & Vassilicos (2014) have examined the evolution of the state of the VGT in a spatially developing turbulent flow generated by a multi-scale space-filling fractal square grid. The study examined the joint p.d.f. between Q and R at three streamwise locations, upstream of the location of peak turbulence intensity (in the 'turbulence production region'), at the location of peak turbulence intensity and downstream of the location of peak turbulence intensity (in the 'turbulence decay region'), for which the flow may be considered to be a fully developed turbulent flow. It is observed that the 'tear drop' shape of the joint p.d.f. gradually unfolds with distance x travelled downstream. In particular sector I is observed to broaden with x at the expense of sector II. Gomes-Fernandes, Ganapathisubramani & Vassilicos (2015) showed that in the near field of the flow past such a fractal square grid the two-point statistics revealed an inverse cascade of turbulent kinetic energy along an attractive axis inclined at some small angle relative to the streamwise direction. Buxton & Ganapathisubramani (2010) showed that sector II is the only sector within the Q-R space that can account for enstrophy attenuation due to the process of vorticity compression, i.e. ω i s ij ω j < 0, which is the inviscid mechanism for such an inverse cascade. In a fully developed flow it has been known since Taylor (1938) that on average ω i s ij ω j > 0, that is vorticity stretching (enstrophy amplification) exceeds vorticity compression on average. The broadening of sector I, which is primarily responsible for enstrophy amplification (Buxton & Ganapathisubramani 2010), at the expense of sector II as x increases is thus consistent with this cascade driven picture. Additionally, the discriminant is known to act as an attractor as R becomes large, encompassing sectors II and III leading to an elongated tail in the Q-R joint p.d.f. (Vieillefosse 1982). This was another feature that was observed to be revealed as x was increased in the study of Gomes-Fernandes et al. (2014). The 'Vieillefosse tail' is one such feature that was able to be accurately predicted by the restricted Euler analysis of Cantwell (1993).
The objective of this paper is to further probe the spatial evolution of the statistical distribution of the Q-R space, and hence the state of the VGT, in a spatially developing inhomogeneous flow. A characteristic feature of many spatially developing turbulent flows is that they contain a significant energy content in coherent dynamics due to, for example, vortex shedding in bluff body flows or Kelvin-Helmholtz instabilities in shear layers. Crucially, in the context of the VGT in turbulent flows, this coherent dynamics is a consequence of large-scale instabilities within the flow that will scale with the largest relevant/global length scales. It is thus not reasonable to consider them statistically homogeneous or isotropic as are the fine scales of turbulence. The presence of such coherent dynamics led to the introduction of a triple decomposition of the form (Hussain & Reynolds 1970)  a decomposition of u may be performed, for example phase-conditioned sampling (Hussain & Reynolds 1970), phase bin averaging (Cantwell & Coles 1983), frequency domain filtering (Brereton & Kodal 1992) or projection of the velocity field onto a modal decomposition of the flow (Perrin et al. 2007;Baj, Bruce & Buxton 2015). For a more detailed exposition on the pros and cons of these various methods the reader is referred to the introduction of Baj et al. (2015). This paper will examine the spatial evolution of the invariants of the VGTs, a ij = ∂u i /∂x j , a φ ij = ∂u φ i /∂x j and a ij = ∂u i /x j as the flow develops in space. The particular flow chosen is that past a high aspect ratio (effectively infinitely long) square cylinder, at a sufficiently high Reynolds number to ensure that the flow is adequately turbulent, which is described in § 2. Such flows have a notable peak in their spectra corresponding to the vortex shedding mode of the cylinder, accounting for a significant energy content in coherent motions. In order to isolate these motions the flow will be subjected to a triple decomposition as laid out in (1.5) and described in § 3. The results and conclusions are presented in § § 4 and 5 respectively.

Data acquisition and processing
The experiments were performed in the water tunnel of the Laboratory for Aero and Hydrodynamics at TU Delft. This facility has a cross section of 600 × 600 mm 2 . Preliminary, planar particle image velocimetry (PIV) experiments revealed that the operating condition was at a free-stream velocity of U ∞ = 0.34 m s −1 with a turbulence intensity of 0.7 %. A high aspect ratio (A = 16) square cylinder of side length D = 32 mm was carefully mounted just downstream of the contraction of the tunnel with great care having been taken to ensure that the cylinder was mounted perpendicularly to the incoming flow. This configuration yielded a Reynolds number based on the cylinder side length and free-stream velocity of Re D = 10 840. Throughout this manuscript a Cartesian coordinate system, x-y-z with corresponding velocity components U-V-W, is adopted with an origin located on the centre of the rear face of the cylinder such that the spanwise extent of the cylinder is −8D z 8D, and the cross-stream extent is −D/2 y D/2 with the rear face of the cylinder identically defined as x = 0. The spanwise-averaged (the flow may be considered as infinite in the spanwise direction since we are at the centre of the cylinder) mean, u, and root-mean-square (r.m.s.), u rms (i.e. before the triple decomposition of (1.5) is applied), streamwise velocity fields are presented in figure 1.
In order to capture the three-dimensional three-component (3D3C) data that were required to compute the statistics of the velocity-gradient tensor tomographic PIV experiments were conducted that imaged the flow immediately downstream of the cylinder in a field of view (FOV) that measured 3.8D × 2D × 0.168D. The flume was homogeneously seeded with polyamide particles of diameter 56 µm which comfortably met the requirement to act as tracers and the FOV was illuminated with the frequency doubled output of a dual cavity/double pulsed Nd:YAG laser (200 mJ pulse −1 output) that was passed through appropriate light cone-forming optics. A variable aperture slit was used to clip the light cone to produce a thickened light sheet (≈7 mm thick) with a steep illumination gradient at the edges. The FOV was imaged with four cameras (LaVision Imager Pro X) with 2048 × 2048 pixel charge-coupled device sensors with 14 bit grey-scale dynamic range, mounted to Scheimpflug adapters and lenses (Nikkor) with a focal length of f = 105 mm. The off-axis viewing angles were deliberately kept small, as were the apertures of the Invariants of the velocity-gradient tensor in a spatially developing flow lenses ( f /16), such that water prisms were not deemed necessary due to the changes in refractive index at the tunnel wall. N = 2000 (×4 cameras) image pairs were captured at an acquisition frequency of 0.76 Hz with an inter-frame time of δt = 1 ms. A convergence study showed that this was a sufficient number of images to converge the computation of the mean fluctuating components of velocity to zero. These images were processed using the multiplicative algebraic reconstruction technique (MART) tomographic PIV algorithm to produce the tomographic volume reconstruction (Elsinga et al. 2006). To improve the quality of this reconstruction the images were pre-processed with a background subtraction (from previously acquired dark images) and Gaussian smoothing using a 3 × 3 pixel filter length and volume self-calibration (Wieneke 2008) was applied. The particle displacement field was obtained from these reconstructed volumes using an iterative cross-correlation technique. The final volume size was 32 × 32 × 32 voxels and an overlap of 75 % between adjacent correlation volumes was used. This translates to a spatial resolution of 0.0486D, with a vector spacing of 0.0122D. This spatial resolution equates to approximately 11η, where η is the Kolmogorov length scale, for x/D 1.5 up to a worst case scenario of ≈17η in the separated shear layers at x/D = 0. Whilst a spatial resolution of ≈3η is generally considered necessary to resolve the smallest, dissipative length scales within a turbulent flow (Worth, Nickels & Swaminathan 2010;Buxton, Laizet & Ganapathisubramani 2011) it has been shown that resolving the characteristic 'tear drop' shape of the joint p.d.f. between Q and R is reliant upon a mix of dissipative and inertial range scales > λ, where λ is the Taylor micro-scale (Buxton 2015). The spatial resolution of the present data easily meet this criteria (it is no worse than 0.4λ and remains better than 0.15λ everywhere other than within the mean recirculation region) and thus this data set may be considered adequately spatially resolved to examine the statistics of the velocity-gradient tensor invariants.

Numerical differentiation and data post-processing
A fourth-order accurate in space, centred finite-difference scheme was used to compute the nine components of the VGT. Only the central five planes of the velocity field were retained due to deteriorating signal-to-noise ratio for the cross-correlations obtained in the volumetric PIV processing towards the edge of the thickened light sheet. Nevertheless, this was satisfactory for the implementation of the fourth-order accurate central-differencing scheme, with only the central plane being retained for  the VGT field. By necessity, a lower-order accurate numerical scheme would have been required to compute the spatial velocity gradients in the outer planes, hence they were discarded. An excellent test of the quality of a 3D3C data set from which the velocity gradients are computed is the scatter on the divergence of the velocity field. For an incompressible flow ∂u i /∂x i = 0. However, the presence of experimental noise on a data set leads to spurious, non-zero divergence.
In order to correct for this non-zero divergence the data were post-processed using the divergence correction scheme of de Silva, Philip & Marusic (2013). This consists of a nonlinear constraint based optimisation that minimally alters the measured velocity field under the constraint of restricting the magnitude of the divergence to a specified level. The tolerable divergence error was specified as |∂ũ i /∂x i | 1 s −1 , in which· denotes the application of the divergence correction scheme. The objective function that is minimised during the optimisation is effectively the ensemble average of the 'turbulent kinetic energy' that is added to the experimentally measured velocity field, i.e.k = 3 i=1 (ũ i − u i ) 2 , which was computed to bek 1/2 /U ∞ = 0.029, which is comparable to the experimental error as estimated by the method of Herpin et al. (2008). The extent of the divergence of the final, corrected velocity field is presented in figure 2 which shows the joint p.d.f. between ∂ũ/∂x and −(∂ṽ/∂y + ∂w/∂z). Henceforth, for convenience, all velocity fields will be computed from the divergence corrected data and thus the· notation is dropped.

Phase averaging and the triple decomposition
Since the data were acquired at a rate that is significantly slower than the characteristic shedding frequency of the cylinder it is necessary to compute the coherent velocity fluctuations, u φ i in (1.5), by computing/approximating a phase average. Whilst Baj et al. (2015) shows that this is inadequate for a multi-scale flow it has proven to be satisfactory for a single cylinder flow (Cantwell & Coles 1983). The determination of the phase of the instantaneous velocity fields was achieved through the use of proper orthogonal decomposition (POD) (Berkooz, Holmes & Lumley 1993). POD is a linear procedure for extracting a basis for modal decomposition of an ensemble of velocity fields, such as that acquired during the tomographic PIV experiments. These modes, Ψ j where 1 j N, are ranked according to their kinetic energy content. The first mode is steady, and corresponds to the base (or time-averaged mean) flow. van Oudheusden et al. (2005) showed that the next two most energetic modes correspond to the vortex shedding of the cylinder and that a low-order model of the flow consisting of the base flow, Ψ 2 and Ψ 3 is practically equivalent to the phase-averaged result. Figure 3 shows the proportion of total kinetic energy in the first 10 modes. It can be seen that after the base flow, mode 1, the majority of the energy is contained within modes 2 and 3 as is expected in the flow downstream of a cylinder (Perrin et al. 2007).
Despite the fidelity of a lower-order model consisting of the base flow (u) and modes Ψ 2 and Ψ 3 we choose to use a phase bin-averaging procedure similar to the approach taken by Perrin et al. (2008), which is shown to more faithfully reproduce the coherent motions. Each POD mode, Ψ j , has a corresponding time varying mode coefficient a j . Since modes Ψ 2 and Ψ 3 are orthogonal to one another a phase angle may be defined from the a 2 -a 3 plane. The scatter plot of a 2 / √ 2ξ 2 versus a 3 / √ 2ξ 3 , in which ξ j is the corresponding eigenvalue to Ψ j (not shown for brevity), shows a pattern scattered around a unit circle with its centre at the origin as expected (van Oudheusden et al. 2005). Each instantaneous velocity field is then assigned a phase angle such that The phase space is discretised into 18 phase bins centred on φ 0m (1 m 18) of size φ = 20 • . All velocity fields for which the phase angle φ ∈ (φ 0m ± φ/2) were ensemble averaged in order to compute u φ i and subsequently u i as below Note that the variable u φ i (m) is discrete since the phase averaging is conducted over bins of size φ = 20 • , however for convenience this discrete notation is dropped for the remainder of the manuscript and the coherent velocity fluctuation is henceforth simply referred to as u φ i . The r.m.s. fields of both u φ and u are presented in figure 4. Both the u φ (x) and u (x) fields were numerically differentiated with the same fourth-order accurate scheme as described in § 2.2 and the divergence of both was computed. Whilst the scatter for both ∂u φ i /∂x i and ∂u i /∂x i was higher than for ∂u i /∂x i (illustrated in figure 2) both were more than acceptable in comparison to previously published data (e.g. Ganapathisubramani, Lakshminarasimhan & Clemens 2007).
The ensemble-averaged turbulent kinetic energy is given by and thus for the triple decomposition to be energy preserving, such that the total turbulent kinetic energy is equal to the sum of that of the coherent and stochastic velocity fluctuation fields, the correlation u φ i u i must be zero. Within the field of view u φ i u i /U ∞ 2 was computed to be 1.4 × 10 −4 , which is zero to within the experimental uncertainty. Nevertheless, this small positive value is a potential explanation for the slightly less strictly observed divergence-free condition in the u φ (x) and u (x) fields, since differentiation is known to amplify any noise present.

Results and discussion
The various fluctuating velocity fields, u (x), u φ (x) and u (x) were differentiated spatially according to the fourth-order accurate scheme described in § 2.2 and the Invariants of the velocity-gradient tensor in a spatially developing flow various invariants of the velocity-gradient tensor were computed. The spatial evolution of the joint p.d.f.s of these invariants was computed along two streamwise traverses, the centreline of the flow (y = 0) and along the shear layers. At a given x location the y profiles of the r.m.s. of u have two local maxima, as indicated in figure 1(b), and thus the shear layers were defined as the locations, y(x), of these maxima. Due to symmetry it was possible to produce statistics that were averaged over both shear layers to aid with convergence. However, due to the typically observed intermittent distribution of the velocity gradients some spatial averaging was performed, in which a square window of size 0.13D × 0.13D was centred on the point of interest (i.e. on the centreline or the shear layer) in order to better converge the statistics. A sensitivity study was conducted to assess the size of this window on the ability to faithfully reproduce the Q-R joint p.d.f. and it was found to be optimal in the sense that it was the minimal window size that generated relatively noise-free statistics that were insensitive to modest increases/decreases in window size. Joint p.d.f.s were then produced of the various Q-R invariants, computed from the {a ij , a φ ij , a ij } fields, at the streamwise locations indicated by the coloured vertical lines of figure 5. In order to make a direct comparison between the joint p.d.f.s computed from the various fields of velocity gradients we wished to choose a contour level that encompasses an equivalent proportion of the overall data available, i.e. (4.1) The choice of Σ is arbitrary, but we wished to choose a contour level defining A that was sufficiently rare to ensure that we captured a broad range of Q-R states but sufficiently common to ensure a reasonably smooth contour from sufficiently converged statistics. We chose a contour level of Σ = 0.683, which is equivalent to the proportion of data bounded by ±σ (one standard deviation) for a univariate Gaussian distribution. The invariants themselves are all normalised by the local value of Q ω , computed from the relevant velocity-gradient field such that the joint p.d.f.s of figure 6(c), for example, are normalised by Q ω calculated from the a φ ij field at the centreline and appropriate x-location etc. as illustrated in figure 5.
It can be seen that the classical 'tear drop' shape for the joint p.d.f. is recovered when calculated from the a ij field. This in itself is a surprising result since the 'tear drop' shape of the Q-R joint p.d.f. is considered to be associated with fully developed turbulent flows. Although there is no clear definition of such a flow they may be generally considered to be close to homogeneous and isotropic (requiring the integral length scale to be smaller than a relevant homogeneity length scale), in the small scales at least, with mean velocity profiles that may be collapsed when scaled by a relevant flow variable such as the centreline velocity for a wake. Evidently these criteria are far from being met in the near-field flow behind a square cylinder. It is clear from inspection of figure 1 that throughout the field of view, but particularly for the cases of x/D 2, the flow is far from homogeneous. Further, Invariants of the velocity-gradient tensor in a spatially developing flow 11 the 'tear drop'-shaped contours are equally visible, and quantitatively similar (with the exception of the x/D = 3.5 cases) on both the centreline and the shear layers in figure 6(a,b). This is despite the fact that the flow physics for both of these regions differ significantly. Figure 1 shows that along the centreline the shape of the Q-R joint p.d.f. contours are similar regardless of whether the data are extracted from the mean circulation region (for the centreline) or downstream of this in the turbulence decay region (x/D = 2). Perversely, x/D = 3.5 is perhaps the location where one might assume that the criteria for a fully developed turbulent flow, for which the 'tear drop'-shaped joint p.d.f. is considered to be a universal feature, are more stringently met than any of the streamwise locations further upstream. Despite this, the joint p.d.f.s from this location are the only ones not to collapse when normalised by Q ω . The clearly defined 'tear drop'-shaped contours of the joint p.d.f.s of figure 6(a,b) are in contrast to the spatially developing fractal square grid flow of (Gomes-Fernandes et al. 2014), in which the 'tear drop' shape was observed to unfold as the flow develops downstream, in particular the 'Vieillefosse tail' and contribution of sector I became increasingly significant whereas the contribution of sector II was reduced.
Comparison of figure 6(a,b) shows that the collapse of the contours of f QR (Q, R) when scaled with the local quantity Q ω is marginally better along the shear layer than the centreline, excluding the data from x/D = 3.5. There is significant energy content in coherent motions in the shear layers, as illustrated in figure 4(a), whereas along the centreline there is virtually none with a significant contribution to the total turbulent kinetic energy from the u i fluctuations. It is thus a surprising result that the collapse of the contours of f QR (Q, R) is better where there is a significant contribution from coherent motions as opposed to primarily the stochastic fluctuations. Nevertheless, the similarity between the contours of figure 6(a,b) is stark, despite the fact that they are computed from regions with very different flow physics to one another. The streamwise location at which the joint p.d.f. most resembles that for fully developed turbulence, with an enhanced contribution from sector I and elongated 'Vieillefosse tail' is at the location of peak turbulence intensity, x/D = 1.104. Figure 6(c,d) shows the equivalent joint p.d.f.s between Q and R computed from the coherent velocity-gradient field, a φ ij , along the centreline and shear layer, respectively whilst (e) and ( f ) show those computed from the stochastic velocity-gradient field, a ij . It is clear that the contours of f QR (Q, R) computed from the a φ ij field do not exhibit anything remotely akin to a 'tear drop' shape, whilst those computed from the a ij field do, and are indeed quantitatively very similar to those computed from the a ij field. If anything, it may be commented that the contours show a better collapse with the local Q ω scaling for the joint p.d.f. computed from the a ij field than from the a ij field. Again, however, the notable exception is the contour extracted from the data at the furthest downstream location, x/D = 3.5, at which we may have expected the flow to best approximate a 'fully developed flow'.
We may thus conclude from figure 6 that in this spatially developing inhomogeneous flow the 'tear drop' shape of the contours of the joint p.d.f. between Q and R is almost entirely due to the stochastic turbulent fluctuations. Even though figure 5 shows that the coherent motions have a non-negligible Q ω , indicative of the average magnitude of the a φ ij tensor, it does not contribute to the kinematics of the overall velocity-gradient tensor through the Q-R joint p.d.f. It thus appears that the 'tear drop' is more ubiquitous than first thought since it appears in flows with significantly varied physics (fully developed/spatially developing/recirculation/separated shear layers etc.) in an observation that is comparable to the 'embarrassment of success'  (Kraichnan 1974). Kolmogorov (1941) predicts the existence of a −5/3 slope of the energy spectrum in homogeneous, isotropic turbulence for which the Reynolds number is sufficiently high to support an inertial range of scales separating the dissipative and large-scale motions. Equation (1.4) shows that the invariant R is the local excess of strain-rate production over enstrophy amplification, ω i s ij ω j . The amplification of enstrophy, through its interaction with the strain-rate field, is the only known inertial mechanism by which enstrophy (and turbulent kinetic energy) are transferred (on average) to smaller scales through the turbulent cascade. Contrastingly, equation (1.3) shows that the invariant Q is representative of the localised topology of the flow, whether rotationally or strain-rate dominated. Numerous studies (e.g. Ruetsch & Maxey 1991) have painted a consistent topological picture of developed, fine-scale turbulence in which sheets of dissipation (high magnitude strain rate) are wrapped around 'worms' of intense enstrophy. Yet figure 6 clearly shows 'tear drop'-shaped distributions of Q and R in a flow that is highly inhomogeneous/anisotropic with wildly varying flow topologies, mirroring the finding of a −5/3 slope in flows at modest Reynolds numbers that are also far from homogeneous/isotropic (e.g. Valente & Vassilicos 2012;Gomes-Fernandes et al. 2015;Laizet, Nedić & Vassilicos 2015). Figure 6 illustrates contours of the joint p.d.f.s between the invariants Q and R at several illustrative downstream locations. However, the nature of the tomographic PIV data that have been acquired allows a closer inspection of the spatial development of the statistical distribution of these invariants. This may be quantified by computing the proportions of the total data located within various sectors of the joint p.d.f. The majority of the total data are locally swirling, i.e. ∆ > 0, and thus for simplicity we focus on the proportion of total data located in sectors I, primarily responsible for enstrophy amplification, and II for which there is mean enstrophy attenuation (Buxton & Ganapathisubramani 2010) as follows Figure 7 shows the spatial evolution of the ratios R I and R II , as defined in (4.2) and (4.3) respectively, from f QR (Q, R) computed from 0.13D × 0.13D square windows as before along the centreline (a,c,e) and the shear layer (b,d,f ). Firstly, it is observed that for a particular velocity-gradient field {a ij , a φ ij , a ij } there is little spatial variation in the ratios R I , R II . This reinforces the observation drawn from figure 6 that despite the rapid changes in flow physics in the near-field flow past a square cylinder the statistical distribution of the VGT invariants remains remarkably constant. For all six panels it is also observed that R I and R II are approximately equal to one another, indicating that the volume of space occupied by enstrophy amplification and enstrophy attenuation are similar, although Buxton & Ganapathisubramani (2010) shows that the magnitude of enstrophy amplification is greater in sector I than that for attenuation in sector II. It is observed that there is an extremely high correlation between the spatial variation of R I and R II computed from the a ij field and the a ij field, along both the shear layers and the centreline, reinforcing the conclusion that we may draw from figure 6, namely that the kinematics of the VGT invariants are driven by the stochastic velocity fluctuations. Whilst the variation Invariants of the velocity-gradient tensor in a spatially developing flow in magnitude may be small for the spatial variation of R I between the centreline and shear layers in figure 7(c,d) it can be seen that they are inversely correlated to one another. Up until x/D = 2 R I increases along the centreline, at the expense of R I decreasing along the shear layer. Since sector I is the sector primarily responsible for enstrophy production this suggests that enstrophy associated with coherent motions is increasingly produced along the centreline at the expense of its production along the shear layers. Evidently, the coherent motions originate from the shear layers and hence the vortex stretching process responsible for the mean cascade of enstrophy acts to homogenise the enstrophy transported in coherent motions across the flow. We may conclude, qualitatively from figure 6 and quantitatively from figure 7, that the statistical behaviour of the invariants of the VGT does not vary significantly in space along either the centreline or the shear layer despite the rapidly changing flow physics. To quantitatively illustrate the spatially developing nature of the flow physics, specific to the VGT itself, figure 8 presents the joint p.d.f.s of the constituent components of the invariant Q; Q ω is the first term on the right-hand side of (1.3) and Q S is the second term on the right-hand side of (1.3). The various joint p.d.f.s are again normalised by the local (and field-specific) values of Q ω and the data are extracted along the centreline and shear layers at the same x locations as with the The topology of the fine-scale turbulent motions can be inferred from such joint p.d.f.s. The −Q s ∼ s ij s ij axis represents points with a locally high rate of dissipation of turbulent kinetic energy which are known to be arranged in sheet-like structures (e.g. Ganapathisubramani, Lakshminarasimhan & Clemens 2008) whereas the Q ω ∼ ω i ω i axis is associated with points with a locally high enstrophy which are arranged in tube-like structures in turbulent flows (e.g. Jiménez et al. 1993). The 45 • line of Invariants of the velocity-gradient tensor in a spatially developing flow 15 −Q s = Q ω represents a balance of enstrophy and dissipation and these points may be somewhat vaguely described as vortex sheets. Nevertheless, it has been shown that a very large proportion of the data for fully developed turbulent mixing layers lies along this line (Soria et al. 1994). Indeed of all the contours of figure 8(e,f ) those extracted from the location of peak turbulence intensity, x/D = 1.104 (the red contours), look the most similar to those computed from fully developed turbulence. Despite the similarity of the flow topology for all of the VGT fields, as illustrated in figure 8, the joint p.d.f.s of figure 6(c,d) are vastly different to the classical 'tear drop' shape for the joint p.d.f. between Q and R.
A particularly clear preference for high enstrophy, tube-like, topology exists in the downstream locations of the a φ ij field, both along the shear layers and the centreline (where admittedly the energy content of u φ (x) is low). This tendency develops with x since the topology of a φ ij shows a tendency for vortex sheets further upstream. In this instance the very different flow physics/topologies at the various downstream locations presented in figure 8(d) do translate into quantitatively different joint p.d.f.s between Q and R in figure 6(d), which are far from 'tear drop' shaped. Contrastingly, it may be observed that (with the exception of x/D = 3.5) the joint p.d.f.s between Q and R computed from the a ij field, figure 6(a,b), collapse onto self-similar 'tear drop' shapes despite the fact that no such collapse is observed in the joint p.d.f.s between Q ω and −Q s in figure 8(a,b). This is true along the centreline of the flow, in which there is initially a recirculation region before spatially developing into a turbulent wake flow, as well as along the shear layers. Along the shear layers there is a clear preference for high enstrophy structures to develop with increasing x, whereas along the centreline there are very different distributions as the flow transitions from a mean recirculation to a vortex street. Nevertheless, all of these regions produce 'tear drop'-shaped joint p.d.f.s between Q and R.
We may exhaustively write out the invariants Q and R, in terms of the components a ij , as follows Q = −(a 2 11 + a 2 22 + a 11 a 22 ) − (a 12 a 21 + a 13 a 31 + a 23 a 32 ) (4.4) R = (a 11 a 22 − a 12 a 21 )(a 11 + a 22 ) + (a 11 a 23 a 32 + a 22 a 13 a 31 ) − (a 12 a 23 a 31 + a 13 a 32 a 21 ). (4.5) We may thus write Q and R as in which Q φ and R are invariants computed from the a φ ij and a ij fields etc. and Q X and R X are cross-terms that involve contributions from both the a φ ij and a ij fields. In general these are defined as Q X = 2(a φ 11 a 11 + a φ 22 a 22 ) + a 11 a φ 22 + a φ 11 a 22 − a 12 a φ 21 − a φ 12 a 21 − a 13 a φ 31 − a φ 13 a 31 − a 23 a φ 32 − a φ 23 a 32 (4.8)  Our particular flow is well approximated as being statistically two-dimensional meaning that there are no gradients of the coherent motions in the z-direction nor a velocity component w φ hence (4.8) and (4.9) may be simplified to Q X 2D = 2(a φ 11 a 11 + a φ 22 a 22 ) + a 11 a φ 22 + a φ 11 a 22 − a 12 a φ 21 − a φ 12 a 21 (4.10) R X 2D = (a 11 a φ 22 + a φ 11 a 22 − a 12 a φ 21 − a φ 12 a 21 )(a φ 11 + a 11 + a φ 22 + a 22 ) + a φ 11 a 23 a 32 − a φ 12 a 23 a 31 − a φ 21 a 13 a 32 + a φ 22 a 13 a 31 . (4.11) Nevertheless, to avoid further assumptions we shall persist with the definitions of Q X and R X from (4.8) and (4.9), cognisant of the fact that terms such as a φ 13 ≈ 0. Figure 9 illustrates the joint p.d.f.s between Q X and R X , as defined in (4.8) and (4.9) respectively, at the various streamwise locations highlighted in figure 5 along the centreline (a) and the shear layer (b). The joint p.d.f.s are computed from the same 0.13D × 0.13D interrogation windows as before, the contour levels are again chosen such that Σ = 0.683 and the axes are normalised by Q ω (x). Whilst the shape of the contours does not exactly resemble the classical 'tear drop' shape of f QR (Q, R), in particular the Vieillefosse tail is absent, some features are present. These include the predominance of states for which ∆ > 0 (swirling), the enhancement of sector I (primarily responsible for enstrophy amplification) and the narrowed sector II (primarily responsible for enstrophy attenuation). Further, it can be seen that all of the contour levels collapse well with each other, particularly from the rear face of the cylinder up until the location of the peak turbulence intensity.
The shape of the joint p.d.f. contours of figure 9 look remarkably similar to f QR (Q, R) extracted from the closest location to the turbulence generating grid of Gomes-Fernandes et al. (2014). It was in this region of the flow that Gomes-Fernandes et al. (2015) observed an inverse cascade of turbulent kinetic energy, along a particular axis, which was partly attributed to the multi-scale nature of the flow generated by their particular space-filling fractal square grid. This region of the flow was termed the turbulence production region since u rms was an increasing function of x. In such a flow the coherent dynamics of the bars of various length f.s between −s ij s jk s ki , the strain self-amplification rate, and ω i s ij ω j , the enstrophy amplification rate along the centreline of the flow computed from the a φ ij field (a) and the a ij field (b). The dashed line marks ω i s ij ω j = −4/3s ij s jk s ki and thus data from below the line correspond to R < 0 whilst the data above the line correspond to R > 0. The contour colours correspond to the streamwise locations depicted in figure 5.
scales interacted with each other, and the background incoherent turbulence. It may thus be hypothesised that the cross-terms, Q X and R X , were important in this region of the flow which perhaps explains the resemblance of f QR (Q, R) to figure 9.
In the transport equations for ω 2 and s ij s ij the inviscid source terms are ω i s ij ω j (in the enstrophy equation) and −s ij s jk s ki and −ω i s ij ω j (in the strain-rate equation).
Clearly ω i s ij ω j is thus a source of enstrophy and a sink for dissipation, and these two terms are the constituents of R as outlined in (1.4). Figure 10 shows joint p.d.f.s, again with contour levels set such that Σ = 0.683, between ω i s ij ω j and −s ij s jk s ki computed along the centreline from the a φ ij field (a) and the a ij field (b). The remarkable collapse of the joint p.d.f.s of figure 6 can perhaps be explained by the self-preserving nature of the statistical distributions of these inviscid source/sink terms, which drive the evolution of enstrophy and dissipation. With the exception of the contour taken immediately downstream of the cylinder for the coherent velocity-gradient tensor field, in which Q ω is small, the contours collapse onto one another. There is a clear qualitative difference in the distributions for both. There is a preference for the production of enstrophy alongside mild (inviscid) production of dissipation in the field of stochastic velocity fluctuations. Intuitively, this would be required for an inhomogeneous flow to evolve into a 'fully developed' flow, in which the fine-scale structure is observed to consist of tubes of intense enstrophy surrounded by more extensive sheets of dissipation. The opposite trend is observed in the field of coherent fluctuations, in which the inviscid production of dissipation is favoured, which again one might expect considering the decreasing significance of the coherent motions with streamwise distance in a turbulent shear flow.

Conclusions
Tomographic PIV experiments were performed in the near field of the turbulent flow around a square cylinder generating a 3D3C data set. The velocity field was decomposed, according to (1.5), into a time invariant base flow, a coherent fluctuation and a residual/stochastic fluctuation. This was conducted by means of phase bin averaging in which the phase angle of each instantaneously acquired 3D3C velocity field was calculated from the time-varying mode coefficients of the second and third modes obtained from a proper orthogonal decomposition of the flow field. These two modes correspond to the vortex shedding downstream of the cylinder. From these three separate velocity fields the invariants of the velocity-gradient tensor were computed and their statistics compared to one another. Since the flow is spatially developing and inhomogeneous these statistics were computed both upstream, including from within the mean recirculation region, and downstream of the location of peak turbulence intensity along the centreline of the flow. Additionally, the statistics were extracted from the shear layer region of the flow; so defined as the local maximum in turbulence intensity at a given downstream, x, location. Joint p.d.f.s between the rotating and straining constituents of the invariant Q, Q ω and Q s respectively, highlighted the very different flow physics and topologies of the various velocity-gradient fields at the various spatial interrogation locations. Despite the differing flow physics and topologies of the various fields of velocity gradients the joint p.d.f.s between the second and third invariants, Q and R resembled the classical 'tear drop' shapes associated with fully developed turbulence when computed from the total and stochastic velocity gradients. In fact the contours of the joint p.d.f.s were found to collapse in a self-similar fashion when appropriately normalised by the local ensemble average of Q ω . Such self-similar 'tear drop' shapes were found along the centreline and the shear layer at all but the farthest downstream locations. Observation of the relative magnitudes of the rotationally dominated sectors I and II showed that these tended to be in balance with one another and only very slowly variant in space, although their values did change according to which velocity-gradient field they were computed from. The statistics of the inviscid source/sink terms in the evolution of the enstrophy and dissipation are observed to also collapse remarkably well with the local value of Q ω , offering a potential explanation for the collapse of the 'tear drop'shaped joint p.d.f.s. As expected the joint p.d.f.s computed from the coherent velocity gradients did not resemble those from fully developed turbulence.
It thus seems that the ubiquitous 'tear drop'-shaped statistical distribution of the invariants of the velocity-gradient tensor is even more ubiquitous than first thought. Not only does it appear in a manner of (homogeneous) fully developed turbulent flows, but this manuscript also shows that it exists in inhomogeneous flows with rapidly spatially varying flow physics. Whilst fully developed flows have a flat spectrum of turbulent kinetic energy at the low wavenumbers the flow past a square cylinder has a clear spike at the vortex shedding wavenumber. From this perspective we should thus consider the 'tear drop'-shaped joint p.d.f. between the invariants of the velocitygradient tensor not as a universal feature of fully developed turbulence (Chacin & Cantwell 2000;Elsinga & Marusic 2010) but as a universal feature of turbulent flows.