Low-order moments of the velocity gradient in homogeneous compressible turbulence

Abstract We derive from first principles analytic relations for the second- and third-order moments of $\boldsymbol{\mathsf{m}}$, the spatial gradient of fluid velocity $\boldsymbol{u}$, $\boldsymbol{\mathsf{m}} = \nabla \boldsymbol{u}$, in compressible turbulence, which generalize known relations in incompressible flows. These relations, although derived for homogeneous flows, hold approximately for a mixing layer. We also discuss how to apply these relations to determine all the second- and third-order moments of the velocity gradient experimentally for isotropic compressible turbulence.


Introduction
In high Reynolds number flows, the velocity field, u, develops very sharp gradients (Frisch 1995;Sreenivasan & Antonia 1997), resulting in extremely large fluctuations of the velocity gradient m, or m ij = ∂ui ∂xj .For this reason, an accurate description of the velocity gradient is essential to understand the small-scale properties of turbulence (Meneveau 2011).In the case of incompressible turbulence, much emphasis has been put on enstrophy, defined as 1 2 ω i ω i , where ω = ∇ × u.Its amplification rate, known as vortex stretching, ω i s ij ω j , where s is the symmetric part of m, or the rate of strain tensor, have been thoroughly studied to investigate the production of small scales in the flow (Tsinober 2009;Buaria et al. 2020).It should be kept in mind that a thorough description of small scales involves the full tensor m, and not just vorticity (Meneveau 2011).
Two remarkable constraints on the second and third moments of m have been established by Betchov (1956) for homogeneous and incompressible flows: tr(m 2 ) = 0 and tr(m 3 ) = 0, (1.1) where denotes ensemble average.In addition, when the flow is isotropic, the identities (1.1) lead to remarkable simplifications, allowing the second and third order moments of the velocity gradient tensors A (2) = m 2 and A (3) = m 3 to be expressed in terms of only one scalar quantity (Pope 2000).This implies that in an isotropic flow, A (2)  and A (3) can be completely determined from the measurement of only one component of the velocity gradient tensor, e.g., using hot-wire probes.Although strictly zero in homogeneous incompressible flows, the values of tr(m 2 ) and tr(m 3 ) remain very small even when spatial inhomogeneity is strong, such as in channel flows (Bradshaw & Perot 1993;Pumir et al. 2016).Effectively, this can be understood as a consequence of the slow variation of the mean flow properties, compared to the very fast variation of turbulence at small scales.Nonetheless, the structure of the velocity gradient tensors, and in particular of A (2) and A (3) , strongly deviate from the isotropic case (Bradshaw & Perot 1993;Vreman & Kuerten 2014;Pumir 2017).
Here, we focus on the tensors A (2) and A (3) in compressible flows.How compressibility affects the structure of the velocity gradient tensor has been studied in homogeneous isotropic flows (Pirozzoli & Grasso 2004;Wang et al. 2012;Fang et al. 2016;Wang et al. 2018), in homogeneous shear flows (Ma & Xiao 2016;Chen et al. 2019), in mixing layers (Vaghel & Madnia 2015) and in boundary layers (Chu et al. 2014).The dynamics of the velocity gradient tensors in compressible flows have been studied by Suman & Girimaji (2009, 2011, 2013) starting from the homogenized Euler equation.
The purpose of this work is to generalize the Betchov relations, Eq. (1.1), to homogeneous compressible flows, which allows us to establish the general structure of A (2)  and A (3) similar to the case of incompressible flows.We validate these relations with direct numerical simulations (DNS) of homogenous isotropic compressible turbulence, and demonstrate that they still approximately hold in flows with strong inhomogeneity, i.e., in a turbulent mixing layer.Whereas only one scalar was sufficient to capture the full structure of A (2) or A (3) in incompressible flows, for compressible turbulence, 2 independent parameters for A (2) and 4 parameters for A (3) are required.Nonetheless, as we show, these relations can be used to construct the full tensors of A (2) and A (3) from stereo-PIV measurements in compressible flows.

