MIXING IN SUPERSONIC TURBULENCE

We study mixing of pollutants in a turbulent flow v(x, t) with a density field ρ(x, t). The mixing process corresponds to the homogenization of the concentration field C(x, t), defined as the ratio of the pollutant density to the local flow density. The evolution of the concentration field is described by the advection–diffusion equation,

$$\begin{equation} \partial _t C + v_i \partial _i C = \frac{1}{\rho } \partial _i (\rho \kappa \partial _i C) + S(\boldsymbol{x},t), \end{equation} \tag{ 1 }$$

where κ denotes the molecular diffusivity, and the term S(x, t) represents the scalar sources. To understand the physics of turbulent mixing, it is crucial to recognize the role of the turbulent velocity. Intuitively, a velocity field redistributes the concentration field and changes its spatial configuration, but does not affect the mass fraction at a given concentration level, since the tracer particles simply follow the flow elements. This is true for any velocity field no matter how complex. Therefore, the velocity field itself does not mix, because true mixing means spatial smearing of pollutants, which can be caused only by molecular diffusion.

The observation that the velocity field is not responsible for true mixing can be formally confirmed from the equation for the variance of the concentration fluctuations, which is usually referred to as the scalar energy in the turbulence literature. To correctly reflect this intuition in compressible flows, one needs to consider the density-weighted variance, $\langle \tilde{\rho } C^2 \rangle,$ where the density-weighting factor $\tilde{\rho }$ is the ratio of ρ to the average flow density, $\bar{\rho }$, and 〈 · · · 〉 denotes an ensemble average, which is equal to the average over the flow domain for statistically homogeneous flows, as in our simulations.

The equation for the density-weighted variance can be derived from Equation (1) with the help of the continuity equation, yielding

$$\begin{equation} \partial _t \langle \tilde{\rho } C^2 \rangle + \partial _i \langle \tilde{\rho } C^2 v_i \rangle = 2 \langle \partial _i (\tilde{\rho } \kappa C \partial _i C) \rangle - 2 \langle \tilde{\rho } \kappa (\partial _i C)^2 \rangle + 2 \langle \tilde{\rho } S C \rangle. \end{equation} \tag{ 2 }$$

The second term on the left-hand side of this equation is an advection term, which corresponds to the spatial transport of concentration fluctuations between different regions. In the statistically homogeneous case, this term vanishes, and without it, Equation (2) does not have an explicit dependence on the velocity field. This proves the intuition that the velocity field does not truly mix. In fact, this is also true in an inhomogeneous flow. While the advection term is not zero in this case, it is a surface term, and thus its integral over the entire flow domain is zero. This means that the advection term conserves the global scalar variance.¹

On the right-hand side of Equation (2), the last term is the source term, which corresponds to the variance increase due to the injection of new pollutants. The other two terms on the right-hand side of this equation arise from the diffusion term in Equation (1). The first of these is again a surface term, which vanishes in the homogeneous case and conserves the global scale variance in the inhomogeneous case. The term $-2\langle \tilde{\rho } \kappa (\partial _i C)^2 \rangle$ represents scalar dissipation. This term is negative definite and thus monotonically reduces the scalar variance, corresponding to homogenization of the scalar fluctuations by molecular diffusion. We denote the dissipation rate as $\bar{\epsilon }_{\rm c} \equiv 2\langle \tilde{\rho } \kappa (\partial _i C)^2 \rangle$ and define a scalar dissipation timescale $\tau _{\rm c} \equiv \langle \tilde{\rho } C^2 \rangle /\bar{\epsilon }_{\rm c}$ to characterize it. The linear dependence on κ in the definition of $\bar{\epsilon }_{\rm c}$ may leave an impression that the scalar dissipation rate strongly depends on κ. However, that is not true because the scalar gradient in the definition also depends on κ. With decreasing κ, larger scalar gradients can develop, which may compensate the decrease in κ. As seen from the physical picture below and in Section 2.2, the dissipation rate and the mixing timescale are essentially independent of the molecular diffusivity.²

The molecular diffusivity κ is tiny in most practical environments, thus the dissipation rate is significant only if the scalar gradient is large, meaning that true mixing occurs only at very small scales. As mentioned earlier, at the injection scale of scalar sources, the molecular diffusion is usually very slow, and a significant mixing rate needs a velocity field to produce small-scale scalar structures, or equivalently, large concentration gradients. In an incompressible turbulent flow, the velocity field stretches, contracts, and folds the scalar field, leading to scalar structures at smaller and smaller scales. The size of the scalar structures keeps decreasing toward the diffusion scale where the molecular diffusion efficiently homogenizes the structures. Since κ is small, a wide scale separation usually exists between the scalar source size and the diffusion scale.

The diffusion scale is essentially the scale at which the diffusion term in Equation (1) becomes larger than the advection term, and thus can be derived by comparing these two terms accounting for their scale dependences. The estimate of the diffusion scale depends on the Schmidt number, Sc, defined as the ratio of viscosity to diffusivity, which determines the diffusion scale relative to the viscous dissipation scale of the velocity field (the scale at which kinetic energy is removed by viscosity). The derivation of the diffusion scale at different Sc values can be found in Monin & Yaglom (1975) and Lesieur (1997). For heavy elements in a neutral ISM, the diffusivity is expected to be smaller than the viscosity due to the larger cross section of heavy elements. However, for relatively light atoms like oxygen, the diffusivity is probably of the same order as the viscosity, and the Schmidt number is close to unity.

A second important dimensionless number for passive scalars is the Peclet number, Pe. It is defined as UL/κ where U and L are characteristic velocity and length scale of the turbulent flow. The Peclet number is related to the Reynolds number, Re, by Pe = ScRe. Therefore, for diffusing species with Sc ≈1, the Peclet number is close to the Reynolds number, ≳10⁶, in typical interstellar clouds.

In summary, a turbulent velocity field acts as a catalyst that enhances the mixing efficiency by feeding molecular diffusion with large-gradient structures. The mixing timescale, τ_c, is thus determined by the rate at which the turbulent velocity field brings the scalar structures down to the diffusion scale. It is essentially independent of κ for the case with Sc ≲1.

The classic picture for the generation of scalar structures at small scales in incompressible turbulence is a cascade similar to the energy cascade of the velocity field (e.g., Lesieur 1997). A constant flux of the scalar energy occurs along the cascade over the convective range, defined as the range of scales between the characteristic length scale of the scalar sources and the diffusion scale. In this range, the scalar fluctuations are regulated primarily by the advecting velocity. The flux feeds the dissipation at the diffusion scale and is thus expected to be equal to the average scalar dissipation rate. Assuming the timescale for a cascade step is of the order of the eddy turnover time at a given scale leads to the following scaling for the concentration difference, δC(l), over a separation l,

$$\begin{equation} \delta C(l)^2 \simeq \bar{\epsilon} _{\rm c} \frac{l}{\delta v(l)}, \end{equation} \tag{ 3 }$$

where δv(l) is the velocity difference over l and $\bar{\epsilon }_{\rm c}$ is the scalar dissipation rate.

In incompressible turbulent flows, we have the Kolmogorov scaling for the velocity difference, δv(l), in the inertial range, i.e., $\delta v(l) \simeq \bar{\epsilon} _{\rm v}^{1/3}l^{1/3}$, where $\bar{\epsilon} _{\rm v}$ is the average dissipation rate of kinetic energy. With this velocity scaling, Equation (3) gives $\delta C(l)^2 \simeq \bar{\epsilon }_{\rm v}^{-1/3} \bar{\epsilon }_{\rm c} l^{2/3}$ in the inertial-convective range (i.e., the intersection of the inertial range for the velocity field and the convective range for the scalar field). This scaling for the scalar difference, usually referred to as the Obukhov–Corrsin scaling law, corresponds to a k^−5/3 scalar spectrum. This spectrum has been confirmed for passive scalars in isotropic incompressible turbulence by both experiments and numerical simulations (Monin & Yaglom 1975; Sreenivasan 1996; Mydlarski & Warhaft 1998; Yeung et al. 2002; Watanabe & Gotoh 2004).

