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.

1. Introduction

In high-Reynolds-number flows, the velocity field, $\boldsymbol {u}$, develops very sharp gradients (Frisch 1995; Sreenivasan & Antonia 1997), resulting in extremely large fluctuations of the velocity gradient $\boldsymbol{\mathsf{m}}$, or ${\mathsf{m}}_{ij} = {\partial u_i}/{\partial x_j}$. 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 $\frac {1}{2} \omega _i \omega _i$, where ${\boldsymbol \omega } = \boldsymbol {\nabla } \times \boldsymbol {u}$. Its amplification rate, known as vortex stretching, $\omega _i {\mathsf{s}}_{ij} \omega _j$, where $\boldsymbol{\mathsf{s}}$ is the symmetric part of $\boldsymbol{\mathsf{m}}$, or the rate of strain tensor, have been thoroughly studied to investigate the production of small scales in the flow (Tsinober 2009; Buaria, Pumir & Bodenschatz 2020). It should be kept in mind that a thorough description of small scales involves the full tensor $\boldsymbol{\mathsf{m}}$, and not just vorticity (Meneveau 2011).

Two remarkable constraints on the second- and third-order moments of $\boldsymbol{\mathsf{m}}$ have been established by Betchov (1956) for homogeneous and incompressible flows:

(1.1a,b) \begin{equation} \langle \mathrm{tr}(\boldsymbol{\mathsf{m}}^2 )\rangle = 0 \quad {\rm and} \quad \langle \mathrm{tr}(\boldsymbol{\mathsf{m}}^3) \rangle = 0,\end{equation}

where $\langle \, \rangle$ denotes ensemble average. In addition, when the flow is isotropic, the identities (1.1a,b) lead to remarkable simplifications, allowing the second- and third-order moments of the velocity-gradient tensors $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{A}}^{(3)}$, defined as ${\mathsf{A}}^{(2)}_{ipjq} = \langle {\mathsf{m}}_{ip} {\mathsf{m}}_{jq} \rangle$ and ${\mathsf{A}}^{(3)}_{ipjqkr} = \langle {\mathsf{m}}_{ip} {\mathsf{m}}_{jq} {\mathsf{m}}_{kr}\rangle$, to be expressed in terms of only one scalar quantity (Pope 2000). This implies that in an isotropic flow, $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{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 $\langle \mathrm {tr}(\boldsymbol{\mathsf{m}}^2 )\rangle$ and $\langle \mathrm {tr}(\boldsymbol{\mathsf{m}}^3) \rangle$ remain very small even when spatial inhomogeneity is strong, such as in channel flows (Bradshaw & Perot 1993; Pumir, Xu & Siggia 2016). Effectively, this can be understood as a consequence of the slow variation of the mean flow properties, compared with the very fast variation of turbulence at small scales. Nonetheless, the structure of the velocity-gradient tensors, and in particular of $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{A}}^{(3)}$, strongly deviate from the isotropic case (Bradshaw & Perot 1993; Vreman & Kuerten 2014; Pumir 2017).

Here, we focus on the tensors $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{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, Wang & Lu 2014). The dynamics of the velocity-gradient tensors in compressible flows has 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 (1.1a,b) to homogeneous compressible flows, which allows us to establish the general structure of $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{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 $\boldsymbol{\mathsf{A}}^{(2)}$ or $\boldsymbol{\mathsf{A}}^{(3)}$ in isotropic incompressible flows, for compressible turbulence, two independent parameters for $\boldsymbol{\mathsf{A}}^{(2)}$ and four parameters for $\boldsymbol{\mathsf{A}}^{(3)}$ are required. Nonetheless, as we show, these relations can be used to construct the full tensors of $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{A}}^{(3)}$ from stereo-PIV measurements in isotropic compressible flows.

2. Second- and third-order relations of compressible turbulence

2.1. Relations for homogeneous compressible flows

Following Betchov (1956), we consider the case when the flow is statistically homogeneous, i.e. $({\partial }/{\partial x_i}) \langle \rangle = 0$. For the second-order moment $\boldsymbol{\mathsf{A}}^{(2)}$, we then have

(2.1)\begin{align} {\mathsf{A}}^{(2)}_{ijji} & = \langle \overline{\boldsymbol{\mathsf{m}}^2} \rangle = \left \langle \frac{\partial u_i}{\partial x_j} \frac{\partial u_j}{\partial x_i} \right\rangle = \frac{\partial}{\partial x_j} \left \langle u_i \frac{\partial u_j}{\partial x_i} \right \rangle - \left\langle u_i \frac{\partial^2 u_j}{\partial x_i \partial x_j } \right\rangle \nonumber\\ & = -\frac{\partial }{\partial x_i} \left\langle u_i \frac{\partial u_j}{\partial x_j } \right\rangle + \left\langle \frac{\partial u_i}{\partial x_i} \frac{\partial u_j}{\partial x_j} \right\rangle = \langle \bar{\boldsymbol{\mathsf{m}}}^2 \rangle = {\mathsf{A}}^{(2)}_{iijj}, \end{align}