Second-and third-order relations of compressible turbulence
Following Betchov (1956), we first consider the second-order moment A (2) : iijj , (2.1) in which we used the notation X = tr(X) and the summation convention, and assumed homogeneity.For the third moments, we notice a general relation between the traces of third-order moments involving the gradients of 3 homogeneous vector fields, Eq. 4.7 in the Appendix, which yields, for the special case where the 3 fields are identical, iijjkk . (2.2) A straightforward consequence of Eqs.(2.1) and (2.2) can be expressed by introducing s ≡ (m + m T )/2 − (m/3)I and w ≡ (m − m T )/2: Equations formally similar to Eqs. (2.1)-(2.4)were derived by Yang et al. (2020) for the perceived velocity gradient tensor based on regular tetrahedra in incompressible homogeneous turbulence (see their Eqs.(3.11), (3.12), and (3.37), (3.38)).

DNS of compressible turbulence
To test the relations presented above in various turbulent flow configurations, we numerically solved the three-dimensional compressible Navier-Stokes equations: in which ρ is the fluid density, P is the pressure, E = 1 2 ρu i u i + P/(γ − 1) is the total energy with γ = 1.4 being the ratio of specific heats, ∂x k δ ij is the viscous stress tensor with the effect of bulk viscosity neglected (Pan & Johnsen 2017), and Q j = κ ∂T ∂xj is the heat flux, where µ is the viscosity coefficient determined by the temperature, T , via Sutherland's law, and κ is the thermal conductivity.The temperature is related to fluid pressure and density via the ideal gas law P = ρRT and the gas constant R. The set of equations (3.1) are solved with high-order finite difference method.Namely, the convection terms are computed by a seventh-order low-dissipative monotonicity-preserving scheme (Fang et al. 2013) in order to capture shock-waves in a compressible flow while preserving the capability of resolving small-scale turbulent structures.The diffusion terms are computed by a sixth-order compact central scheme (Lele 1992) with a domain decoupling scheme for parallel computation (Fang et al. 2019).The time-integration is performed by a three-step third-order total variation diminishing Runge-Kutta method (Gottlieb & Shu 1998).The flow solver used in the present study is ASTR, an in-house code previously tested in DNS of various compressible turbulent flows with and without shock-waves (Fang et al. 2013(Fang et al. , 2014(Fang et al. , 2015(Fang et al. , 2020)).
We first discuss the results of decaying isotropic compressible turbulence.The computational domain is a (2π) 3 cube with periodic boundary conditions in all three directions.The initial flow is a divergence-free random velocity field with a sharply decaying spectrum: E(k) = Ak 4 e −2k 2 /k 2 0 , which peaks at k = k 0 = 4, and the value of A determines the kinetic energy at time t = 0.The density, pressure and temperature are all initialized to constant values, similar to the "IC4" run of Samtaney et al. (2001).The initial turbulent Mach number based on the root-mean-square velocity u ′ = √ u i u i Here, we use the initial large-eddy-turnover time τ 0 = ( ∞ 0 E(k)/kdk)/u ′ to normalize time when quantifying the flow evolution.The computational domain is discretized with a 512 3 uniform grid, leading to a resolution down to the dissipation scale over the entire simulation duration: the ratio of the Kolmogorov scale η over the grid size ∆ grows from η/∆ = 0.98 at t = 0 to η/∆ = 2.88 at t = 20.
Figure 1a shows that the turbulent Mach number M t , the turbulent kinetic energy K = 1 2 ρu i u i , and the Reynolds number R λ , all decay monotonically with time.Note that even at time t/τ 0 > 10, M a t 1, the flow is still highly compressible with a large number of spatially distributed shocklets.The skewness of the longitudinal velocity derivative, S = (∂u i /∂x i ) 3 / (∂u i /∂x i ) 2 3/2 shown in Fig. 1b, however, develops a sharp negative peak at t/τ 0 ≈ 0.35 and then gradually returns to a value fluctuating around −2, consistent with Samtaney et al. (2001).The larger magnitude of S in the later stage compared with S ≈ −0.6 in Samtaney et al. (2001) is most likely due to the higher grid resolution in this work.In fact, Wang et al. (2012) demonstrated in DNS of forced compressible turbulence that for a given value of R λ and M t , the magnitude of S increases when the grid is refined.We note here that fully resolving the shocklets is not feasible as their sizes are out-of-reach of the current DNS with a fixed grid.The peak of S is a consequence of compressibility.In the inset of Fig. 1b, we show the evolution of the turbulence dissipation rate, separated into the solenoidal (enstrophy) part ε s ≡ µω i ω i and the dilational part ε d ≡ 4 3 µm 2 .The solenoidal dissipation ε s grows first due to the build-up of turbulent structures from the initial random field, reaches a peak at t/τ 0 ≈ 0.8, then decreases gradually due to the decay of the turbulent fluctuation.The dilational dissipation rate, ε d , represents the contribution of compressibility to dissipation, which is exactly zero in incompressible turbulence.With our choice of solenoidal initial condition, ε d starts at zero and first grows rapidly, as shocklets are forming (Samtaney et al. 2001).The peak of ε d is reached at t/τ 0 ≈ 0.65, slightly earlier than ε s .In our simulation, ε d contributes to approximately half the energy dissipation at the peak, but the solenoidal part remains the major cause of dissipation at later times.
We now discuss the structure of the velocity gradient correlations in this flow.Figure 2 shows the evolution of the l.h.s. and r.h.s. of Eq. ( 2.1) and (2.2), the identities for the invariants of A (2) and A (3) , respectively.The ratios between the two sides of those equations, shown in the insets, are exactly 1, as predicted for homogeneous flows.
Our analysis predicts, for homogeneous and isotropic turbulence, further relations among the components of A (2) and A (3) .Figure 3a shows that A (2) 1122 and A (2) 1221 are equal during the entire simulation, as predicted by Eq. (2.6).In Fig. 3b, the ratio of the components A (2) 1212 /A (2) 1111 also lies within 1/3 and 2 as predicted.It starts at 2 since the flow is initially solenoidal.The rapid drop to values lower than 1 is concurrent with the rise of dilational dissipation ε d .The minimal possible value of 1/3 for A Figure 4a shows the evolution of the invariants of A (3) after the early stage when the flow rapidly adjusts in response to the initial divergence-free condition.Remarkably, the magnitude of w 2 m is very small compared with the others, consistent with the observation that in compressible homogeneous turbulence, vorticity and dilation are nearly uncorrelated (Erlebacher & Sarkar 1993;Wang et al. 2012).This, in view of Eq. (2.18), provides an additional relation: 3a 2 ≈ 3a 3 + 2a 4 − 2a 5 , which then leaves only 3 independent quantities for a complete determination of A (3) .Plugging Eq. (2.14) to (2.17) into this relation leads to Figure 4b This result may be helpful to reconstruct complete expression ipjqkr from planar PIV or 2C LDV data, instead of requiring stereo-PIV or 3C LDV.
Next we perform a DNS of a planar compressible mixing layer, formed by two comoving free streams with Mach numbers M a 1 = U 1 /c 1 = 7.5 and M a 2 = U 2 /c 2 = 1.5, respectively, in which U 1 and U 2 are the mean velocities in the upper and lower free streams, and c 1 and c 2 are the speeds of sound in the two streams.Similar with the setup of Li & Jaberi (2011), at the inflow plane x = 0, the mean stream velocity profile U (0, y, z) is specified to be U (0, y, z) = 1 2 [U 1 + U 2 + (U 1 − U 2 ) tanh (2y/δ ω0 )] , with the inflow vorticity thickness δ ω0 = 1.The convective Mach number is M a c = (U 1 −U 2 )/(c 1 + c 2 ) = 3 and the Reynolds number is Re c = ρ 1 (U 1 − U 2 )δ ω0 /µ 1 = 3500.The mean temperature profile at the inlet is given by the Crocco-Busemann law with a uniform mean pressure.The effective size of the computational domain is 450δ ω0 × 100δ ω0 × 16δ ω0 in the x, y and z directions, respectively.The domain is discretized using a mesh of 3050 × 400 × 128 nodes uniformly distributed in the x and z directions and stretched in the y direction with higher resolution in the centre of the mixing layer.Near the outflow plane, a sponge layer with a highly stretched mesh is added to damp fluctuations near the boundary.Random velocity fluctuations are superposed on the mean profile at the inlet plane to trigger turbulence, which develops downstream, forming a large number of shocklets in both upper and lower parts of the mixing layer.When it is not too close to the inlet plane, the momentum thickness, θ, grows linearly downstream and the mean velocity profiles are self-similar, i.e., U (x, y) at different x locations collapse to U (ξ) with ξ = (y − y c )/θ, where y c is the centre of the mixing layer.The result is validated by checking the balance of turbulence kinetic energy budget (not shown).
Figure 5 shows the ratio between the l.h.s. and the r.h.s. of Eq. (2.1) and (2.2) in the center region of the mixing layer (−6 y/θ 6) at the downstream locations xδ ω0 = 350 (dashed lines) and x/δ ω0 = 400 (solid lines).Although the flow is not homogeneous, the ratios m 2 / m 2 and m 3 /( 3 2 m 2 m − 1 2 m 3 ) are very close to unity, so Eqs.(2.1) and (2.2) are still approximately valid.Fig. 6 shows the profiles of various third-order invariants m 3 , w 2 m , s 2 m , s 3 and wsw .Despite the inhomogeneity and anisotropy, the vorticity-dilation correlation w 2 m remains very small compared to other invariants, which could help to obtain more relations among invariants.