On the other hand, the scalar structure function and spectrum have not been studied in supersonic turbulence. Therefore, it is unknown whether Equation (3) originally proposed for passive scalar mixing in incompressible turbulence is valid in supersonic turbulence. The scaling of the velocity difference in supersonic turbulence has been shown to be significantly steeper than the Kolmogorov scaling (e.g., Kritsuk et al. 2007a). Thus, if Equation (3) holds for turbulent flows at any Mach number, one would expect the scalar structure functions to become flatter with increasing Mach number. We will measure the scaling exponents, ζ_v and ζ_c, of the second-order velocity structure functions $(\langle \delta v(l)^2 \rangle \propto l^{\zeta _{\rm v}})$ and the scalar structure functions $(\langle \delta C(l)^2 \rangle \propto l^{\zeta _{\rm c}})$ from our simulations at different Mach numbers. Equation (3) predicts the following relation between the two exponents,

$$\begin{equation} \zeta _{\rm c} \simeq 1- \zeta _{\rm v}/2, \end{equation} \tag{ 4 }$$

which will be tested against our simulation data.

Equation (3) also suggests that, if the scalar field and the flow are driven at similar scales, the dissipation timescales for scalar fluctuations and for kinetic energy will be similar. They are both determined by the turnover time of the largest eddies at the driving/source scale, which is of the order of the dynamical timescale, because the cascade becomes progressively faster at smaller scales, and the most time is spent at large scales. Numerical simulations confirmed that these two timescales are similar in incompressible turbulence (see, e.g., Donzis et al. 2005 and references therein) and that the scalar dissipation timescale is about half the energy dissipation timescale. Furthermore, in the supersonic case, the energy dissipation rate has been found to be similar to that in incompressible turbulence, i.e., close to the dynamical timescale (Stone et al. 1998; Mac Low 1999; Padoan & Norlund 1999), implying an energy cascade similar to that in incompressible turbulence. However, in the supersonic case, the dissipation timescale for forced scalars has never been computed.

Our simulations also give us an opportunity to explore the dissipation of kinetic energy in supersonic turbulence. The momentum equation in a driven turbulent flow is given by

$$\begin{equation} \partial _t v_i + v_j \partial _j v_i = - \frac{\partial _i p}{\rho } + \frac{1}{\rho } \partial _j (\sigma _{ij}) + f_i (\boldsymbol{x}, t), \end{equation} \tag{ 5 }$$

where p(x, t) is the pressure, f(x, t) is the driving force, and σ_ij is the viscous stress tensor. For an ideal gas, the bulk viscosity is negligible and the viscous tensor $\sigma _{ij} = \frac{1}{2} \rho \nu (\partial _i v_j +\partial _j v_i - \frac{2}{3} \partial _k v_k \delta _{ij})$ where ν is the kinematic viscosity.

From the momentum equation, Equation (5), and the continuity equation, we can derive an equation for the average kinetic energy per unit mass, $\langle \frac{1}{2} \tilde{\rho } v^2 \rangle$,

$$\begin{eqnarray} \partial _t \left\langle \frac{1}{2} \tilde{\rho } v^2 \right\rangle &=& \frac{1}{\bar{\rho }} \left\langle p \partial _i v_i \right\rangle - \frac{1}{2} \left\langle \tilde{\rho } \nu \left(\partial _i v_j +\partial _j v_i - \frac{2}{3} \partial _k v_k \delta _{ij} \right)^2 \right\rangle \nonumber\\ &&+ \langle \tilde{\rho } f_i v_i \rangle. \end{eqnarray} \tag{ 6 }$$

The three terms on the right-hand side of Equation (6) represent the pdV work, the viscous dissipation of kinetic energy, and the energy injection from the driving force, respectively.

The pdV term, which vanishes in incompressible flows, is a two-way energy exchange between kinetic energy and thermal energy. Although the conversion by pdV work between the two energy forms is reversible, the exchange rates in the two directions can be asymmetric. In fact, the term preferentially converts kinetic energy to heat due to an anticorrelation between the density (or pressure) and the velocity divergence. Consider, for example, an expanding region and a converging region with the same |∇ · v|. In the converging region, the density is typically larger, and thus the rate at which the pdV work converts kinetic energy to heat is faster than the conversion from thermal to kinetic energy in the corresponding expanding region. Therefore, in supersonic turbulence, the pdV work provides another mechanism for kinetic energy loss in addition to viscous dissipation. We denote the rate of kinetic energy loss by this mechanism as $\bar{\epsilon }_{\rm p} \equiv - \langle \frac{p}{\bar{\rho }} \partial _i v_i \rangle$, and point out that this effect has not been explicitly considered in previous studies.

Viscous dissipation, corresponding to the second term on the right-hand side of Equation (6), is an irreversible conversion of kinetic energy to heat. If ρν is constant, which is true if the flow temperature is roughly constant, the overall average kinetic energy dissipation rate in the flow can be written as $\bar{\epsilon }_{\nu } = \tilde{\rho }\nu [\langle (\nabla \times \boldsymbol{v})^2 \rangle + \langle \frac{4}{3} (\nabla \cdot \boldsymbol{v})^2 \rangle ]$ using integration by parts. We define a viscous dissipation timescale, τ_diss, as $\frac{1}{2} \langle \tilde{\rho } v^2 \rangle /\bar{\epsilon }_{\nu },$ and point out that this timescale for kinetic energy loss purely by viscous dissipation, to our knowledge, has also never been evaluated in supersonic turbulence.

In fact, the only timescale for kinetic energy loss estimated in earlier studies is the overall dissipation timescale τ_v that characterizes the total rate for the kinetic energy loss by both pdV work and viscous dissipation, $\tau _{\rm v} = \frac{1}{2} \langle \tilde{\rho } v^2 \rangle /(\bar{\epsilon }_{\nu } + \bar{\epsilon }_{\rm p})$ (Stone et al. 1998; Mac Low 1999; Padoan & Norlund 1999; Lemaster & Stone 2009). For example, the timescale for kinetic energy loss, t_diss defined in Lemaster & Stone (2009) corresponds to the timescale τ_v defined here. Below we will evaluate both τ_diss and τ_v, separating out the relative importance of kinetic energy loss by pdV work and viscous dissipation as a function of Mach number.

To drive turbulence in our simulation, we made use of the FLASH3 "Stir" package (Benzi et al. 2008), which we heavily modified for our purposes. The flow velocity was driven by the forcing term, f(x, t) in Equation (5). The forcing amplitude was adjusted to achieve different Mach numbers. In our 512³ simulations, we drove the turbulence with solenoidal modes only, i.e., ∇ · f = 0 and took f to be a Gaussian random vector with an exponential temporal correlation. A finite correlation timescale is adopted in the forcing scheme (e.g., Padoan & Nordlund 1999). This is different from studies using an infinitesimal correlation timescale with independent forcing at each time step (e.g., Lemaster & Stone 2009) or an infinite correlation timescale with a fixed driving force (e.g., Kritsuk et al. 2007a). The driving scheme can be summarized as $\langle f_i(\boldsymbol{k}, t)f_j(\boldsymbol{k}, t^{\prime }) \rangle = \mathcal {P}_{\rm f}(k) (\delta _{ij} - \frac{k_i k_j}{k^2}) \exp (-\frac{t-t^{\prime }}{t_{\rm f}}),$ where t_f is the forcing correlation timescale. The forcing wave numbers are chosen to be in the range 1 ⩽ k/2π ⩽ 3 in units in which the domain size is 1, and each independent wave vector in this range is given the same forcing energy per unit time, which results in a quite flat forcing spectrum. This pattern was used in all the simulated flows, so that the shape of the forcing spectrum was identical for all Mach numbers. We also tried several 256³ simulations with a different driving scheme where we set 1/3 of the forcing energy to be in compressible modes. The driving force for these flows was otherwise the same as for the 512³ simulations.

To characterize the large-scale properties of the flow, we define a forcing length scale, L_f, from the forcing spectrum, $L_{\rm f} = \int \!\frac{2 \pi }{k} \mathcal {P}_{\rm f}(k) d\boldsymbol{k}/\!\int\!\mathcal {P}_{\rm f}(k) d\boldsymbol{k}$. For the forcing spectrum chosen here, we find L_f = 0.46 in units in which the domain size is 1. Lemaster & Stone (2009) adopted a forcing spectrum, $\mathcal {P}_{\rm f}(k) \propto k^6 \exp (- 8 k /k_{\rm pk})$, that strongly peaks at a wave number k_pk. Inserting this spectrum into our definition for L_f, we obtain that L_f is equal to 2π/k_pk, which is exactly the same as the characteristic flow length scale used in that study. For each simulation in our study, we define a dynamical time, τ_dyn ≡ L_f/v_rms, where v_rms is the three-dimensional density-weighted rms velocity.