in which we used the notation $\bar {\boldsymbol{\mathsf{X}}} = \mathrm {tr}(\boldsymbol{\mathsf{X}})$ and the summation convention. For the third-order moments $\boldsymbol{\mathsf{A}}^{(3)}$, we use a general relation, derived in Appendix A (A7) among the gradients of three vector fields in homogeneous compressible turbulence. In the special case where the three fields are identical, (A7) reduces to the following relation for $\boldsymbol{\mathsf{A}}^{(3)}$:

(2.2) \begin{equation} {\mathsf{A}}^{(3)}_{ijjkki} = \langle \overline{\boldsymbol{\mathsf{m}}^3} \rangle = \tfrac{3}{2} \langle \bar{\boldsymbol{\mathsf{m}}} \overline{\boldsymbol{\mathsf{m}}^2}\rangle - \tfrac{1}{2} \langle \bar{\boldsymbol{\mathsf{m}}}^3 \rangle = \tfrac{3}{2} {\mathsf{A}}^{(3)}_{iijkkj} - \tfrac{1}{2} {\mathsf{A}}^{(3)}_{iijjkk}.\end{equation}

A straightforward consequence of (2.1) and (2.2) can be expressed by introducing $\boldsymbol{\mathsf{s}} \equiv (\boldsymbol{\mathsf{m}} + \boldsymbol{\mathsf{m}}^T)/2 - (\bar {\boldsymbol{\mathsf{m}}}/3) \boldsymbol{\mathsf{I}}$ and $\boldsymbol{\mathsf{w}} \equiv (\boldsymbol{\mathsf{m}} - \boldsymbol{\mathsf{m}}^T)/2$:

(2.3)$$\begin{gather} \langle \overline{\boldsymbol{\mathsf{s}}^2} \rangle ={-} \langle \overline{\boldsymbol{\mathsf{w}}^2} \rangle + \tfrac{2}{3} \langle \bar{\boldsymbol{\mathsf{m}}}^2 \rangle , \end{gather}$$

(2.4)$$\begin{gather}\langle \overline{\boldsymbol{\mathsf{s}}^3} \rangle ={-} 3 \langle \overline{\boldsymbol{\mathsf{wsw}}} \rangle + \tfrac{1}{2} \langle \bar{\boldsymbol{\mathsf{m}}} \overline{\boldsymbol{\mathsf{m}}^2} \rangle - \tfrac{5}{18} \langle \bar{\boldsymbol{\mathsf{m}}}^3 \rangle . \end{gather}$$

Equations formally similar to (2.1)–(2.4) were derived by Yang, Pumir & Xu (2020) for the perceived velocity-gradient tensor based on regular tetrahedra in incompressible homogeneous turbulence (see their (3.11), (3.12), and (3.37), (3.38)).

2.2. Isotropic flows

In the restricted case of isotropic flows, $\boldsymbol{\mathsf{A}}^{(2)}$ is expressible in terms of the Kronecker $\delta$-tensor ($\delta_{ij} = 1$ if $i=j$, and 0 otherwise), as (Pope 2000):

(2.5)\begin{equation} {\mathsf{A}}^{(2)}_{ijkl} = \alpha \delta_{ij} \delta_{kl} + \beta \delta_{ik} \delta_{jl} + \gamma \delta_{il} \delta_{jk},\end{equation}

where $\alpha$, $\beta$ and $\gamma$ are three scalar quantities. With this notation, $\langle \overline {\boldsymbol{\mathsf{m}}^2} \rangle = 9 \alpha + 3 \beta + 3 \gamma$ and $\langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle = 3 \alpha + 3 \beta + 9 \gamma$, so (2.1) implies that $\alpha = \gamma$, which means that only two quantities are necessary to fully determine the second-order tensor $\boldsymbol{\mathsf{A}}^{(2)}$. These can be determined in an experiment measuring ${\partial u_1}/{\partial x_1}$ and ${\partial u_1}/{\partial x_2}$ by, for example, planar PIV, or 2-component laser doppler velocimetry (2C LDV) or hot-wire anemometry with cross-wires using Taylor's frozen turbulence hypothesis, and noticing that ${\mathsf{A}}^{(2)}_{1111} = \alpha + \beta + \gamma$, ${\mathsf{A}}^{(2)}_{1122} = \alpha$ and ${\mathsf{A}}^{(2)}_{1212} = \beta$. Interestingly, $\alpha = \gamma$ implies that

