1. Introduction
A complex turbulent flow resulting from the interaction of jets has significant engineering applications, including cooling systems (Geers et al. Reference Geers, Hanjalić and Tummers2006; Caliskan et al. Reference Caliskan, Baskaya and Calisir2014), combustors (Burr & Ken Reference Burr and Yu2019; Yang et al. Reference Yang, Wu, Zeng, Wang and Zhao2021) and air-conditioning systems (Medaouar et al. Reference Medaouar, Loukarfi, Braikia, Khelil and Naji2019). In environmental flows, similar interactions are observed in plumes from stacked chimneys (Bornoff & Mokhtarzadeh-Dehghan Reference Bornoff and Mokhtarzadeh-Dehghan2001) and exhausts from air-cooled condenser arrays (Liu et al. Reference Liu, Duan and Zhao2009). High-speed jet interactions are also relevant to supersonic aircraft, where jet-induced noise is a critical factor in design and development (Raman & Taghavi Reference Raman and Taghavi1996). These practical challenges have motivated extensive research on jet interactions in various configurations, as summarised in the literature review by Boussoufi et al. (Reference Boussoufi, Sabeur-Bendehina, Ouadha, Morsli and El Ganaoui2017).
The simplest jet interaction occurs in confluent jets issued from two identical nozzles. This flow configuration, involving both planar and circular jets, has been extensively studied (Miller & Comings Reference Miller and Comings1960; Lin & Sheu Reference Lin and Sheu1990; Anderson & Spall Reference Anderson and Spall2001; Manohar et al. Reference Manohar, Sundararajan, Ramjee and Kumar2004; Meslem et al. Reference Meslem, Nastase and Allard2010; Hao et al. Reference Hao, Tian and Zhou2021). These studies have identified distinct regimes in jet interaction, including converging, merging and combined regions, highlighting the flow evolution during the interaction. Initially, the separated jets converge towards the centreline between the two jets. For planar jets, a recirculation zone forms between the jets in the converging region. Strong jet interaction occurs in the merging region, where the centreline velocity between the jets increases. Further downstream lies the combined region, where the centreline mean velocity decays similarly to a single jet. Multiple confluent jets issued from nozzle arrays into free space have also been studied, though less extensively than dual jets (Yimer et al. Reference Yimer, Becker and Grandmaison1996; Meslem et al. Reference Meslem, Nastase and Allard2010; Ghahremanian et al. Reference Ghahremanian, Svensson, Tummers and Moshfegh2014; Svensson et al. Reference Svensson, Rohdin and Moshfegh2016). These studies reveal entrainment of surrounding fluid into the freely evolving multiple jets, which contributes to the decay of the mean streamwise velocity. The nozzle shape and configuration influence flow development in the near field but have minimal impact on the far field, where the combined jet exhibits characteristics similar to a single jet. Research on supersonic jet interactions has primarily focused on jet noise (Coltrin et al. Reference Coltrin, Blotter, Maynes and Gee2013, Reference Coltrin, Maynes, Blotter and Gee2014). Additionally, numerous studies have investigated impinging multiple jets due to their applications in heating and cooling devices (Thielen et al. Reference Thielen, Jonker and Hanjalić2003; Geers et al. Reference Geers, Tummers and Hanjalić2005).
Confined multiple confluent jets have also been studied extensively (Teunissen Reference Teunissen1975; Tatsumi et al. Reference Tatsumi, Tanaka, Woodfield and Nakabe2010; Tan et al. Reference Tan, Xu, Qi and Ni2023; Mori et al. Reference Mori, Watanabe and Nagata2024). When a large number of nozzles are arranged in a square array, jet interaction leads to the formation of decaying homogeneous isotropic turbulence (HIT). A similar formation of decaying HIT has been reported for unconfined square jet arrays in the early stage of the combined region, where the flow near the centreline becomes nearly homogeneous in the cross-streamwise plane (Ghahremanian et al. Reference Ghahremanian, Svensson, Tummers and Moshfegh2014). As detailed below, the behaviour of decaying HIT strongly depends on the turbulence generation mechanism. Accordingly, some previous studies have conducted experiments with multiple confluent jets to examine HIT (Tan et al. Reference Tan, Xu, Qi and Ni2023; Mori et al. Reference Mori, Watanabe and Nagata2024). The aforementioned studies primarily focus on continuous jets. However, the formation of HIT has also been observed in studies of unsteady jet interactions, such as pulse-operated jets (Bellani & Variano Reference Bellani and Variano2014; Carter et al. Reference Carter, Petersen, Amili and Coletti2016; Pérez-Alvarado et al. Reference Pérez-Alvarado, Mydlarski and Gaskin2016) and synthetic jets (Hwang & Eaton Reference Hwang and Eaton2004; Variano et al. Reference Variano, Bodenschatz and Cowen2004; Yamamoto et al. Reference Yamamoto, Watanabe and Nagata2022b ). Some of these studies explore the collision of multiple jets from several arrays. The turbulence generation mechanisms differ between continuous and unsteady jets: the former is associated with turbulent kinetic energy production due to mean shear, while the latter is influenced by temporal variations in jet velocity.
HIT is one of the most fundamental issues in turbulence research (Davidson Reference Davidson2004). The assumption of statistical homogeneity and isotropy simplifies the problem, facilitating the development of statistical theories. Most theories have been formulated for the statistics of small-scale motions, which are believed to exhibit some universal properties. Fewer theories target large scales, despite their dominant role in various turbulent phenomena, such as turbulent mixing. Notable exceptions include the theories of Saffman and Batchelor turbulence, which describe HIT with three-dimensional energy spectra
$E(k)= L k^2/4\pi ^2$
and
$I k^4/24\pi ^2$
, respectively, with a small wavenumber
$k$
(Batchelor & Proudman Reference Batchelor and Proudman1956; Saffman Reference Saffman1967; Davidson Reference Davidson2004). Saffman turbulence is characterised by linear momentum, whose conservation results in the Saffman integral
$L$
being constant. In scenarios where
$L=0$
, the Loitsyansky integral
$I$
related to angular momentum remains approximately constant when long-range correlations are negligible, leading to the Batchelor theory. These invariants constrain root-mean-squared (r.m.s.) velocity fluctuations
$U$
and integral length scale
$L$
during decay, where
$U^2 L^3=\mathrm {const.}$
for Saffman turbulence and
$U^2 L^5=\mathrm {const.}$
for Batchelor turbulence, provided that the large scales evolve in a self-similar manner. The variation of the non-dimensional dissipation rate
$C_\varepsilon =\varepsilon /(U^3/L)$
in decaying HIT is often described by
$C_\varepsilon \sim (t-t_0)^{n_C}$
, where
$\varepsilon$
is the dissipation rate of turbulent kinetic energy (TKE) per unit mass and
$t_0$
is a virtual origin of decay (Krogstad & Davidson Reference Krogstad and Davidson2010). A constant
$C_\varepsilon$
at high Reynolds numbers implies
$n_C=0$
. Non-equilibrium dissipation scaling or low-Reynolds-number effects result in
$n_C\gt 0$
(Meldi & Sagaut Reference Meldi and Sagaut2013; Kitamura et al. Reference Kitamura, Nagata, Sakai, Sasoh, Terashima, Saito and Harasaki2014; Valente & Vassilicos Reference Valente and Vassilicos2014). The evolution of
$U^2$
and
$L$
can be expressed as
\begin{align} U^2 &\sim (t-t_0)^{-10(1+n_C)/7}, \,\,\, L \sim (t-t_0)^{2(1+n_C)/7} \,\,\,\text {for Batchelor turbulence}, \end{align}
\begin{align} U^2 &\sim (t-t_0)^{-6(1+n_C)/5}, \,\,\, L \sim (t-t_0)^{2(1+n_C)/5}\,\,\,\text {for Saffman turbulence}, \end{align}
where (1.1) is referred to as the Kolmogorov decay law (Kolmogorov Reference Kolmogorov1941; Krogstad & Davidson Reference Krogstad and Davidson2010).
Experiments of HIT have been conducted to validate theories and models of turbulence, including Saffman and Batchelor theories. Grid turbulence is often studied as an approximation to HIT, employing various grids designed to generate HIT with desired properties. These grids vary in their geometries, including shapes and configurations of grid bars. Some of grid turbulence experiments demonstrate that the decay law follows that of Saffman turbulence, for various types of grid (Praud et al. Reference Praud, Fincham and Sommeria2005; Krogstad & Davidson Reference Krogstad and Davidson2010, Reference Krogstad and Davidson2011; Kitamura et al. Reference Kitamura, Nagata, Sakai, Sasoh, Terashima, Saito and Harasaki2014; Sinhuber et al. Reference Sinhuber, Bodenschatz and Bewley2015; Watanabe et al. Reference Watanabe, Zheng and Nagata2022). However, the observation of the Kolmogorov decay law has been limited to an early study of grid turbulence at extremely low Reynolds numbers (Batchelor & Townsend Reference Batchelor and Townsend1948). This evidence is also considered weak because it relies solely on the decay of
$U$
, rather than on the invariance of
$U^2L^5$
or the variations of both
$U$
and
$L$
. Although the theories discussed above predict the decay laws for turbulence with
$ E(k) \sim k^2$
or
$ E(k) \sim k^4$
, they do not explain under what conditions turbulence with such spectra are generated, such as how to produce turbulence with or without net linear momentum. Whether the decay of HIT follows the laws of Saffman turbulence, Batchelor turbulence or neither is expected to depend on the method of turbulence generation. Therefore, it is essential to examine decaying HIT generated by various methods and compare the results with the theoretical decay laws. It is important to note that numerical studies on Saffman and Batchelor turbulence typically examine the decay of isotropic turbulence, where the spectrum is artificially initialised to follow
$ E(k) \sim k^2$
or
$ E(k) \sim k^4$
. The formation of Saffman or Batchelor turbulence relates to the transition from laminar to turbulent flow, such as turbulence generated by grids. Such studies are primarily limited to grid turbulence, as discussed earlier. While multiple jet interactions are known to generate HIT, their relevance to the decay laws of Saffman and Batchelor turbulence has been scarcely explored.
The present study investigates the decay of nearly HIT generated by multiple jet interactions. Direct numerical simulation (DNS) is performed for temporally evolving circular jets. Temporal shear flows are commonly considered in fundamental turbulence studies as approximations of spatial flows, such as jets (da Silva & Pereira Reference da Silva and Pereira2008; Sadeghi et al. Reference Sadeghi, Oberlack and Gauding2018; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021a
), wakes (Diamessis et al. Reference Diamessis, Spedding and Domaradzki2011; Watanabe et al. Reference Watanabe, Riley, de Bruyn Kops, Diamessis and Zhou2016a
; Zecchetto & da Silva Reference Zecchetto and da Silva2021), mixing layers (Gampert et al. Reference Gampert, Boschung, Hennig, Gauding and Peters2014; Watanabe et al. Reference Watanabe, da Silva, Sakai, Nagata and Hayase2016b
; Blakeley et al. Reference Blakeley, Olson and Riley2022) and boundary layers (Kozul et al. Reference Kozul, Chung and Monty2016; Watanabe et al. Reference Watanabe and Nagata2018; Zhang et al. Reference Zhang, Watanabe and Nagata2018). Temporal simulations employ periodic boundary conditions in the streamwise direction, allowing the flow to evolve over time rather than spatially. Previous studies have demonstrated that the turbulence generation process due to mean shear is similar for both temporal and spatial flows. Using this approach, DNS of grid turbulence has been conducted as a temporal wake of a grid, showing good agreement with spatially evolving grid turbulence observed in experiments (Watanabe & Nagata Reference Watanabe and Nagata2018; Watanabe et al. Reference Watanabe, Zheng and Nagata2022). The present study adopts the same methodology for jet interactions. Temporal simulations are computationally less expensive than spatial simulations, enabling the use of a very large computational domain where turbulence decay is unaffected by boundary conditions. Additionally, periodic boundary conditions in the streamwise direction allow direct evaluation of the three-dimensional energy spectrum, facilitating comparison with
$ E(k) \sim k^2$
and
$ E(k) \sim k^4$
. Consequently, this study focuses on temporally evolving jets. DNS is conducted for both subsonic and supersonic jets. This DNS demonstrates that turbulence arising from jet interactions behaves similarly to Batchelor turbulence. In the supersonic case, fluid compressibility potentially affects turbulence decay, with such effects prominent in both transitional and decay regimes. A supersonic jet generates shock waves during the initial transitional regime. These waves propagate into the far field, influencing small-scale flow features, such as velocity gradient statistics (Yamamoto et al. Reference Yamamoto, Ishida, Watanabe and Nagata2022a
). The present DNS also aims to explore how the modulation of small scales by these waves affects the large-scale flow features described in the theories of Batchelor and Saffman turbulence.
Previous DNS studies on Batchelor and Saffman turbulence have primarily focused on validating the respective theories (Ishida et al. Reference Ishida, Davidson and Kaneda2006; Anas et al. Reference Anas, Joshi and Verma2020; Davidson et al. Reference Davidson, Okamoto and Kaneda2012). These studies typically employ DNS with initial conditions specifically designed so that the initial energy spectrum
$ E(k)$
follows
$ k^4$
or
$ k^2$
. The decay of turbulence with
$ E(k) \sim k^4$
or
$ k^2$
is then compared against the theoretical predictions. Such investigations require the turbulence to be incompressible, homogeneous and strictly isotropic, as assumed in the theories. In contrast, the present study aims to examine the decay properties of turbulence generated by jet interaction. The present DNS does not impose a spectral shape of
$ E(k) \sim k^4$
or
$ k^2$
in the initial condition. Instead, these spectral shapes may emerge during the turbulent transition from the jet interaction, as the initial condition is laminar. In this sense, the present DNS differs from earlier studies of Batchelor and Saffman turbulence. It provides insights into how Batchelor or Saffman turbulence can develop from a laminar state, rather than focusing solely on whether the theoretical decay laws of Batchelor and Saffman turbulence are recovered in turbulence with pre-imposed spectral shapes of
$ E(k) \sim k^4$
or
$ k^2$
.
The paper is organised as follows. Section 2 describes the numerical details of temporally evolving jets. Section 3 presents the DNS results, including the formation of nearly HIT from the jet interaction, compressiblity effects on small-scale characteristics of turbulence and decay behaviour. Finaly, the findings from the DNS are summerised in § 4.
Figure 1. A schematic of the initial field for the DNS of temporally evolving jet interactions. The initial jet regions are depicted in grey. The number of jets shown is illustrative and does not match the actual value of 40. The figure also displays the coordinate system
$(Y, Z)$
, which is used for statistical evaluations.
2. Direct numerical simulations for the interaction of temporally evolving jets
2.1. Temporally evolving jets
DNS is used to study temporally evolving multiple jets, applying methodologies of temporal jets and grid turbulence (da Silva & Pereira Reference da Silva and Pereira2008; Sadeghi et al. Reference Sadeghi, Oberlack and Gauding2018; Watanabe & Nagata Reference Watanabe and Nagata2018; Watanabe et al. Reference Watanabe, Zheng and Nagata2022). Figure 1 illustrates a schematic of the initial field used in the present DNS. The number of jets shown is illustrative and does not correspond to the actual number in the simulations. Temporal jets, with an initial diameter of
$ D$
, are positioned at an equal spacing of
$ S$
in the
$ y$
- and
$ z$
-directions. The initial velocity profile of each jet matches that of temporal round jets studied previously (Pineau & Bogey Reference Pineau and Bogey2020; Takahashi et al. Reference Takahashi, Fukui, Tsujimoto, Ando and Shakouchi2023), with its mathematical expression provided below. The temporal jets evolve over time within a triply periodic domain, becoming turbulent due to instabilities associated with their mean velocity profiles. Interactions among the jets result in the formation of nearly HIT. As all jets are statistically equivalent, statistics are evaluated by spatial averages in the streamwise direction and ensemble averages across all jets, expressed as functions of the local jet coordinate
$(Y, Z)$
centred at each jet, as shown in figure 1. The average is denoted by
$ \langle f \rangle$
, fluctuations by
$ f' = f - \langle f \rangle$
and the r.m.s. value by
$ f_{rms} = \langle f'^2 \rangle ^{1/2}$
. Two simulations are performed, one for subsonic and the other for supersonic jet velocities.
The governing equations are the compressible Navier–Stokes equations and the equation of state for an ideal gas, expressed as follows:
\begin{align} \frac {\partial \rho }{\partial t} +\frac {\partial \rho u_{j}}{\partial x_{j}} &=0, \end{align}
\begin{align} \frac {\partial \rho u_{i}}{\partial t} +\frac {\partial \rho u_{i}u_{j}}{\partial x_{j}} &=-\frac {\partial P}{\partial x_{i}} +\frac {\partial \tau _{ij}}{\partial x_{j}}, \end{align}
\begin{align} c_v\frac {\partial \rho T}{\partial t} +c_v\frac {\partial \rho Tu_{j}}{\partial x_{j}} &=-P\frac {\partial u_{j}}{\partial x_{j}} + \frac {\partial }{\partial x_{j}}\left (\kappa \frac {\partial T}{\partial x_{j}}\right ) +\tau _{ij} \frac {\partial u_{i}}{\partial x_{j}}, \end{align}
\begin{align} P&=\rho R T. \end{align}
The viscous stress tensor
$\tau _{ij}$
is written as
\begin{align} \tau _{ij}= \mu \left ( \frac {\partial u_{i}}{\partial x_{j}} +\frac {\partial u_{j}}{\partial x_{i}} -\frac {2}{3}\delta _{ij}\frac {\partial u_{k}}{\partial x_{k}} \right ). \end{align}
Here,
$\rho$
is the density,
$u_i$
is the velocity,
$P$
is the pressure,
$T$
is the temperature,
$\mu$
is the viscosity,
$\kappa$
is the thermal conductivity and
$\delta _{ij}$
is the Kronecker delta. The simulations assume air as the working gas, characterised by a gas constant
$R=287$
J (kg
$\cdot$
K)−1, a specific heat ratio
$\gamma = 1.4$
and a Prandtl number
$Pr = 0.71$
. The viscosity coefficient
$\mu$
is determined using Sutherland’s law.
The mean streamwise velocity of each jet, centred at coordinates
$(y_C, z_C)$
, is defined as
$\langle u \rangle = 0.5U_{J} + 0.5U_{J}\mathrm {tanh} [(D - 2|r|)/4\theta _{J} ]$
with
$ r(y,z) = \sqrt {(y-y_C)^2 + (z-z_C)^2}$
, where
$ U_{J}$
is the initial jet velocity,
$ D$
is the diameter and
$ \theta _J = 0.03D$
represents the shear layer thickness at the jet edge. Hereafter, the subscript
$ J$
refers to quantities related to the initial jet. The initial mean velocity is uniform in the streamwise direction. This methodology follows the DNS studies of a single temporal jet (da Silva & Pereira Reference da Silva and Pereira2008). The combined mean velocity profiles of all jets establish the initial mean velocity field. Weak spatially correlated random noise with a scale of
$0.25D$
and r.m.s. values of
$0.02U_J$
is superimposed on the mean velocity inside the jets to trigger the shear instability (Nagata et al. Reference Nagata, Watanabe and Nagata2018). The r.m.s. value is close to those measured at a jet nozzle in several studies (Ghahremanian et al. Reference Ghahremanian, Svensson, Tummers and Moshfegh2014; Watanabe et al. Reference Watanabe, Sakai, Nagata, Ito and Hayase2014). The initial conditions are set with constant temperature
$ T=300$
K and pressure
$ P=1.013 \times 10^5$
Pa. The flow is characterised by the jet Reynolds number and Mach number, given as
$ Re_J = \rho _J U_J D / \mu _J$
and
$ M_J = U_J / a_J$
, where
$ a_J = \sqrt {\gamma R T_J}$
is the speed of sound in the jet.
2.2. Numerical methods and parameters
The simulations are conducted for
$M_J = 0.6$
and 1.6 with
$Re_J = 3000$
, referred to as M06 and M16, respectively. The corresponding convective Mach numbers, defined as
$ M_C = M_J / 2$
, are 0.3 and 0.8, respectively. For
$ M_C = 0.3$
, the initial shear instability is expected to be minimally affected by compressibility, while for
$ M_C = 0.8$
, the transition to turbulence is delayed, as demonstrated below (Freund et al. Reference Freund, Lele and Moin2000; Pantano & Sarkar Reference Pantano and Sarkar2002). The normalised jet spacing
$S/D$
is set to
$1.5$
and the number of jets,
$N_J^2$
, is
$40^2$
. We have confirmed through additional DNS with a smaller jet number that the formaton of nearly HIT is minimally affected by
$S/D$
between 1.125 and 2.0, which mainly influences the time required for turbulence development. Figure 2(
$a$
) shows the initial field with the local Mach number distribution.
Figure 2. Local Mach number
$M=\sqrt {u^2+v^2+w^2}/a$
on a
$y$
–
$z$
plane for M16 at (
$a$
)
$t/T_J=0$
, (
$b$
)
$t/T_J=8$
and (
$c$
)
$t/T_J=13$
, where
$a$
is the local speed of sound.
The governing equations are solved with an inhouse DNS code based on finite difference schemes, whose details can be found from Yamamoto et al. (Reference Yamamoto, Ishida, Watanabe and Nagata2022a
). Spatial discretisation is achieved through a hybrid scheme that combines a sixth-order central difference scheme and a fifth-order weighted essentially non-oscillatory scheme used together with the advection upstream splitting method (Jiang & Shu Reference Jiang and Shu1996; San & Kara Reference San and Kara2015). The former is applied to smooth regions, while the latter is used for highly compressive regions, defined as
$ \Theta \lt -5\Theta _{rms}$
, where
$ \Theta _{rms}$
represents the r.m.s. fluctuations of the dilatation
$ \Theta =\partial u_j /\partial x_j$
. The equations are temporally integrated using a third-order total-variation-diminishing Runge–Kutta scheme (Gottlieb & Shu Reference Gottlieb and Shu1998).
The computational domain is cubic with a side length of
$ l = N_J S$
. The reference time scale of the jet is defined as
$ T_J = D / U_J$
. Here, the present study defines the time scale using the jet diameter
$ D$
rather than the spacing
$ S$
, as the development of turbulence through jet interaction is characterised by
$ D$
(Tan et al. Reference Tan, Xu, Qi and Ni2023). Time is advanced until
$ t = 5000T_J$
. Statistics are evaluated as functions of time and
$(Y,Z)$
with streamwise averages and ensemble averages across different jets. The time step is set with a constant Courant–Friedrichs–Lewy number of 0.3. As the Kolmogorov scale increases over time, the grid resolution is adjusted accordingly. The numbers of grid points for time intervals
$ 0 \leqslant t/T_J \leqslant 40$
,
$ 40 \leqslant t/T_J \leqslant 100$
,
$ 100 \leqslant t/T_J \leqslant 520$
and
$ 520 \leqslant t/T_J \leqslant 5000$
are denoted by
$ N_1^3$
,
$ N_2^3$
,
$ N_3^3$
and
$ N_4^3$
, respectively. For M06 and M16, the grid resolutions
$ (N_1, N_2, N_3, N_4)$
are
$ (2592, 1728, 1152, 576)$
and
$ (3456, 2592, 1152, 576)$
, respectively. Changes in grid settings are managed using a third-order Lagrange polynomial interpolation. For
$M_J=1.6$
, another DNS was conducted using a single grid setting of
$N_1$
up to
$t/T_J=60$
. By comparing turbulence statistics from this set-up with others, we have confirmed that grid coarsening does not affect the flow evolution. The Kolmogorov scale is defined as
$\eta =(\langle \mu \rangle /\langle \rho \rangle )^{3/4}\varepsilon ^{-1/4}$
with the kinetic energy dissipation rate
$\varepsilon =\langle \tau _{ij}{\mathsf{\textit S}}_{ij}\rangle /\langle \rho \rangle$
, where
${\mathsf{\textit S}}_{ij}$
is the rate-of-strain tensor. The spatial resolution is maintained at less than 1.1 times the Kolmogorov scale for M06 and less than 0.6 times for M16 after nearly HIT forms. A higher spatial resolution for the supersonic case is required to capture shocklets generated with large velocity fluctuations at an early time. This resolution for M16 was determined by the grid sensitivity test for isotropic turbulence with a high turbulent Mach number, where shocklets are generated by turbulent motion (Watanabe et al. Reference Watanabe, Tanaka and Nagata2021).
The integral scales of turbulence increase with time. At the end of the simulations, the longitudinal integral scales of the streamwise velocity, evaluated by integrating the auto-correlation function, are approximately
$ 4D$
in M06 and
$ 3D$
in M16. The side length of the cubic computational domain,
$ l = 60D$
, is approximately 15 and 20 times larger than the integral scales in M06 and M16, respectively. A previous DNS study of decaying isotropic turbulence has suggested that the decay begins to be influenced by confinement effects when the domain size
$ l$
is less than approximately three times the integral scale (Anas et al. Reference Anas, Joshi and Verma2020). The present DNS employs a domain size significantly larger than this threshold, ensuring that confinement effects are negligible.
3. Results and discussions
3.1. Development of nearly homogeneous isotropic turbulence
We first examine the development of turbulence arising from temporal jet interactions. Figure 2 visualises the local Mach number on a
$ y$
–
$ z$
plane at
$ t/T_J = 0$
, 8 and 13 in M16. The initial condition, shown in figure 2(
$ a$
), is characterised by jets, each exhibiting a circular distribution of high streamwise velocity. From the initial state, each jet transitions into a turbulent state in figure 2(
$b$
). The imprint of the jets remains visible at
$ t/T_J = 8$
as a non-uniform velocity distribution. However, as the jets interact with each other, the velocity of the initial jets becomes less pronounced by
$ t/T_J = 13$
.
Figure 3. Temporal variations of (
$a$
) mean streamwise velocity
$\langle u\rangle$
and (
$b$
) r.m.s. fluctuations of streamwise velocity
$u_{rms}$
at
$(Y/S,Z/S)=(0,0)$
,
$(0.5,0)$
and
$(0.5,0.5)$
for M06 and M16.
Figure 3(
$ a$
) shows the temporal variations of the mean streamwise velocity,
$ \langle u \rangle$
, at three locations. As illustrated in figure 1,
$ (Y/S, Z/S) = (0, 0)$
corresponds to the jet centre, while
$ (Y/S, Z/S) = (0.5, 0)$
and
$ (0.5, 0.5)$
are located between the jets. At
$ (Y/S, Z/S) = (0, 0)$
, the potential core regime of the jet is characterised by a constant value of
$ \langle u \rangle / U_J = 1$
, observed for
$ t/T_J \lesssim 5$
in M06 and
$ t/T_J \lesssim 7$
in M16. The turbulence development becomes slower as
$M_J$
increases, as expected from the delay in the turbulent transition of a supersonic jet (Bogdanoff Reference Bogdanoff1983; Nagata et al. Reference Nagata, Watanabe and Nagata2018). As the turbulent jets evolve, the mixing of high-speed jet fluid with low-speed ambient fluid progresses, resulting in a spatially uniform mean velocity distribution. This uniform mean flow is established at
$ t/T_J \approx 25$
. Ghahremanian et al. (Reference Ghahremanian, Svensson, Tummers and Moshfegh2014) conducted near-field measurements for
$ 6 \times 6$
unconfined low-speed jets and observed that a nearly uniform mean velocity in a cross-sectional plane forms at a non-dimensionalised streamwise distance of approximately
$ x/D \approx 20$
.
Figure 3(
$ b$
) shows the r.m.s. fluctuations of the streamwise velocity,
$ u_{rms}$
, at the same three locations. An initial rise in
$u_{rms}$
reflects the turbulence production associated with the non-uniform mean velocity. Following its peak,
$u_{rms}$
starts to decay. As
$ u_{rms}$
decays from its peak, it becomes independent of the positions
$(Y, Z)$
. The uniform distribution of
$\langle u \rangle$
and
$u_{rms}$
is achieved before
$t/T_J \approx 25$
for both
$M_J$
values. The evolution of the mean velocity and r.m.s. velocity aligns well with the experiments on nearly HIT generated by
$ 6 \times 6$
unconfined low-speed jets (Ghahremanian et al. Reference Ghahremanian, Svensson, Tummers and Moshfegh2014). In the present DNS, the flow becomes statistically independent of positions on the the cross-sectional plane. Hereafter, the results at the jet centre are presented.
Figure 4. Temporal variations of the ratio between streamwise and lateral r.m.s. velocity fluctuations
$u_{rms}/v_{rms}$
at
$(Y,Z)=(0,0)$
.
Figure 4 shows the variations of the ratio between
$u_{rms}$
and
$v_{rms}$
. Once the mean streamwise velocity estabilishes the uniform distribution at
$t/T_J\approx 25$
,
$u_{rms}/v_{rms}\approx 1.1$
hardly varies with time. Thus, nearly HIT has developed by the interaction of temporal jets. This value of
$u_{rms}/v_{rms}\approx 1.1$
is consistent with observations in grid turbulence (Krogstad & Davidson Reference Krogstad and Davidson2010; Isaza et al. Reference Isaza, Salazar and Warhaft2014; Kitamura et al. Reference Kitamura, Nagata, Sakai, Sasoh, Terashima, Saito and Harasaki2014). In addition, experiments on both confined and unconfined multiple jets have observed
$ u_{rms} / v_{rms} \approx 1.1 {-} 1.2$
(Ghahremanian et al. Reference Ghahremanian, Svensson, Tummers and Moshfegh2014; Mori et al. Reference Mori, Watanabe and Nagata2024).
Figure 5. Temporal variations of (
$a$
) the turbulent Reynolds number
$Re_\lambda$
and (
$b$
) turbulent Mach number
$M_T$
at
$(Y,Z)=(0,0)$
.
Figure 5(
$ a$
) shows the turbulent Reynolds number,
$ Re_\lambda = \langle \rho \rangle u_{rms} \lambda _x / \langle \mu \rangle$
, calculated using the streamwise velocity. The Taylor microscale,
$ \lambda _x$
, is defined as
$ \lambda _x = u_{rms} / (\partial u / \partial x)_{rms}$
. The Reynolds number decays from its peaks, which are 89 in M06 and 107 in M16. The range of
$ Re_\lambda$
remains within the low-
$ Re_\lambda$
regime, where turbulence decay is influenced by viscous effects, as observed in incompressible turbulence (Meldi & Sagaut Reference Meldi and Sagaut2013).
Figure 5(
$ b$
) presents the turbulent Mach number,
$ M_T = \sqrt {u_{rms}^2 + v_{rms}^2 + w_{rms}^2} / \sqrt {\gamma R \langle T \rangle }$
, where
$\sqrt {\gamma R \langle T \rangle }$
is the speed of sound based on the mean temperature. As the mean temperature remains nearly constant over time, except during the initial transitional regime, the variation of
$ M_T$
closely follows that of
$ u_{rms}$
. The maximum values of
$ M_T$
are 0.19 and 0.53 for M06 and M16, respectively, indicating that compressibility effects on turbulence are not negligible for M16. For a single jet at the same Mach number of 1.6, it has been shown that the temporally evolving supersonic jet initially generates pressure waves resembling spherical shock waves, characterised by strong compression followed by blunt expansion (Nagata et al. Reference Nagata, Watanabe and Nagata2018). In the case of statistically steady isotropic turbulence subjected to solenoidal linear forcing, the statistical properties at small scales, such as the
$ Re_\lambda$
dependence of velocity derivative skewness and flatness, are not affected by compressibility effects even at
$ M_T = 0.3$
(Watanabe et al. Reference Watanabe, Tanaka and Nagata2021). In M06, the maximum value of
$ M_T$
remains below 0.2, indicating negligible compressibility effects due to turbulent velocity fluctuations, as will also be discussed below.
Figure 6. (
$a$
) Temporal variations of the r.m.s. fluctuations of density (
$ \rho _{rms}$
) and pressure (
$ P_{rms}$
), normalised by the mean density and mean pressure, respectively. (
$b$
) Dependence of the normalised r.m.s. pressure fluctuation on the turbulent Mach number (
$ M_T$
). All results are evaluated at
$ (Y, Z) = (0, 0)$
.
3.2. Compressibity effects
Figure 6(
$a$
) presents the r.m.s. values of density and pressure fluctuations, denoted
$\rho _{rms}$
and
$P_{rms}$
, respectively, at the jet centre. The r.m.s. fluctuations are normalised by the mean density
$\langle \rho \rangle$
or mean pressure
$\langle P \rangle$
. The maximum values of
$\rho _{rms}$
and
$P_{rms}$
in M06 are less than 2 % of the mean density and pressure, respectively. During the decay of nearly HIT in M06,
$\rho _{rms}/\langle \rho \rangle$
and
$P_{rms}/\langle P \rangle$
are as small as 0.1 %, indicating that compressibility effects appear insignificant for this case. However,
$\rho _{rms}/\langle \rho \rangle$
and
$P_{rms}/\langle P \rangle$
increase up to approximately 10 % in M16. These highest fluctuations are observed in an early time, and both density and pressure fluctuations decay with time. As shown below, the shock waves are generated in an early time in M16, influencing turbulence at a later time by the long-time propagation.
Figure 6(
$b$
) plots
$ P_{rms} / \langle P \rangle$
as a function of the turbulent Mach number
$ M_T$
. The DNS data are shown for the time period corresponding to the decay of
$ M_T$
. In statistically steady isotropic turbulence subjected to solenoidal forcing (Wang et al. Reference Wang, Gotoh and Watanabe2017; Watanabe et al. Reference Watanabe, Tanaka and Nagata2021), the normalised r.m.s. pressure fluctuations follow the relation:
\begin{align} P_{rms}/\langle P \rangle =A\gamma M_T^2/3, \end{align}
with
$A\approx 1$
, as originally derived by Donzis & Jagannathan (Reference Donzis and Jagannathan2013). Figure 6(
$b$
) also illustrates this relation with
$ A = 1$
. The turbulence generated by jet interaction in each case follows this relation at large
$ M_T$
, which corresponds to the early period of the decay. Thus, (3.1) is valid for a short time after turbulence is generated by the jet interaction. Equation (3.1) provides a time-local relationship between pressure fluctuations and velocity fluctuations, assuming that pressure fluctuations at a given time are related to fluid motions at the same time. However, in decaying compressible turbulence, pressure waves, such as shock waves generated at an early time, decay more slowly than velocity fluctuations. As
$ M_T$
decreases with time, (3.1) begins to underestimate
$ P_{rms} / \langle P \rangle$
because the pressure fluctuations at time
$ t = t_0$
can be partially influenced by turbulent motions at earlier times (
$ t \lt t_0$
), when
$ M_T$
was higher than at
$ t = t_0$
. A similar non-local compressibility effect due to pressure wave propagation has been reported in the context of velocity gradient statistics (Yamamoto et al. Reference Yamamoto, Watanabe and Nagata2022b
).
Figure 7. (
$a$
) Distribution of dilatation
$\Theta =\partial u_i /\partial x_i$
on a
$y$
–
$z$
plane at
$t/T_J=80$
for M16. (
$b$
) Probability density function of
$\Theta$
normalised by its r.m.s. value
$\Theta _{rms}$
. The black line indicates the Gaussian distribution.
The supersonic jets in M16 generate shock waves, which propagate in the periodic domain, even though
$M_T$
decays with time. Figure 7(
$a$
) visualises dilatation
$\Theta =\partial u_i/\partial x_i= -(1/\rho )D\rho /Dt$
normalised by its r.m.s. value
$\Theta _{rms}$
in M16. The visualised instance is
$t/T_J=80$
, at which the turbulent Mach number is smaller than 0.1. Shock waves are identified as thin layers with large negative
$\Theta$
, which indicate strong fluid compression. Figure 7(
$b$
) presents the probability density function (p.d.f.) of
$\Theta /\Theta _{rms}$
from
$t/T_J=38$
to 340. The presence of shock waves affects the p.d.f. of
$\Theta$
, which is negatively skewed under the influence of shock waves because thin shock waves with large negative
$\Theta$
occupy a small fraction of the flow. The p.d.f. follows the Gaussian distribution in M06, where no shock waves propagate in turbulence. Negatively skewed distributions are observed in M06. The shape of the p.d.f. in M16 does not change with time, indicating that the shock waves propagate for long time, affecting the late-time behaviour of turbulence even at low
$M_T$
.
In the case of a single jet, the evolution of various velocity statistics in temporal simulations shows similarities to those observed in a spatial jet (da Silva & Pereira Reference da Silva and Pereira2008; van Reeuwijk & Holzner Reference van Reeuwijk and Holzner2014; Hayashi et al. Reference Hayashi, Watanabe and Nagata2021b ). However, the interaction of temporal supersonic jets exhibits notable differences compared with the wind tunnel experiments. In the experiments, a counter pressure gradient, observed before a homogeneous mean flow fully develops due to jet interaction, causes a deceleration of the mean flow, associated with a weak compression (Mori et al. Reference Mori, Watanabe and Nagata2024). This deceleration is not observed in the temporal jets. Moreover, shock wave propagation in the wind tunnel is influenced by the test section walls, which cause the waves to reflect. This wave reflection does not occur in DNS due to periodic boundary conditions. However, wave propagation from the near-nozzle region towards the decay region in the wind tunnel bears similarity to that from early to late times in temporal simulations. Both spatial and temporal jets produce turbulence through the mean shear associated with the non-uniform mean velocity of the jets. The cross-sectional distribution of mean shear is virtually identical for both spatial and temporal jets.
Figure 8. (
$a$
) Dependence of the velocity derivative skewness,
$ S_{\partial u / \partial x}$
, on the turbulent Reynolds number,
$ Re_\lambda$
. The DNS results are taken at
$(Y, Z) = (0, 0)$
after
$ Re_\lambda$
reaches its peak. (
$b$
) Dependence of the velocity derivative flatness,
$ F_{\partial u / \partial x}$
, on
$ Re_\lambda$
. The DNS results are taken at
$(Y, Z) = (0, 0)$
after
$ Re_\lambda$
reaches its peak. The figure includes data from studies of incompressible turbulence (Kuo & Corrsin Reference Kuo and Corrsin1971; Antonia & Chambers Reference Antonia and Chambers1980; Van Atta & Antonia Reference Van Atta and Antonia1980; Kerr Reference Kerr1985; Mydlarski & Warhaft Reference Mydlarski and Warhaft1996; Sreenivasan & Antonia Reference Sreenivasan and Antonia1997; Mi et al. Reference Mi, Xu and Zhou2013; Watanabe & Nagata Reference Watanabe and Nagata2018; Watanabe et al. Reference Watanabe, Zhang and Nagata2019; Zheng et al. Reference Zheng, Nagata and Watanabe2021).
Figure 8 plots the skewness,
$ S_{\partial u / \partial x} = \langle (\partial u' / \partial x)^3 \rangle / \langle (\partial u' / \partial x)^2 \rangle ^{3/2}$
, and flatness,
$ F_{\partial u / \partial x} = \langle (\partial u' / \partial x)^4 \rangle / \langle (\partial u' / \partial x)^2 \rangle ^2$
, of the longitudinal velocity gradient
$ \partial u / \partial x$
as functions of
$ Re_\lambda$
. The results are shown after
$ Re_\lambda$
reaches its peak and are compared with previous studies of incompressible turbulence. The flatness values in both M06 and M16 align well with those of incompressible turbulence at comparable
$ Re_\lambda$
, decreasing as
$ Re_\lambda$
declines with turbulence decay. Experiments and DNS have shown that statistically steady compressible turbulence tends to exhibit higher flatness values than incompressible turbulence (Donzis & John Reference Donzis and John2020; Watanabe et al. Reference Watanabe, Tanaka and Nagata2021; Yamamoto et al. Reference Yamamoto, Ishida, Watanabe and Nagata2022a
). However, this trend is not observed in the present DNS of decaying turbulence. The dependence on Mach number is more pronounced for the skewness. It is known that skewness plotted against
$ Re_\lambda$
shows more flow dependence than flatness, even in incompressible turbulence (Zheng et al. Reference Zheng, Nagata and Watanabe2021). In M06,
$-S_{\partial u / \partial x}$
decreases significantly for
$ Re_\lambda \lesssim 10$
, consistent with grid turbulence data reported by Zheng et al. (Reference Zheng, Nagata and Watanabe2021). In M16, this decrease occurs at larger
$ Re_\lambda$
and
$-S_{\partial u / \partial x}$
tends to be smaller than in M06. This reduction in
$-S_{\partial u / \partial x}$
due to compressibility effects has been reported in other DNS studies of decaying compressible turbulence, although it has received limited attention (Samtaney et al. Reference Samtaney, Pullin and Kosović2001; Li et al. Reference Li, Fu, Ma and Liang2010). While Samtaney et al. (Reference Samtaney, Pullin and Kosović2001) also found that decaying compressible turbulence with high initial
$ M_T$
has smaller
$-S_{\partial u / \partial x}$
values than low-
$ M_T$
cases, the difference is less significant than in the present DNS. Samtaney et al. (Reference Samtaney, Pullin and Kosović2001) conducted DNS over a period of less than 10 times the integral time scale. In contrast, the present DNS covers a much longer duration, over 2000 times the integral time scale at the time of maximum
$ Re_\lambda$
, which might explain the larger deviation of
$-S_{\partial u / \partial x}$
in M16 from incompressible flow data.
Figure 9. Temporal variations of
$u_{rms}^2L_u^3$
and
$u_{rms}^2L_u^5$
in (
$a$
) M06 and (
$b$
) M16.
Figure 10. (
$a$
) Dependence of the non-dimensional energy dissipation rate,
$ C_\varepsilon$
, on the turbulent Reynolds number,
$ Re_\lambda$
. (
$b$
) Temporal variations of
$ C_\varepsilon$
. The data are taken for the period following the peak of
$ Re_\lambda$
.
3.3. Decay properties and energy spectra
The decay laws of (1.1) and (1.2) are based on the invariance of
$ u_{rms}^2L_u^5$
and
$ u_{rms}^2L_u^3$
during the decay, respectively, under the assumption of self-similarity at large scales. Here,
$ L_u$
denotes the longitudinal integral scale of
$ u$
. At a given time,
$ L_u$
is evaluated using the auto-correlation function of
$ u$
,
$ f_u(r_x, Y, Z) = \langle u'(x, Y, Z) u'(x + r_x, Y, Z) \rangle / u_{rms}^2(Y, Z)$
, as
$ L_u(Y, Z) = \int _0^{r_{x0}} f_u(r_x, Y, Z) {\rm{d}}r_x$
, where
$ r_{x0}$
is the
$ r_x$
value at which
$ f_u$
first crosses zero. Figure 9 shows the temporal variations of
$ u_{rms}^2L_u^5$
and
$ u_{rms}^2L_u^3$
for both Mach number cases. Nearly HIT develops by approximately
$ t/T_J = 25$
. However, both
$ u_{rms}^2L_u^5$
and
$ u_{rms}^2L_u^3$
vary with time even after turbulence generation. In the present DNS,
$ u_{rms}^2L_u^3$
decreases consistently over time until
$ t/T_J = 5000$
for both cases. Conversely,
$ u_{rms}^2L_u^5$
becomes nearly time-independent at later times, particularly after
$ t/T_J \approx 2000$
for M06 and
$ t/T_J \approx 500$
for M16. The variations of
$ u_{rms}^2 L_u^5$
observed in the present DNS are consistent with findings from previous studies on decaying isotropic turbulence initialised with an energy spectrum
$ E(k) \sim k^4$
. The DNS results of Ishida et al. (Reference Ishida, Davidson and Kaneda2006) indicate that in such cases,
$ u_{rms}^2 L_u^5$
increases over time before reaching a constant value. In addition, grid turbulence exhibits a similar behaviour (Watanabe & Nagata Reference Watanabe and Nagata2018):
$ u_{rms}^2L_u^3$
becomes time-independent only after the flow has evolved for a sufficiently long time in a nearly homogeneous and isotropic state. For
$ t/T_J \geqslant 2000$
, the least squares method yields
$ u_{rms}^2L_u^5 / U_J^2D^5 \sim (t/T_J)^{0.038}$
and
$ u_{rms}^2L_u^3 / U_J^2D^3 \sim (t/T_J)^{-0.561}$
for M06. Similarly, the results for M16 yield
$ u_{rms}^2L_u^5 / U_J^2D^5 \sim (t/T_J)^{0.021}$
and
$ u_{rms}^2L_u^3 / U_J^2D^3 \sim (t/T_J)^{-0.533}$
. Thus,
$ u_{rms}^2L_u^5$
exhibits a weaker time dependence than
$ u_{rms}^2L_u^3$
, suggesting that the decay law of Batchelor turbulence is valid.
Figure 10 plots
$C_\varepsilon = \varepsilon /(u_{rms}^3/L_u)$
as a function of
$Re_\lambda$
or time. Initially,
$C_\varepsilon$
exhibits an increase with time, following the scaling
$C_\varepsilon \sim Re_\lambda ^{-1}$
of non-equilibrium turbulence (Valente & Vassilicos Reference Valente and Vassilicos2014). This scaling near the nozzles has also been observed in the near-field of grid turbulence (Valente & Vassilicos Reference Valente and Vassilicos2014; Mora et al. Reference Mora, Muñiz Pladellorens, Riera Turró, Lagauzere and Obligado2019) as well as the near field of turbulence generated by the jet interaction (Mori et al. Reference Mori, Watanabe and Nagata2024). In M06, a gradual increase in
$C_\varepsilon$
continues beyond the non-equilibrium phase due to the low-
$Re_\lambda$
effects. Conversely, in M16,
$C_\varepsilon$
declines after the non-equilibrium phase is over. Shock wave propagation, from an early time with large
$M_T$
in M16, amplifies velocity fluctuations and reduces integral scales, potentially contributing to the observed decrease in
$C_\varepsilon$
(Kitamura et al. Reference Kitamura, Nagata, Sakai, Sasoh and Ito2017; Tanaka et al. Reference Tanaka, Watanabe, Nagata, Sasoh, Sakai and Hayase2018). Such deviations in small-scale statistics from incompressible values are also documented for inhomogeneous or decaying compressible turbulence (Samtaney et al. Reference Samtaney, Pullin and Kosović2001; Xinliang et al. Reference Xinliang, Dexun and Yanwen2002; Yamamoto et al. Reference Yamamoto, Ishida, Watanabe and Nagata2022a
).
Figure 11. Temporal variations of
$ u_{rms}^2$
,
$ L_u$
and
$ C_\varepsilon$
over time: (
$a$
) M06; (
$b$
) M16. Thin black solid lines indicate power laws with exponents determined using the least squares method. Dashed and dash-dotted line lines represent the power laws for
$ u_{rms}^2$
and
$ L_u$
with the exponents for Batchelor and Saffman turbulence.
The decay properties are investigated during the phase with nearly constant
$ u_{rms}^2 L_u^5$
to compare the variations of
$ u_{rms}$
and
$ L_u$
with theoretical predictions. The examined time intervals are
$ 2000 \leqslant t/T_J \leqslant 5000$
for M06 and
$ 800 \leqslant t/T_J \leqslant 5000$
for M16. Additionally, a sensitivity test for the fitting range is provided below. The present study considers the following power laws:
\begin{align} u_{rms}^2\sim (t-t_0)^n,\quad L_u\sim (t-t_0)^{n_L},\quad C_\varepsilon \sim (t-t_0)^{n_C}, \end{align}
where
$ t_0$
is the virtual origin. The value of
$ t_0$
is determined from the decay of
$ u_{rms}^2$
. Following the methodology in a previous study, a linear least squares method is applied with a predetermined
$ t_0$
that minimises the fitting error for
$ u_{rms}^2$
(Watanabe et al. Reference Watanabe, Zheng and Nagata2022). Specifically, the fitting for
$ (t - t_0, u_{rms}^2)$
using the linear least squares method is iteratively performed by varying
$ t_0$
in the range of
$ -300T_J \leqslant t_0 \leqslant 300T_J$
with an increment of
$ 0.1T_J$
. The
$ t_0$
value that minimises the error is selected and used to determine the power law exponents
$ (n, n_L, n_C)$
. This approach has been shown to yield decay exponents comparable to those obtained using a nonlinear least squares method (Levenberg–Marquardt algorithm), which directly determines
$ (t_0, n)$
from the data
$ (t, u_{rms}^2)$
(Watanabe et al. Reference Watanabe, Zheng and Nagata2022). The values of
$ t_0$
determined here are also used in the subsequent analysis. A least squares method is applied to analyse
$ L_u$
and
$ C_\varepsilon$
, thereby determining the power law exponents
$ n_L$
and
$ n_C$
. Figure 11 presents the variations of
$ u_{rms}^2$
,
$ L_u$
and
$ C_\varepsilon$
with
$ t - t_0$
, compared with the power laws (thin black solid lines) using the exponents evaluated by the described method. The variations of these quantities are well approximated by the power laws. When confinement effects due to a finite computational domain become significant, the decay of
$ u_{rms}^2$
accelerates (Anas et al. Reference Anas, Joshi and Verma2020). In the present study, such confinement-induced acceleration is not observed in the decay of
$ u_{rms}^2$
.
Table 1 summarises the virtual origin
$ t_0$
and the power law exponents (
$ n$
,
$ n_L$
,
$ n_C$
), comparing them with the theoretical predictions for Saffman and Batchelor turbulence. The theoretical values from (1.1) and (1.2) are evaluated using the DNS-derived values of
$ n_C$
and are denoted as
$ ( n^{(B)}, n_L^{(B)} )$
for Batchelor turbulence, and
$ ( n^{(S)}, n_L^{(S)} )$
for Saffman turbulence. The behaviour of
$ u_{rms}^2L_u^5$
and
$ u_{rms}^2L_u^3$
indicates that the decay aligns with the Batchelor turbulence theory. For both Mach number cases,
$ n$
and
$ n_L$
closely match the values predicted for Batchelor turbulence. The table also includes the relative differences between the exponents and the theoretical values, defined as
$\Delta n^{(\alpha )} = (n - n^{(\alpha )})/n^{(\alpha )}$
and
$\Delta n_L^{(\alpha )} = (n_L - n_L^{(\alpha )})/n_L^{(\alpha )}$
, where
$ \alpha = B$
(Batchelor) or
$ S$
(Saffman). For Batchelor turbulence, the relative difference in exponents is less than 4 %, while for Saffman turbulence, the differences are approximately 15%–27 % for
$ u_{rms}^2$
and
$ L_u$
. As also visualised in figure 11, the theory for Saffman turbulence predicts a slower decay of
$u_{rms}^2$
and a faster growth of
$ L_u$
than the DNS results. Although the exponents
$ (n, n_L)$
in both M06 and M16 are evaluated during the period with nearly constant
$ u_{rms}^2 L_u^5$
,
$ (n, n_L)$
differ between these two cases. This difference is attributed to the variation of
$ C_\varepsilon$
. Thus, compressibility effects influence the decay laws through
$ C_\varepsilon$
. These findings suggest that the turbulence decay in the jet interaction follows the Batchelor turbulence theory. The Appendix provides comparisons of the decay exponents with the theories of Batchelor and Saffman turbulence using experimental data available in the literature (Tan et al. Reference Tan, Xu, Qi and Ni2023; Mori et al. Reference Mori, Watanabe and Nagata2024). The exponents evaluated in these studies are also consistent with the predictions of Batchelor turbulence. However, these experiments provided measurement results over a limited range of streamwise positions. Further experimental investigations are encouraged to comprehensively assess the decay of turbulence generated by jet interactions.
Table 1. Power law exponents for
$ u_{rms}^2$
,
$ L_u$
and
$ C_\varepsilon$
, denoted as
$ n$
,
$ n_L$
and
$ n_C$
, respectively, and the virtual origin
$ t_0$
. The theoretical values of Batchelor and Saffman turbulence for
$ n$
and
$ n_L$
are denoted by superscripts
$ (B)$
and
$ (S)$
, respectively. The relative differences between the DNS results and the theoretical predictions are calculated as
$\Delta n^{(\alpha )} = (n - n^{(\alpha )})/n^{(\alpha )}$
and
$ \Delta n_L^{(\alpha )} = (n_L - n_L^{(\alpha )})/n_L^{(\alpha )}$
, where
$ \alpha = B$
for Batchelor turbulence and
$ \alpha = S$
for Saffman turbulence.
Figure 12 examines the dependence of the decay exponent
$ n$
of
$ u_{rms}^2$
on the fitting range
$ t_S \leqslant t \leqslant t_E$
, with varying
$ t_S$
or
$ t_E$
. The decay exponents only slightly vary between
$ -1.554$
and
$ -1.549$
for M06 in figure 12(
$a$
,
$b$
) and between
$ -1.472$
and
$ -1.464$
for M16 in figure 12(
$c$
,
$d$
), depending on the choice of
$ t_S$
and
$ t_E$
. This demonstrates that the selection of the fitting range hardly affects the evaluation of the decay exponent. Additionally, different virtual origins from those determined above were tested. Changing the virtual origin by 10% resulted in a change of the decay exponent
$ n$
by less than 2 %. This weak dependence on the virtual origin is partly attributed to the fitting range encompassing much larger
$ t$
values compared with
$ t_0$
, a phenomenon also observed in grid turbulence (Watanabe et al. Reference Watanabe, Zheng and Nagata2022).
Figure 12. Dependence of the decay exponent
$ n$
of
$ u_{rms}^2$
on the fitting range
$ t_S \leqslant t \leqslant t_E$
: (
$a$
)
$ t_S$
dependence with
$ t_E = 5000T_J$
fixed in M06; (
$b$
)
$ t_E$
dependence with
$ t_S = 2000T_J$
fixed in M06; (
$c$
)
$ t_S$
dependence with
$ t_E = 5000T_J$
fixed in M16; (
$d$
)
$ t_E$
dependence with
$ t_S = 800T_J$
fixed in M16.
Figure 13. Temporal evolutions of the three-dimensional energy spectrum,
$ E(k)$
, for (
$a$
) M06 and (
$b$
) M16. The spectra are premultiplied by
$k^{-4}$
or
$ k^{-2}$
to facilitate comparison with
$ E(k) \sim k^4$
of Batchelor turbulence and
$ E(k) \sim k^2$
of Saffman turbulence.
Figure 13 presents the three-dimensional energy spectrum
$E(k)$
during the decay. The calculation of
$E(k)$
employs the shell averaging method in wavenumber space (Valente et al. Reference Valente, da Silva and Pinho2016). For comparisons with
$ E(k) \sim k^4$
and
$ k^2$
,
$ E(k)$
is multiplied by
$ k^{-4}$
or
$ k^{-2}$
so that these power laws appear as horizontal lines. For both Mach numbers,
$k^{-4}E$
tends towards a constant value and
$k^{-2}E$
decreases following
$k^{-2}E\sim k^2$
as
$ k$
becomes small, indicating that the energy spectrum at small wavenumbers follows a
$ k^4$
scaling rather than
$ k^2$
. This behaviour confirms the formation of Batchelor turbulence. The spectral slope within the
$ k^4$
range remains consistent over time, with small variations related to statistical convergence. The formation of
$ E(k) \sim k^4$
occurs earlier than the invariance of
$ u_{rms}^2 L_u^5$
observed in figure 9. It should be noted that
$ E(k) \sim k^4$
is not a sufficient condition for the constancy of
$ u_{rms}^2 L_u^5$
in Batchelor turbulence. In DNS of decaying isotropic turbulence, where the flow is initialised with a prescribed spectral shape of
$ E(k) \sim k^4$
,
$ u_{rms}^2 L_u^5$
becomes time-independent only after the turbulence evolves over a certain period (Ishida et al. Reference Ishida, Davidson and Kaneda2006). This behaviour is consistent with the turbulence generated by jet interactions observed in the present study.
4. Conclusion
DNS is performed to investigate the interaction of temporally evolving circular jets with jet Mach numbers
$M_J=0.6$
and 1.6, and a jet Reynolds number of 3000. For both cases, the jet interaction leads to the formation of nearly HIT, with the ratio of streamwise to transverse r.m.s. velocity fluctuations around 1.1, although the turbulent transition is delayed in the higher
$ M_J$
case. The r.m.s. pressure fluctuations,
$ P_{rms}$
, normalised by the mean pressure, decrease as the turbulent Mach number decreases. Their relationship shortly after turbulence generation is well described by (3.1), which is often applied to statistically steady isotropic turbulence, although it underestimates
$ P_{rms}$
after the turbulence has sufficiently decayed. In the
$M_J=1.6$
case, the turbulent Mach number becomes sufficiently high to generate shock waves, confirmed by the negatively skewed probability distribution of dilatation. These shock waves propagate through the turbulence over a long period and continue to affect the flow at later times, even after the turbulent Mach number decreases. Consistent with previous studies of decaying compressible turbulence (Samtaney et al. Reference Samtaney, Pullin and Kosović2001; Li et al. Reference Li, Fu, Ma and Liang2010), the magnitude of the velocity derivative skewness is smaller for
$M_J=1.6$
than for
$M_J=0.6$
. Additionally, the non-dimensional dissipation rate
$C_\varepsilon$
is lower for
$M_J=1.6$
. Thus, the statistical properties of small-scale motions are altered by compressibility effects. The late-time decay of turbulence is consistent with Batchelor turbulence theory. For the r.m.s. velocity fluctuations
$ u_{rms}$
and the integral scale
$ L_u$
of the streamwise velocity,
$ u_{rms}^2 L_u^5$
remains nearly constant over time, while
$ u_{rms}^2 L_u^3$
continuously decreases, indicating that the turbulence follows Batchelor turbulence theory. The power-law exponents of
$ u_{rms}^2$
and
$ L_u$
are also consistent with the theoretical predictions for Batchelor turbulence, with corrections due to the temporal variation of
$ C_\varepsilon$
. The three-dimensional energy spectrum further supports this, exhibiting
$ E(k) \sim k^4$
for small
$k$
assumed in Batchelor turbulence. This
$ k^4$
spectral shape develops even before
$ u_{rms}^2 L_u^5$
becomes nearly constant. These behaviours, consistent with Batchelor turbulence, are observed even for
$M_J=1.6$
, where small-scale turbulence characteristics deviate from incompressible flows due to compressibility effects. The modulations of small-scale motions, such as shock wave propagation, appear to have little influence on the decay laws determined by the large-scale flow characteristics. The present numerical simulations provide new insights and encourage further experimental investigations into Batchelor turbulence generated by jet interactions.
Acknowledgements
DNS were performed using the high-performance computing systems at the Japan Agency for Marine-Earth Science and Technology and Nagoya University. This work was supported by Collaborative Research Project on Computer Science with High-Performance Computing in Nagoya University.
Funding
This work was supported by JSPS KAKENHI Grant Nos. JP22K03903 and JP22H01398.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Appendix. Experimental studies of decaying turbulence generated by jet interaction
The interaction of jets has also been investigated in early studies of turbulence generated by jet grids, which consist of a conventional grid combined with a nozzle array (Gad-el Hak & Corrsin Reference Gad-el Hak and Corrsin1974). In these set-ups, jets are injected from the grid into a coflow. However, in most experiments, the jets operate at very low flow rates, resulting in turbulence generation dominated by the grid itself. An exception is the study by Tan et al. (Reference Tan, Xu, Qi and Ni2023), which examined jet grids with a very weak coflow, where turbulence is primarily generated by the jets. Similarly, Mori et al. (Reference Mori, Watanabe and Nagata2024) reported experiments on turbulence generated by the interaction of numerous supersonic jets, which decay within a confined test section. In this appendix, we analyse the measurement results from these studies. Fundamental velocity statistics related to homogeneity and isotropy were reported in the original studies (Tan et al. Reference Tan, Xu, Qi and Ni2023; Mori et al. Reference Mori, Watanabe and Nagata2024; Watanabe et al. Reference Watanabe, Mori, Ishizawa and Nagata2024). The mean velocity generated by continuous jet interactions becomes homogeneous within the test section. These studies examined parallel jets issued into a test section, the end of which functions as an outlet. Consequently, the mean velocity profile is unlikely to be influenced by large-scale circulations, which may play a significant role when continuous jets are issued into a closed vessel. Other experiments on HIT generated by jet interaction have been conducted using unsteady-jet generators (Hwang & Eaton Reference Hwang and Eaton2004; Variano et al. Reference Variano, Bodenschatz and Cowen2004; Bellani & Variano Reference Bellani and Variano2014; Carter et al. Reference Carter, Petersen, Amili and Coletti2016; Pérez-Alvarado et al. Reference Pérez-Alvarado, Mydlarski and Gaskin2016; Yamamoto et al. Reference Yamamoto, Watanabe and Nagata2022b ). However, the interaction of unsteady jets generates turbulence differently from continuous jets. In continuous jets, temporal velocity fluctuations are initially produced by shear instability arising from the mean shear. In contrast, unsteady jets generate velocity fluctuations directly through their operational variations. Since the decay of HIT strongly depends on the turbulence generation mechanism, we focus on the interaction of continuous jets.
First, we examine the velocity data measured in the vertical octagonal non-corrosive stirred energetic turbulence (V-ONSET) facility reported by Tan et al. (Reference Tan, Xu, Qi and Ni2023). Velocity measurements using particle image velocimetry were conducted for a nozzle spacing of
$ S = 40$
mm and a nozzle diameter of
$D=5$
mm. A square grid with jet holes was installed in a water tunnel with an octagonal cross-section test section. Jets were issued from the grid at a velocity of 5.5 ms−1 within a coflow at 0.27 ms−1. The measurements covered a streamwise distance
$ x$
normalised by the jet diameter
$ D$
within
$ 53 \leqslant x/D \leqslant 165$
. As observed in the present DNS, their measurements confirmed the generation of nearly HIT decaying within the test section, where a power-law decay was observed for
$ u_{rms}$
. The integral scale, evaluated with the auto-correlation function, was not available in this study. Therefore, the discussion focuses solely on the decay of
$ u_{rms}^2$
. The parameters in the scaling
$ u_{rms}^2 \sim (x - x_0)^{n}$
were determined with the same method as in the present DNS, using a linear least squares method with a predetermined virtual origin
$ x_0$
to minimise the fitting error. Their data yielded
$ x_0 = 2.5D$
and
$ n = -2.2$
, with the nozzle diameter
$ D$
. In decaying HIT, the relation
$ u_{rms}^2 \sim (x - x_0)^{n}$
implies
$ \varepsilon \sim (x - x_0)^{(n-1)}$
, where
$ \varepsilon$
is the dissipation rate. The turbulent Reynolds number,
$ Re_\lambda = u_{rms} \lambda / \nu$
, defined with the Taylor microscale
$ \lambda = u_{rms} \sqrt {15\nu / \varepsilon }$
, evolves as
$ Re_\lambda \sim (x - x_0)^{(n+1)/2}$
. The measurements were conducted in the near-nozzle region of
$ 53 \leqslant x/D \leqslant 165$
. In turbulence generated by jet interaction in this range,
$ C_\varepsilon$
increases as turbulence decays, following the non-equilibrium dissipation scaling (Mori et al. Reference Mori, Watanabe and Nagata2024). In this non-equilibrium region, where
$ C_\varepsilon \sim Re_\lambda ^{-1}$
, it follows that
$ C_\varepsilon \sim (x - x_0)^{-(n+1)/2}$
. The power law decay of
$ u_{rms}^2$
with
$ n = -2.2$
predicts
$ C_\varepsilon \sim (x - x_0)^{0.6}$
. With this exponent for
$C_\varepsilon$
, (1.1) and (1.2) predict
$ n = -2.3$
for Batchelor turbulence and
$ n = -1.9$
for Saffman turbulence. The experimentally derived exponent
$ n = -2.2$
closely aligns with the Batchelor turbulence theory. Although this analysis assumes non-equilibrium turbulence, which is not explicitly verified in their experimental data, the jet grid used by Tan et al. (Reference Tan, Xu, Qi and Ni2023) likely generates turbulence consistent with the theory of Batchelor turbulence.
Additionally, the velocity data measured in the multiple-jet wind tunnel described by Mori et al. (Reference Mori, Watanabe and Nagata2024) are examined here. This dataset was also analysed by Watanabe et al. (2024). These papers include various velocity statistics, comparisons with other turbulent flows and uncertainty tests. This facility generates decaying, nearly HIT through the interaction of
$6\times 6$
parallel supersonic jets within a test section of 1 m length and a cross-section of
$ 0.01 \times 0.01 \, \text {m}^2$
. The nozzle diameter is
$D=4.31$
mm and the nozzle spacing is
$S=12$
mm. Velocity measurements using particle tracking velocimetry (PIV) were conducted for ideally expanded, supersonic parallel jets with a jet Mach number of
$M_J = 1.36$
and a jet Reynolds number of
$Re_J = 1.9 \times 10^5$
. The measurements were conducted for
$25\leqslant x/D \leqslant 180$
. The experiments in the nearly HIT region were conducted repeatedly until approximately 700 statistically independent vector profiles were collected. Velocity statistics are evaluated by taking ensemble averages. It was observed that nearly HIT forms before
$x/D \approx 80$
. As the integral scale grows with the decay of turbulence, the test-section side wall causes confinement effects on turbulence. Beyond
$x/D\approx 140$
, the confinement effects become prevalent, halting the increase in the integral scale and accelerating the decay of velocity fluctuations, as also observed in decaying HIT with confinement effects (Skrbek & Stalp Reference Skrbek and Stalp2000; Morize & Moisy Reference Morize and Moisy2006). Consequently, decay properties are analysed within the range of
$81 \leqslant x/D \leqslant 139$
. Because turbulence decays along the streamwise direction, the spatial longitudinal auto-correlation function of
$u$
is assessed using two statistically different points from PIV data, complicating interpretation. Therefore, the present study focuses on the statistics of
$v$
. The longitudinal auto-correlation function of
$v$
is evaluated using ensemble and spatial averages to mitigate statistical errors. Within each PIV measurement area, the correlation function is computed at eight streamwise locations, spaced equidistantly. Spatial averages for
$f_v$
are calculated over a length of approximately
$2S$
in the
$x$
direction and over the full height of the measurement area in the
$y$
direction. The integral of
$f_v$
provides the integral scale
$L_v$
.
Figure 14. Streamwise variations of velocity variance
$v_{rms}^2$
, integral scale
$L_v$
and non-dimensional dissipation rate
$C_\varepsilon$
. The error bars indicate possible statistical errors estimated with reduced samples. The velocity variance is normalised by the uniform mean velocity in the decay region,
$U_0$
.
The energy dissipation rate
$\varepsilon$
was evaluated from the decay of TKE,
$k_T=(u_{rms}^2+2v_{rms})/2$
with the TKE transport equation in decaying HIT, and the non-dimensional dissipation rate was assessed as
$C_\varepsilon =\varepsilon /(v_{rms}^3/L_v)$
, as detailed by Mori et al. (Reference Mori, Watanabe and Nagata2024). Here, the streamwise gradient of
$k_T$
,
$\partial k_T/\partial x$
, is evaluated by fitting a power law to
$k_T(x)$
in a dimensional form and is used to assess
$\varepsilon$
with the TKE equation. The evaluation
$\varepsilon$
is validated by comparing the inertial range statistics normalised
$\varepsilon$
and kinematic viscosity with other turbulent flows (Mori et al. Reference Mori, Watanabe and Nagata2024; Watanabe et al. 2024). Figure 14 depicts the variations of
$v_{rms}^2$
,
$L_v$
and
$C_\varepsilon =\varepsilon /(v_{rms}^3/L_v)$
as functions of
$(x-x_0)/D$
, where
$ x_0$
is the virtual origin of the TKE decay. Possible statistical errors are evaluated by dividing staistical samples randomly into two datasets, whose root-mean-squared differences in the statistical quantities are shown as error bars, as detailed by Mori et al. (Reference Mori, Watanabe and Nagata2024). The virtual origin is determined by applying the same method used to the present DNS data. A least squares method is employed to determine the power-law exponents for
$ v_{rms}^2 \sim (x-x_0)^{n}$
,
$ L_v \sim (x-x_0)^{n_L}$
and
$C_\varepsilon \sim (x-x_0)^{n_C}$
, resulting in
$n=-1.70$
,
$ n_L=0.36$
and
$n_C=0.16$
. The streamwise increase in
$ C_\varepsilon$
is most likely attributable to non-equilibrium turbulence, as shown in the plot of
$(Re_\lambda , C_\varepsilon )$
by Mori et al. (Reference Mori, Watanabe and Nagata2024) and by the present DNS. With
$n_C=0.16$
, the theories of Batchelor and Saffman turbulence, (1.1)and (1.2), predict
$(n, n_L) = (-1.66, 0.33)$
and
$ (-1.39, 0.47)$
, respectively. These exponents indicate
$v_{rms}^2 L_v^5\sim (x-x_0)^{0.10}$
and
$v_{rms}^2 L_v^3\sim (x-x_0)^{-0.62}$
. As observed in the present DNS,
$v_{rms}^2 L_v^5$
weakly depends on the streamwise position, while
$v_{rms}^2 L_v^3$
continuously declines. Therefore, the decay of turbulence generated by jet interaction aligns with the Kolmogorov decay law of Batchelor turbulence.
These two experimental studies have provided measurement results for a limited range of streamwise locations, although they do not contradict the far-field decay observed in the present DNS. The decay properties are likely influenced by how turbulence is generated and depend on various parameters such as jet Reynolds and Mach numbers, as well as jet nozzle arrangements. The findings of this study highlight the need for further experimental investigations into jet interactions.