We simulated turbulent flows at six different Mach numbers ranging from transonic to highly supersonic. At each Mach number, the simulation covered more than 10 dynamical times. After an initial phase of steady increase, the rms velocity saturated in a few dynamical times, and a statistically steady state was reached. In the steady state, the rms velocity was not exactly constant, but had small variations with time, with an (rms) amplitude of about 2%–5%, depending on the Mach number.

Following Lemaster & Stone (2009), we used the density-weighted (three-dimensional) rms velocity in our definition for the Mach number, i.e., $M= \langle \tilde{\rho } v^2 \rangle ^{1/2}/c_{\rm s}$ where c_s is the sound speed. For each flow, we averaged the density-weighted rms velocity over about 10 dynamical times, and the average Mach numbers in the six flows were 0.9, 1.4, 2.1, 3.0, 4.6, and 6.1. If the volume-weighted rms velocity was used, the Mach number was only slightly larger (by about 7% except for the M = 0.9 case where the two are essentially equal). Because all the flows were driven with the same forcing pattern, L_f was the same at all Mach numbers, and the dynamical timescale τ_dyn = L_f/v_rms was proportional to 1/M.

We studied both forced scalars and decaying scalars. For the forced case, we continuously added pollutants in the flow by the source term S(x, t) in the advection–diffusion equation (Equation (1)). Similar to the velocity driving, S was taken to be a Gaussian variable³ with a spectrum given by $\langle S(\boldsymbol{k}, t) S(\boldsymbol{k}, t^{\prime }) \rangle = \mathcal {P}_{\rm s}(k) \exp (-\frac{t-t^{\prime }}{t_{\rm s}})$. We considered several different forcing schemes for $\mathcal {P}_{\rm s}(k)$. In the first scheme, $\mathcal {P}_{\rm s}(k)$ was chosen to have the same shape as the velocity forcing spectrum, i.e., $\mathcal {P}_{\rm s}(k) \propto \mathcal {P}_{\rm f}(k)$. With exactly the same forcing pattern as the velocity, this scheme allows for a direct comparison of the dissipation timescale of the scalar variance to that of the kinetic energy. As discussed below, this reveals important differences in the physics of scalar and kinetic energy dissipation, especially concerning the role of compressible modes and shocks.

Like the rms velocity, the scalar variance exhibited temporal variations in the steady state. In order to reduce their amplitude, we used three independent scalars, each of which has the same source spectrum with identical scalar and flow forcing. Averaging over the three independent scalars resulted in a much smaller amplitude (in the range of 2%–5%) in the temporal variations of the rms concentration than that of each individual scalar.

To test the dependence of the mixing timescales on the characteristic length scale of the scalar sources, we considered three other schemes in which the scalars were forced at larger wave numbers. In these cases, we chose the source wave numbers to be 3.5 ⩽ k/2π ⩽ 4.5,7.5 ⩽ k/2π ⩽ 8.5, and 16.2 ⩽ k/2π ⩽ 16.5, respectively. Again, each independent mode in these ranges was given the same forcing power. Similar to the length scale, L_f, of the flow driving force, we defined a source scale, $L_{\rm s} = \int \!\frac{2 \pi }{k} \mathcal {P}_{\rm s}(k) d\boldsymbol{k}/\!\int\!\mathcal{P}_{\rm s}(k) d\boldsymbol{k}$. The source length scales are 0.25, 0.124, and 0.061 the size of the computation box for the three cases, respectively.

In our simulations, we take the scalar source term to have zero mean and do not consider the evolution of the mean concentration, which is simply determined by the rate at which pollutants are injected. By definition, the concentration is positive and bounded in the range 0 ⩽ C ⩽ 1. The negative concentration values in our simulations should be understood as relative to the mean concentration. The amplitude of the scalar fluctuations in our simulations is arbitrary (due to the arbitrarily chosen amplitude for the source term). This does not affect the study of mixing physics, because the advection–diffusion equation is linear with the concentration.

In all the simulations, we took both t_f and t_s to be 1/4 of the sound crossing timescale, defined as the box size divided by the sound speed. We also tried different values of t_f and t_s (but with a resolution of 256³), and found that the results were not affected by the specific choice of the forcing/source correlation timescales.

For the decaying case, we set an initial configuration for the scalar field when the turbulent flows were well developed, and let it evolve with no continuing sources, i.e., subsequently setting S(x, t) to zero. For each Mach number, we studied four decaying scalars with different initial spectra. The initial scalar spectra were taken to be the same as the source spectra chosen for forced cases, i.e., the wave number ranges were set to be 1 ⩽ k/2π ⩽ 3, 3.5 ⩽ k/2π ⩽ 4.5,7.5 ⩽ k/2π ⩽ 8.5, and 16.2 ⩽ k/2π ⩽ 16.5, respectively.

Although the main motivation for considering scalars forced at different scales is theoretical for a full understanding of mixing, they are also of practical interest. For example, if the main energy source for the ISM turbulence is SN explosions, then, for mixing of new metals from SNe, the source length scale would be same as the flow driving scale. On the other hand, the scalar injection scale would be different from the flow forcing scale if the pollutant sources and the flow energy sources are different. An example is self-enrichment in star-forming clouds, where the scalar source length scale may be smaller than the flow driving scale if the turbulent energy is mainly from the cascade of the interstellar turbulence at larger scales. The decaying case may be applicable to mixing in galaxies with an episodic star formation history: the mixing process for elements produced by core-collapse SNe in between star bursts in these galaxies would be a decaying case.

Finally, we point out that neither the viscous term, $\frac{1}{\rho }\partial _j \sigma _{ij}$, nor the diffusion term, $\frac{1}{\rho }\partial _i (\rho \kappa \partial _i C)$, was explicitly included in our simulations. Rather, kinetic energy dissipation and scalar dissipation were both controlled by numerical diffusion. Therefore, both the viscous dissipation scale and the diffusion scale are expected to be around the resolution scale. The effective viscosity and diffusivity are thus similar, with the effective Schimdt number, Sc ≈1, and the Peclet number about equal to the Reynolds number. For the case with identical scalar and flow forcing, the convective scale range coincides with the inertial range, while the scalars forced at larger wave numbers has a smaller convective range than the inertial range.

To calculate the scalar dissipation timescale, τ_c, we need the scalar dissipation rate, _c, which could not be directly computed in our simulations because the dissipation was controlled by numerical diffusion. However, in the case of forced scalars, after the scalar fluctuations reach a statistically stationary state, we have a balance between the dissipation term and the source term in Equation (2), i.e., $\epsilon _{\rm c} = \epsilon _{\rm s}\equiv \langle \tilde{\rho } S C \rangle$. The computation of the source term is straightforward, and we calculated the scalar dissipation timescale as $\tau _{\rm c} =\langle \tilde{\rho } C^2 \rangle /\langle \tilde{\rho } S C \rangle$. For the scheme with identical scalar and flow forcing, we used the average of the source term over the three independent scalars.

Similarly, the viscous dissipation of kinetic energy was obtained using energy balance, $\bar{\epsilon }_{\rm f} -\bar{\epsilon }_{\rm p} - \bar{\epsilon }_{\nu } =0$, in the statistically stationary state. The energy injection by the driving force per unit time, $\bar{\epsilon }_{\rm f}$, and the contribution from pdV work, $\bar{\epsilon }_{\rm p}$ were calculated by $\langle \tilde{\rho } f_i v_i \rangle$ and $\bar{\rho }^{-1} \langle p \partial _i v_i \rangle$ in Equation (6), respectively. In the stationary state, this allowed us to evaluate the viscous dissipation timescale as $\tau _{\rm diss} = \langle \frac{1}{2}\tilde{\rho } v^2 \rangle / (\bar{\epsilon }_{\rm f} -\bar{\epsilon }_{\rm p})$, and the timescale for the overall kinetic energy loss including the contribution from the pdV work as $\tau _{\rm v} = \langle \frac{1}{2}\tilde{\rho } v^2 \rangle /\bar{\epsilon }_{\rm f}$.