(2.6)\begin{equation} \left \langle \frac{\partial u_1}{\partial x_1}\frac{\partial u_2}{\partial x_2}\right \rangle = {\mathsf{A}}^{(2)}_{1122} = {\mathsf{A}}^{(2)}_{1221} = \left \langle \frac{\partial u_1}{\partial x_2} \frac{\partial u_2}{\partial x_1}\right \rangle.\end{equation}

Additionally, we notice that $\langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle = 9 \alpha + 3 \beta + 3 \gamma = 3(4\alpha +\beta ) \geqslant 0$, thus $4\alpha + \beta \geqslant 0$, and $\langle \overline {\boldsymbol{\mathsf{w}}^2} \rangle = \frac {1}{2} ( \langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle - \langle \overline{\boldsymbol{\mathsf{m}}\boldsymbol{\mathsf{m}}^{\boldsymbol{T}}}\rangle) =3(\alpha - \beta ) \leqslant 0$, thus $\alpha \leqslant \beta$. These two inequalities constrain the ratio of components of $\boldsymbol{\mathsf{A}}^{(2)}$:

(2.7)\begin{equation} \frac{1}{3} = \frac{\beta}{2\beta + \beta} \leqslant \frac{{\mathsf{A}}^{(2)}_{1212}}{{\mathsf{A}}^{(2)}_{1111}} = \frac{\beta}{2\alpha+\beta} = 2 \times \frac{\beta}{(4\alpha + \beta)+\beta} \leqslant 2 .\end{equation}

The two limiting equalities arise when $(4\alpha + \beta ) = 0$, i.e. when the flow is incompressible, or when $\alpha = \beta$, i.e. when the flow is irrotational.

In isotropic flows, the third-order tensor $\boldsymbol{\mathsf{A}}^{(3)}$ can be expressed in terms of the Kronecker $\delta$ tensor, and of 5 scalars, $a_1\ldots a_5$ as (Pope 2000)

(2.8)\begin{align} {\mathsf{A}}^{(3)}_{ipjqkr} &= a_1 \delta_{ip} \delta_{jq} \delta_{kr} + a_2 ( \delta_{ip} \delta_{jk} \delta_{qr} + \delta_{jq} \delta_{ik} \delta_{pr} + \delta_{kr} \delta_{ij} \delta_{pq}) \nonumber\\ &\quad + a_3 ( \delta_{ip} \delta_{jr} \delta_{qk} + \delta_{jq} \delta_{ir} \delta_{pk} + \delta_{kr} \delta_{iq} \delta_{pj})+ a_4 (\delta_{iq} \delta_{pk} \delta_{jr} + \delta_{ir} \delta_{pj} \delta_{qk}) \nonumber\\ &\quad + a_5 (\delta_{ij} \delta_{pk} \delta_{qr} + \delta_{ij} \delta_{qk} \delta_{pr} + \delta_{ik} \delta_{pj} \delta_{qr} + \delta_{ik} \delta_{rj} \delta_{pq} + \delta_{jk} \delta_{qi} \delta_{pr} + \delta_{jk} \delta_{ri} \delta_{pq}), \end{align}

from which it is easy to obtain following expressions for the invariants of $\boldsymbol{\mathsf{m}}$:

(2.9)$$\begin{gather} \langle \bar{\boldsymbol{\mathsf{m}}}^3 \rangle = 27 a_1 + 27 a_2 + 27 a_3 + 6 a_4 + 18 a_5, \end{gather}$$

(2.10)$$\begin{gather}\langle \overline{\boldsymbol{\mathsf{m}}^3} \rangle = 3 a_1 + 9 a_2 + 27 a_3 + 30 a_4 + 36 a_5, \end{gather}$$

(2.11)$$\begin{gather}\langle \bar{\boldsymbol{\mathsf{m}}} \overline{\boldsymbol{\mathsf{m}}^2} \rangle = 9 a_1 + 15 a_2 + 33 a_3 + 18 a_4 + 30 a_5, \end{gather}$$

(2.12)$$\begin{gather}\langle \bar{\boldsymbol{\mathsf{m}}} \overline{\boldsymbol{\mathsf{m}}\boldsymbol{\mathsf{m}}^{\boldsymbol{T}}} \rangle = 9 a_1 + 33 a_2 + 15 a_3 + 6 a_4 + 42 a_5, \end{gather}$$

(2.13)$$\begin{gather}\langle \overline{\boldsymbol{\mathsf{m}}^2\boldsymbol{\mathsf{m}}^{\boldsymbol{T}}} \rangle = 3 a_1 + 21 a_2 + 15 a_3 + 12 a_4 + 54 a_5 . \end{gather}$$

