1. Introduction
Universal scaling remains a central pursuit in turbulence research, with homogeneous isotropic turbulence serving as a canonical framework that retains the essential dynamics while avoiding complications from inhomogeneity and boundaries. In incompressible flows, the Taylor-microscale Reynolds number
$\textit{Re}_{\!\lambda }$
characterises the turbulence. In compressible turbulence, however, interactions among vortical, acoustic and entropic modes (Kovásznay Reference Kovásznay1953; Chu & Kovásznay Reference Chu and Kovásznay1958) render
$\textit{Re}_{\!\lambda }$
alone insufficient. For fixed thermophysical parameters, including the specific heat ratio
$\gamma$
and Prandtl number
$\textit{Pr}$
, conventional theory typically predicts normalised statistics as functions of
$\textit{Re}_{\!\lambda }$
and the turbulent Mach number
$\textit{Ma}_t$
. In compressible homogeneous isotropic turbulence, the two governing parameters are defined as
\begin{equation} \textit{Re}_{\!\lambda } = \frac { u_{\textit{rms}} \lambda \langle \rho \rangle }{\langle \mu \rangle }, \,\,\,\textit{Ma}_{t} = \frac {\sqrt {3} u_{\textit{rms}}}{\langle \sqrt {\gamma \textit{RT}} \rangle }, \end{equation}
where
$\rho$
is the density,
$\mu$
is the shear viscosity,
$R$
is the specific gas constant,
$T$
is the temperature,
$u_{\textit{rms}} = \sqrt { \langle u^2 \rangle / 3 }$
is the root-mean-square (r.m.s.) velocity,
$\lambda$
is the Taylor microscale
$\lambda$
and angular brackets are suitably defined ensemble averages. The Helmholtz decomposition (Helmholtz Reference Helmholtz1867) uniquely separates the velocity field into two orthogonal components, namely
$\boldsymbol{u}=\boldsymbol{u}_s+\boldsymbol{u}_d$
; the divergence-free solenoidal velocity
$\boldsymbol{u}_s$
(i.e.
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_s=0$
) is associated with the vortical mode, and the irrotational dilatational velocity
$\boldsymbol{u}_d$
(i.e.
$\boldsymbol {\nabla }\boldsymbol{\times }\boldsymbol{u}_d=\boldsymbol{0}$
) is related to the change rate of the volume element. Donzis & John (Reference Donzis and John2020) demonstrated that self-similar scaling in compressible turbulence requires the consideration of dilatational motion. And they introduced the ratio of the r.m.s. dilatational velocity
$u_{\textit{d,rms}}=\sqrt { \langle u_d^2 \rangle / 3 }$
to the r.m.s. solenoidal velocity
$u_{\textit{s,rms}}=\sqrt { \langle u_s^2 \rangle / 3 }$
. This velocity-ratio framework provides a more universal scaling and has since been extended to mean shear, thermal, reactive (Donzis & John Reference Donzis and John2020; John & Donzis Reference John and Donzis2024), magnetohydrodynamic (Li et al. Reference Li, Yang, Matthaeus, Jiang, Wan and Chen2024) and bulk-viscosity flows.
Pressure may also be decomposed into solenoidal and dilatational components, i.e.
$p = p_s + p_d$
(Donzis & John Reference Donzis and John2020). Unlike the turbulent kinetic energy (TKE) per unit mass
$\langle u^2\rangle /2 = \langle u_s^2\rangle /2 + \langle u_d^2\rangle /2$
, whose solenoidal–dilatational split is exact, the square of the r.m.s. pressure
$p_{\textit{rms}}^2 = p_{\textit{s,rms}}^2 + p_{\textit{d,rms}}^2 + 2 r_p p_{\textit{s,rms}} p_{\textit{d,rms}}$
has an additional contribution from the correlation coefficient
$r_p$
between
$p_s$
and
$p_d$
. Dilatational pressure
$p_d$
reflects an acoustic dynamics with a characteristic time scale
$t\sim 1/(c_s| \boldsymbol{\nabla }|)$
corresponding to the dilatational length scale
$1/|\boldsymbol{\nabla }|$
(
$\boldsymbol {\nabla}$
is the nabla operator of the spatial derivative), and the speed of sound
$c_s=\sqrt {\gamma \langle p\rangle /\langle \rho \rangle }$
. From mass conservation, the r.m.s. dilatational pressure satisfies
$p_{\textit{d,rms}}/\langle p\rangle \sim \gamma \textit{Ma}_d$
(Donzis & John Reference Donzis and John2020), where
$\textit{Ma}_d = u_{\textit{d,rms}} \textit{Ma}_t / u_{\textit{rms}}$
is the dilatational turbulent Mach number. The pressure Poisson equation further yields
$p_{\textit{s,rms}} \sim \langle \rho \rangle u_{\textit{s,rms}}^2$
in the incompressible limit. Within the pressure decomposition framework, scaling with respect to
$\textit{Ma}_d$
becomes clear: when
$p_{s}$
dominates,
$p_{\textit{rms}}/\langle p\rangle \sim \gamma \textit{Ma}_d^2 u_{\textit{s,rms}}^2/u_{\textit{d,rms}}^2$
, whereas dilatational dominance gives
$p_{\textit{rms}}/\langle p\rangle \sim \gamma \textit{Ma}_d$
, the so-called dilatationally dominated
$p$
-equipartition (DDE) regime (Donzis & John Reference Donzis and John2020).
Modelling the dissipation in compressible turbulence requires understanding of the pressure-dilatation term
$\langle p\vartheta \rangle$
and dilatational dissipation
$\epsilon _d = (\chi +4/3)\langle \mu \vartheta ^2 \rangle /\langle \rho \rangle$
, where
$\vartheta =\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}$
is the dilatation, and
$\chi =\mu _V/\mu$
is the bulk-to-shear viscosity ratio. Conventional pseudo-sound theory (Ristorcelli Reference Ristorcelli1995, Reference Ristorcelli1997; Ristorcelli & Blaisdell Reference Ristorcelli and Blaisdell1997) forecasts both terms to scale with
$\textit{Ma}_t^2$
. Within this framework,
$\langle p\vartheta \rangle$
, which characterises the interconversion between TKE and internal energy, is assumed to scale with the product of
$\textit{Ma}_t^2$
and the solenoidal dissipation
$\epsilon _s = \langle \mu \omega ^2 \rangle /\langle \rho \rangle$
in isotropic turbulence (Sarkar Reference Sarkar1992; Ristorcelli Reference Ristorcelli1995, Reference Ristorcelli1997; Ristorcelli & Blaisdell Reference Ristorcelli and Blaisdell1997), where
$\boldsymbol{\omega } = \boldsymbol{\nabla }\boldsymbol{\times }\boldsymbol{u}$
denotes the vorticity. Similarly,
$\epsilon _d$
satisfies
$\epsilon _d/\epsilon _s\sim \textit{Ma}_t^2+O(\textit{Ma}_t^4)$
(Sarkar et al. Reference Sarkar, Erlebacher, Hussaini and Kreiss1991; Ristorcelli Reference Ristorcelli1997; Ristorcelli & Blaisdell Reference Ristorcelli and Blaisdell1997). A more recent approach (Donzis & John Reference Donzis and John2020) relates
$\epsilon _d/\epsilon _s$
to the velocity ratio
$u_{\textit{d,rms}}/u_{\textit{s,rms}}$
, with direct numerical simulations (DNSs) (Donzis & John Reference Donzis and John2020; John & Donzis Reference John and Donzis2024) supporting the scaling
$\epsilon _d/\epsilon _s\sim u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
. This result extends the
$\textit{Ma}_t$
-based conventional theory and offers new insight into compressible turbulence modelling. However, the dependence of
$\langle p\vartheta \rangle$
on the velocity ratio remains unexplored.
High-order moments of the velocity field exhibit unique phenomena associated with
$u_{\textit{d,rms}}/u_{\textit{s,rms}}$
for compressible turbulence. The skewness
$S_u$
and flatness
$F_u$
, i.e.
\begin{equation} S_u = {\left \langle \left [(\partial u_x/\partial x)^3 + (\partial u_y/\partial y)^3 + (\partial u_z/\partial z)^3\right ]/3\right \rangle \over \left \langle \left [(\partial u_x/\partial x)^2 + (\partial u_y/\partial y)^2 + (\partial u_z/\partial z)^2\right ]/3\right \rangle ^{3/2}}, \end{equation}
\begin{equation} F_u = {\left \langle \left [(\partial u_x/\partial x)^4 + (\partial u_y/\partial y)^4 + (\partial u_z/\partial z)^4\right ]/3\right \rangle \over \left \langle \left [(\partial u_x/\partial x)^2 + (\partial u_y/\partial y)^2 + (\partial u_z/\partial z)^2\right ]/3\right \rangle ^2}, \end{equation}
which are closely related to vortex stretching and intermittency, deviate markedly from the canonical values observed in fully developed incompressible turbulence once
$\textit{Ma}_t u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
exceeds a critical threshold, as a consequence of the emergence of intense shocklets (Samtaney, Pullin & Kosović Reference Samtaney, Pullin and Kosović2001; Wang et al. Reference Wang, Shi, Wang, Xiao, He and Chen2011). When the so-called,
$S_u$
divergence, occurs, John & Donzis (Reference John and Donzis2024) conjectured that it signifies the onset of a regime where dilatational motions begin to influence the turbulent energy cascade.
The dynamics of dilatational kinetic energy depends on several mechanisms, with bulk viscosity
$\mu _V$
exerting a particularly strong influence. The bulk viscosity of a gas is strongly correlated with its internal degrees of freedom; for example,
$\chi$
at
$300\rm\,K$
is
$0.73$
for
$N_2$
,
$7.4$
for water vapour,
$30$
for
$H_2$
and as high as
$3000$
for
$\rm CO_2$
(Cramer Reference Cramer2012). As a thermophysical property,
$\mu _V$
enhances the dissipation of dilatational motions and thereby mitigates flow compressibility (Pan & Johnsen Reference Pan and Johnsen2017). Increasing
$\mu _V$
leads the velocity ratio to scale as
$\textit{Ma}_t^2$
up to turbulent Mach numbers of
$0.6$
(Chen et al. Reference Chen, Wang, Wang, Wan, Li and Chen2019). Sakurai & Ishihara (Reference Sakurai and Ishihara2023, Reference Sakurai and Ishihara2024) showed that the dilatational energy spectrum exhibits a
$k^{-3}$
power-law scaling at high wavenumbers, which is attributed to the presence of shocklets. As
$\mu _V$
increases, both the amplitude of the high-wavenumber dilatational spectrum and the associated
$k^{-3}$
scaling are progressively suppressed. Despite its relevance, the robustness of established turbulence-scaling behaviours under varying bulk viscosity has not been systematically assessed. To address this issue, we conduct numerical simulations to investigate how bulk viscosity affects the universality of key scaling relations. We examine the scaling of
$p_{\textit{rms}}$
across
$\mu _V$
, propose
$\chi $
-dependent modifications for
$\langle p\vartheta \rangle$
and
$\epsilon _d$
and show that the threshold associated with the onset of the ‘
$S_u$
divergence’ increases with
$\chi$
. These results highlight the regulatory role of bulk viscosity in modulating the evolution of compressible turbulence.
2. Numerical details
Numerical simulations of the Navier–Stokes–Fourier equations are performed using the hybrid method developed by Wang et al. (Reference Wang, Wang, Xiao, Shi and Chen2010), which combines a seventh-order weighted essentially non-oscillatory (WENO) scheme with an eighth-order compact finite-difference scheme. The governing equations follow the macroscopic formulation of Qi et al. (Reference Qi, Wang, Guo and Chen2023), as detailed in the supplemental material is available at https://doi.org/10.1017/jfm.2026.11466, and the numerical method solves the conservative form of the total energy equation. Here, we emphasise that the results presented here do not constitute a DNS in the strict sense, since strong shocks with a thickness smaller than the grid spacing are not directly resolved but are effectively modelled through the inherent numerical subgrid dissipation of the WENO scheme. The viscous stress in the momentum equation is
\begin{equation} \boldsymbol{\sigma } = \mu \big [ \boldsymbol{\nabla }\boldsymbol{u}+( \boldsymbol{\nabla }\boldsymbol{u})^T + ( \chi - 2/3 )\vartheta \boldsymbol{I} \big ], \end{equation}
where
$\mu$
satisfies Sutherland,s law (Sutherland Reference Sutherland1893), a superscript
$T$
denotes the transpose for the velocity gradient and
$\boldsymbol{I}$
is the identity tensor. To mitigate the impact of dilatational forcing (Petersen & Livescu Reference Petersen and Livescu2010), solenoidal large-scale forcing is applied by maintaining prescribed kinetic energies in the two lowest-wavenumber shells, consistent with the
$k^{-5/3}$
spectrum. The specific heat ratio and Prandtl number are set to
$\gamma = 1.4$
and
$\textit{Pr} =0.7$
, respectively.
Table 1 in the supplementary material summarises the key parameters and statistical quantities for the simulations, which cover the parameter ranges of
$\textit{Re}_{\!\lambda } \in (99,181)$
,
$\textit{Ma}_t \in (0.39, 1.04)$
and
$\chi \in [0,30]$
. The Taylor microscale
$\lambda$
and Kolmogorov scale
$\eta$
, along with the two parts of the Favre-averaged TKE,
$K_s$
and
$K_d$
, the r.m.s. vorticity
$\omega _{\textit{rms}}$
and dilation
$\vartheta _{\textit{rms}}$
, the integral time scale
$L_I$
, eddy-turnover time
$T_e$
and kinetic energy spectrum
$E_K$
are defined as
\begin{eqnarray} \lambda &=& \frac { \sqrt {3} u_{\textit{rms}} }{ \sqrt { \left \langle({\partial u_x/\partial x} )^2 + ({\partial u_y/\partial y})^2 + ({\partial u_z/\partial z})^2 \right \rangle } } , \,\eta = \left [ \frac { \langle \mu /\rho \rangle ^3 }{ \varepsilon } \right ]^{1/4} , \nonumber \\[6pt]K_s &=& {\big\langle \rho u_s^2 \big\rangle \over 2\langle \rho \rangle } , \,K_d = {\big\langle \rho u_d^2 \big\rangle \over 2\langle \rho \rangle } , \,\omega _{\textit{rms}} = \sqrt { \langle \boldsymbol{\omega } \boldsymbol{\cdot }\boldsymbol{\omega } \rangle } , \,\vartheta _{\textit{rms}} = \sqrt { \langle \vartheta ^2 \rangle } , \nonumber \\[6pt]L_I &=& { \pi \over 2u_{\textit{rms}}^2 }\int _{0}^{\infty } {E_K(k) \over k} {\rm d}k , \,T_e = { L_I \over u_{\textit{rms}} } , \,E_K(k,t) = \sum \limits _{k-{1\over 2}\lt |\boldsymbol{k}|\leqslant k+{1\over 2}} {|\hat {\boldsymbol{u}}(\boldsymbol{k},t)|^2 \over 2}, \end{eqnarray}
with
$\hat {\boldsymbol{u}}$
denoting the Fourier-transformed velocity in the wavenumber space
$\boldsymbol{k}$
. Stationary ensemble statistics are obtained from spatial averages with additional temporal averaging over more than
$28 T_e$
. All simulations satisfy
$k_{max}\eta \gt 2.27$
, ensuring convergence of high-order velocity-gradient moments (Wang, Gotoh & Watanabe Reference Wang, Gotoh and Watanabe2017).
3. Results
3.1. Favre-averaged velocity decomposition
(a) Scaling of
$K_d/K_s$
vs.
$u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
and (b)
$\textit{Ma}_t-u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
diagram. Data are represented as the `mean
$\pm$
standard error, in all figures.
In compressible turbulence modelling, the TKE is typically defined using the Favre (density-weighted) average. We therefore begin by analysing the solenoidal and dilatational Favre-averaged components,
$K_s$
and
$K_d$
. Although this decomposition is not orthogonal
$\langle \rho \boldsymbol{u}_s\boldsymbol{\cdot }\boldsymbol{u}_c\rangle \neq 0$
, figure 1(a) shows that
$K_d/K_s$
collapses onto
$u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
for all cases, demonstrating that the Favre-averaged TKE ratio provides a reliable proxy for the velocity ratio under the present density-variation conditions. Besides, the inset in figure 1(a) shows that approximating
$K$
by
$K_s+K_d$
yields deviations within
$5\,\%$
, with the discrepancy increasing at higher velocity ratios.
The
$\textit{Ma}_t-u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
diagram in figure 1(b) includes two established criteria. The `DDE’ line (
$\mathcal{D} = \textit{Ma}_t {u_{\textit{s,rms}}^2 / (u_{\textit{rms}} u_{\textit{d,rms}})} \gt 0.5$
) (Donzis & John Reference Donzis and John2020) indicates that all cases lie in the solenoidal-pressure-dominated regime, with solenoidal contributions strengthened at larger bulk viscosity. The `
$S_u-div$
’ line (
$\textit{Ma}_t u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2 \gt 3\times 10^{-2}$
) (John & Donzis Reference John and Donzis2024) suggests that several simulations satisfy the empirical onset condition. Regarding the scaling of the turbulent Mach number, Donzis & John (Reference Donzis and John2020) proposed that solenoidally forced isotropic turbulence may follow
$\textit{Ma}_t \sim C (u_{\textit{d,rms}}/u_{\textit{s,rms}})^\alpha$
, with the same exponent
$\alpha$
, while most of the cases they examined assume a zero bulk viscosity. By contrast, our results demonstrate a pronounced dependence of
$\alpha$
on the bulk viscosity. The exponent decreases systematically with increasing
$\chi$
and is well described by the empirical relation
$\alpha \approx (0.460\chi +4.57)/(\chi +5.99)$
.
3.2. Pressure
Scaling of (a)
$p_{\textit{rms}}^2/\langle p\rangle ^2$
with
$\textit{Ma}_t$
and (b)
$p_{\textit{rms}}^2/( \langle p\rangle \gamma \textit{Ma}_d )^2$
with
$\mathcal{D}$
.
Figure 2(a) shows the variation of the normalised r.m.s. pressure
$p_{\textit{rms}}^2/\langle p\rangle ^2$
with
$\textit{Ma}_t$
. Classical pseudo-sound theory (Ristorcelli Reference Ristorcelli1997; Donzis & Jagannathan Reference Donzis and Jagannathan2013; Wang et al. Reference Wang, Gotoh and Watanabe2017) predicts a quartic scaling
$p_{\textit{rms}}^2/\langle p\rangle ^2\sim \textit{Ma}_t^4$
, yet the measured exponents are consistently smaller than four. As
$\chi$
increases, however, the scaling gradually approaches the quartic law. This trend can be attributed to the underlying pressure-generation mechanisms: pseudo-sound theory derives from the pressure Poisson equation at the incompressible limit, which reflects primarily solenoidal contributions. Elevated bulk viscosity suppresses dilatational motions, thereby enhancing the relative weight of solenoidal pressure and driving the overall behaviour closer to the pseudo-sound prediction. A unified representation of r.m.s. pressure was proposed by Donzis & John (Reference Donzis and John2020), through a decomposition into solenoidal and dilatational components
\begin{equation} {p_{\textit{rms}}^2\over ( \langle p\rangle \gamma \textit{Ma}_d )^2} \approx C_1\, \textit{Ma}_d^2 {u_{\textit{s,rms}}^4 \over u_{\textit{d,rms}}^4} + C_2 = C_1 \mathcal{D}^2 + C_2 , \end{equation}
where
$\mathcal{D} = \textit{Ma}_d {u_{\textit{s,rms}}^2 /u_{\textit{d,rms}}^2} = \textit{Ma}_t {u_{\textit{s,rms}}^2 / (u_{\textit{rms}} u_{\textit{d,rms}})}$
is the reciprocal of that defined in Donzis & John (Reference Donzis and John2020). and the first term captures the solenoidal contribution while the second represents dilatational pressure. The detailed derivations of (3.1)–(3.3) are provided in the supplementary material. Figure 2(b) indicates that the normalised r.m.s. pressure follows the
$\mathcal{D}^2$
scaling reasonably well with
$C_1\approx 0.08$
and
$C_2\approx 0$
, and the agreement improves markedly with increasing
$\chi$
, further confirming that the pressure field becomes progressively solenoidal-dominated when bulk-viscosity effects are strong.
3.3. Dissipation
Scaling of (a)
$\langle p\vartheta \rangle /(\langle \rho \rangle \epsilon _s \textit{Ma}_t^2)$
with
$\textit{Ma}_t$
and (b)
$\langle p\vartheta \rangle \textit{Ma}_d ({ u_{\textit{s,rms}} / u_{\textit{d,rms}} })^{3}/(\langle \rho \rangle \epsilon _s)$
with
$\mathcal{D}$
.
We first examine the pressure-dilatation term
$\langle p\vartheta \rangle$
, representing the net work of pressure on a fluid element. In homogeneous isotropic turbulence, positive
$\vartheta$
converts internal energy to TKE, while negative
$\vartheta$
transfers TKE to internal energy. Figure 3(a) shows that the classical pseudo-sound scaling
$\langle p\vartheta \rangle /\langle \rho \rangle \epsilon _s \sim \textit{Ma}_t^2$
(Sarkar Reference Sarkar1992) is not supported; the data display a stronger Mach-number dependence. Notably, with zero bulk viscosity (
$\chi =0$
), the transfer is strictly from TKE to internal energy, whereas for
$\chi \geqslant 5$
internal-to-TKE transfer appears at higher Mach numbers.
Based on the velocity ratio, an alternative model is constructed for
$\langle p\vartheta \rangle$
. Following the methodology of Zeman (Reference Zeman1991), we rewrite the term as
$\langle p\vartheta \rangle \approx - \left \langle (p_s + p_d) {\partial p_d / \partial t} \right \rangle / \gamma \langle p\rangle$
and simplify using the acoustic time scale
$t\sim 1/(c_s|\boldsymbol{\nabla }|)$
, the solenoidal dissipation
$\epsilon _s \approx u_{\textit{s,rms}}^3/\mathcal{L}$
in the equilibrium turbulence (Kolmogorov Reference Kolmogorov1941a
,
Reference Kolmogorovb
; Rubinstein & Clark Reference Rubinstein and Clark2017) and the aforementioned model for
$p_{\textit{rms}}$
(Donzis & John Reference Donzis and John2020), leading to the relations
\begin{eqnarray} \langle p\vartheta \rangle \approx - { c_s| \boldsymbol{\nabla }| \over \gamma \langle p\rangle } \big ( r_p p_{\textit{s,rms}} p_{\textit{d,rms}} + p_{\textit{d,rms}}^2 \big) = \langle \rho \rangle c_s^3 |\boldsymbol{\nabla }| \textit{Ma}_d^2 \left ( C_3 \mathcal{D} + C_4 \right )\! , \end{eqnarray}
\begin{eqnarray} {\langle p\vartheta \rangle \over \langle \rho \rangle \epsilon _s} \approx \mathcal{L} | \boldsymbol{\nabla }| { c_s^3 \over u_{\textit{s,rms}}^3} \textit{Ma}_d^2 \left ( C_3 \mathcal{D} + C_4 \right ) = \mathcal{L} | \boldsymbol{\nabla }| \left ({ u_{\textit{d,rms}} \over u_{\textit{s,rms}} }\right )^{3} { C_3 \mathcal{D} + C_4 \over \textit{Ma}_d} , \end{eqnarray}
where
$\mathcal{L}$
is the energy-containing scale with the undetermined coefficients
$C_3$
and
$C_4$
. As shown in figure 3(b), the model (represented by the line of the same colour) performs poorly for
$\chi =0$
, consistent with the same limitations of the
$p_{\textit{rms}}$
models (figure 2
b), but agreement improves systematically as bulk viscosity increases.
$\mathcal{L} |\boldsymbol {\nabla }|$
may vary with
$\textit{Re}_\lambda$
and
$\textit{Ma}_t$
, warranting further study.
Scaling of (a)
$\vartheta _{\textit{rms}}/\omega _{\textit{rms}}$
and (b)
$\epsilon _d/\epsilon _s$
with
$u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
.
Figure 4(a) shows the dependence of
$\vartheta _{\textit{rms}}/\omega _{\textit{rms}}$
on
$u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
. For a given bulk viscosity, the data follow a clear scaling
$\vartheta _{\textit{rms}}/\omega _{\textit{rms}}\sim (u_{\textit{d,rms}}/u_{\textit{s,rms}})^\beta$
, with the exponent
$\beta$
decreasing as
$\chi$
increases, indicating non-zero bulk viscosity alters the dilatational motions rather than simply rescaling for flows without bulk viscosity. This trend suggests that the ratio
$\epsilon _d/\epsilon _s$
should likewise depend on the velocity ratio. Figure 4(b) confirms this, showing that
$\epsilon _d/\epsilon _s$
collapses onto an empirical relation of the form
\begin{align} {\epsilon _d\over \epsilon _s} \approx C_5 \left ( u_{\textit{d,rms}}\over u_{\textit{s,rms}}\right )^{2C_6}, \, \lg {C_5} = {0.927 \chi + 18.6 \over \chi + 10.2} , \, C_6 = {0.694 \chi + 18.9 \over \chi + 9.46} , \, \chi \in [0,30] . \end{align}
Notably, the exponent
$C_6$
for all simulated cases consistently exceeds the empirical value of unity proposed by Donzis & John (Reference Donzis and John2020) and John & Donzis (Reference John and Donzis2024), and our coefficients may be influenced by the finite dataset and should therefore be regarded as provisional. Increasing
$\chi$
substantially amplifies the dilatational-dissipation contribution, demonstrating that
$\epsilon _d$
cannot be neglected in compressible flows with large bulk viscosity.
3.4. High-order velocity statistics
Finally, we examine higher-order statistics of the velocity gradient, namely the skewness
$S_u$
and flatness
$F_u$
. The former reflects the vortex-stretching dynamics (Champagne Reference Champagne1978), whereas the latter characterises the turbulent intermittency (Pope Reference Pope2000). In incompressible turbulence, both quantities exhibit only weak
$\textit{Re}_\lambda$
dependence, with empirical scalings (Ishihara et al. Reference Ishihara, Kaneda, Yokokawa, Itakura and Uno2007):
$-S_u=(0.32\pm 0.02)\textit{Re}_\lambda ^{0.11\pm 0.01}$
and
$F_u=(1.14\pm 0.19)\textit{Re}_\lambda ^{0.34\pm 0.03}$
. In compressible flows, however, previous studies (Donzis & John Reference Donzis and John2020) suggest that, once the condition
$\textit{Ma}_t u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2 \gt 3\times 10^{-2}$
is satisfied, shocklets induce extreme values of skewness and flatness (Wang et al. Reference Wang, Shi, Wang, Xiao, He and Chen2011), signalling their growing influence on turbulence statistics.
Skewness of the velocity gradient
$S_u$
with (a)
$\textit{Ma}_t$
and (b)
$\textit{Ma}_t u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
. The dotted line is the average value
$0.460$
for the cases before
$S_u$
-divergence.
Flatness of the velocity gradient
$F_u$
with (a)
$\textit{Ma}_t$
and (b)
$\textit{Ma}_t u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
. The dotted line is the average value
$4.95$
for the cases before
$S_u$
-divergence.
Figure 5 shows the plot of
$S_u$
. For most cases without `
$S_u$
-divergence’, the mean skewness is approximately
$-0.460$
, consistent with incompressible scaling at the present
$\textit{Re}_\lambda$
. For the different
$\chi$
cases simulated here, the comparison suggests that the critical threshold for `
$S_u$
-divergence’ may increase as
$\chi$
rises, and for
$\chi =30$
, no divergence is observed despite the phase diagram in figure 1(b). Flatness,
$F_u$
(figure 6), shows a similar trend, with the mean value
$4.95$
for non-divergent cases. Increasing bulk viscosity suppresses shocklet-induced intermittency, indicating that sufficiently large
$\chi$
restores near-incompressible statistics even at relatively high Mach numbers.
4. Discussion and conclusions
We conduct numerical simulations of forced homogeneous isotropic turbulence over varying bulk viscosities, Reynolds numbers and Mach numbers, and examine the scaling of Favre-averaged TKE (via Helmholtz decomposition), r.m.s. pressure, pressure dilatation, dilatational dissipation and higher-order velocity-gradient moments. New models are developed for pressure dilatation and for the bulk-viscosity dependence of dilatational dissipation.
Although the solenoidal and dilatational components of the Favre-averaged TKE are not strictly orthogonal, their ratio is well approximated by the squared ratio of the respective r.m.s. velocities, as confirmed numerically. The r.m.s. pressure departs from classical pseudo-sound predictions (Ristorcelli Reference Ristorcelli1997) but tends towards that limit as bulk viscosity increases. It remains close to the Donzis r.m.s. pressure model (Donzis & John Reference Donzis and John2020), with the solenoidal contribution increasingly dominant at larger bulk viscosity.
For compressible dissipation, we analyse pressure dilatation and dilatational dissipation. Pseudo-sound theory is found to be inaccurate across Mach numbers. Without bulk viscosity, pressure dilatation transfers energy mainly from TKE to mean internal energy; finite bulk viscosity can reverse this transfer at high Mach numbers. Extending Zeman's (Reference Zeman1991) approach and the Donzis pressure model (Donzis & John Reference Donzis and John2020), we propose a new pressure-dilatation model based on the velocity ratio, with improved accuracy at large bulk viscosity. For dilatational dissipation, we show that the scaling exponent in the Donzis velocity-ratio model is not universal but varies systematically with bulk viscosity, involving an amended exponent and prefactor.
We further evaluate the velocity-gradient skewness and flatness. Donzis & John (Reference Donzis and John2020) identified the shocklet-induced `
$S_u$
-divergence’ as governed by the parameter
$\textit{Ma}_t u_{\textit{d,rms}}^2/u_{\textit{s,rms}}^2$
. The divergence emerges once this parameter exceeds a critical threshold, but this threshold may increase with bulk viscosity; at sufficiently large bulk viscosity `
$S_u$
-divergence’ may be fully suppressed even when the empirical value is exceeded.
Overall, the velocity-ratio scaling framework offers a promising route toward universal laws for compressible turbulence. By incorporating bulk-viscosity effects, the present study refines this framework and delivers new modelling for TKE decomposition, pressure dilatation and dilatational dissipation. A natural extension is the development of separate transport equations for the solenoidal and dilatational components of TKE.
Supplementary material.
Supplementary material is available at https:/doi.org/10.1017/jfm.2026.11466.
Funding.
This work has been supported by the Chinese Academy of Sciences Project for Young Scientists in Basic Research (YSBR-107), the National Natural Science Foundation of China (NSFC award numbers 91852205, 42075071), NSFC Basic Science Center Program (Award number 11988102), Department of Science and Technology of Guangdong Province (2023B1212060001) and Guangdong-Hong Kong-Macao Joint Laboratory for Data-Driven Fluid Mechanics and Engineering Applications (2020B1212030001).
Declaration of interests.
The authors report no conflicts of interest.
Data availability statement.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