The kinetic energy loss by both the pdV work and the viscous dissipation was stored as thermal energy. We point out that the isothermal equation of state adopted in our simulations was used to imitate a roughly constant temperature due to the effect of efficient radiative cooling in interstellar clouds. In nature, the heat converted from kinetic energy would be radiated away.

We find that, if the forcing correlation timescale, t_f, was much smaller than the typical timescale for the kinetic energy loss to thermal energy, τ_v, the energy injection rate goes like $\bar{\epsilon }_{\rm f} \simeq t_{\rm f} \int \mathcal {P}_{\rm f}(k) d{\boldsymbol{k}}$. On the other hand, if t_f is larger than τ_v, as is the case for our choice of τ_f, $\bar{\epsilon }_{\rm f} \propto \tau _{\rm {v}} \int \mathcal {P}_f(k) d{\boldsymbol{k}}$ because only the forcing energy within a timescale τ_v remains as kinetic energy in the flow. The same applies to the variance input from scalar sources, i.e., $\bar{\epsilon }_{\rm s} \simeq t_{\rm s} \int \mathcal {P}_{\rm s} (k) d{\boldsymbol{k}}$ if t_s ≪ τ_c or $\simeq \tau _{\rm c} \int \mathcal {P}_{\rm s} (k) d{\boldsymbol{k}}$ for t_s ≫ τ_c.

Similar to the velocity and scalar variances, the other quantities needed for the calculation of timescales also show temporal variations to different degrees. These quantities include the kinetic energy injection by the driving force, the pdV work, and the scalar variance input by the sources. In our calculations, we used the temporal average of all these quantities over about 10 dynamical times. Comparing with results from 256³ simulations for the corresponding flows, we find that the timescales have converged at our resolution of 512³ zones (see also Lemaster & Stone 2009).

In the left panels of Figure 1, we show the X-component of the velocity field on a slice (X–Y plane) of the simulation grid at one snapshot for two Mach numbers, M = 0.9 and M = 6.1. The snapshots are taken at t = 8τ_dyn and 5.6τ_dyn for M = 0.9 and M = 6.1, respectively. The concentration field shown here is from a scalar forced at the flow driving scale. The velocity field exhibits extremely different features in these two limits. In the M = 0.9 case, the flow contains a large number of vortices, which persist down to small scales. In the M = 6.1 case, on the other hand, the field appears to be dominated by large-scale expansions and compressions, and the most prominent structures are shock fronts. However, we point out that, although not clearly visible in Figure 1, a significant fraction of the kinetic energy is contained in solenoidal modes even at M = 6.1. For example, strong shears and vortices are usually found behind narrow postshock regions. As we will show later, the solenoidal modes play important roles in both energy dissipation and mixing at all Mach numbers.

Zoom In Zoom Out Reset image size

In Table 1, we list the overall properties of the flow as a function of Mach number. In the second column, we give the ratio of the density-weighted divergence variance and vorticity variance, $\langle \tilde{\rho } (\nabla \cdot \boldsymbol{v})^2 \rangle /\langle \tilde{\rho } (\nabla \times \boldsymbol{v})^2 \rangle$. This ratio characterizes the contribution of compressible modes to viscous dissipation relative to solenoidal modes, and it is smaller than 1 even for M = 6.1, suggesting that the majority of kinetic energy is dissipated by solenoidal modes even in highly supersonic flows (Kritsuk et al. 2007a). However, we point out that this ratio is only a rough estimate for the relative contribution from the two modes, as it is not clear how the kinetic energy is dissipated by numerical diffusion by the high-order PPM solver. In the third column, we also give the same ratios, but with no density weighting. For both cases, the ratio increases with the Mach number at small M. The density-weighted ratio saturates at M ≈ 3, while the volume-weighted one saturates at a little larger M, and the density-weighted ratios are all smaller because of an anticorrelation between the density and divergence.

Table 1. Flow Properties at Different Mach Numbers

M	$\langle \tilde{\rho } (\nabla \cdot \boldsymbol{v})^2 \rangle / \langle \tilde{\rho } (\nabla \times \boldsymbol{v})^2 \rangle$	〈(∇ · v)2〉/〈(∇ × v)2〉	fcomp	fcomp(k/2π ⩾ 4)	fpdV
0.9	0.06	0.06	0.05	0.07	0.15
1.4	0.19	0.20	0.10	0.16	0.30
2.1	0.41	0.49	0.15	0.25	0.35
3.0	0.52	0.75	0.17	0.30	0.30
4.6	0.55	0.88	0.18	0.33	0.23
6.1	0.54	0.86	0.19	0.32	0.17

Download table as:  ASCIITypeset image

The fourth column lists the fraction of kinetic energy in compressible modes. We computed this fraction by decomposing the energy spectrum into a solenoidal part and a potential part and comparing the kinetic energy contained in each part. Although the flow is driven by solenoidal modes only, a finite fraction of kinetic energy appears in the compressible modes, which are converted from the solenoidal modes. The fraction increases from 5% at M = 0.9 to 20% at M = 6.1 in our simulations. We note that these values are not universal, but depend on the ratio of the compressible modes to the solenoidal modes in the driving term, as discussed below.

We are particularly interested in the velocity fluctuations in the inertial range, which regulate the scalar cascade to the diffusion scale. For scalars forced at the flow driving scale, the cascade starts at the same scale as the energy cascade, i.e., at k/2π ≈ 4. Thus, we calculated the kinetic energy contained in the compressible modes and the solenoidal modes in the inertial range by integrating the potential and solenoidal spectra over all wave numbers k/2π ⩾ 4, which also includes the negligible contribution from the dissipation range. The fifth column in Table 1 gives the fraction of kinetic energy in compressible modes at k/2π ⩾ 4, which increases at small M and saturates at about 1/3 for M ≳ 3. The fraction of 1/3 implies an energy equipartition between the solenoidal modes and the potential modes in the inertial range, which is reached at M ≈ 3. We find the slope of the compressible spectrum is close to that of the solenoidal spectrum in the inertial range, meaning that a equipartition is established at each wave number in the inertial range for M ≈ 3. We find this equipartition is universal regardless the compressible fraction in the flow driving (see below).

Note that the fractions given in Columns 4 and 5 do not account for density weighting, which cannot be perfectly applied to the fraction of kinetic energy in compressible modes at large M. Although the velocity field can be decomposed into a solenoidal part, v_s, and a potential part, v_p, the density-weighed ratio $\langle \tilde{\rho } v_{\rm p}^2 \rangle /(\langle \tilde{\rho } v_{\rm p}^2 \rangle + \langle \tilde{\rho } v_{\rm s}^2 \rangle)$ cannot be interpreted as the energy fraction at high Mach numbers. This is because the denominator is not equal to the total energy, which is given by $\langle \tilde{\rho } v_{\rm p}^2 \rangle + 2\langle \tilde{\rho } \boldsymbol{v}_{\rm p} \cdot \boldsymbol{v}_{\rm s} \rangle + \langle \tilde{\rho } v_{\rm s}^2 \rangle$. In fact, we find that the density-weighted correlation term $\langle \tilde{\rho } \boldsymbol{v}_{\rm p} \cdot \boldsymbol{v}_{\rm s} \rangle$ is not zero. Rather, it is negative, decreases with increasing M, and becomes significant at M ≳ 2–3. Therefore, density weighting for the fraction of energy in the compressible modes cannot be well defined at large Mach numbers. However, there is no such problem for the volume-weighted fraction because it can be shown that 〈v_p · v_s〉 is exactly zero in the homogeneous case, and we thus only show the volume-weighted fraction in Table 1.

The last column in Table 1 lists the ratio of the kinetic energy loss by pdV work to the total kinetic energy loss. This will be discussed in detail in Section 4.3.1.

To test whether and how the flow properties listed in Table 1 depend on the flow driving, we performed a set of (256³) simulations covering a similar range of Mach numbers but with 1/3 forcing energy contained in the compressible modes. We found that, with this different driving scheme, most quantities listed in Table 1 remain essentially unchanged, except the overall fraction of energy in the compressible modes (Column 4). This is expected since the overall compressible fraction depends on the driving scales, which contain most of the kinetic energy, while the rest of the quantities depend on the statistics at smaller scales, which should be driving independent.