In the incompressible case, with $\bar {\boldsymbol{\mathsf{m}}} = 0$ and $\langle \overline {\boldsymbol{\mathsf{m}}^3 } \rangle = 0$ (Betchov 1956), the left-hand sides of (2.9)–(2.12) are all zero, which provides four constraints to express $a_1, \ldots, a_5$ in terms of one scalar quantity. For compressible turbulence, only one relation is obtained by plugging equations (2.9)–(2.11) into (2.2), which leads to $a_1 = 3 a_3 - 2 a_4$. Thus four independent scalars are needed to completely determine $\boldsymbol{\mathsf{A}}^{(3)}$. Their experimental determination would require techniques such as stereo-PIV or 3-component Laser Doppler Velocimetry (3C LDV) with frozen turbulence hypothesis, giving access to spatial derivatives of the third velocity component normal to the plane of imaging, e.g. ${\mathsf{A}}^{(3)}_{113232}$. With these components, and using (2.8), one can determine four independent scalars, say, $a_2, \ldots, a_5$, as

(2.14)$$\begin{align} a_2 &= {\mathsf{A}}^{(3)}_{113232}, \end{align}$$

(2.15)$$\begin{align}a_3 &= \tfrac{1}{6} {\mathsf{A}}^{(3)}_{111111} - \tfrac{1}{2} {\mathsf{A}}^{(3)}_{212111}, \end{align}$$

(2.16)$$\begin{align}a_4 &= \tfrac{1}{3} {\mathsf{A}}^{(3)}_{111111} - {\mathsf{A}}^{(3)}_{212111} - \tfrac{1}{2} {\mathsf{A}}^{(3)}_{111122} + \tfrac{1}{2} {\mathsf{A}}^{(3)}_{113232}, \end{align}$$

(2.17)$$\begin{align}a_5 &= \tfrac{1}{2} {\mathsf{A}}^{(3)}_{212111} - \tfrac{1}{2} {\mathsf{A}}^{(3)}_{113232}. \end{align}$$

The other invariants of $\boldsymbol{\mathsf{m}}$, e.g. $\langle \overline {\boldsymbol{\mathsf{s}}^3} \rangle$, $\langle \overline {\boldsymbol{\mathsf{wsw}}} \rangle$, $\langle \bar {\boldsymbol{\mathsf{m}}} \overline {\boldsymbol{\mathsf{s}}^2} \rangle$ and $\langle \bar {\boldsymbol{\mathsf{m}}} \overline {\boldsymbol{\mathsf{w}}^2} \rangle$, can also be represented by these scalars. In particular, we note that

(2.18)\begin{equation} \langle \bar{\boldsymbol{\mathsf{m}}} \overline{\boldsymbol{\mathsf{w}}^2} \rangle ={-}9 a_2 + 9 a_3 + 6 a_4 -6 a_5 . \end{equation}

3. DNS of compressible turbulence

To test the relations presented above, we numerically solved the three-dimensional compressible Navier–Stokes equations in various turbulent flow configurations:

(3.1)\begin{equation} \left.\begin{array}{c} \dfrac{\partial \rho}{\partial t}+\dfrac{\partial \rho u_i}{\partial x_i}=0, \\ \dfrac{\partial \rho u_i}{\partial t}+\dfrac{\partial }{\partial x_j} ( \rho u_i u_j + P\delta_{ij})-\dfrac{\partial}{\partial x_j} \sigma_{ij}=0, \\ \dfrac{\partial E}{\partial t}+\dfrac{\partial }{\partial x_j} [(E + P) u_j ]-\dfrac{\partial}{\partial x_j} (\sigma_{ij} u_i-Q_j )=0, \end{array}\right\} \end{equation}

in which $\rho$ is the fluid density, $P$ is the pressure, $E=\frac {1}{2}\rho u_iu_i+P/(\gamma -1)$ is the total energy with $\gamma =1.4$ being the ratio of specific heats, $\sigma _{ij}=\mu ( m_{ij} + m_{ji} - \frac {2}{3}m_{kk} \delta _{ij} )$ is the viscous stress tensor with the effect of bulk viscosity neglected (Pan & Johnsen 2017), and $Q_j=-\kappa ({\partial T}/{\partial x_j})$ is the heat flux, the temperature $T$ is related to fluid pressure and density via the ideal gas law $P=\rho R T$ and the gas constant $R$, $\mu$ is the viscosity depending on the temperature via Sutherland's law $\mu = \mu _{ref} ({T}/{T_{ref}})^{{3}/{2}}(({T_{ref}+T_s})/({T+T_s}))$ with $\mu _{ref}$, $T_{ref}$ and $T_s$ being constants, and the thermal conductivity $\kappa$ is determined from the viscosity by a constant Prandtl number $Pr = \mu C_p / \kappa = 0.72$, with $C_p = {\gamma }/({\gamma - 1})R$ being the specific heat at constant pressure. The set of (3.1) are solved with a high-order finite difference method. Namely, the convection terms are computed by a seventh-order low-dissipative monotonicity-preserving scheme (Fang, Li & Lu 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 open-source code previously tested in DNS of various compressible turbulent flows with and without shock waves (Fang et al. 2013, 2014, 2015, 2020).

