# Baroclinic interaction of forced shock waves with random thermal gradients

Joaquim P. Jossy  Prateek Gupta<EMAIL_ADDRESS>https://web.iitd.ac.in/~prgupta/  Department of Applied Mechanics, Indian
Institute of Technology Delhi, New Delhi 110016, India

###### Abstract

Density gradients aligned at an angle to pressure gradients result in
baroclinic torque in fluid flows, generating vorticity. In this work, we study
the vorticity generated by the baroclinic torque exerted by the interaction of
pressure jumps across random two-dimensional shock waves with density
gradients. A field of random two-dimensional shock waves has acoustic spectral
energy scaling as $\widehat{E}_{k}\sim\varepsilon^{2/3}\ell^{-1/3}k^{-2}$
where $k$ is the wavenumber, $\varepsilon$ is the energy dissipation, and
$\ell$ is the integral length scale of the field. Since the acoustic energy is
broadband, pressure and velocity gradients exist in a wide range of length
scales. We study the interaction of these broadband gradients with isobaric
thermal gradients localized at a length scale in the spectral space. We show
that the method of generating shock waves or injection of wave energy in the
system governs the baroclinic interactions. For stochastically forced shock
waves, baroclinic terms are negligible. Broadband vorticity with energy at
least two orders of magnitude smaller is generated due to continuous variation
in curvature of shock waves caused by stochastic forcing. On the other hand,
shock waves maintained by energy rescaling result in the generation of
coherent vorticity. We also discuss the relative magnitude of the baroclinic
torque generated due to total density gradients compared to the one generated
due to non-isentropic density gradients within the shock waves interacting
with the pressure gradients.

††preprint: Preprint submitted to Physics of fluids on February 28, 2023

##  I Introduction

Finite amplitude or nonlinear acoustic waves exist in various engineering and
scientific applications such as aerospace Goldstein and Wundrow (1990);
Liepmann and Roshko (2001); Baars and Tinney (2013); Bonciolini, Boujo, and
Noiray (2017); Gupta, Lodato, and Scalo (2017), nondestructive testing
Breazeale and Thompson (1963); Lissenden (2021), astrophysics Sagdeev (1979),
nuclear energy Brouillette (2002), and medicine Lukes  _et al._ (2016).
Nonlinear effects such as acoustic streaming Lighthill (1978); Shimizu and
Sugimoto (2010) and wave steepening Whitham (2011); Gupta and Scalo (2018);
Thirani, Gupta, and Scalo (2020) govern the propagation of such finite
amplitude nonlinear acoustic waves in a compressible fluid. Due to steepening,
nonlinear acoustic waves form propagating shock waves Whitham (2011); Gupta
and Scalo (2018). Propagation, dissipation, and dispersion of these shock
waves can be altered by the thermodynamic properties of the medium in which
they propagate. Furthermore, interaction of shock waves with inhomogeneities
in the medium may generate additional hydrodynamic quantities of interest such
as vorticity. In this work, we study the propagation of forced shock waves in
an inhomogeneously heated medium in two dimensions.

Ellermeier (1993) showed that non-uniform cross-section and density
stratification affect the nonlinear distortion of planar waves. Tyagi and
Sujith (2003) found that entropy gradients adversely affect the nonlinear wave
steepening and shock formation. Prasad (2006) studied one-dimensional shock
waves in slowly varying isentropic one-dimensional flows. Using multiple
scales analysis, Prasad (2006) examined the effect of non-uniform flow in a
duct on evolution and transport of acoustic power in the duct and found that
the mean flow has an impact on the residence time of shock waves in the duct
due to viscous dissipation inside the shock waves. In this work, we focus on
the propagation of shock waves in a quiescent thermally inhomogeneous medium
in two dimensions. Due to the misalignment of pressure gradients across the
shock waves and density gradients, baroclinic torque generates vorticity. Such
misalignment may exist due to the thermoviscous dissipation inside the shock
waves or externally imposed thermal gradients. The vorticity caused by such
baroclinic interaction can be used to enhance fuel-air mixing in air-breathing
engines Budzinsky, Zukoski, and Marble (1992); Andreopoulos, Agui, and
Briassulis (2000); Romagnosi  _et al._ (2011); Tian  _et al._ (2017); Wong
_et al._ (2022) or active flow-control Singh  _et al._ (2019). We perform
shock-resolved direct numerical simulations (DNS) of two-dimensional randomly
generated shock waves and study their interaction with background thermal
gradients. We focus on the length scales and magnitude of vorticity generated
due to baroclinic interaction between pressure gradients across the shock
waves and background density stratification generated by thermal gradients.

In baroclinic flows, fluid density is a function of pressure, temperature, and
composition or concentration of the dissolved constituents Kundu, Cohen, and
Dowling (2015). This dependency can create a misalignment in the direction of
the density and pressure gradients, unlike barotropic flows in which the
density and pressure gradients are parallel. The resulting torque due to the
misalignment of density and pressure gradients is called baroclinic torque and
may result in increase or decrease of the vorticity and circulation in the
fluid, depending on the existing vorticity field. In inertially confined
fusion reactors, baroclinic torque results in Richtmyer-Meshkov instability
Mikaelian (1985); Brouillette (2002) which causes mixing of the fuel and
capsule material limiting the efficiency of the reactor. In supersonic
combustion chambers, shock-flame interaction results in vorticity generated
due to the resulting baroclinic torque which is argued as the primary cause of
deflagration to detonation transition Yang and Radulescu (2021). Additionally,
baroclinic production is attributed as the primary mechanism of enhanced
mixing downstream of an oblique shock wave as streamwise vorticity passes
through it Yu  _et al._ (2020); Wei  _et al._ (2022). Usually, only
stratification or density gradient is considered while analyzing the
baroclinic torque. However, in a multi-dimensional field of shock waves
propagating in a homogeneous medium, the density gradient is not necessarily
parallel to the pressure gradient due to local entropy generation.
Consequently, the interaction of pressure gradients with density gradients in
a two-dimensional or three-dimensional field of shock waves is possible
without prescribed or background stratification as well. In this work, we
analyze the interaction of two-dimensional field of random shock waves with
random background thermal gradients using two-dimensional shock-resolved
direct numerical simulation (DNS) of fully compressible Navier-Stokes
equations. Such background thermal gradients result in the stratification of
the fluid in which the shock waves propagate. We generate these shock waves
using two contrasting forcing methods. In one method, we force the momentum
equations using a stochastic process, also used to generate equilibrium
turbulence Eswaran and Pope (1988); Jagannathan and Donzis (2016). In the
other method, we rescale the density-weighted wave spectral energy in the
system in a band of small wavenumbers (large length scales). In the latter
forcing method, energy loss to smaller length scales due to acoustic spectral
energy cascade Gupta, Lodato, and Scalo (2017); Gupta and Scalo (2018) is
compensated at each time step of the simulation. Furthermore, the relative
phasing of shock waves is maintained since only amplitudes are rescaled. Using
numerical simulations, we show that the vorticity generated due to different
forcing methods is at different length scales. For random forcing, the
vorticity and enstrophy are broadband, while they are coherent for energy
rescaling. For energy rescaling, we discuss the enstrophy budgets to elucidate
the mechanism of coherent vorticity generation in the system. To the best of
our knowledge, this is the first study in which the interaction of random
shock waves with a inhomogeneously stratified medium has been studied.

In Section II we discuss the governing equations of our simulations and the
theoretical derivation of wave spectral energy equation using second-order
nonlinear perturbation equations. Since vorticity is expected to be generated
due to baroclinic interactions, we also discuss the decomposition of flow
fields in rotational and dilatational components to isolate the effects on the
spectral energy dynamics of waves. In Section III we discuss the numerical
setup and problem formulation. Particularly, we outline the two different
forcing methods used to maintain equilibrium with shock waves without any
preference of direction of propagation. We also show that for simulations with
no background thermal gradients, density wave spectral energy
$\widehat{E}^{f}_{wk}$ scales as $\varepsilon^{2/3}\ell^{-1/3}k^{-2}$ (see
Section III for definitions). In Section IV we discuss the length scales of
the two-dimensional field of random shock waves generated due to forcing, the
magnitude and length scales of vorticity generated, and the enstrophy budgets
before concluding in Section V.

##  II Flow field decomposition and spectral energy transfer

In this section, we discuss the theoretical background and the derivation of
nonlinear governing equations for perturbations in density, velocity, and
pressure field from fully compressible Navier-Stokes equations for an ideal
gas. Using these nonlinear perturbation equations, we derive the expressions
for the flux of spectral wave energy and spectral dissipation. To this end, we
discuss the possible decomposition of fields into the wave and the vortical
components.

Dimensionless fully compressible Navier-Stokes equations for an ideal gas are
given by,

$\displaystyle\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bm{u})=0,$ (1a)
$\displaystyle\frac{\partial\bm{u}}{\partial
t}+\bm{u}\cdot\nabla\bm{u}+\frac{\nabla p}{\rho}=\frac{1}{\rho
Re_{ac}}\left(\nabla\cdot(2\mu\bm{S})\right)+$ $\displaystyle\frac{1}{\rho
Re_{ac}}\left(\nabla\cdot\left(\kappa-\frac{2\mu}{3}\right)\bm{D}\right)+\bm{F},$
(1b) $\displaystyle\frac{\partial p}{\partial t}+\bm{u}\cdot\nabla p+\gamma
p\nabla\cdot\bm{u}=\frac{1}{Re_{ac}\hskip
2.27621pt\mathrm{Pr}}\nabla\cdot(\alpha\nabla T)+$
$\displaystyle\frac{\gamma-1}{Re_{ac}}\left(2\mu\bm{S}:\bm{S}+\left(\kappa-\frac{2\mu}{3}\right)\bm{D}:\bm{D}\right),$
(1c)

combined with the dimensionless ideal gas equation of state,

$\gamma p=\rho T.$ (2)

In Eqs. (1b) and (1c), $\bm{S}$ and $\bm{D}$ are the strain rate and
dilatation tensors, respectively, given by,

$\bm{S}=\frac{1}{2}\left(\nabla\bm{u}+\nabla\bm{u}^{T}\right),~{}\bm{D}=\nabla\cdot\bm{u}\bm{I}.$
(3)

In Eq. (1b), $\bm{F}$, is the external force vector which is used to maintain
random shocks in a two-dimensional domain and $\kappa$, $\mu$, $\alpha$,
$Re_{ac}$, and Pr are the dimensionless bulk viscosity, dimensionless dynamic
viscosity, dimensionless thermal conductivity, acoustic Reynolds number, and
Prandtl number, respectively. As discussed in detail in next section, the
external force vector $\bm{F}$ may represent the stochastic forcing or the
energy rescaling. Furthermore, to obtain the dimensionless Eqs. (1a), (1b),
and (1c), following scales have been used,

$[\rho]=\rho_{m},~{}[p]=\gamma
p_{m},~{}[T]=T_{m},~{}[\bm{u}]=(c_{m},c_{m}),~{}[x,y]=(L,L).$ (4)

Subscript $()_{m}$ in above relations denotes characteristic values
corresponding to the quiescent medium. Velocity scale $c_{m}$ denotes the
speed of sound at these characteristic values and $L$ corresponds to a length
scale. Using scales in Eq. (4), $Re_{ac}$ is defined as,

$Re_{ac}=\frac{\rho_{m}c_{m}L}{\mu_{m}},$ (5)