The change in the overall compressible fraction is slight for M ≲ 2, only a few percent larger than the case with solenoidal driving, but there is a significant increase in the fraction at M ≳ 3. With the compressible driving scheme, the fraction saturates at ≈24% for M ≳ 3, considerably larger than ≈ 18% given in Table 1. We note that the overall compressible fraction in the flow turns out to be smaller than that (1/3) in the driving force at all Mach numbers, indicating a preferential conversion from compressible modes to solenoidal modes at the driving scales. This is different from the case with solenoidal driving where a fraction of kinetic energy in solenoidal modes is converted to compressible modes.

Of particular interest is that, in the inertial range, the equipartition between solenoidal and compressible modes is also found in the M ≳ 3 flows that are forced with 1/3 driving energy in compressible modes, suggesting that this equipartition at high Mach numbers is universal. On the other hand, no equipartition between compressible and solenoidal modes exists at the driving scales even in the case where the compressible fraction in the driving force is set to the equipartition value of 1/3.

In the right panels of Figure 1, we show the concentration on slices taken from simulations with Mach numbers M = 0.9 and M = 6.1. As was the case for the velocity field, there are large differences between the two runs. At M = 0.9, the scalar field shows numerous small-scale structures, which are clearly produced by stretching and shear in vortices. There are many structures with sharp concentration contrasts, which are usually referred to as "cliffs and ramps" in the literature for passive scalar in incompressible turbulence (e.g., Shraiman & Siggia 2000). These cliffs and ramps are produced by turbulent stretching, which rearranges the scalar gradients that initially exist only at large scales and brings fluid elements with very different concentration levels next to each other. The "cliffs and ramps" have large scalar gradients, and are thus strongly dissipative.

On the other hand, the scalar field at M = 6.1 primarily reflects expansions and compressions in the flow, similar to the velocity field at this Mach number. Also like the velocity case, the visual impression of the figure is biased toward compressible modes, and does not sufficiently reflect the effect of solenoidal modes on the scalar field, which we will show is dominant for mixing.

An interesting property of this scalar field is the existence of scalar "edges" that coincide with the locations of velocity shocks. Unlike the cliff and ramp structures in the M = 0.9 case, these edges are produced not by stretching, but rather by compression by shocks, which squeeze the concentration profile and bring spatially well-separated regions with different concentrations close to each other. In this way, the length scale of the concentration profile is reduced by a factor about equal to the density jump across the shock, and the concentration gradient is amplified by the same factor.

Consider for example a shock with a density jump factor of 20 in our simulation. Squeezing by this shock could bring two fluid elements in the preshock region with a distance of 20 times resolution scales to a separation of one computation cell. If the two elements have quite different concentrations, then, when they pass the shock, one would see a significant concentration change across the shock front, i.e., a scalar edge. However, if the fluid elements in the preshock region have the same concentration, i.e., if the preshock scalar profile is uniform, no scalar edge would form. This explains why scalar edges are not found at every velocity shock front.

Gradient amplification is not limited to shock fronts, but instead occurs in any region that has been compressed. For example, in postshock regions, the scalar field is also squeezed by compression, and thus large gradients are expected behind shocks. In fact, in Figure 1, we see there is usually a second sharp concentration change right behind a scalar edge, which may correspond to the amplified gradient in the postshock region. However, extended regions with strong gradient amplification are not observed behind the edges, perhaps because postshock regions are narrow. Also these regions are subject to rapid mixing because the scalar gradients there are amplified by both compression and strong shears and vortices, which usually exist behind shocks.

As compression events give rise to larger scalar gradients, they increase the mixing efficiency, which depends on the generation of large-gradient structures. On the other hand, the compressible modes also include expansions, which reduce the scalar gradients, and thus decrease the mixing rate. Thus, it is not clear if the presence of compressible modes leads to a net increase of the mixing efficiency. This will be addressed in more detail in Section 5, where we also compare the efficiency of the compressible modes relative to the solenoidal modes in amplifying the scalar gradients.

Note that there is a fundamental difference between scalar edges and velocity shocks. Although the scalar edges in Figure 1 appear to be discontinuous, they are actually continuous. This can be shown by considering the concentration change across a shock front in the limit of infinitesimal shock thickness. From the formation mechanism of the scalar edges, the concentration change over a shock front can be roughly estimated by the product of three factors: the preshock gradient, the density jump factor (i.e., the squeezing factor), and the shock thickness. Therefore, when the shock thickness approaches zero (corresponding to decreasing viscosity), the concentration change across the shock front would also approach zero, since the density jump factor would remain constant and finite as the shock width decreases to infinitesimal. This means that the concentration is essentially continuous, unlike the velocity shocks, which are exact discontinuities in the limit of infinitesimal viscosity. This argument is confirmed by simulations for one-dimensional supersonic flows at very high resolutions.

The continuity of the scalar field across a shock of infinitesimal width is also expected from the conservation of the tracer mass. Mass conservation implies that the tracer density has the same jump as the flow density across a shock, and thus the concentration, which is the ratio of the two densities, would be continuous across the front.

The reason the edges appear to be discontinuous in our simulation is due to the limited resolution. A finite and significant shock width in our simulations means that the materials across the shock front could come from well-separated regions (especially for strong shocks), making a considerable scalar change across the shock front possible. With increasing resolution, we would expect the scalar structures around the shocks appear to be more continuous, while the shocks would clearly remain discontinuous. If a shock had an infinitesimal width, its effect on the scalar field is just to amplify the scalar gradient in the postshock region.

The difference between the scalar edges and the velocity shocks has an interesting implication. Every shock is intrinsically a strong dissipative structure for kinetic energy, but it is not true that every shock produces a strong dissipative scalar structure, considering that a scalar edge is essentially not a discontinuity. In fact, most shocks cannot produce scalar structures that dissipate scalar fluctuations at the same level as shocks for energy dissipation. To produce such a structure, a shock needs to reduce the length scale of a scalar structure close to the diffusion scale. Whether a shock can give rise to this scale reduction depends on the shock strength, which determines the squeezing factor across the shock. We will show that strong shocks capable of such a scale reduction are very rare in Section 5, and this point is responsible for the difference in the roles of shocks in kinetic energy dissipation and in scalar dissipation.

4.3.1. Kinetic Energy Dissipation

In Figure 2, we plot the timescales defined in Sections 3.1 and 3.2 as a function of Mach number, normalized to the dynamical time, τ_dyn. The behavior of the normalized timescales as a function of M is expected to essentially reflect the contribution from compressible modes relative to solenoidal modes.

Zoom In Zoom Out Reset image size

The dotted line in Figure 2 shows the normalized timescale for the kinetic energy loss by viscous dissipation, τ_diss/τ_dyn. This timescale steadily decreases as M increases from 0.9 to 3, and shocks become stronger and more frequent. Shocks are strong dissipative structures for kinetic energy, and thus compressible modes provide an extra fast channel for the viscous dissipation of kinetic energy, which does not exist in incompressible flows. The steady decrease in the τ_diss/τ_dyn curve with M for M ≲ 3 parallels the increase in the ratio of the divergence variance to the vorticity variance listed in Table 1 with and without density weighting. The ratio becomes roughly constant for M ≳ 3, corresponding to the slower decrease of τ_diss/τ_dyn in that range of M. This suggests that the faster viscous dissipation at larger M may be primarily due to a larger contribution from the compressible modes. Overall the density-weighted divergence to vorticity variance ratio increases by about 50% as the Mach number increases from M = 0.9 to 6.1, generally consistent with the decrease in the viscous dissipation timescale by a factor of 1.8.

The dashed curve in Figure 2 is τ_v/τ_dyn, the normalized timescale for the overall kinetic energy loss by both pdV work and viscous dissipation. This timescale is significantly smaller than τ_diss because the pdV work gives a considerable contribution to the kinetic energy loss in the chosen range of M.

The fraction of kinetic energy that is converted to heat by pdV work is given in the sixth column of Table 1, which shows that the fraction peaks at the intermediate Mach number, M ≃ 2. At low Mach numbers, increasing M allows more compressible modes to appear and larger |∇ · v| to be available, leading to a rapid increase in the rate of pdV work. On the other hand, the importance of pdV work in converting kinetic energy into heat decreases with M for M ≳ 2. At these large Mach numbers, the amplitude of the divergence scales linearly with M, corresponding to the saturation of the fraction of energy in compressible modes (see Table 1) and the increase in the amplitude of pressure fluctuations is slower than ∝M. Therefore, the conversion rate by pdV work does not increase faster than ∝M². Since the viscous dissipation rate increases as ∝v³_rms, this means that the fraction of energy loss by pdV work decreases at large M.