3.1. Compressible homogeneous and isotropic turbulence

We first discuss the results of decaying isotropic compressible turbulence. The computational domain is a $(2{\rm \pi} )^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 ) =A k^4\, \textrm {e}^{-2k^2/k_{0}^{2} }$, 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, Pullin & Kosovic (2001). The initial turbulent Mach number based on the root-mean-square velocity $u'=\sqrt {{u_i}{u_i}}$ is $Ma_t= u' / \langle {c}\rangle =2.0$, where $c= \sqrt {\gamma R T}$ is the speed of sound. The initial Reynolds number based on the Taylor microscale $\lambda =u' / \langle ( {\partial {u_1}}/{\partial x_1} ) ^2 + ( {\partial {u_2}}/{\partial x_2} ) ^2 + ( {\partial {u_3}}/{\partial x_3} ) ^2 \rangle ^{{1}/{2}}$ is $R_\lambda ={\langle {\rho }\rangle u'\lambda }/({\sqrt {3}\langle {\mu }\rangle})=450$. Here, we use the initial large-eddy-turnover time $\tau _0 = (\int _0^\infty E(k)/k \, \textrm {d}k)/{u'}^{3}$ to normalize time when quantifying the flow evolution. The computational domain is discretized with a $512^3$ uniform grid. Once the turbulent regime is established, for $t/\tau _0 \gtrsim 2.5$, the Kolmogorov scale $\eta$ continuously increases, from a value comparable to the grid size $\varDelta$, up to $\eta /\varDelta = 2.88$ at the end of the run, therefore ensuring that the flow is well resolved down to the dissipation scale. Similarly, the integral scale of the turbulence, $L = u'^3/(\varepsilon /\langle \rho \rangle )$ also varies rapidly in the initial stage and increases slowly with time, but remains less than $1/5$ of the simulation box size throughout the entire run. We stress that, because of the peculiar divergence-free initial velocity field, strong compression develops in the initial stage ($t/\tau _0 \lesssim 1$) of the DNS run before a turbulent regime sets in.

Figure 1(a) shows that the turbulent Mach number $Ma_t$, the turbulent kinetic energy $K=\frac {1}{2}\langle \rho {u_i}{u_i}\rangle$ and the Reynolds number $R_\lambda$, all decay monotonically with time. Specifically, $R_\lambda$ is less than $\sim 50$ for $t/\tau _0 \ge 1$. Note that even at time $t/\tau _0 > 10$, $Ma_t \lesssim 1$, the flow is still highly compressible with a large number of spatially distributed shocklets. The skewness of the longitudinal velocity derivative, $S = \langle ( {\partial {u_1}}/{\partial x_1} )^3 + ( {\partial {u_2}}/{\partial x_2} )^3 + ( {\partial {u_3}}/{\partial x_3} )^3 \rangle / ( \langle ( {\partial {u_1}}/{\partial x_1} )^2 \rangle ^{{3}/{2}} +\langle ( {\partial {u_2}}/{\partial x_2} )^2 \rangle ^{{3}/{2}} +\langle ( {\partial {u_3}}/{\partial x_3} )^2 \rangle ^{{3}/{2}} )$ shown in figure 1(b), however, develops a sharp negative peak at $t/\tau _0 \approx 0.35$ before returning to an approximately constant value, consistent with Samtaney et al. (2001). The sharp peak reflects the transient from the initial divergence-free velocity field. Although not obviously related to any experimental situation, this regime reflects some interesting properties of the dynamics. The larger value of $-S \approx 2$ in the later stage compared with $-S \approx 0.6$ in Samtaney et al. (2001) is most likely due to the larger Mach number in this work. As noticed by Donzis & John (2020), the skewness depends on the product of $Ma_t$, multiplied by the ratio $\delta ^2 \equiv \langle \varepsilon _d \rangle /\langle \varepsilon _s \rangle$. Quantitatively, we find that for $t /\tau _0 \ge 5$, the product $M_t \times \delta ^2 \approx 0.2$ and $-S$ fluctuates around $1.8$, consistent with the results of Donzis & John (2020).

Figure 1. (a) Evolution of the turbulent Mach number $Ma_t$, the turbulence kinetic energy normalized with its initial value $K(t)/K(0)$ and the Reynolds number $R_\lambda$ (shown in inset). (b) The evolution of the skewness of the longitudinal velocity derivative $S$. (c) The solenoidal and the dilatational parts of the energy dissipation rates $\varepsilon _s$ and $\varepsilon _d$ (solid lines), together with their approximations using averaged viscosity $\langle \mu \rangle \langle { \omega _{i}\omega _{i}} \rangle$ and $\frac {4}{3} \langle \mu \rangle \langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle$ (dashed lines). The inset shows the ratios of the approximations to the true values. (d) Comparison of the invariants of $\boldsymbol{\mathsf{A}}^{(2)}$ (full lines) involved in $\varepsilon _s$ and $\varepsilon _d$ ($\varepsilon _s > \varepsilon _d$) and the values determined from the approximate expressions given by (3.2a,b) (dashed lines).

In figure 1(c), we show the evolution of the turbulence dissipation rate, separated into the solenoidal (enstrophy) part $\varepsilon _s \equiv \langle {\mu \omega _{i}\omega _{i}} \rangle$ and the dilatational part $\varepsilon _d\equiv \frac {4}{3} \langle \mu \bar {\boldsymbol{\mathsf{m}}}^2 \rangle$. The solenoidal dissipation $\varepsilon _s$ grows first due to the build-up of turbulent structures from the initial random field, reaches a peak at $t/\tau _0 \approx 0.8$, then decreases gradually due to the decay of the turbulent fluctuation. The dilatational dissipation rate, $\varepsilon _d$, represents the contribution of compressibility to dissipation, which is exactly zero in incompressible turbulence. With our choice of solenoidal initial condition, $\varepsilon _d$ starts at zero and first grows rapidly, as shocklets are forming (Samtaney et al. 2001). The peak of $\varepsilon _d$ is reached at $t/\tau _0 \approx 0.65$, slightly earlier than $\varepsilon _s$. In our simulation, $\varepsilon _d$ contributes to approximately half the energy dissipation at the peak, but the solenoidal part remains the major cause of dissipation at later times. In figure 1(c), the dashed lines show the dissipations approximated with the averaged viscosity, $\varepsilon _s \approx \langle \mu \rangle \langle { \omega _{i}\omega _{i}} \rangle$ and $\varepsilon _d \approx \frac {4}{3} \langle \mu \rangle \langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle$, which have been shown to give values with good accuracy for steady compressible channel flows (Huang, Coleman & Bradshaw 1995). In our DNS, judging from the ratios shown in the inset, this approximation also works well except at the initial stage when the flow is adjusting from the initial conditions, before a turbulent regime is reached. Using this approximation, the energy dissipations $\varepsilon _s$ and $\varepsilon _d$ can be determined from two components of velocity derivatives by taking advantage of the isotropic expression (2.5) for the full tensor $\boldsymbol{\mathsf{A}}^{(2)}$. In particular, note that we have shown that ${\mathsf{A}}^{(2)}_{1111} = 2 \alpha + \beta$ and ${\mathsf{A}}^{(2)}_{1212} = \beta$, and that $\langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle = 12 \alpha + 3 \beta$ and $\langle \overline {\boldsymbol{\mathsf{w}}^2} \rangle = 3(\alpha - \beta )$, from which $\varepsilon _s$ and $\varepsilon _d$ can be expressed as

(3.2a,b)\begin{equation} \varepsilon_s \approx 3 \langle \mu\rangle ( 3 {\mathsf{A}}^{(2)}_{1212} - {\mathsf{A}}^{(2)}_{1111}) \quad \textrm{and} \quad \varepsilon_d \approx 4 \langle \mu\rangle ( 2{\mathsf{A}}^{(2)}_{1111} - {\mathsf{A}}^{(2)}_{1212}) .\end{equation}

Figure 1(d), demonstrates that the values obtained with (3.2a,b), shown by the dashed lines, compare very well with the numerical values of $\varepsilon _s$ and $\varepsilon _d$ (note that $\varepsilon _s > \varepsilon _d$). We interpret the small deviation as a result of a small residual anisotropy in our DNS.

We now discuss the structure of the velocity-gradient correlations in this flow. Figure 2 shows the evolution of the left-hand side and right-hand side of (2.1) and (2.2), the identities for the invariants of $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{A}}^{(3)}$, respectively. The ratios between the two sides of those equations, are exactly 1, as predicted for homogeneous flows.