where $\mu_{m}$ is dynamic viscosity at $T_{m}$. Throughout this work, we
analyse results from numerical simulations of dimensionless equations. Hence,
values of scales in Eqs. (4) are irrelevant and only values of dimensionless
numbers are specified.

In this work, we study the interaction of random nonlinear acoustic fields
(weak shock waves Gupta and Scalo (2018), also see results in Section III.3)
with random thermal gradients. To this end, we perform DNS of Eqs. (1a)-(1c)
(see Section III) and extract relevant physical quantities using the second
order nonlinear governing equations for perturbations in density, velocity,
and pressure. To derive these equations, we consider a quiescent base state
with isobaric thermal gradients. The fields of dependent variables
$\bm{u},p,\rho,$ and $T$ can be written as a decomposition of the base field
and the random perturbations denoting the random nonlinear acoustic fields,

$\bm{u}=\bm{u}^{\prime},~{}p=\frac{1}{\gamma}+p^{\prime},~{}\rho=f(x,y)+\rho^{\prime},$
(6)

and

$T=\frac{1}{f(x,y)}+T^{\prime}.$ (7)

Field $1/f(x,y)$ corresponds to the base state thermal gradients (which may be
interpreted as base state stratification as well). For a homogeneous medium,
$f(x,y)=1$. As we discuss in Section III, we initialize these thermal
gradients within bands of wavenumbers which generates an inhomogeneous medium
in which the random shock waves are forced and maintained at statistical
steady state. .

To understand energy cascade and interaction with background thermal
gradients, we first derive nonlinear governing equations in $p^{\prime}$ and
$\bm{u}^{\prime}$. Since the thermoviscous terms can only have second order
contribution at the leading order (see Gupta and Scalo Gupta and Scalo (2018)
for a detailed discussion), we obtain the following dimensionless nonlinear
equations correct up to second order,

$\displaystyle\frac{\partial\bm{u}^{\prime}}{\partial
t}+\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}+\frac{\nabla
p^{\prime}}{f}-\frac{p^{\prime}\nabla p^{\prime}}{f^{2}}=$
$\displaystyle\frac{1}{fRe_{ac}}\left(\nabla\cdot\left(2\mu\bm{S}^{\prime}\right)+\nabla\cdot\left(\left(\kappa-\frac{2\mu}{3}\right)\bm{D}^{\prime}\right)\right),$
(8a) $\displaystyle\frac{\partial p^{\prime}}{\partial
t}+\bm{u}^{\prime}\cdot\nabla p^{\prime}+\nabla\cdot u^{\prime}+\gamma
p^{\prime}\nabla\cdot\bm{u}^{\prime}=\frac{1}{Re_{ac}\hskip
2.27621pt\mathrm{Pr}}\nabla\cdot(\alpha\nabla T^{\prime}).$ (8b)

In Eq. (8b), temperature perturbation $T^{\prime}$ is related to pressure and
density perturbations ($p^{\prime}$ and $\rho^{\prime}$, respectively) via the
dimensionless ideal gas equation of state (see Eq. (2)),

$\gamma
p^{\prime}=\left(f+\rho^{\prime}\right)\left(\frac{1}{f}+T^{\prime}\right)=1+\frac{\rho^{\prime}}{f}+fT^{\prime}+\rho^{\prime}T^{\prime}.$
(9)

Density perturbation $\rho^{\prime}$ is governed by,

$\frac{\partial\rho^{\prime}}{\partial
t}+\nabla\cdot\left(f\bm{u}^{\prime}\right)+\nabla\cdot\left(\rho^{\prime}\bm{u}^{\prime}\right)=0.$
(10)

A spectral energy equation is possible for perturbation energy only if a true
energy corollary can be derived for the governing system of equations Gupta
and Scalo (2018). To formulate the spectral energy equation from Eqs. (8a) and
(8b), we highlight that the linearized governing equations for inviscid
acoustic waves in a stratified medium,

$\displaystyle\frac{\partial\bm{u}^{\prime}}{\partial t}+\frac{\nabla
p^{\prime}}{f}=0~{}\mathrm{and}$ (11) $\displaystyle\frac{\partial
p^{\prime}}{\partial t}+\nabla\cdot\bm{u}^{\prime}=0,$ (12)

do not admit a modal decomposition for non-uniform $f$. For a homogeneous
background medium ($f(x,y)=1$ everywhere), the modal decomposition of the
linearized Eqs. (11) and (12) is given by,