This trend is also expected from another perspective. The viscous dissipation rate by compressible modes is a second-order function of the velocity divergence and increases with M faster than the pdV work, which scales linearly with the divergence. Thus, as M increases above 2, the fraction of kinetic energy loss by pdV work starts to decrease as the viscous dissipation rate by compressible modes becomes a significant fraction of the total viscous dissipation rate.

The increase in the fraction of kinetic energy loss by pdV work for M ≲ 2 makes the τ_v/τ_dyn curve a little steeper than the curve for τ_diss/τ_dyn, the normalized timescale for viscous dissipation only. On the other hand, the τ_v/τ_dyn curve is flatter than the τ_diss/τ_dyn curve for M ≳ 3 due to the decrease in the pdV work fraction, which happens to cancel the weak increase in the viscous dissipation rate. The overall kinetic energy loss rate in highly supersonic flows (M > 3) is faster than in the transonic case with M = 0.9 also by a factor of ≈1.8.

Our simulations with 1/3 forcing energy in compressible modes give essentially the same results for the energy dissipation timescales as those found from pure solenoidal driving. Clearly, the energy dissipation occurs at small scales and basically depends on the cascade in the inertial range. As discussed earlier, the statistical measures at small scales listed in Table 1 (i.e., except the fourth column) are the same in the two different driving schemes, and thus the timescales for kinetic energy dissipation are expected to be also driving independent.

We note that, although our forcing scheme is very different from that in Lemaster & Stone (2009), the timescale τ_v from our simulations is in good agreement with the results given in their Table 1 (the column for t_diss/t_f).

4.3.2. Scalar Dissipation

4.3.3. Scale Dependence of the Mixing Timescale

In Figure 3, we show the dependence of τ_c on the characteristic length scale of scalar sources. The three curves are for three scalars forced in different wave number ranges: 3.5 ⩽ k/2π ⩽ 4.5, 7.5 ⩽ k/2π ⩽ 8.5, and 16.2 ⩽ k/2π ⩽ 16.5. The three mixing timescales are normalized to that of the scalar forced in the same wave number range as the flow driving. The corresponding mixing timescales normalized to the dynamical time are given in Table 2.

Zoom In Zoom Out Reset image size

Table 2. Energy Dissipation and Mixing Timescales

M	τdiss/τdyn	τv/τdyn	τc/τdyn
			1 ⩽ k/2π ⩽ 3	3.5 ⩽ k/2π ⩽ 4.5	7.5 ⩽ k/2π ⩽ 8.5	16.2 ⩽ k/2π ⩽ 16.5
0.9	1.84	1.56	0.85	0.57	0.37	0.23
1.4	1.76	1.23	0.91	0.63	0.41	0.27
2.1	1.55	1.01	1.00	0.73	0.50	0.33
3.0	1.26	0.88	1.03	0.77	0.55	0.37
4.6	1.10	0.85	0.99	0.75	0.55	0.37
6.1	1.03	0.86	0.97	0.74	0.54	0.38

Download table as:  ASCIITypeset image

The results in Figure 3 suggest that the mixing timescale is proportional to the turnover time of eddies at the length scale of the scalar sources. At all Mach numbers, the mixing timescales decrease with decreasing source length scale. This is because the turnover timescale in smaller eddies is shorter, and thus, starting from a smaller source scale, the cascade to the diffusion scale is faster. Figure 3 also shows that, with increasing Mach numbers, the scale dependence becomes weaker, which corresponds to the steepening of the velocity structure functions (see Section 4.4). At larger Mach number, the velocity field has relatively less power at small scales, and thus the increase in the scalar cascade toward small scales becomes slower, because the decrease in the eddy turnover time is slower.

Quantitatively, we find that the scaling of mixing timescales with the source size is consistent with that of the eddy turnover time with the eddy size. We consider the two forced scalars with source length scales L_s = 0.124 and 0.061, which lie in the inertial range. Roughly fitting the mixing timescales for these two scalars with a power law, τ_c ∝ L^α_s, gives α = 0.64, 0.62, 0.60, 0.57, 0.54, 0.52 for Mach numbers from 0.9 to 6.1. As we will see below, these numbers are quite close to the scaling exponents of the turnover time, l/δv(l), as a function of the eddy size l, using the inertial-range scaling of δv(l). We thus conclude that the mixing timescale is essentially given by the turnover time of eddies at the scalar injection scale.

In Figure 4, we show the second-order structure functions for the velocity and scalar fluctuations. Note that these results are from volume-weighted averages, although using density weighting would be more consistent with the dissipation timescales for the density-weighted velocity and scalar variances studied above. However, unlike one-point statistics, the choice of an appropriate density-weighting scheme is not trivial for two-point statistical measures, and will be pursued in future studies. Nevertheless, the volume-weighed structure functions and spectra shown here uncover interesting differences between the velocity and scalar structures.

Zoom In Zoom Out Reset image size

The curves in the left panel of Figure 4 are the longitudinal velocity structure functions, $S_{\rm ll} (l) \equiv \langle [(\boldsymbol{v}(\boldsymbol{x} +\boldsymbol{l})- \boldsymbol{v}(\boldsymbol{x})) \cdot \frac{\boldsymbol{l}}{l}]^2\rangle$, at different Mach numbers. The separation l is in units of the resolution scale. At M = 0.9, an extended inertial range (about a decade) exists already at the resolution of 512³. In this range, the structure function is well fit by a power law with a slope of 0.67, consistent with the Kolmogorov scaling for incompressible turbulence.

At larger Mach numbers, the structure function is significantly steeper, as found in previous studies (e.g., Kritsuk et al. 2007a). We find that the scale range with good power-law fits decreases with M, such that at M = 6.1 there is barely any range where the structure function is well fit by a straight line in the log plot. This decrease of the power-law range may be caused by stronger numerical viscosity in the FLASH code to stabilize shocks at larger Mach numbers (see discussions below). To estimate the best-fit slopes of the structure functions at high Mach numbers, we choose a small range of l from 20 to 70, far from both the dissipation range and the driving scale. In this range, the exponents, ζ_v, are 0.67, 0.73, 0.83, 0.89, 0.92, and 0.95 for Mach numbers from 0.9 to 0.61.

The structure functions steepen with M for two reasons. First, the number of shocks and the average shock intensity both increase with increasing Mach number, and shocks tend to give rise to a structure function that scales linearly with separation. This is why at large Mach numbers the scaling exponent approaches 1, which is sometimes referred to as the Burgers scaling (e.g., Kritsuk et al. 2007a). Second, as discussed earlier, the pressure term preferentially converts kinetic energy into thermal energy. Thus along the energy cascade some kinetic energy may be lost by the pdV work, resulting in smaller velocities toward small scales, and thus a steeper structure function. This effect may peak at M ≈ 2 where the fraction of kinetic energy loss by pdV work is maximum. These two reasons are interrelated as pdV work converts kinetic energy to heat across shock fronts.

The second-order structure function for the scalar field S_c(l) ≡ 〈(C(x + l) − C(x))²〉 is shown in the right panel of Figure 4. Its behavior as a function of the Mach number is very different from the velocity structure function. The scalar structure function first flattens as the Mach number increases from 0.9 to 2. In this range of Mach numbers, we find an extended inertial-convective range for the scalar structure function, and the scaling exponents can be accurately measured. The scaling exponent ζ_c is 0.67 at M = 0.9, in agreement with the Obukohov–Corrsin scaling, and it decreases to 0.62 and 0.59 at M = 1.4 and 2.1, respectively.

These values are in good agreement with the prediction from the cascade picture, represented by Equation (4). As found earlier, the scaling exponents for the velocity structure function are ζ_v = 0.73 and 0.83 at M = 1.4 and M = 2.1. Therefore, Equation (4) would predict ζ_c = 0.64 and 0.59 at these two Mach numbers, which are very close to the values measured from the simulation data. This confirms the validity of the cascade picture for scalars in flows with M ≲ 2. Physically, the flattening of the scalar structure function with M occurs because, as the scaling of δv(l) steepens at larger Mach numbers, there is relatively less velocity power at small scales, and thus the cascade of scalar fluctuations becomes slower, leading to a flatter scalar structure function and scalar spectrum.