Figure 2. Evolution of invariants of $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{A}}^{(3)}$ in compressible homogeneous isotropic decaying turbulence. (a) Second-order invariants, $\langle \overline {\boldsymbol{\mathsf{m}}^2} \rangle /\langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle$ and $\langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle / \langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle$. (b) Third-order invariants $\langle \overline {\boldsymbol{\mathsf{m}}^3} \rangle /\langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle ^{3/2}$ and $(\frac {3}{2}\langle \overline {\boldsymbol{\mathsf{m}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle - \frac {1}{2} \langle \bar {\boldsymbol{\mathsf{m}}}^3 \rangle ) / \langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle ^{3/2}$.

Our analysis predicts, for homogeneous and isotropic turbulence, further relations among the components of $\boldsymbol{\mathsf{A}}^{(2)}$ and $\boldsymbol{\mathsf{A}}^{(3)}$. Figure 3(a) shows that ${\mathsf{A}}^{(2)}_{1122}$ and ${\mathsf{A}}^{(2)}_{1221}$ are equal during the entire simulation, as predicted by (2.6). In figure 3(b), the ratio of the components ${\mathsf{A}}^{(2)}_{1212} / {\mathsf{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 dilatational dissipation $\varepsilon _d$. The minimal possible value of $1/3$ for ${\mathsf{A}}^{(2)}_{1212} / {\mathsf{A}}^{(2)}_{1111}$ corresponds to a compressible irrotational flow. The observed value of ${\mathsf{A}}^{(2)}_{1212} / {\mathsf{A}}^{(2)}_{1111}$, close to $1.5$ at later times, indicates instead the prevalence of the solenoidal part.

Figure 3. (a) Time evolution of ${\mathsf{A}}^{(2)}_{1122}$ and ${\mathsf{A}}^{(2)}_{1221}$, which should be equal in isotropic compressible turbulence. (b) The ratio of ${\mathsf{A}}^{(2)}_{1212} / {\mathsf{A}}^{(2)}_{1111}$, which should lie between $1/3$ and $2$ in isotropic turbulence.

Figure 4(a) shows the evolution of the invariants of $\boldsymbol{\mathsf{A}}^{(3)}$ after the early stage when the flow rapidly adjusts in response to the initial divergence-free condition. Remarkably, the magnitude of $\langle \overline {\boldsymbol{\mathsf{w}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle$ is very small compared with the others, consistent with the observation that in compressible homogeneous turbulence, vorticity and dilatation are nearly uncorrelated (Erlebacher & Sarkar 1993; Wang et al. 2012). This, in view of (2.18), provides an additional relation: $3 a_2 \approx 3 a_3 + 2 a_4 - 2 a_5$, which then leaves only three independent quantities for a complete determination of $\boldsymbol{\mathsf{A}}^{(3)}$. Plugging (2.14)–(2.17) into this relation leads to

(3.3)\begin{equation} \tfrac{7}{6} {\mathsf{A}}^{(3)}_{111111} \approx \tfrac{9}{2} {\mathsf{A}}^{(3)}_{121211} + {\mathsf{A}}^{(3)}_{111122} + {\mathsf{A}}^{(3)}_{113232}. \end{equation}

Figure 4(b) shows the left-hand side and the right-hand side of (3.3) and their ratio, demonstrating that (3.3), obtained from the simplification of $\langle \overline {\boldsymbol{\mathsf{w}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle \approx 0$, approximately holds even in the early stage, with a deviation of less than $20\,\%$. This result may be helpful to reconstruct the complete isotropic expression of ${\mathsf{A}}^{(3)}_{ipjqkr}$ from planar PIV or 2C LDV data, instead of requiring stereo-PIV or 3C LDV.

Figure 4. (a) Evolution of the invariants $\langle \bar {\boldsymbol{\mathsf{m}}}^3 \rangle$, $\langle \overline {\boldsymbol{\mathsf{w}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle$, $\langle \overline {\boldsymbol{\mathsf{s}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle$, $\langle \overline {\boldsymbol{\mathsf{wsw}}} \rangle$ and $\langle \overline {\boldsymbol{\mathsf{s}}^3} \rangle$, all normalized by $\langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle ^{3/2}$. (b) The values of the left-hand side and right-hand side of (3.3), normalized by $\langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle ^{3/2}$, and their ratio (shown in inset).

3.2. Compressible mixing layer

Next we perform a DNS of a planar compressible mixing layer, formed by two co-moving free streams with Mach numbers $Ma_1=U_1/c_1=7.5$ and $Ma_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 set-up 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) =\frac {1}{2}[U_1+U_2 +(U_1-U_2)\tanh (2y/\delta _{\omega 0} )],$ with the inflow vorticity thickness ${\delta _{\omega 0}} = 1$. The speeds of sound in the incoming upper and lower streams are equal, $c_1 = c_2$. The convective Mach number is $Ma_c = (U_1-U_2)/(c_1+c_2)=3$ and the Reynolds number is $Re_c =\rho _1 (U_1-U_2){\delta _{\omega 0}}/\mu _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{\delta _{\omega 0}}\times 100{\delta _{\omega 0}}\times 16{\delta _{\omega 0}}$ in the $x$, $y$ and $z$ directions, respectively. The momentum thickness of the mixing layer grows nearly linearly with downstream distance and reaches approximately $2 \delta _{\omega0}$ at $x /\delta _{\omega0} = 450$. Therefore, the simulation domain size of $100 \delta _{\omega0}$ in the $y$ direction is large enough to cover the entire mixing layer. The domain is discretized using a mesh of $3050\times 400\times 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, $\theta$, grows linearly downstream and the mean velocity profiles are self-similar, i.e. $U(x,y)$ at different $x$ locations collapse to $U(\xi )$ with $\xi =(y-y_c)/\theta$, 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 left-hand side and the right-hand side of (2.1) and (2.2) in the centre region of the mixing layer $(-6 \leqslant y/\theta \leqslant 6)$ at the downstream locations $x/\delta _{\omega 0} = 350$ (dashed lines) and $x / \delta _{\omega 0}=400$ (solid lines). Although the flow is not homogeneous, the ratios $\langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle / \langle \overline {\boldsymbol{\mathsf{m}}^2} \rangle$ and $\langle \overline {\boldsymbol{\mathsf{m}}^3} \rangle /(\frac {3}{2} \langle \overline {\boldsymbol{\mathsf{m}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle - \frac {1}{2} \langle \bar {\boldsymbol{\mathsf{m}}}^3 \rangle )$, shown in the insets, are very close to unity, so (2.1) and (2.2) are still approximately valid. Figure 6 shows the profiles of various third-order invariants $\langle \bar {\boldsymbol{\mathsf{m}}}^3 \rangle$, $\langle \overline {\boldsymbol{\mathsf{w}}^2}\bar {\boldsymbol{\mathsf{m}}} \rangle$, $\langle \overline {\boldsymbol{\mathsf{s}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle$, $\langle \overline {\boldsymbol{\mathsf{s}}^3} \rangle$ and $\langle \overline {\boldsymbol{\mathsf{wsw}}} \rangle$. Despite the inhomogeneity and anisotropy, the vorticity-dilatation correlation $\langle \overline {\boldsymbol{\mathsf{w}}^2}\bar {\boldsymbol{\mathsf{m}}} \rangle$ remains very small compared with other invariants, which could help to obtain more relations among invariants.

Figure 5. Approximate validity of the homogeneous relations in the compressible mixing layer. (a) Normalized second-order invariants, $\langle \overline {\boldsymbol{\mathsf{m}}^2} \rangle / \langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle$, $\langle \bar {\boldsymbol{\mathsf{m}}}^2 \rangle / \langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle$, and their ratio (shown in the inset). (b) Normalized third-order invariants, $\langle \overline {\boldsymbol{\mathsf{m}}^3} \rangle / \langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle ^{3/2}$, $(\frac {3}{2}\langle \overline {\boldsymbol{\mathsf{m}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle - \frac {1}{2} \langle \bar {\boldsymbol{\mathsf{m}}}^3 \rangle ) / \langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle ^{3/2}$, and their ratio (shown in the inset). In both plots, dashed lines correspond to $x/\delta _{\omega 0} = 350$ and solid lines to $x/\delta _{\omega 0} = 400$.

Figure 6. (a) Profiles of the invariants $\langle \bar {\boldsymbol{\mathsf{m}}}^3 \rangle$, $\langle \overline {\boldsymbol{\mathsf{w}}^2}\bar {\boldsymbol{\mathsf{m}}} \rangle$, $\langle \overline {\boldsymbol{\mathsf{s}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle$, $\frac {1}{3}\langle \overline {\boldsymbol{\mathsf{s}}^3}\rangle$ and $\langle \overline {\boldsymbol{\mathsf{wsw}}} \rangle$ in the mixing layer, all normalized by $\langle \overline {\boldsymbol{\mathsf{mm}}^{\boldsymbol{T}}} \rangle ^{3/2}$. (b) Relative ratios $\langle \overline {\boldsymbol{\mathsf{m}}^3} \rangle$/$\langle \overline {\boldsymbol{\mathsf{s}}^3} \rangle$, $\langle \overline {\boldsymbol{\mathsf{w}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle$/$\langle \overline {\boldsymbol{\mathsf{s}}^3} \rangle$, $\langle \overline {\boldsymbol{\mathsf{s}}^2} \bar {\boldsymbol{\mathsf{m}}} \rangle$/$\langle \overline {\boldsymbol{\mathsf{s}}^3} \rangle$ and $\langle \overline {\boldsymbol{\mathsf{wsw}}} \rangle$/$\langle \overline {\boldsymbol{\mathsf{s}}^3} \rangle$. In both plots, dashed lines correspond to $x/\delta _{\omega 0} = 350$ and solid lines to $x/\delta _{\omega 0} = 400$.

4. 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, for example, to determine separately the solenoidal and the dilatational energy dissipation rates from only two velocity derivatives ${\partial u_1}/{\partial x_1}$ and ${\partial u_2}/{\partial x_1}$ in isotropic compressible turbulence. In the future, it would be interesting to investigate the structure of the velocity gradient implied by these relations, as done by Betchov (1956) for incompressible turbulence.