$\left(u^{\prime},v^{\prime},p^{\prime}\right)^{T}=\sum_{n}\widehat{\bm{\phi}}_{\pm}e^{(i(k_{nx}x+k_{ny}y\pm\omega_{n}t)}+\sum_{n}\widehat{\bm{\phi}}_{0}e^{(i(k_{nx}x+k_{ny}y)}.$
(13)

where $u^{\prime}$ and $v^{\prime}$ are components of velocity field
perturbation $\bm{u}^{\prime}$. Modes $\widehat{\bm{\phi}}_{0}$ and
$\widehat{\bm{\phi}}_{\pm}$ are given by,

$\displaystyle\widehat{\bm{\phi}}_{0}=\frac{1}{|\bm{k}|}\left(-k_{ny},k_{nx},0\right)^{T},$
(14)
$\displaystyle\widehat{\bm{\phi}}_{\pm}=\frac{1}{|\bm{k}|\sqrt{2}}\left(k_{nx},k_{ny},\pm|\bm{k}|\right)^{T},$
(15)

where $\widehat{\bm{\phi}}_{0}$ represents the vortical mode which has no
contributions from the pressure perturbations and $\widehat{\bm{\phi}}_{\pm}$
represent the acoustic modes. Since
$\widehat{\bm{\phi}}_{0},\widehat{\bm{\phi}}_{\pm}$ are orthogonal to each
other, any two-dimensional compressible flow field can be projected on
$\widehat{\bm{\phi}}_{0},\widehat{\bm{\phi}}_{\pm}$. For second order
nonlinear Eqs. (8a)- (8b), writing the instantaneous solution as an expansion
of $\widehat{\bm{\phi}}_{0},\widehat{\bm{\phi}}_{\pm}$, we obtain,

$\displaystyle\left(u^{\prime},v^{\prime},p^{\prime}\right)^{T}=$
$\displaystyle\sum_{n}\Big{(}\widehat{\Omega}(\bm{k}_{n},t)\widehat{\bm{\phi}}_{0}$
$\displaystyle+\widehat{R}_{+}(\bm{k}_{n},t)\widehat{\bm{\phi}}_{+}+\widehat{R}_{-}(\bm{k}_{n},t)\widehat{\bm{\phi}}_{-}\Big{)}e^{i(k_{nx}x+k_{ny}y)}.$
(16)

Here $\widehat{\Omega}(\bm{k}_{n},t)$ corresponds to the vortical contribution
of the perturbation field and $\widehat{R}_{+}(\bm{k}_{n},t),$ and
$\widehat{R}_{-}(\bm{k}_{n},t)$ correspond to the acoustic contribution.
Equation (16) yields a spectral energy density definition (since the modes
$\widehat{\bm{\phi}}_{0},\widehat{\bm{\phi}}_{\pm}$ are orthogonal) as,

$\widehat{E}_{\bm{k}}(t)=\frac{1}{2}|\widehat{\Omega}(\bm{k}_{n},t)|^{2}+\frac{1}{2}|\widehat{R}_{+}(\bm{k}_{n},t)|^{2}+\frac{1}{2}|\widehat{R}_{-}(\bm{k}_{n},t)|^{2}.$
(17)

Since there is no such modal decomposition possible for waves in an
inhomogeneous medium, one may follow Miura and Kida (1995) and use the spatial
Fourier modes of density-weighted velocity perturbations
$\bm{w}^{\prime}=\sqrt{f}\bm{u}^{\prime}$ to define the spectral energy as,

$\displaystyle\widehat{E}_{\bm{k}}$
$\displaystyle=\frac{1}{2}\left(\widehat{\bm{w}}^{*}_{\bm{k}}\cdot\widehat{\bm{w}}_{\bm{k}}+\widehat{p}^{*}_{\bm{k}}\widehat{p}_{\bm{k}}\right)$
$\displaystyle=\underbrace{\frac{1}{2}\left(\widehat{\bm{w}}^{*}_{v\bm{k}}\cdot\widehat{\bm{w}}_{v\bm{k}}\right)}_{\widehat{E}^{f}_{v\bm{k}}}+\underbrace{\frac{1}{2}\left(\widehat{\bm{w}}^{*}_{w\bm{k}}\cdot\widehat{\bm{w}}_{w\bm{k}}+\widehat{p}^{*}_{\bm{k}}\widehat{p}_{\bm{k}}\right)}_{\widehat{E}^{f}_{w\bm{k}}},$
(18)

where $()^{*}$ denotes complex conjugate, $\widehat{\bm{w}}_{\bm{k}}$ is the
Fourier transform of $\bm{w}^{\prime}$ and
$\bm{w}_{v{\bm{k}}},\bm{w}_{w{\bm{k}}}$ are obtained using the Helhmholtz
decomposition in the Fourier space,

$\widehat{\bm{w}}_{\bm{k}}=\widehat{\bm{w}}_{v{\bm{k}}}+\widehat{\bm{w}}_{w{\bm{k}}},\\\
$ (19a)
$i\bm{k}\cdot\widehat{\bm{w}}_{\bm{k}}=i\bm{k}\cdot\widehat{\bm{w}}_{w{\bm{k}}},~{}\mathrm{and}\\\
$ (19b)
$i\bm{k}\times\widehat{\bm{w}}_{\bm{k}}=i\bm{k}\times\widehat{\bm{w}}_{v{\bm{k}}}.$
(19c)

Energy can be decomposed as $\widehat{E}^{f}_{w\bm{k}}$ and
$\widehat{E}^{f}_{v\bm{k}}$ where the superscript $()^{f}$ denotes energy
evaluated using Fourier modes of $\bm{w}^{\prime}$. However, as we show in
further sections, $E^{f}_{v\bm{k}}$ is not a good choice for quantifying the
actual vortical energy in the system corresponding to the local rotation rate
of the fluid. Hence, we reserve $\widehat{E}_{v\bm{k}}$ for the _actual
vortical energy_ defined as,

$\widehat{E}_{v\bm{k}}=\frac{1}{2}\widehat{\bm{u}}^{*}_{v\bm{k}}\cdot\widehat{\bm{u}}_{v\bm{k}},$
(20)

where $\widehat{\bm{u}}_{v\bm{k}}$ satisfies,

$i\bm{k}\times\widehat{\bm{u}}_{v\bm{k}}=i\bm{k}\times\widehat{\bm{u}}_{\bm{k}}.$
(21)

In sections below, we use $\widehat{E}^{f}_{wk}$ to represent the spectral
energy in waves. In the physical space, perturbation energy equation can be
written as,

$\displaystyle\frac{\partial}{\partial
t}\left(\frac{f|\bm{u}^{\prime}|^{2}}{2}+\frac{p^{\prime
2}}{2}\right)+\nabla\cdot\left(\bm{u}^{\prime}p^{\prime}\right)+f\bm{u}^{\prime}\cdot\left(\bm{u}^{\prime}\cdot\nabla\bm{u}^{\prime}\right)+$
$\displaystyle\gamma p^{\prime
2}\nabla\cdot\bm{u}^{\prime}=f\bm{u}^{\prime}\cdot D_{u}+p^{\prime}D_{p},$
(22)

where $D_{u}$ and $D_{p}$ correspond to the thermoviscous dissipation terms
from Eq. (8a) and Eq. (8b), respectively. Equation (22) is obtained by taking
the dot product of Eq. (8a) with $f\bm{u}^{\prime}$ and Eq. (8b) with
$p^{\prime}$.

As shown by Gupta and Scalo (2018), the $\gamma p^{\prime
2}\nabla\cdot\bm{u}^{\prime}$ term can be further decomposed into an infinite
series expansion to yield an exact energy corollary for planar one-dimensional
nonlinear acoustic waves. However, such a decomposition is hard in higher
dimensions even for waves propagating in a homogeneous medium. Consequently we
choose to work with the remaining terms for obtaining the spectral energy
cascade terms. In spectral space, the energy equation can be represented as,

$\frac{d\widehat{E}_{\bm{k}}}{dt}+\widehat{\mathcal{T}}_{\bm{k}}=\widehat{\mathcal{D}}_{\bm{k}},$
(23)

where $\widehat{\mathcal{T}}_{\bm{k}}$ is the spectral energy transfer
function and $\widehat{\mathcal{D}}_{\bm{k}}$ is the spectral dissipation
function. From numerical simulations, we use discretized spectral quantities
($\widehat{E}_{\bm{k}},\widehat{\mathcal{T}}_{\bm{k}},\widehat{\mathcal{D}}_{\bm{k}}$
and derived quantities) to investigate spectral energy dynamics using the
corresponding binned quantities as,

$\widehat{\Phi}_{k}=\sum^{|\bm{k}|<k+\Delta k/2}_{|\bm{k}|>k-\Delta
k/2}\Phi_{\bm{k}}=b(\widehat{\Phi}_{\bm{k}}),$ (24)

where $\Phi$ may denote spectral energy, spectral energy transfer, spectral
dissipation or derived quantities. Binned spectral energy transfer and
spectral dissipation are obtained by multiplying Eq. (8a) with $\sqrt{f}$ and
then taking the dot-product with $\widehat{w}^{*}_{wk}$ in the spectral space
(and similarly multiplying Eq. (8b) with $\widehat{p}^{*}_{k}$ in the spectral
space) as,

$\displaystyle\widehat{\mathcal{T}}_{k}=b\Bigg{(}\mathrm{Re}\left[\widehat{\bm{w}}^{*}_{\bm{k}}\cdot\left(\widehat{\frac{\nabla
p^{\prime}}{\sqrt{f}}}\right)_{\bm{k}}+\widehat{p}^{*}_{\bm{k}}\widehat{\nabla\cdot\bm{u}^{\prime}}_{\bm{k}}\right]+$
$\displaystyle\mathrm{Re}\left[\widehat{\bm{w}}^{*}_{\bm{k}}\cdot\left(\widehat{\bm{w}\cdot\nabla\bm{u}^{\prime}}\right)_{\bm{k}}-\widehat{\bm{w}}^{*}_{\bm{k}}\cdot\left(\widehat{\frac{p^{\prime}\nabla
p^{\prime}}{\sqrt{f}}}\right)_{\bm{k}}\right]+$
$\displaystyle\mathrm{Re}\left[\gamma\widehat{p}^{*}_{\bm{k}}\left(\widehat{p^{\prime}\nabla\cdot\bm{u}^{\prime}}\right)_{\bm{k}}+\widehat{p}^{*}_{\bm{k}}\left(\widehat{\bm{u}^{\prime}\cdot\nabla
p^{\prime}}\right)_{\bm{k}}\right]\Bigg{)},$ (25)

and

$\begin{split}\widehat{\mathcal{D}}_{k}=b\Bigg{(}\mathrm{Re}\left[\widehat{\bm{w}}^{*}_{\bm{k}}\cdot\left(\widehat{\frac{\mu}{\sqrt{f}\hskip
2.27621ptRe_{ac}}\left(\nabla\cdot(\nabla u^{{}^{\prime}}+\nabla u^{\prime
T})\right)}\right)_{\bm{k}}\right]+\\\
\mathrm{Re}\left[\widehat{\bm{w}}^{*}_{\bm{k}}\cdot\left(\widehat{\frac{1}{\sqrt{f}\hskip
2.27621ptRe_{ac}}\nabla\cdot\left(\kappa-\frac{2\mu}{3}(\nabla\cdot
u^{\prime}))\right)}\right)_{\bm{k}}\right]+\\\
\mathrm{Re}\left[\widehat{\bm{p}}^{*}_{\bm{k}}\cdot\left(\widehat{\frac{1}{Re_{ac}\hskip
2.27621pt\mathrm{Pr}}\nabla\cdot(\alpha\nabla
T^{\prime})}\right)_{\bm{k}}\right]\Bigg{)},\end{split}$ (26)

respectively. For uniform background temperature and density ($f=1$), Eq. (25)
yields the transfer function for a nonlinear acoustic system in a homogeneous
medium,

$\displaystyle\widehat{\mathcal{T}}_{k}=b\Bigg{(}\mathrm{Re}\left[\widehat{\bm{u}}^{*}_{\bm{k}}\cdot\left(\widehat{\bm{u}\cdot\nabla\bm{u}^{\prime}}\right)_{\bm{k}}-\widehat{\bm{u}}^{*}_{\bm{k}}\cdot\left(\widehat{p^{\prime}\nabla
p^{\prime}}\right)_{\bm{k}}\right]+$
$\displaystyle\mathrm{Re}\left[\gamma\widehat{p}^{*}_{\bm{k}}\left(\widehat{p^{\prime}\nabla\cdot\bm{u}^{\prime}}\right)_{\bm{k}}+\widehat{p}^{*}_{\bm{k}}\left(\widehat{\bm{u}^{\prime}\cdot\nabla
p^{\prime}}\right)_{\bm{k}}\right]\Bigg{)}.$ (27)

Throughout this work, we quantify the results from numerical simulations of
fully nonlinear governing Eqs. (1a)-(1c) (discussed in Section III) using
second order nonlinear Eqs. (8a) and (8b) and corresponding spectral energy
cascade quantities derived.

##  III Numerical Simulations

In this section we discuss the computational setup for our numerical
simulations of two-dimensional fully compressible dimensionless Navier-Stokes
Eqs. (1a)-(1c). We solve Eqs. (1a)-(1c) numerically using the fourth order
Runge-Kutta time stepping and standard Fourier pseudospectral method Boyd
(2001) combined with $2/3$ deliasing for numerical stability and slab-
decomposition for parallel distributed computing using MPI Mortensen and
Langtangen (2016). For all the simulations presented in this work, we consider
$3072\times 3072$ Fourier modes ensuring shock-resolved DNS (see Section III.3
and IV.1 for details). Smooth variation of binned spectral energy, spectral
energy transfer function, and cumulative spectral dissipation at higher
wavenumbers also highlight the shock-resolved nature of the simulations. For
all simulations, we consider constant viscosity ($\mu=1$) and thermal
conductivity ($\alpha=1$). We show in Appendix A that temperature dependent
viscosity and thermal conductivity have no significant affect on shock wave
field. Since the shock waves in our simulations are weak, temperature
fluctuations and hence viscosity fluctuations are very small compared to their
respective base state values to have any significant impact. Furthermore, for
simplicity, we ignore bulk viscosity in our simulations ($\kappa=0$).

At time $t=0$, we prescribe quiescent conditions $\bm{u}=0$ along with
isobaric thermal gradients. To generate a thermal gradient at a wavevector
$\bm{k}$ randomly, we define the initial temperature field $1/f(x,y)$ in Eq.
(7) as,

$\frac{1}{f(x,y)}=1+\epsilon_{T}T^{b}(x,y),$ (28)

where $T^{b}(x,y)$ is defined as a random periodic field for a wavenumber
vector $\bm{k}$ in the Fourier space as,

$\widehat{T}^{b}_{\bm{k}}=\exp(2\pi i\theta_{\bm{k}}),$ (29)

where $\theta_{\bm{k}}$ is drawn from a normal distribution over $(0,1]$ for
every wavenumber vector $\bm{k}$. Since these thermal gradients are assumed to
be isobaric, the initial density field is set to $f(x,y)$ ensuring no pressure
perturbations due to thermal initialization. To study the interaction of these
thermal gradients with random shock waves, we consider two types of forcings
namely random forcing and energy rescaling as discussed below.

###  III.1 Random forcing

(a)(b)

Figure 1: (a) Instantaneous contours of $-\nabla\cdot\bm{w}^{\prime}$ showing
the two-dimensional field of shock waves generated by stochastic forcing and
(b) instantaneous contours of $\nabla\times\bm{u}^{\prime}$ showing the
vorticity field generated due to these shock waves in a inhomogeneous medium
(for case 6d in Table 1) at different times.

To generate random shocks, we force acoustic waves stochastically using the
Uhlenbeck-Ornstein (UO) process Eswaran and Pope (1988), typically used for
generating homogeneous isotropic turbulence in a three-dimensional box. The
forcing $\bm{F}$ in Eq. (1b) is defined in the Fourier space using four
independent solutions of a UO process $a_{i,\bm{k}}(t)$ for $i=1,2,3,$ and $4$
(one complex number in each direction). Following conditions are imposed on
each of the $a_{i}$,

$\displaystyle\left\langle a_{i,\bm{k}}(t)\right\rangle=0,$ (30)
$\displaystyle\left\langle{a}_{i,\bm{k}}(t){a}_{j,\bm{k}}^{*}(t+s)\right\rangle=2\sigma^{2}\delta_{ij}\exp(-s/T_{L}),$
(31)

where $\left\langle\cdot\right\rangle$ denotes the ensemble average, $()^{*}$
denotes the complex conjugate, $\sigma^{2}$ is the variance, and $T_{L}$ is
the forcing time-scale. We restrict $\bm{F}$ to a band of wavenumbers
$0<|\bm{k}|<k_{F}$. Spectral wave energy cascades from the forcing band to the
higher wavenumbers generating shock waves (see Fig. 1). In this work, we
choose $k_{F}=5$ ensuring energy injection in large length scales Yeung,
Sreenivasan, and Pope (2018); Bell  _et al._ (2022). The total rate of energy
addition $\epsilon$ in the system depends on the variance $\sigma^{2}$, time-
scale $T_{L}$, and number of wavenumber vectors within the forcing band
$N_{F}$ as,

$\displaystyle\epsilon=4N_{F}T_{L}\sigma^{2}.$ (32)

Combining independent processes $a_{1},a_{2},a_{3}$, and $a_{4}$, the forcing
vector in the Fourier space can be obtained as,

$\widehat{\bm{a}}_{\bm{k}}=\left(a_{1,\bm{k}}+ia_{2,\bm{k}},a_{3,\bm{k}}+ia_{4,\bm{k}}\right)^{T}.$
(33)

To force only the wave field, we remove any vortical component from
$\bm{a}_{\bm{k}}$ yielding the forcing $\bm{F}$ introduced in Eq. (1b) in the
Fourier space as,

$\displaystyle\widehat{\bm{F}}_{\bm{k}}(t)=\frac{(\widehat{\bm{a}}_{\bm{k}}(t)\cdot\bm{k})}{|\bm{k}|^{2}}\bm{k}.$
(34)

Consequently, Eq. (32) yields an upper bound on the rate of energy addition
rather than holding exactly. Even though we remove the vortical component from
the forcing, randomly forced acoustic wave turbulence exhibits broadband
vortical energy (see Fig. 2) because the forcing $\bm{F}$ (chosen from a
stochastic process) may point in a different direction at different times.
Consequently, acoustic density gradient $\nabla\rho(t_{1})$ generated by
$\bm{F}(t_{1})$ can interact with $\bm{F}(t_{2})$ at time $t_{2}>t_{1}$, when
$\bm{F}(t_{2})$ is not parallel to $\nabla\rho(t_{1})$, thus changing the
curvature in the shocks generating vorticity. Figure 1 shows the two-
dimensional field of random shock waves generated due to the stochastic
forcing and the vorticity generated by the shock waves interacting with the
background thermal gradients ($k_{T}=60$) and the forcing vector $\bm{F}$.

Figure 2: Wave energy $\widehat{E}_{wk}$, and actual vortical energy
$\widehat{E}_{vk}$ spectra of stochastically forced simulations for case 6d in
Table 1. Stochastic forcing generates broadband vorticity with highest actual
vortical energy in forced wavenumbers due to continuously changing curvature
of shock waves.

###  III.2 Density weighted energy rescaling

Since the forcing parameter $\widehat{\bm{a}}_{\bm{k}}$ in Eq. (33) is a
solution of a Langevin equation, a randomly forced wave turbulence can be
compared to the Langevin thermostat in the molecular simulations J Evans and P
Morriss (2007). At every time-step, components of $\widehat{\bm{a}}_{k}$
satisfy the Langevin equation independently. Hence, phase of the acoustic
waves generated by stochastic forcing is random. Consequently, vorticity
generated as a result of the baroclinic torque due to shock waves propagating
in one direction is quickly nullified by shock waves propagating in opposite
direction at some later time. Hence, vorticity is not accumulated in the
domain at any particular length scale due to baroclinic torque. However, due
to interaction of the forcing vector $\bm{F}$ with the shock waves, a
broadband vorticity with higher energy at forced length scales is generated.

Based on this reasoning, we use an _energy rescaling_ forcing, which is
inspired by the velocity rescaling thermostat in molecular simulations J Evans
and P Morriss (2007) and the deterministic forcing used in hydrodynamic
turbulence Ghosal  _et al._ (1995). For one set of simulations with different
length scales of background thermal gradients, we restart the simulations by
disabling random forcing $\bm{F}$ after statistical stationarity is achieved
approximately. In simulations with energy rescaling, density-weighted wave
spectral energy $\widehat{E}^{f}_{wk}$ is rescaled in a band of wavenumbers
$\mathcal{K}$ by forcing the contribution of smaller bins of width $\Delta k$
in $\mathcal{K}$ uniform. For instance, let a band of wavenumbers
$\mathcal{K}$ defined as $\mathcal{K}=\\{\bm{k}:k_{1}<|\bm{k}|<k_{2}\\}$ be
composed of partitions
$\mathcal{K}_{i}=\\{\bm{k}:k_{i}-1/2<|\bm{k}|<k_{i}+1/2\\}$ where $k_{i}$ is
between $k_{1}$ and $k_{2}$. We rescale the variables
$\widehat{\bm{w}}_{w\mathcal{K}_{i}}$ and $\widehat{p}_{\mathcal{K}_{i}}$ to
$\tilde{\widehat{\bm{w}}}_{w\mathcal{K}_{i}}$ and
$\tilde{\widehat{p}}_{\mathcal{K}_{i}}$, respectively,

$\displaystyle\tilde{\widehat{\bm{w}}}_{w\mathcal{K}_{i}}=\epsilon_{A}\frac{\alpha_{k_{i}}\widehat{\bm{w}}_{w\mathcal{K}_{i}}}{\sqrt{\frac{1}{2}\sum_{j}\alpha^{2}_{k_{j}}\left(|\widehat{\bm{w}}_{w\mathcal{K}_{j}}|^{2}+|\widehat{p}_{\mathcal{K}_{j}}|^{2}\right)}},$
(35)
$\displaystyle\tilde{\widehat{p}}_{\mathcal{K}_{i}}=\epsilon_{A}\frac{\alpha_{k_{i}}\widehat{p}_{\mathcal{K}_{i}}}{\sqrt{\frac{1}{2}\sum_{j}\alpha^{2}_{k_{j}}\left(|\widehat{\bm{w}}_{w\mathcal{K}_{j}}|^{2}+|\widehat{p}_{w\mathcal{K}_{j}}|^{2}\right)}},$
(36)

where the weight coefficients $\alpha_{k_{i}}$ are chosen as,

$\alpha_{k_{i}}=\frac{1}{\sqrt{\frac{N_{i}}{2}\sum_{\mathcal{K}_{i}}\left(|\widehat{\bm{w}}_{w\mathcal{K}_{i}}|^{2}+|\widehat{p}_{\mathcal{K}_{i}}|^{2}\right)}},$
(37)

where $N_{i}$ is the number of wavenumber vectors $\bm{k}$ within the
partition $\mathcal{K}_{i}$. Thus, all the partitions $\mathcal{K}_{i}$
contribute equally to the total wave energy in $\mathcal{K}$. The total wave
energy in $\mathcal{K}$ is equal to $\epsilon^{2}_{A}$ after every rescaling
operation. Rescaling of wave components of $\bm{w}^{\prime}$ in the Fourier
space results in rescaling of the amplitudes in the real space, maintaining
the propagation direction of the waves at each time step. Furthermore, since
we are forcing the density-weighted wave components $\bm{w}^{\prime}_{w}$,
energy rescaling also results in external injection of vorticity. However, as
we show in Section IV.3, the baroclinic interaction between the shock waves
and the background thermal gradients are also significant thus resulting in
_coherent_ vorticity, unlike random forcing which results in broadband
vorticity.

###  III.3 Simulation parameters

Figure 3: Time series of $E^{f}_{w}$ for all the simulation cases in Table 1
and the corresponding baseline cases in Table 2. All simulations are run with
stochastic forcing till $t=100$. Only simulations 1d - 6d are restarted with
the density-weighted energy rescaling (marked with r in Table 1). Due to
forcing, spatially averaged pressure and temperature rise in the domain. The
shown time series are obtained using pressure perturbations which are
calculated by subtracting the spatial average of pressure from total pressure.

We define our simulation parameter space as the wavenumber of thermal
perturbations $k_{T}$ and the rate of energy addition $\varepsilon$. We choose
$k_{T}$ values between $10$ and $60$ (increments of $10$) and
$\varepsilon/4N_{f}$ values as $1.5\times 10^{-5}$, $2.0\times 10^{-5}$, and
$2.5\times 10^{-5}$ (see Table 1). We run these simulations till $t=100$ and
then restart with energy rescaling enabled and stochastic forcing disabled to
isolate the coherent vorticity generated due to the interaction between the
shock waves and background thermal gradients. Since the parameter
$\varepsilon$ does not influence rescaled runs, we only restart simulations
with $\varepsilon=1.5\times 10^{-5}$ with energy rescaling (1dr-6dr in Table
1). For each $\varepsilon$ value, we also run a baseline case with no thermal
gradients ($f=1$). Figure 3 shows time series of density-weighted wave energy
${E}^{f}_{w}$ for all simulations in Tables 1 and 2. For simulations with
stochastic forcing, spatially averaged pressure and temperature rise due to
constant dissipation of energy caused by shock waves in the system. In Fig. 3,
spatially averaged pressure is removed from $E^{f}_{w}$, thus showing the true
perturbation energy in the system. For all cases with random gradients,
thermal gradients are defined for a range of wavenumber vectors
$k_{T}-\frac{1}{2}<|\bm{k}|<k_{T}+\frac{1}{2}$.

Table 1: Simulation parameter space considered in this study and the corresponding simulation names. We run all the simulations using random forcing till $t=100$ and then restart the simulation cases marked with r with energy rescaling to isolate the coherent baroclinic vorticity. The rate of energy injection parameter $\varepsilon$ loses its significance in rescaled forcing simulations since there is no random forcing. Other fixed parameters and their values are : $\varepsilon_{A}=0.15,~{}\varepsilon_{T}=0.1,~{}Re_{ac}=2500,$ and$~{}\mathrm{Pr}=0.72,~{}$.  |  | $k_{T}$
---|---|---
10 | 20 | 30 | 40 | 50 | 60
$\epsilon/4N_{F}$ | $1.5\times 10^{-5}$ | 1dr | 2dr | 3dr | 4dr | 5dr | 6dr
$2.0\times 10^{-5}$ | 1e | 2e | 3e | 4e | 5e | 6e
$2.5\times 10^{-5}$ | 1f | 2f | 3f | 4f | 5f | 6f
Table 2: Baseline simulation cases for $f=1$.  $\epsilon/4N_{F}\times 10^{5}$ | $0.25$ | $0.5$ | $1.0$ | 1.5 | 2.0 | 2.5
---|---|---|---|---|---|---
| 0a | 0b | 0c | 0d | 0e | 0f

The acoustic Reynolds number $Re_{ac}$ and the rate of energy injection
$\epsilon$ (equal to rate of energy dissipation at stationarity) are related
as Gupta and Scalo (2018); Augier, Mohanan, and Lindborg (2019),

$Re_{ac}\sim\frac{1}{\eta}(\epsilon\ell)^{-1/3},$ (38)

where $\eta$ is the smallest length scale (Kolmogorov length scale in
hydrodynamic turbulence or shock thickness) and $\ell$ is the integral length
scale (average distance between two shocks). Since 3072 Fourier modes are used
in each direction with $2/3$ dealiasing, maximum wavenumber captured in the
simulations is $k_{max}=1024$. Hence, for shock-resolved simulations,
$k_{max}\eta>1$ must hold (see Section IV.1). We choose $Re_{ac}=2500$ in our
simulations such that $\eta>1/1024$ for all the values of $\varepsilon$ and
the field of shock waves is well resolved (see Table 3). For all the values of
$\varepsilon$, the shock waves in the field are weak shock waves (see Fig.
4a).

As the rate of energy injection $\epsilon$ increases, stronger shocks are
generated in the domain. We calculate the mean Mach number $\langle M\rangle$
using the entropy perturbation at each point in the simulation domain and
calculating the weighted mean of the shock Mach number. At each point, we
evaluate the dimensionless entropy perturbation $s^{\prime}$ using,

$s^{\prime}=\frac{1}{\gamma-1}\ln\left(\gamma
p^{\prime}\right)-\frac{\gamma}{\gamma-1}\ln\left(\frac{\rho^{\prime}}{f}\right),$
(39)

using which, we approximate the entropy jump at each point as,

$\Delta s^{\prime}=\Bigg{|}\frac{\partial s^{\prime}}{\partial
x}\Bigg{|}dx+\Bigg{|}\frac{\partial s^{\prime}}{\partial y}\Bigg{|}dy.$ (40)

Using $\Delta s^{\prime}$, we calculate the Mach number field assuming weak
shocks as,

$M=\sqrt{1+\left(\frac{3\Delta
s^{\prime}(\gamma+1)^{2}}{2\gamma}\right)^{1/3}}.$ (41)

(a)(b)

Figure 4: (a) Time averaged normalized histogram (PDF) of Mach number, (b)
variation of weighted mean Mach number $\langle M\rangle$ with forcing
parameter $\epsilon$ for baseline cases 0a-0f (see Table 2).

In Fig. 4 we show the normalized histogram or probability distribution
function (PDF) of Mach number $M$ for all the baseline simulations. Since the
shock waves are very thin, most of the domain exhibits Mach number very close
to 1 ($s^{\prime}\approx 0$). To obtain a single numerical value
representative of the Mach number of shock waves, we compute weighted average
of the Mach number $\langle M\rangle$, such that contribution of points where
$s^{\prime}<25\%$ of the maximum entropy perturbation is set to 0. Figure 4(b)
shows the variation of weighted mean Mach number $\langle M\rangle$ of shock
waves in the domain with the forcing energy injection rate $\epsilon$. Since
the shocks become stronger for higher energy injection rate for same forcing
wavenumber band, the smallest length scale of fields decreases (see Eq. (38))
and the mean Mach number of shocks increases. Even for the highest injection
rate considered in this study (case 0f) the shocks are weak ($\langle
M\rangle<1.2$). Thus the governing equations correct up to second order in
perturbations (Eqs. (8a)-(8b)) are sufficient to study the cascade of spectral
wave energy. To elucidate the generation of vorticity, we consider the fully
nonlinear Eq. (1b) (see Section IV.2). Below, we discuss the features of field
of shock waves generated using stochastic forcing and the rescaled forcing in
detail. We study the length scales of the field of shock waves following which
we focus on vorticity generated due to the two forcing schemes and finally
discuss the enstrophy budgets.

##  IV Results and Discussion

In this section, we discuss the interaction of stochastically forced shocks
and shocks maintained by energy rescaling with the background thermal
gradients. We show that broadband vorticity is generated due to stochastically
forced shocks irrespective of the presence of the background thermal
gradients. Shocks maintained by energy rescaling interact with the background
thermal gradients and generate coherent vorticity, with eddies of the same
length scale as that of the thermal gradient. We also show that the background
thermal gradients do not significantly affect the length scales of the shocks
and spectral energy cascade in shocks.

###  IV.1 Length scales

(a)(b)(c)

Figure 5: (a) Dimensionless scaled density-weighted wave energy spectra
$\widehat{E}^{f}_{wk}$, (b) spectral flux of wave energy $\widehat{\Pi}_{wk}$,
and (c) cumulative dissipation $\widehat{\mathcal{D}}_{wk}$ for baseline cases
(0b-0f in Table 2). The wave energy spectra $\widehat{E}^{f}_{wk}$ scale with
dissipation $\varepsilon$ and integral length scale $\ell$ as in Eq. (42).

Due to shock waves, the wave energy spectra decay as $k^{-2}$ in the
wavenumber range where spectral energy flux is high. Gupta and Scalo (2018)
derived such scaling for acoustic wave energy spectra for decaying one-
dimensional acoustic wave turbulence in a homogeneous medium as,
$\widehat{E}_{wk}\sim\varepsilon^{2/3}\ell^{-1/3}k^{-2}$
($\widehat{E}^{f}_{wk}=\widehat{E}_{wk}$ for $f=1$). Augier et al. Augier,
Mohanan, and Lindborg (2019) derived a similar scaling for stochastically
forced isotropic shocks in shallow water equations. In Fig. 5 we show the
normalized density-weighted wave energy spectra, spectral flux, and cumulative
spectral dissipation for baseline cases 0a-0f (see Table 2). The wave energy
spectra scale as $k^{-2}$ with a similar dependence on dissipation and
integral length scale as for decaying one-dimensional acoustic wave turbulence
Gupta and Scalo (2018) and two-dimensional shallow water wave turbulence
Augier, Mohanan, and Lindborg (2019), indicating that the scaling relation,

$\frac{\widehat{E}^{f}_{wk}}{\ell^{-5/3}\epsilon^{2/3}}\sim\frac{1}{k^{2}\ell^{2}},$
(42)

holds where the spectral energy flux is high ($10\lesssim k\ell\lesssim 300$).
It is not possible to term the range of scales where the spectral energy flux
is high as _inviscid range_ or _inertial range_ because the dissipation exists
at all scales in shock waves propagating with a physical diffusion term. In
Appendix A we show the effect of hyperdiffusion (widely used in two-
dimensional wave turbulence literature) and how it modifies the distribution
of spectral energy cascade flux $\widehat{\Pi}_{k}$ in shock waves,
particularly concentrating the dissipation at small scales. Since the wave
energy spectra decay as $k^{-2}$ in a field of shock waves, dissipation is
active at all length scales. Consequently, cascade of wave energy spectra in
compressible flows may be regarded as a variable flux energy transfer Verma
(2019, 2021). The scaling in Eq. (42) also indicates that the shocks generated
by stochastic forcing in all the simulations are isotropic Augier, Mohanan,
and Lindborg (2019).

(a)(b)(c)

Figure 6: (a) Dimensionless scaled density-weighted wave energy spectra
$\widehat{E}^{f}_{wk}$, (b) spectral flux of wave energy $\widehat{\Pi}_{wk}$,
and (c) cumulative dissipation $\widehat{\mathcal{D}}_{wk}$ for inhomogeneous
cases (1d-6d in Table 1 and the baseline case (0d in Table 2). The density-
weighted wave energy spectra are almost identical for all the cases and follow
the scaling in Eq. (42). The Taylor microscale $\lambda$ for all the cases is
almost identical highlighting that background thermal gradients and length
scales of these gradients do not affect the shock waves significantly.

The integral length scale $\ell$ quantifies the large length scales at which
energy is injected in the system. Since we choose $k_{F}=5$, we use
$\ell=\frac{1}{5/2}$ in Fig. 5 for scaling. We identify the Taylor microscale
as Gupta and Scalo (2018)

$\lambda\sim\sqrt{\frac{\sum_{k}E_{k}}{\sum_{k}k^{2}E_{k}}}.$ (43)

Assuming the spectral energy $\widehat{E}^{f}_{wk}$ constant from
$0<k<2/\ell$, we obtain the relation,

$\lambda\sim\sqrt{\eta\ell},$ (44)

which indicates that the Taylor microscale is in the middle of the integral
length scale and the Kolmogorov length scale, similar to one-dimensional
acoustic wave turbulence and unlike classical hydrodynamic turbulence in which
$\lambda\sim\eta^{2/3}\ell^{1/3}$ and is biased towards the smallest length
scale indicating that the dissipation is higher at smaller length scales Pope
(2000). Hence, scaling in Eq. (44) also shows that the dissipation occurs
uniformly across all length scales smaller than the integral length scale,
particularly in the range of wavenumbers for which $\widehat{\Pi}_{wk}$ is
high. Figure 6 shows the wave energy spectra $\widehat{E}^{f}_{wk}$, flux of
spectral wave energy $\widehat{\Pi}_{wk}$, and cumulative dissipation of
spectral wave energy $\widehat{\mathcal{D}}_{wk}$ for simulation cases 0d-6d.
Background thermal gradients do not affect the wave energy spectra and the
length scales associated with the cascade of energy. Table 3 shows the time-
averaged values of the Taylor microscale $\langle\lambda\rangle_{t}$
calculated using Eq. (43) and the Kolmogorov length scale
$\langle\eta\rangle_{t}$ calculated using Eq. (38) for cases 0d-6d (see Table
1). Both the length scales change insignificantly as the length scales of
thermal gradients change. Moreover, for all cases $k_{max}\eta$ remains
greater than 1 confirming the shock-resolved nature of simulations. As we
discuss further, baroclinic interactions with background thermal gradients
result in change of vortical length scales generated in the domain.

Table 3: Time-averaged Taylor microscale $\langle\lambda\rangle_{t}$ and Kolmogorov scale $\langle\eta\rangle_{t}$.  Case | 0d | 1d | 2d | 3d | 4d | 5d | 6d
---|---|---|---|---|---|---|---
$\langle\lambda\rangle_{t}$ | 0.087 | 0.088 | 0.083 | 0.089 | 0.084 | 0.086 | 0.088
$\langle\eta\rangle_{t}$ | 0.0030 | 0.0027 | 0.0028 | 0.0028 | 0.0028 | 0.0030 | 0.0028

###  IV.2 Vorticity

Due to spectral energy cascade, pressure gradients exist at a wide range of
length scales (from largest length scale $\ell$ to the smallest length scale
$\eta$). These pressure gradients can interact with the density gradients in
the domain to generate local baroclinic rotation of fluid, quantified by
vorticity. Motivated by the decomposition in Eqs. (19c), we identify vorticity
$\bm{\Omega}$ and density-weighted vorticity $\bm{\Omega}_{f}$ as,

$\bm{\Omega}=\nabla\times\bm{u},~{}~{}\bm{\Omega}_{f}=\nabla\times\left(\sqrt{\rho}\bm{u}\right).$
(45)

Taking curl of Eq. (1b) (with constant $\mu$ and $\kappa=0$), we obtain the
governing equation for $\bm{\Omega}$ (in 2D) as,

$\displaystyle\frac{\partial\bm{\Omega}}{\partial t}$
$\displaystyle=\frac{\nabla\rho\times\nabla
p}{\rho^{2}}+\frac{\mu\nabla^{2}\bm{\Omega}}{\rho
Re_{ac}}+\frac{\mu}{Re_{ac}}\frac{\nabla\rho\times[\nabla\times\bm{\Omega}]}{\rho^{2}}$
$\displaystyle-\left(\frac{4\mu}{3Re_{ac}}\right)\frac{\nabla\rho\times\nabla(\nabla\cdot\bm{u})}{\rho^{2}}-\nabla\cdot(\bm{u}\bm{\Omega})+\bm{F}^{r}_{\Omega}.$
(46)

Figure 7: Vortical energy $\widehat{E}_{vk}$ and density-weighted vortical
energy $\widehat{E}^{f}_{vk}$ for homogeneous (baseline) simulation 0d (see
Table 2) and simulation with thermal gradients 6d (see Table 1) with
stochastically forced acoustic waves. While $\widehat{E}^{f}_{vk}$ peaks near
$k_{T}=60$ for 6d, $\widehat{E}_{vk}$ exhibits no such peak, thus highlighting
that the rotational energy $\widehat{E}_{vk}$ is broadband when shocks are
forced stochastically.

(a)(b)

Figure 8: Time series of (a) density-weighted wave energy $E^{f}_{w}$ and (b)
Vortical energy $E_{v}$ for simulations 1dr-6dr before and after energy
rescaling is enabled. Overall vortical energy in the system decreases and the
density-weighted wave energy quickly saturates near $\epsilon^{2}_{A}$.

Equation (46) holds without any perturbation decomposition of the field
variables. First term on the right hand side of Eq. (46) is the baroclinic
term. As we show in the next section, interaction of density gradients with
pressure gradients is the primary source of vorticity (and enstrophy)
generation in our simulations. These density gradients include both the
background thermal inhomogeneity and the density gradients generated in the
shock waves. The second term is the viscous diffusion of vorticity. The next
two terms represent the interaction of density gradients with viscous stresses
in the domain. In a field of shock waves, bulk viscous stress terms
(proportional to $\nabla\left(\nabla\cdot\bm{u}\right)$) are significantly
higher than shear stresses (proportional to $\nabla\times\Omega$) because
velocity dilatation is orders of magnitude higher than the rotation rate (c.f.
Figs. 1 and 10). Hence, we expect the third term in Eq. (46) to be
significantly smaller in magnitude compared to the fourth term. Furthermore,
an interaction of density gradients with bulk viscous stresses (captured by
the fourth term) does not necessarily denote dissipation of vorticity and
could also participate in generation of vorticity, as noted by Kida and Orszag
(1990). The next term represents advection of vorticity. The last term
$\bm{F}^{r}_{\Omega}$ denotes any forcing of vorticity due to density-weighted
wave energy rescaling. Since it is difficult to write the energy rescaling as
a right hand side in the momentum equation. Hence, we do not simplify
$\bm{F}^{r}_{\Omega}$ further. Since the curl of stochastic forcing $\bm{F}$
is set to zero, there is no stochastic forcing of vorticity $\bm{\Omega}$. In
the next section, we derive the enstrophy budget terms using Eq. (46).
Similarly, a governing equation for
$\bm{\Omega}_{f}=\nabla\times\left(\sqrt{\rho}\bm{u}\right)$ can be derived
as,

$\displaystyle\frac{\partial\bm{\Omega}_{f}}{\partial t}$
$\displaystyle=\frac{\nabla\sqrt{\rho}\times\nabla
p}{\rho}+\frac{\mu\nabla^{2}\bm{\Omega}}{\sqrt{\rho}Re_{ac}}+\frac{\mu}{Re_{ac}}\frac{\nabla\sqrt{\rho}\times[\nabla\times\bm{\Omega}]}{\rho}$
$\displaystyle-\left(\frac{4\mu}{3Re_{ac}}\right)\frac{\nabla\sqrt{\rho}\times\nabla(\nabla\cdot\bm{u})}{\rho}-\nabla\times(\sqrt{\rho}\bm{u}\cdot\nabla\bm{u})$
$\displaystyle+\nabla\sqrt{\rho}\times\bm{F}.$ (47)

Figure 9: Vortical energy $\widehat{E}_{vk}$ and density-weighted vortical
energy $\widehat{E}^{f}_{vk}$ for homogeneous (baseline) simulation 0d (see
Table 2) and simulation with thermal gradients 6d (see Table 1) with energy
rescaling.

Since energy rescaling only rescales energy in density-weighted wave modes,
density-weighted vorticity $\bm{\Omega}^{f}$ is not affected by energy
rescaling. However, the final term in Eq. (47) is due to random forcing
interacting with the density gradients. Figure 7 shows the spectra of
$\widehat{E}_{vk}$, which quantifies the vorticity $\bm{\Omega}$ in the domain
and $\widehat{E}^{f}_{vk}$, which quantifies the density-weighted vorticity
$\bm{\Omega}_{f}$ for a stochastically forced homogeneous simulation (0d in
Table 2) and a stochastically forced inhomogeneous simulation (6d in Table 1).
The random forcing results in generation of a broadband vorticity
$\bm{\Omega}$. The peak in $\widehat{E}^{f}_{vk}$ near the wavenumber $k_{T}$
occurs due to $\sqrt{f}$ in the definition of $\widehat{E}^{f}_{vk}$, hence
falsely indicating coherent vorticity in the simulation. Consequently, in our
further analysis of vorticity in this section and enstrophy budgets in the
next section, we consider $\widehat{E}_{vk}$, $\bm{\Omega}$, and the
corresponding velocity field $\bm{u}_{v}$ (c.f. Eq. (21)).

As shown in Fig. 7, vorticity is not localized at any particular length scale
in both the homogeneous and inhomogeneous simulations. The broadband vorticity
is generated due to changing shock curvature by the action of stochastic
forcing. This also indicates that if the curvature of the shock waves is
maintained, then the interaction with the thermal gradients must generate
coherent vorticity at length scale at which thermal gradients are
concentrated. To this end, we consider the energy rescaling runs which we
restart from the stochastically forced runs at $t=100$. During restart, we
remove the vortical velocity $\bm{u}_{v}$ (c.f. Eq. (21)) and put $\bm{F}=0$.
During the simulations, the density-weighted spectral energy
$\widehat{E}^{f}_{wk}$ is rescaled as per Eqs. (35)- (37) for
$\varepsilon_{A}=0.15$.

(a)(b)

Figure 10: (a) Instantaneous contours of $-\nabla\cdot\bm{w}^{\prime}$ showing
the two-dimensional field of shock waves generated by rescaled forcing and (b)
instantaneous contours of $\nabla\times\bm{u}^{\prime}$ showing the coherent
vorticity field generated due to the interaction shock waves with the
inhomogeneous background medium (for case 6d in Table 1) at different times.

Figure 8 shows the time series of total density-weighted wave energy
$E^{f}_{w}$ and vortical energy $E_{v}$ before and after energy rescaling
restart. The wave energy readjusts to the energy rescaling and saturates
around $\epsilon^{2}_{A}$ (see Fig. 8(a)). Since we remove any vortical
component $\bm{u}_{v}$ from the velocity field, the vortical energy drops but
immediately recovers to a small but finite value (approximately $1\times
10^{-5}$ in Fig. 8b) indicating that the vortical velocity field is atleast
one order smaller than the wave velocity field. Figure 9 shows the spectra of
$\widehat{E}_{vk}$ and $\widehat{E}^{f}_{vk}$ for simulations with energy
rescaling in a homogeneous medium (0dr in Table 2) and a heterogeneous medium
with thermal gradients at $k_{T}=60$ (6dr in Table 1). The vortical energy
$\widehat{E}^{f}_{vk}$ is narrow-band and is localized around the length scale
corresponding to $k_{T}$ for the inhomogeneous simulation, indicating coherent
vorticity generation. In simulations with energy rescaling, spectral wave
energy $\widehat{E}^{f}_{wk}$ is rescaled (restored) so as to maintain the
total wave energy at large length scales constant in time. Recoursing to
concepts in hydrodynamic turbulence, energy rescaling ensures _permanence of
large eddies_ in time. Consequently, the direction in which shock waves are
propagating (phasing of shock waves) is never altered by energy rescaling. As
a result, there is no broadband vorticity generated in the domain due to
changing directions and curvature of shock waves. Instead, a coherent
vorticity at the same length scale corresponding to the thermal gradient is
generated. Figure 10(a) shows the snapshots of two-dimensional field of shock
waves generated by the energy rescaling in an inhomogeneous medium with
thermal gradients at $k_{T}=60$ (6dr in Table 1). Figure 10(b) shows the
coherent vorticity generated due to these shock waves.

In the next section, we discuss the budgets of enstrophy defined as
$Q=|\bm{\Omega}|^{2}/2$ indetifying the primary mechanism for generation of
coherent vorticity due to shock waves driven by energy rescaling.

Figure 11: Enstrophy spectra $\widehat{Q}_{k}$ for simulations with
inhomogeneous thermal gradients 1dr \- 6dr (see Table 1) with energy
rescaling.

###  IV.3 Enstrophy Budget

Generation of coherent vorticity is highlighted in enstrophy spectra (see Fig.
11) defined as,

$\widehat{Q}_{k}=\frac{1}{2}\big{|}\widehat{\bm{\Omega}}_{k}\big{|}^{2}.$ (48)

Using Eq. (46), we obtain the governing equation for enstrophy spectra as,

$\displaystyle\frac{d\widehat{Q}_{k}}{dt}=\widehat{B}_{k}-\widehat{D}_{k}-\widehat{A}_{k}+\widehat{F}^{r}_{k},$
(49)

where

$\displaystyle\widehat{B}_{k}$
$\displaystyle=\mathrm{Re}\left[\widehat{\bm{\Omega}}^{*}_{k}\cdot\widehat{\left(\frac{\nabla\rho\times\nabla
p}{\rho^{2}}\right)}_{k}\right],$ (50) $\displaystyle\widehat{D}_{k}$
$\displaystyle=-\frac{\mu}{Re_{ac}}\mathrm{Re}\Bigg{[}\widehat{\bm{\Omega}}^{*}_{k}\cdot\Bigg{(}\frac{\nabla^{2}\bm{\Omega}}{\rho}+\frac{\nabla\rho\times\nabla\times\bm{\Omega}}{\rho^{2}}$
$\displaystyle-\frac{4}{3}\frac{\nabla\rho\times\nabla(\nabla\cdot\bm{u})}{\rho^{2}}\Bigg{)}_{k}\Bigg{]},~{}$
(51) $\displaystyle\widehat{A}_{k}$
$\displaystyle=\mathrm{Re}\left[\widehat{\bm{\Omega}}^{*}_{k}\cdot\left(\widehat{\nabla\cdot\left(\bm{u}\bm{\Omega}\right)}\right)_{k}\right],~{}$
(52)

and $\widehat{F}^{r}_{k}$ denotes the enstrophy generation due to energy
rescaling (term corresponding to $\bm{F}^{r}_{\Omega}$ in Eq. (46)).

(a)(b)

Figure 12: (a) Time series of total enstrophy $\sum_{k}\widehat{Q}_{k}$,
baroclinic term $\sum_{k}\widehat{B}_{k}$, nonlinear baroclinic term
$\sum_{k}\widehat{B}^{NL}_{k}$, dissipation $\sum_{k}\widehat{D}_{k}$, and the
advection $\sum_{k}\widehat{A}_{k}$ term for simulation 6dr. The baroclinic
production oscillates about a positive mean value (averaged in time). (b) Time
averaged total enstrophy $\langle Q\rangle_{t}$, baroclinic term $\langle
B\rangle_{t}$, nonlinear baroclinic term $\langle B^{NL}\rangle_{t}$ negative
of dissipation $\langle D\rangle_{t}$, and advection $\langle A\rangle_{t}$
terms. Magnitude of $\langle D\rangle_{t}$ is higher than $\langle
B\rangle_{t}$. Difference between the two is accounted for by the energy
injection due to energy rescaling $\widehat{F}^{r}_{k}$. Legends in both (a)
and (b) are same.

(a)(b)(c)

Figure 13: (a) Spectra of baroclinic term $\widehat{B}_{k}$, negative of
nonlinear baroclinic term $\widehat{B}^{NL}_{k}$, and the negative of (c)
dissipation term $-\widehat{D}_{k}$ in enstrophy Eq. (49).

First term on the right of Eq. (49), $\widehat{B}_{k}$, denotes the baroclinic
generation of enstrophy and captures the effect of the baroclinic torque on
enstrophy. The next three terms (combined in $\widehat{D}_{k}$) appear due to
the viscous forces in the momentum equation. Not all three terms corresponding
to the viscous forces necessarily destroy enstrophy. The first term of these
three terms which is proportional to the Laplacian of vorticity is always non-
positive and has primary contribution in viscous dissipation of enstrophy. The
second term and the third term of these three terms can be positive as well as
negative. The next term ($\widehat{A}_{k}$) arises due to convection term in
the momentum equation and has very small contribution in changing enstrophy of
the system. Lastly, $\widehat{F}^{r}_{k}$ is due to density-weighted energy
rescaling. As discussed in Section IV.2, a simplified form of energy rescaling
as a right hand side term in governing equations is difficult so we do not
simplify $\widehat{F}^{r}_{k}$.

Figure 11 shows the spectra of enstrophy for simulations with energy
rescaling. As noted for vortical energy spectra $\widehat{E}^{f}_{vk}$ in Fig.
9, the vorticity peaks at the same wavenumber at which background thermal
gradients exist ($k_{T}$). Consequently, the enstrophy spectra also peak at
$k_{T}$. Such localized (in the spectral space) enstrophy is generated due to
both the baroclinic production term $\widehat{B}_{k}$ as well as the energy
rescaling term $\widehat{F}^{r}_{k}$. Figure 12(a) shows the time series of
right hand side terms of Eq. (49) and total enstrophy for simulation 6dr. The
baroclinic term oscillates about a positive mean value while the advection
terms have negligible contribution to the change of enstrophy in time. Figure
12(b) shows the time averaged values of the enstrophy, the full baroclinic
production term and the nonlinear baroclinic production term, negative of the
dissipation, and the advection terms for simulations 1dr-6dr
($\langle\cdot\rangle_{t}$ denoting time average). As the wavenumber $k_{T}$
increases, the enstrophy in the system increases. Accordingly, the baroclinic
generation and the dissipation also increase. Enstrophy $\widehat{Q}_{k}$ is
statistically stationary since the vortical energy $E_{v}$ is statistically
stationary as shown in Fig. 8(b). Consequently, $\widehat{F}^{r}_{k}$ accounts
for the difference between the two since the L.H.S of Eq. (49) vanishes.
Density gradients in curved shock waves need not be parallel to the pressure
gradients due to thermoviscous generation of entropy within the shock waves.
Furthermore, due to imposed thermal inhomogeneity, shock waves get reflected
and refracted continuously in the domain. A combination of these two effects
can be seen by decomposing the baroclinic term $\widehat{B}_{k}$ as,

$\widehat{B}_{k}=\mathrm{Re}\left[\widehat{\bm{\Omega}}^{*}_{k}\cdot\widehat{\left(\frac{\nabla
f\times\nabla
p}{\rho^{2}}\right)}_{k}\right]+\underbrace{\mathrm{Re}\left[\widehat{\bm{\Omega}}^{*}_{k}\cdot\widehat{\left(\frac{\nabla(\rho-f)\times\nabla
p}{\rho^{2}}\right)}_{k}\right]}_{\widehat{B}^{NL}_{k}},$ (53)

where $\widehat{B}^{NL}_{k}$ captures the nonlinear baroclinic interactions
between pressure gradients and density gradients. Figure 12(a) shows the time
series of the nonlinear baroclinic term for simulation 6dr and Fig. 12(b),
time average of the nonlinear baroclinic interaction term is shown. The
overall contribution of the nonlinear baroclinic term to enstrophy is very
small. However, as shown in Fig. 13(b), the nonlinear baroclinic interaction
is local in the wavenumber space and negative in value, indicating its
contribution in destruction of enstrophy. Additionally, the nonlinear
baroclinic term peaks near $k_{T}$, indicating that any nonlinear baroclinic
interaction is amplified due to the background thermal gradients. Similarly,
the dissipation term $\widehat{D}_{k}$ is also highly negative and localized
near $k_{T}$. On the other hand, the total baroclinic term $\widehat{B}_{k}$
has a nonlocal contribution indicating the participation of a wide range of
length scales in generation of enstrophy. As shown in Fig. 12(a) and 12(b),
the overall contribution of $\widehat{B}_{k}$ is positive and large.
Furthermore, the difference between $\widehat{B}_{k}$ and $\widehat{D}_{k}$ in
Fig. 12(b) represents $\widehat{F}^{r}_{k}$ in Eq. (49). In Appendix B, we
discuss the relative magnitudes of dissipation terms in Eq. (51).

##  V Conclusions

We have studied forced two-dimensional shock waves in a periodic domain with
an inhomogeneous background temperature field using numerical simulations. We
showed that random forcing generates an isotropic field of shock waves and the
density-weighted wave energy spectra scale with wavenumber as $k^{-2}$ as in
Eq. (42). The flux of energy spectra does not remain constant in the cascade
range of wavenumbers if physical diffusion is considered indicating the
spectral energy transfer to be of variable flux type Verma (2021). We have
also shown that dissipation exists in all length scales, unlike hydrodynamic
turbulence, where dissipation is localized in smaller length scales. Such
energy spectra cascade characteristics remain similar when the field of shock
waves propagates in an inhomogeneous medium with randomized thermal gradients.

While the interaction of shock waves with background thermal gradients results
in a baroclinic torque generating vorticity, the type of energy injection in
shock waves governs the evolution of this vorticity field in the domain. We
showed that for stochastically forced shock waves, no coherent vorticity is
generated, and only a broadband vorticity (distributed over all length scales)
exists. For stochastic forcing, the pressure gradient field in the domain is
also a stochastic process. Since the direction of the pressure gradient is
also stochastic, time-averaged baroclinic torque vanishes since the leading
order term is linearly dependent on the pressure gradient. Consequently, only
broadband vorticity is observed due to changing direction of pressure
gradients. The broadband vorticity is high in the forced length scales and
decays for smaller length scales, further highlighting that it is generated
due to the stochastic forcing. To investigate the generation of coherent
vorticity, we have introduced a new forcing based on wave energy rescaling. In
this forcing, we rescale the energy in the density-weighted wave modes at each
time step. The cascade is thus driven by waves that do not change phase and
propagation direction continuously. Consequently, coherent vorticity is
generated in the domain at the same length scale as that of the thermal
gradients. We have also shown that density-weighted variables introduced by
Miura and Kida (1995) are not suitable for studying rotational and
dilatational flow fields separately in compressible flows. Though in the
current study we have investigated the baroclinic interaction of forced weak
shock waves, we expect the baroclinic vorticity production to increase for
stronger shock waves due to higher pressure gradients. However, such
simulations are beyond the scope of the current study. For all the simulations
we used constant viscosity and thermal conductivity, however, as shown in
Appendix A, temperature dependent viscosity and conductivity have negligible
affects on the energy spectra of shock waves. Consequently, the baroclinic
interaction is not expected to be affected by the temperature dependent
viscosity and conductivity values in the weak-shock regime.

In practice, the generation of random weak shock waves may be due to intense
unsteady combustion events or thermoacoustic instabilities. Spatially
distributed unsteady combustion events are better represented by stochastic
forcing due to randomly generated pressure waves Dowling and Mahmoudi (2015),
while instability-driven shock waves are better represented by forcing based
on energy rescaling since only particular modes exhibit such instabilities
Gupta, Lodato, and Scalo (2017) depending on the geometry of the combustion
chamber. Besides the practical relevance in understanding unsteady shock
induced mixing, we also highlighted that the method of injection of energy in
the so-called _wave-turbulent_ Nazarenko (2011) systems in which nonlinear
waves interact with each other is a richer problem compared to the
hydrodynamic turbulent systems. In hydrodynamic turbulence, there is no wave
velocity of propagation of eddies due to which phasing is never considered in
hydrodynamic turbulence forcing Rosales and Meneveau (2005). However, phasing
or phase difference between pressure and velocity in both space and time is an
important parameter governing the acoustic wave turbulence and interaction
with background medium. Accordingly, future research directions include a
detailed investigation of random shock wave forcing and the possibility of
sustaining hydrodynamic turbulence via forced weak shock waves in three
dimensions.

###### Acknowledgements.

We acknowledge the financial support received from Science and Engineering
Research Board (SERB), Government of India under Grant No. SRG/2022/000728. We
also thank IIT Delhi HPC facility for computational resources.

##  Appendix A Effects of hyperdiffusion and variable viscosity

(a)(b)(c)

Figure 14: (a) Dimensionless scaled density-weighted wave energy spectra
$\widehat{E}^{f}_{wk}$, (b) spectral flux of wave energy $\widehat{\Pi}_{wk}$,
and (c) cumulative dissipation $\widehat{\mathcal{D}}_{wk}$ for a baseline
case (0d in Table 2), compared against a simulation with viscosity varying
with the Sutherland law in Eq. (54) (0s) and a simulation with hyperdiffusion
(0h). For simulation 0s, all parameters are kept same as 0d simulation except
for viscosity and thermal conductivity. Prandtl number $\mathrm{Pr}$ is fixed
at 0.72. For simulation 0h, hyperdiffusion is used with $r=2$ and
$\nu=5.44\times 10^{-10}$ (c.f. Eqs. (55) and (56)).

Throughout our simulations, we consider the dimensionless viscosity constant
($\mu=1$). However, for high thermal perturbations, changes in viscosity due
to inhomogeneous temperature may be significant. To confirm that such thermal
dependence of viscosity has an insignificant impact on the spectral wave
energy cascade in our simulations, we also perform a simulation with
dimensionless viscosity governed by the Sutherland’s law

$\mu=T^{3/2}\left(\frac{T_{m}+S_{ref}}{TT_{m}+S_{ref}}\right),$ (54)

where $T_{m}=273.15~{}\mathrm{K}$ and is the reference temeperature scale (see
Eq. (4)), $S_{ref}=110.4~{}\mathrm{K}$, and $T$ is the dimensionless
temperature. To obtain a wide range of length scales in a given resolution,
hyperdiffusion is used in wave-turbulence studies Augier, Mohanan, and
Lindborg (2019); Farge and Sadourny (1989). Such hyperdiffusion modifies all
the viscous stress tensor and thermal conduction terms in the momentum and
pressure governing Eqs. (1b) and (1c) as,

$\displaystyle\frac{\partial\bm{u}}{\partial
t}+\bm{u}\cdot\nabla\bm{u}=-\frac{\nabla p}{\rho}+\nu\nabla^{2r}\bm{u},$ (55)
$\displaystyle\frac{\partial p}{\partial t}+\bm{u}\cdot\nabla p+\gamma
p\nabla\cdot\bm{u}=\nu\nabla^{2r}p,$ (56)

where $r$ and $\nu$ are parameters of hyperdiffusivity. To assert the validity
of constant dynamic viscosity, we compare the energy spectra of simulation 0d
in Table 2, with two simulations, one with viscosity varying with local
temperature as per Eq. (54) (simulation 0s) and one with hyperdiffusion with
$r=2$ and $\nu=5.44\times 10^{-10}$ (simulation 0h) in Fig. 14.

As shown in Fig. 14 temperature dependent viscosity does not affect the energy
spectra significantly. For simulation with hyperdiffusion (0h), the energy
spectra extend with a $-2$ slope for larger range of wavenumbers. However, the
flux of spectral energy $\widehat{\Pi}$ and cumulative spectral dissipation
$\widehat{\mathcal{D}}_{k}$ are modified significantly. The flux is
approximately constant in the spectral transfer range and the dissipation is
accumulated in higher wavenumbers. However, due to physical viscous terms, the
spectra flux decays smoothly and dissipation exists in all wavenumbers.
Furthermore, even in a hypothetical simulations with infinite resolution,
since the energy spectra will decay as $k^{-2}$, the dissipation will exist in
all scales (since $\widehat{\mathcal{D}}_{k}\sim k^{2}\widehat{E}_{k}$).
Hence, replacing physical dissipation terms with a hyperdiffusion term tends
to modify the fundamental nature of spectral energy transfer. Since the
variation of viscosity has negligible effect on energy spectra and
hyperdiffusion modifies the fundamental nature of spectral energy transfer in
nonlinear acoustic wave turbulence, we choose dimensionless viscosity $\mu=1$
in all our simulations.

##  Appendix B Enstrophy dissipation

(a)(b)(c)

Figure 15: (a) Spectra of dissipation due to density variation interacting
with velocity divergence, (b) curl of vorticity, and (c) Laplacian of
vorticity as derived in Eq. (51).

Three terms contribute to enstrophy dissipation as derived in Eq. (51). In
Fig. 15 we show the spectra of all three terms contributing to the dissipation
of enstrophy in simulations 1dr-6dr (see Table 1). Similar to vorticity and
enstrophy, all the dissipation terms are localized around $k_{T}$ signifying
coherent vorticity in the domain. Figure 15(a) shows spectra of dissipation
due to misalignment between density gradient and bulk viscous forces
(quantified by $\nabla\left(\nabla\cdot\bm{u}\right)$. Figure 15b shows the
spectra of dissipation due to misalignment between density gradient and a
combination of bulk and shear viscous forces (since
$\nabla\times\Omega=\nabla\left(\nabla\cdot\bm{u}\right)-\nabla^{2}\bm{u}$).
Together, the terms shown in Figs. 15(a) and 15(b) maybe considered as viscous
counterparts of baroclinic torque in which misalignment of density gradients
with viscous forces contributes to dissipation. We highlight that these two
terms are not necessarily dissipative in nature Kida and Orszag (1990). Figure
15(c) shows the spectra of the Laplacian dissipation of vorticity and is
necessarily negative. The Laplacian dissipation is highest in rescaled
simulations denoting highly localized (in spectral space) vorticity. Since the
vorticity and production of enstrophy are highly sensitive to the method of
forcing random shock waves (stochastic forcing compared with rescaled forcing
in this work), we expect the enstrophy dissipation to also change with the
forcing methodology.

## References

  * Goldstein and Wundrow (1990) M. Goldstein and D. Wundrow, “Spatial evolution of nonlinear acoustic mode instabilities on hypersonic boundary layers,” Journal of Fluid Mechanics 219, 585–607 (1990).
  * Liepmann and Roshko (2001) H. W. Liepmann and A. Roshko,  _Elements of gasdynamics_ (Courier Corporation, 2001).
  * Baars and Tinney (2013) W. Baars and C. Tinney, “Transient wall pressures in an overexpanded and large area ratio nozzle,” Experiments in fluids 54, 1–17 (2013).
  * Bonciolini, Boujo, and Noiray (2017) G. Bonciolini, E. Boujo, and N. Noiray, “Output-only parameter identification of a colored-noise-driven van-der-pol oscillator: thermoacoustic instabilities as an example,” Physical Review E 95, 062217 (2017).
  * Gupta, Lodato, and Scalo (2017) P. Gupta, G. Lodato, and C. Scalo, “Spectral energy cascade in thermoacoustic shock waves,” Journal of Fluid Mechanics 831, 358–393 (2017).
  * Breazeale and Thompson (1963) M. Breazeale and D. Thompson, “Finite-amplitude ultrasonic waves in aluminum,” Applied Physics Letters 3, 77–78 (1963).
  * Lissenden (2021) C. J. Lissenden, “Nonlinear ultrasonic guided waves—principles for nondestructive evaluation,” Journal of Applied Physics 129, 021101 (2021).
  * Sagdeev (1979) R. Z. Sagdeev, “The 1976 oppenheimer lectures: Critical problems in plasma astrophysics. i. turbulence and nonlinear waves,” Reviews of Modern Physics 51, 1 (1979).
  * Brouillette (2002) M. Brouillette, “The richtmyer-meshkov instability,” Annual Review of Fluid Mechanics 34, 445–468 (2002).
  * Lukes  _et al._ (2016) P. Lukes, F. Fernández, J. Gutiérrez-Aceves, E. Fernández, U. Alvarez, P. Sunka, and A. Loske, “Tandem shock waves in medicine and biology: a review of potential applications and successes,” Shock waves 26, 1–23 (2016).
  * Lighthill (1978) J. Lighthill, “Acoustic streaming,” Journal of sound and vibration 61, 391–418 (1978).
  * Shimizu and Sugimoto (2010) D. Shimizu and N. Sugimoto, “Numerical study of thermoacoustic taconis oscillations,” Journal of Applied Physics 107, 034910 (2010).
  * Whitham (2011) G. B. Whitham,  _Linear and nonlinear waves_ (John Wiley & Sons, 2011).
  * Gupta and Scalo (2018) P. Gupta and C. Scalo, “Spectral energy cascade and decay in nonlinear acoustic waves,” Physical Review E 98, 033117 (2018).
  * Thirani, Gupta, and Scalo (2020) S. Thirani, P. Gupta, and C. Scalo, “Knudsen number effects on the nonlinear acoustic spectral energy cascade,” Physical Review E 101, 023101 (2020).
  * Ellermeier (1993) W. Ellermeier, “Nonlinear acoustics in non-uniform infinite and finite layers,” Journal of Fluid Mechanics 257, 183–200 (1993).
  * Tyagi and Sujith (2003) M. Tyagi and R. I. Sujith, “Nonlinear distortion of travelling waves in variable-area ducts with entropy gradients,” Journal of Fluid Mechanics 492, 1–22 (2003).
  * Prasad (2006) D. Prasad, “Weakly nonlinear shock propagation in slowly varying one-dimensional flows,” Physics of Fluids 18, 036101 (2006).
  * Budzinsky, Zukoski, and Marble (1992) J. Budzinsky, E. Zukoski, and F. Marble, “Rayleigh scattering measurements of shock enhanced mixing,” in  _28th Joint Propulsion Conference and Exhibit_ (1992) p. 3546.
  * Andreopoulos, Agui, and Briassulis (2000) Y. Andreopoulos, J. H. Agui, and G. Briassulis, “Shock wave—turbulence interactions,” Annual Review of Fluid Mechanics 32, 309–345 (2000).
  * Romagnosi  _et al._ (2011) L. Romagnosi, A. Ingenito, D. Cecere, G. Eugenio, and C. Bruno, “The role of the baroclinic term in supersonic fuel/air mixing enhancement,” in  _49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition_ (2011) p. 401.
  * Tian  _et al._ (2017) Y. Tian, F. A. Jaberi, Z. Li, and D. Livescu, “Numerical study of variable density turbulence interaction with a normal shock wave,” Journal of Fluid Mechanics 829, 551–588 (2017).
  * Wong  _et al._ (2022) M. L. Wong, J. R. Baltzer, D. Livescu, and S. K. Lele, “Analysis of second moments and their budgets for richtmyer-meshkov instability and variable-density turbulence induced by reshock,” Physical Review Fluids 7, 044602 (2022).
  * Singh  _et al._ (2019) B. Singh, L. K. Rajendran, P. Gupta, C. Scalo, P. P. Vlachos, and S. P. Bane, “Experimental and numerical study of flow induced by nanosecond repetitively pulsed discharges,” in  _AIAA Scitech 2019 Forum_ (2019) p. 0740.
  * Kundu, Cohen, and Dowling (2015) P. K. Kundu, I. M. Cohen, and D. R. Dowling,  _Fluid mechanics_ (Academic press, 2015).
  * Mikaelian (1985) K. O. Mikaelian, “Richtmyer-meshkov instabilities in stratified fluids,” Physical Review A 31, 410 (1985).
  * Yang and Radulescu (2021) H. Yang and M. I. Radulescu, “Dynamics of cellular flame deformation after a head-on interaction with a shock wave: reactive richtmyer–meshkov instability,” Journal of Fluid Mechanics 923, A36 (2021).
  * Yu  _et al._ (2020) B. Yu, M. He, B. Zhang, and H. Liu, “Two-stage growth mode for lift-off mechanism in oblique shock-wave/jet interaction,” Physics of Fluids 32, 116105 (2020).
  * Wei  _et al._ (2022) F. Wei, R. Yang, W. Liu, Y. Zhao, Q. Wang, and M. Sun, “Flow structures of strong interaction between an oblique shock wave and a supersonic streamwise vortex,” Physics of Fluids 34, 106102 (2022).
  * Eswaran and Pope (1988) V. Eswaran and S. B. Pope, “An examination of forcing in direct numerical simulations of turbulence,” Computers & Fluids 16, 257–278 (1988).
  * Jagannathan and Donzis (2016) S. Jagannathan and D. A. Donzis, “Reynolds and mach number scaling in solenoidally-forced compressible turbulence using high-resolution direct numerical simulations,” Journal of Fluid Mechanics 789, 669–707 (2016).
  * Miura and Kida (1995) H. Miura and S. Kida, “Acoustic energy exchange in compressible turbulence,” Physics of Fluids 7, 1732–1742 (1995).
  * Boyd (2001) J. P. Boyd,  _Chebyshev and Fourier spectral methods_ (Courier Corporation, 2001).
  * Mortensen and Langtangen (2016) M. Mortensen and H. P. Langtangen, “High performance python for direct numerical simulations of turbulent flows,” Computer Physics Communications 203, 53–65 (2016).
  * Yeung, Sreenivasan, and Pope (2018) P. Yeung, K. Sreenivasan, and S. Pope, “Effects of finite spatial and temporal resolution in direct numerical simulations of incompressible isotropic turbulence,” Physical Review Fluids 3, 064603 (2018).
  * Bell  _et al._ (2022) J. B. Bell, A. Nonaka, A. L. Garcia, and G. Eyink, “Thermal fluctuations in the dissipation range of homogeneous isotropic turbulence,” Journal of Fluid Mechanics 939, A12 (2022).
  * J Evans and P Morriss (2007) D. J Evans and G. P Morriss,  _Statistical mechanics of nonequilbrium liquids_ (ANU Press, 2007).
  * Ghosal  _et al._ (1995) S. Ghosal, T. S. Lund, P. Moin, and K. Akselvoll, “A dynamic localization model for large-eddy simulation of turbulent flows,” Journal of fluid mechanics 286, 229–255 (1995).
  * Augier, Mohanan, and Lindborg (2019) P. Augier, A. V. Mohanan, and E. Lindborg, “Shallow water wave turbulence,” Journal of Fluid Mechanics 874, 1169–1196 (2019).
  * Verma (2019) M. K. Verma,  _Energy transfers in fluid flows: multiscale and spectral perspectives_ (Cambridge University Press, 2019).
  * Verma (2021) M. K. Verma, “Variable energy flux in turbulence,” Journal of Physics A: Mathematical and Theoretical 55, 013002 (2021).
  * Pope (2000) S. B. Pope,  _Turbulent flows_ (Cambridge university press, 2000).
  * Kida and Orszag (1990) S. Kida and S. A. Orszag, “Enstrophy budget in decaying compressible turbulence,” Journal of scientific computing 5, 1–34 (1990).
  * Dowling and Mahmoudi (2015) A. P. Dowling and Y. Mahmoudi, “Combustion noise,” Proceedings of the Combustion Institute 35, 65–100 (2015).
  * Nazarenko (2011) S. Nazarenko,  _Wave turbulence_ , Vol. 825 (Springer Science & Business Media, 2011).
  * Rosales and Meneveau (2005) C. Rosales and C. Meneveau, “Linear forcing in numerical simulations of isotropic turbulence: Physical space implementations and convergence properties,” Physics of fluids 17, 095106 (2005).
  * Farge and Sadourny (1989) M. Farge and R. Sadourny, “Wave-vortex dynamics in rotating shallow water,” Journal of Fluid Mechanics 206, 433–462 (1989).