On the other hand, Equation (4) is no longer valid at Mach numbers larger than 2 where the scalar structure function starts to become steeper. The scaling exponents are ζ_c = 0.60, 0.67, and 0.72 for M = 3, 4.6, and 6.1, respectively. As in the velocity case, the inertial-convective range becomes smaller at larger Mach numbers, and these exponents are measured in the range of l∈ [20–70]. More accurate measurement for these large Mach numbers requires higher resolution.

The reason for the steepening of the scalar structure functions is probably different from that for the velocity structure functions. First, for scalars there is no counterpart for the effect of pdV work in the velocity case, which, as discussed above, can cause kinetic energy loss along the cascade and make the velocity scaling steeper. Instead, the scalar energy is exactly conserved along the scalar cascade. Second, although at large Mach numbers scalar "edges" can contribute to the steepening of the scalar structure function, it is not clear how much contribution these edges can make. Unlike the velocity shocks, they are not discontinuities (see Section 4.2), and thus their contribution to steepening of the scalar structure functions may be smaller than that of shocks to that of the velocity structure functions.

We think that the more important reason responsible for the steepening of the scalar structure function for M ≃ 3 is related to strong compressions and expansions in these flows. At large Mach numbers, strong compressions concentrate most of the mass in small postshock regions, while most of the volume is occupied by expanding regions. In the expansion regions, the typical length scale of the scalar field increases as the scalar follows the expansion. On the other hand, the scalar length scale is reduced in compressed regions. Since the scalar structure function considered here is volume weighted, the effect of producing large structures by expansions dominates over the opposite effect by compressions. This leads to more power in scalar fluctuations at large scales, and thus a steeper volume-weighted structure function. This tends to counteract the flattening of the scalar structure function discussed above, and starts to dominate at M ≳ 3 when the density fluctuations become very strong. We speculate that this effect, which may also contribute to the steepening of the volume-weighted velocity structure function, could be removed by introducing an appropriate density weighting scheme. However, while several previous studies have shown the importance of density weighting in interpreting the statistics of supersonic turbulence (e.g., Kritsuk et al. 2007a; Pan & Padoan 2009; Pan et al. 2009) it is still not clear what would be an appropriate density weighting factor for two-point statistical measures.⁵ The topic of density weighting for structure functions is out of the scope of the current paper and will be explored separately in a future study.

To summarize these results, we find that for M ≲ 2 the scaling exponents for the scalar and velocity structure functions agree well with Equation (4), supporting the cascade picture for the scalar fluctuations at these Mach numbers. On the other hand, the steepening of the scalar structure function with the Mach number for M > 2 is in contradiction to the prediction by Equation (4). However, this does not mean that the cascade picture is invalid for scalars in high Mach number flows, but rather suggests the need for an appropriate density-weighting scheme when studying their statistics. In fact, the validity of the cascade picture at all Mach numbers is supported by our results for the mixing timescale and its scale dependence.

Figure 5 shows the power spectra of both the velocity and scalar fields, again employing volume-weighted averages. At the current resolution, there is barely any inertial range in the spectra for all Mach numbers. At small Mach numbers, there exist strong bumps, named bottlenecks (Falkovich 1994; Kritsuk 2007a), close to the dissipation range (see the compensated spectra in the insets). The bottlenecks appear to be weaker at larger Mach numbers. This suggests that the numerical scheme in the FLASH code probably applies a larger numerical viscosity to stabilize shocks at larger Mach numbers. We find that essentially no bottlenecks exist at M ≳ 3, which is different from simulations with less diffusive numerical codes (e.g., Kritsuk et al. 2007a; Lemaster & Stone 2009). As mentioned above, strong numerical diffusion at high Mach numbers may be responsible for the reduction in the inertial scale range with good power-law fits for the structure functions shown in Figure 4. However, as discussed below, we find evidence that strong numerical diffusion in the FLASH code at large M does not affect our conclusions on the behavior of the power spectra and structure functions in the inertial range.

Zoom In Zoom Out Reset image size

Due to the very limited inertial range in the spectra in Figure 5 at all Mach numbers, an accurate measurement of the spectra slope is not possible. However, the overall trend in the power spectra as a function of Mach number is the same as that in the structure functions. As can be seen from the insets, the velocity spectra keep steepening with increasing M, while the scalar spectra flatten at first and start to be steeper at M ≳ 3. Also similar to the structure functions, the scalar spectrum is significantly flatter than the velocity spectrum at all Mach numbers except the M = 0.9 case where the slopes are about equal. Future higher-resolution simulations are needed to quantify these results more accurately.

Finally, we discuss the potential effects of strong numerical diffusion in the FLASH code at large Mach numbers. We will show that our conclusions above are actually not affected. Due to the limited inertial range at current resolution, strong numerical diffusion may give rise to overestimates in the slopes of the power spectra and structure functions. Consider the structure functions for example. A large numerical viscosity would reduce the inertial range, and increase the chance that some of the scales in the range we chose to measure the structure functions would be affected by numerical diffusion. This may overestimate the slope of the structure function. Therefore, one may suspect that the steepening of the structure function in the inertial range found in our simulations is due to larger numerical viscosity at larger Mach numbers. We find that is not the case by comparing our results at M = 6 to those from simulations by Kristuk et al. (2007) using less diffusive codes, who find prominent bottlenecks in their velocity spectra.

Kritsuk et al. (2007a) found a slope of 0.95 for the volume-weighted longitudinal structure function in the inertial range of a M = 6 flow. This slope is exactly the same as that measured in the scale range of [20–70] resolution scales in our simulations at the same M. This suggests that the scale range we chose is essentially independent of bottlenecks. As seen in Figure 7 of Kritsuk et al. (2007a), bottlenecks may contaminate the velocity spectra only below ∼30 resolution scales. At those small scales, i.e., below 20–30 times the resolution scale, the velocity structure function in our M = 6 flow is steeper than that in Kristuk et al. (2007), consistent with stronger numerical diffusion in the FLASH code at large M. In conclusion, above ∼20 resolution scales, the slopes of the structure functions are not affected by bottlenecks, and the steepening trend at large M found in our simulations is realistic. This conclusion also applies to the scalar structure functions since the numerical diffusion for scalars in the FLASH code is essentially the same as the numerical viscosity.

Furthermore, a comparison of the velocity spectrum in our M = 6 flow to that in Kritsuk et al. (2007a) shows that the stronger numerical viscosity in our simulations makes the spectrum slightly steeper (by a few percent). However, this uncertainty of a few percent is smaller than the slope change from M = 3 to M = 6.1 found in our simulations, suggesting that the steepening of the velocity spectrum at large M is also realistic.

In this subsection, we study the probability distribution function (PDF) of forced scalars in compressible turbulence. Here, we first calculated the one-point density-weighted PDF, P(C, t), in each snapshot by $P(C, t) = \frac{1}{V} \int _V \delta (C-C(\boldsymbol{x},t))\tilde{\rho }(\boldsymbol{x},t) d\boldsymbol{x}$ where V is the volume of the simulation box. Then for each Mach number we obtained the average PDF, P(C), over 10 snapshots covering about 10 dynamical times. The results for Mach numbers 0.9 and 6.1 are shown in Figure 6, where all the PDFs are normalized to have unity variance. Note that the PDFs shown here are not the PDFs of the scalar differences between two points at different separations, which will not be considered in the present paper.

Zoom In Zoom Out Reset image size

The four lines in each panel of this figure correspond to different forcing schemes for the scalars. We find that, although the scalar source term is set to be Gaussian in all our simulations, the scalar PDF is Gaussian only if the scalar is forced at the flow driving scale. This agrees with the Gaussian scalar PDF found in Watanabe & Gotoh (2004) for incompressible turbulence where the scalar source was also injected at the flow driving scale. For scalars forced at smaller scales, the PDFs show broad non-Gaussian tails, which become broader with decreasing source length scale at all Mach numbers.