Concluding Remarks
In summary, we derived exact relations among invariants of the moments of velocity gradients in compressible homogeneous turbulence and verified these relations by DNS of decaying compressible turbulence.Interestingly, these relations, derived under homogeneity assumptions, hold approximately in a compressible mixing layer.We also devised approaches to determine the full tensor with a minimal set of measurements in experiments.These relations could help, e.g., to determine separately the solenoidal and the dilational energy dissipation rates from only two velocity derivatives ∂u1 ∂x1 and ∂u2 ∂x1 .In the future, it would be interesting to investigate the structure of the velocity gradient implied by these relations, as done by Betchov (1956) in which h a = ∇a etc.For divergence-free fields, this yields ∇a∇b∇c = 0.An equivalent form of this special case has been shown in Appendix D of Eyink (2006).

Figure 1 .
Figure 1.(a) Evolution of the turbulent Mach number Mt, the turbulence kinetic energy normalized with its initial value K(t)/K(0), and the Reynolds number R λ (shown in inset).(b) The evolution of the skewness of the longitudinal velocity derivative S. The inset shows the energy dissipation rates: the solenoidal part εs and the dilational part ε d .isM a t = u ′ / c =2.0, where c = √ γRT is the speed of sound.The initial Reynolds number based on the Taylor microscale λ = u ′ / ( ∂ui ∂xi ) 2 1/2 is R λ = ρ u ′ λ
later times, indicates instead the prevalence of the solenoidal part.

Figure 4 .
Figure 4. (a) Evolution of the invariants m 3 , w 2 m , s 2 m , and wsw , all normalized by s 3 .(b) The values of the l.h.s. and r.h.s. of Eq. (3.2) and their ratio (shown in inset).

Figure 5 .Figure 6 .
Figure 5. Approximate validity of the homogeneous relations in compressible mixing layer.(a) Normalized second-order invariants, m 2 / mm T , m 2 / mm T , and their ratio (shown in the inset).(b) Normalized third-order invariants, m 3 / m 2 m T , (32 m 2 m − 1 2 m 3 )/ m 2 m T , and their ratio (shown in the inset).In both plots, dashed lines correspond to x/δω0 = 350 and solid lines to x/δω0 = 400.
for incompressible turbulence.