This behavior is likely to be related to the PDF of the velocity difference as a function of separation. The velocity PDF at the driving scale is nearly Gaussian, and the PDF for the scalar forced at this scale is nearly Gaussian as well. On the other hand, the PDF of the velocity difference across a small separation is known to have broad tails, which become broader as the separation decreases, a phenomenon known as intermittency (e.g., Frisch 1995). Since the scalar cascade is controlled by the velocity field, such small-scale non-Gaussian velocity structures are likely to leave a signature on the PDF tails of scalars forced at small scales. Furthermore, the small-scale velocity structures are known to become more intermittent at larger Mach numbers (see, e.g., Padoan et al. 2004; Kritsuk et al. 2007b). Again, this may be directly responsible for the fact that the PDF tails for scalars forced at small scales broaden with increasing Mach number, as can be seen by comparing the two panels in Figure 6.

We point out that, although the non-Gaussian velocity structures contribute to the broad tails in the one-point PDFs of scalars forced at small scales, it is not clear whether the non-Gaussianity in the velocity structures is the only cause for non-Gaussian tails in the scalar PDFs. In other words, it is possible that the scalar PDFs would exhibit non-Gaussian tails even if the flow velocity were Gaussian at all scales, and that non-Gaussian velocity structures only enhance an already existing effect. In fact, previous studies of mixing in incompressible turbulence show that at small scales the scalar difference PDF is non-Gaussian even if the advecting velocity field is exactly Gaussian (see, e.g., Shraiman & Siggia 2000), meaning that the non-Gaussianity in the scalar difference statistics is an intrinsic feature of mixing itself. Whether this is also true for the one-point PDFs of scalars forced at small scales is a question for future studies.

Finally, we consider the evolution of decaying scalar fields with no continuous sources (i.e., pollutants added to the flow only once). As mentioned earlier, we included four decaying scalars in our simulations, each of which has a different initial spectrum (or length scale). Figure 7 shows the decay of the density-weighted concentration variance for the four scalars, focusing on the results for M = 0.9 and M = 6.1. The results for other Mach numbers are qualitatively similar. The variance decay as a function of time and its dependence on the initial scalar scale is similar to that found by Eswaran & Pope (1988) for the incompressible case (see their Figure 10). Each curve in Figure 7 has a short initial period where it is flat and smooth. During this period, the cascade develops scalar structures toward smaller scales, and true mixing by molecular diffusion (numerical diffusion in our simulations) simply waits for the structures to reach the diffusion scale. Once structures of the size of the diffusion scale appear, the variance starts to decay steadily.

Zoom In Zoom Out Reset image size

We find that the concentration variance can be approximately fit by a single exponential function for decaying scalars with initial length scale equal to the flow driving scale (see 1 ⩽ k/2π ⩽ 3 curves in Figure 7). This is true for any Mach number. Fitting the variance decay of these scalars by exponentials shows that the decaying timescale is in the range of 0.5–0.6τ_dyn for M ≈ 1–6. The change of the decay timescale with M is similar to that for the forced cases, i.e., it increases with M at small M and becomes constant at M ≳ 3. We note that these timescales are smaller than those for the corresponding forced cases (see Table 2). The reason for this is that, without continuous sources at large scales, the spectrum of a decaying scalar is flatter than that of a forced scalar, and, with relatively more fluctuations at smaller scales, the mixing process is faster.

For scalars with initial scales smaller than the flow driving scale, the variance evolution cannot be fit by a single exponential function (with the 3.5 ⩽ k/2π ⩽ 4.5 curve for M = 6 being the only exception). The decay generally consists of two phases, an early fast phase and a slower phase that starts at a few dynamical timescales. In the fast decay phase, the timescale for the variance decay decreases with decreasing initial scalar scale, consistent with the scale dependence of the mixing timescale in the forced case. For example, roughly fitting exponentials to the fast phases of the bottom three curves in the M = 6.1 panel of Figure 7, we find that the timescales, from top to bottom, are respectively, 0.56τ_dyn, 0.4τ_dyn, and 0.3τ_dyn. Again, these timescales for the variance decay are smaller than the corresponding timescales for the forced cases. In the fast phase, the scalar fluctuations are significantly homogenized, with the variance reduced by a factor of 100–1000 in a few dynamical timescales.

After that, a slower phase starts, which corresponds to the evolution of the scalar spectrum. We find that if the initial scalar spectrum is within the inertial range of the flow, the spectrum spreads out in both directions, i.e., toward smaller scales, due to the cascade, and also toward larger scales. The spread to larger scales is due to the fact that interactions of scalar fluctuations at a given scale with velocity fluctuations at larger scales can produce scalar fluctuations at larger scales. We refer to this effect as the backward transfer. In a few dynamical timescales, the scalar power at small scales is significantly erased by mixing, and the fluctuation power starts to be dominated by structures at larger scales caused by the backward transfer. This results in an increase in the characteristic length scale for the scalar fluctuations, leading to slower mixing rate at later times. The backward transfer is expected to stop when the typical scalar length scale becomes close to the driving scale of the flow, L_f, beyond which there are no considerable velocity power. Therefore, in the later phase, the scalar length scale is always close to the flow driving scale, and the decay time is approximately given by the turnover time of the eddies at the flow driving scale. This is why at late time the decay rates for all scalars tend to become similar.

This picture also explains why the variance of a decaying scalar can be approximately fit by a single exponential function if its initial length scale is close to the flow forcing scale. In that case, the characteristic length scale of the scalar is fixed around that flow driving scale and does not change with time.

In summary, the timescale for the variance decay is smaller by 30%–40% than that obtained in the corresponding forced case. It is about 0.5–0.6 dynamical timescale if the initial scalar length scale is close to the flow driving scale. The dependence of the decay timescale on the initial length scale and on the Mach number for decaying scalars is generally consistent with that in the forced case.

Figure 8 shows the evolution of the density-weighted scalar probability distribution, P(C, t), for the scalar whose initial spectrum is in the range 1 ⩽ k/2π ⩽ 3, i.e., the same as the flow forcing spectrum. The initial PDFs are approximately Gaussian. The PDFs become narrower with time, indicating homogenization of the scalar fluctuations. We point out that Figure 8 is only meant to illustrate the general problem of how the PDF of a decaying scalar evolves. In practical applications, the initial PDF is typically not Gaussian, but instead is likely to be bimodal with two peaks, representing the unmixed pollutants and unpolluted flow, respectively (see, e.g., Pan 2008).

Zoom In Zoom Out Reset image size

In Figure 8, we see that at later times exponential tails develop in the PDFs. Note that in the corresponding forced cases, i.e., for scalars forced at the flow driving scale, the concentration PDF shown in Figure 6 is Gaussian (see discussion below). A likely origin for the fat tails in the PDFs of decaying scalars shown in Figure 8 at the later time is that mixing of concentration levels at the PDF tails is more difficult than that for the central part. Significantly moving the tails toward the central part requires contact and mixing between fluid elements with extreme concentrations corresponding to opposite tails of the scalar PDF, i.e., one from the high tail and the other from the low tail. Clearly, these events are much rarer than those for the mixing of the central part of the PDF. This suggests that the tails are more persistent, and this persistency tends to give fat tails.

Exponential tails also exist in the PDFs of decaying scalars with smaller initial length scales. Comparing the PDFs of different decaying scalars, we find that, at similar variances (which correspond to different evolution times for different scalars since their decaying timescales are different), the PDF tails are broader for scalars with smaller initial length scales. In other words, the exponential tails develop faster for decaying scalars with smaller initial length scales. The breadth of the tails also increases with increasing Mach numbers. These trends are consistent with those found for forced scalars. Similarly, this trend is related to the non-Gaussian velocity structures at small scales.

The exponential tails for the decaying scalars with initial scale at the flow driving scale raise the question why the PDFs of scalars continuously forced at the flow driving scale are Gaussian, considering that scalar sources injected at early times have "decayed" and would thus give exponential contributions to the PDFs at the current time. The answer relies on the fact that the main contribution to the PDF of a forced scalar at a given time is from sources injected within about a mixing timescale since the older sources have already been significantly mixed. We find that a decaying scalar with initial scale at the flow driving scale develops considerable exponential tails only after more than a mixing timescale. This explains why the PDFs in the corresponding forced cases are Gaussian. On the other hand, for scalars forced at smaller scales, recent scalar sources have larger contributions to the tails of the PDFs at the current time because fat tails develop faster for sources injected at smaller scales. This is consistent with exponential tails for scalars forced at small scales as found in Figure 6.