2. Numerical simulations

Three-dimensional compressible turbulence for an ideal gas is simulated in a periodic box of $(2{\rm\pi})^{3}$ by using a hybrid scheme (Wang et al. Reference Wang, Wang, Xiao, Shi and Chen2010). Namely, for the convective term an eighth-order compact scheme is used in the smooth region away from the shocklets and a seventh-order weighted essentially non-oscillatory scheme is used near the shocklets. For time advancing we employ an explicit low-storage second-order Runge–Kutta method. The flow is driven by a large-scale force which can be applied to both solenoidal and compressive components of the velocity field. More details about the numerical method are reported in Wang et al. (Reference Wang, Yang, Shi, Xiao, He and Chen2013a ). Four cases are simulated with different driving forces and turbulent Mach numbers $M_{t}=U_{\mathit{rms}}/c$ , where $U_{\mathit{rms}}$ is the root-mean-square (r.m.s.) value of the velocity magnitude and $c$ is the mean speed of sound. The resolution is fixed at $512^{3}$ . The parameters for the simulations are listed in table 1. Cases S045, S075 and S104 are driven by purely solenoidal forces. They have different values of $M_{t}$ and similar Taylor microscale Reynolds numbers $R_{{\it\lambda}}$ . Case C069 is driven by a combination of compressive and solenoidal forces with a ratio of 4:1.

Table 1. Parameters of the numerical simulations. Columns from left to right: case name, the ratio of the solenoidal force to the compressible force, the turbulent Mach number, the Taylor microscale Reynolds number, the viscous time scale and the integral time scale. In the case names the leading letter indicates purely solenoidal driving (S) or the driving force dominated by the compressive component (C), and the following three digits indicate the Mach number.

A Lagrangian module has been developed for particle tracking (Yang et al. Reference Yang, Wang, Shi, Xiao, He and Chen2013). For each case, after the flow reaches the statistically steady state, a million passive tracers are seeded uniformly and advected by the background flow. The flow information at the exact locations of particles is obtained by a trilinear interpolation of the values at surrounding grid points. The particle locations are then updated by using the same time-marching scheme as the Eulerian solver. The trilinear interpolation is adequate in our simulations because (1) our mesh size is slightly smaller than the Kolmogorov scale and (2) linear interpolation may avoid any unphysical oscillation near shocklets where the flow quantities change abruptly. The flow then evolves for a time period of approximately $11T^{E}$ . Here, $T^{E}=\sqrt{3}L/u_{\mathit{rms}}$ is the large-eddy turnover time, where $L$ is the integral length scale. During the last $8T^{E}$ , two sets of tracer data are stored. One set of data covers the whole $8T^{E}$ period with a sampling rate of approximately $0.1{\it\tau}_{{\it\eta}}$ , where ${\it\tau}_{{\it\eta}}=({\it\nu}/{\it\epsilon})^{1/2}$ denotes the Kolmogorov time scale, with ${\it\nu}$ being the mean kinematic viscosity and ${\it\epsilon}$ the energy-dissipation rate per unit mass. Another set of data consists of several segments, all of which last a time period of $0.4{\it\tau}_{{\it\eta}}$ and are separated by $T^{E}/3$ from one another. Within each segment the sampling rate is approximately $0.01{\it\tau}_{{\it\eta}}$ . From the former dataset we calculate the Lagrangian statistics for large time scales, which are usually averaged over a time period of approximately $5T^{E}$ and include in total $4\times 10^{7}$ data points. The statistics of small time scales are calculated from the latter dataset, which covers a time period of approximately $3T^{E}$ and contains $10^{7}$ data points.

Figure 1. Contours of $\log ({\it\rho}/{\it\rho}_{0})$ on the midplane of the domain for (a) case S045, (b) case S075, (c) case S104 and (d) case C069. Due to the small density variation in case S045, the colour map in (a) has half of the value range compared with the others.

In figure 1 we show the contours of the logarithm of density $\log ({\it\rho}/{\it\rho}_{0})$ on one slice, where ${\it\rho}_{0}$ is the mean density. For case S045 with solenoidal driving at the smallest Mach number, the density is non-uniform but no distinct shocklet is observed. As the Mach number increases, the inhomogeneity of the density is enhanced and the shocklets become stronger. These shocklets are still localized in space, i.e. the size of the shocklets is much smaller than the domain size. For case C069 with hybrid driving, the morphology of the shock structures is totally different from those with solenoidal force. Large-scale shocks develop and extend through the entire domain. The flow field is divided into blocks by these large shock structures, and within each block the variation of density is relatively weak.

3. Lagrangian statistics

Different force schemes have profound influences on the behaviour of the LSFs, as illustrated in figure 2, where we plot the LSFs up to the sixth order. For case S045, the LSFs are quite similar to those of the incompressible turbulence, since at such a low Mach number the compressibility effect is negligible, which has also been confirmed by simulations with different numerical schemes (Benzi et al. Reference Benzi, Biferale, Fisher, Kadanoff, Lamb and Toschi2008). As $M_{t}$ increases, the LSFs at ${\it\tau}>10{\it\tau}_{{\it\eta}}$ show no visible changes. However, a new self-similar region emerges below the Kolmogorov time scale and expands as $M_{t}$ increases. For case C069, the LSFs exhibit a self-similar region in the range $0.2{\it\tau}_{{\it\eta}}<{\it\tau}<8{\it\tau}_{{\it\eta}}$ and saturate for larger ${\it\tau}$ . The shape is very similar to those at higher $M_{t}$ (Konstandin et al. Reference Konstandin, Federrath, Klessen and Schmidt2012). All of these results indicate that for the Mach number considered here a purely solenoidal driving only has notable effects on the LSFs at short time scale, while a compressive driving can change the behaviour of the LSFs from the sub-Kolmogorov time scale to much larger scales.

Figure 2. Even-order LSFs for (a) case S045, (b) case S075, (c) case S104 and (d) case C069: green, $p=2$ ; blue, $p=4$ ; red, $p=6$ . The error bar is estimated by the standard deviation of the temporary fluctuation. The dark grey shade marks the range where the relative scaling exponent ${\it\zeta}_{p}$ is calculated and the light grey shade marks the range where ${\it\xi}_{p}$ is calculated.

We performed an extended self-similarity (ESS) analysis (Benzi et al. Reference Benzi, Ciliberto, Tripiccione, Baudet, Massaioli and Succi1993) and calculated the local slope ${\it\chi}_{p}=\text{d}\log (S_{p}({\it\tau}))/\text{d}\log (S_{2}({\it\tau}))$ . In figure 3 we plot the curves of ${\it\chi}_{p}$ for $p=4$ and 6. All of the curves approach the dimensional non-intermittent value $p/2$ as ${\it\tau}$ decreases to zero. This confirms that our Lagrangian simulation has sufficient time resolution. For case S045, the curve shares similar characteristics to those for incompressible turbulence (Arnèodo et al. Reference Arnèodo, Benzi, Berg, Biferale, Bodenschatz, Busse, Calzavarini, Castaing, Cencini, Chevillard, Fisher, Grauer, Homann, Lamb, Lanotte, Lévèque, Lüthi, Mann, Mordant, Müller, Ott, Ouellette, Pinton, Pope, Roux, Toschi, Xu and Yeung2008; Biferale et al. Reference Biferale, Bodenschatz, Cencini, Lanotte, Ouellette, Toschi and Xu2008; Toschi & Bodenschatz Reference Toschi and Bodenschatz2009). A dip appears at approximately ${\it\tau}\approx 3{\it\tau}_{{\it\eta}}$ , which is caused by small-scale vortex filaments (Biferale et al. Reference Biferale, Boffetta, Celani, Lanotte and Toschi2005; Bec et al. Reference Bec, Biferale, Cencini, Lanotte and Toschi2006). We refer to this as the vortex dip. Interestingly, a new dip emerges at sub-Kolmogorov time scale and becomes deeper when $M_{t}$ increases, as shown by the curves of cases S075 and S104. Meanwhile the vortex dip becomes shallower and eventually invisible in the curve of case S104. Thus, a new type of intermittency develops at small time scale when the compressibility is strong enough. For case C069, the situation is totally different. As ${\it\tau}$ increases from small time scale, ${\it\chi}_{p}$ quickly drops to near unity and stays at this low value over a relatively wide range of time scales, and then slowly increases before the LSFs saturate. This new dip, starting from sub-Kolmogorov time scale, will be referred to as the shocklet dip. We will show below that it can be related to shocklets.

Figure 3. Extended self-similarity local slope ${\it\chi}_{p}$ for $p=4$ and 6: green, case S045; blue, case S075; red, case S104; grey, case C069. The error bar is estimated by the standard deviation of the temporary fluctuation and is only shown at every other data point for clarity. For case C069, ${\it\chi}_{p}$ is computed in the range of ${\it\tau}<12{\it\tau}_{{\it\eta}}$ before the LSFs saturate.

Our results indicate that two different types of self-similarities may exist at different time scales. Thus, two relative scaling exponents can be extracted for both self-similar regions. For case S045, only one relative scaling exponent ${\it\zeta}_{p}$ is calculated by a linear fitting of $(\log (S_{2}({\it\tau}))\log (S_{p}({\it\tau})))$ at $8<{\it\tau}/{\it\tau}_{{\it\eta}}<T_{max}=20$ with $T_{max}$ comparable to $T^{E}$ . This region corresponds to the plateau at the scales above the vortex dip in the ${\it\chi}_{p}$ curves, as shown in figure 3. Another relative scaling exponent ${\it\xi}_{p}$ is calculated by a linear fitting at $0.2<{\it\tau}/{\it\tau}_{{\it\eta}}<0.7$ for case S075 and $0.2<{\it\tau}/{\it\tau}_{{\it\eta}}<1.0$ for case S104. For case C069, we only calculate ${\it\xi}_{p}$ over the region $0.2<{\it\tau}/{\it\tau}_{{\it\eta}}<8.0$ . These regions are marked by grey shade in figure 2. The exponents are shown in figure 4. One may expect that ${\it\zeta}_{p}$ will recover the value of the incompressible flow as $M_{t}$ decreases to zero. This is exactly the case in our results, as shown in figure 4, where ${\it\zeta}_{p}$ from our simulations is compared with the values reported in Benzi et al. (Reference Benzi, Biferale, Fisher, Lamb and Toschi2010). For cases S045 and S075 with lower $M_{t}$ , ${\it\zeta}_{p}$ takes a value very close to that of the incompressible case, although our simulations have smaller Reynolds numbers. For higher $M_{t}$ such as case S104, ${\it\zeta}_{p}$ is considerably higher than the incompressible value at high order. This indicates that compressibility starts to affect the overall turbulent statistics over a wide range of time scales.

Figure 4. Extended self-similarity relative scaling exponents up to $p=8$ for ${\it\zeta}_{p}$ (red upward triangles) and ${\it\xi}_{p}$ (blue downward triangles) for (a) case S045, (b) case S075, (c) case S104 and (d) case C069. The grey circles in the plots for cases S045, S075 and S104 mark the values for the incompressible case with $Re_{{\it\lambda}}\sim 600$ (Benzi et al. Reference Benzi, Biferale, Fisher, Lamb and Toschi2010). For cases S045 and C069, only the applicable exponents are shown.

The relative exponent ${\it\xi}_{p}$ calculated in the new self-similar region shows totally different behaviour compared with ${\it\zeta}_{p}$ . It is slightly larger than ${\it\zeta}_{p}$ for $p=1$ and smaller than ${\it\zeta}_{p}$ for $p>2$ , indicating a much stronger intermittency. It saturates when $p>4$ for cases S075 and S104. The saturated exponent ${\it\xi}_{\infty }$ , which is calculated by averaging ${\it\xi}_{p}$ over $p=5$ –8, is approximately 1.47 for case S075 and 1.19 for case S104. For case C069, ${\it\xi}_{p}$ does not saturate completely but increases very slowly when $p>4$ . We define ${\it\xi}_{\infty }={\it\xi}_{8}\approx 1.19$ for case C069. Based on multi-fractal theory, the saturation of the exponent suggests that the tails of the probability density function (PDF) $P({\it\delta}u,{\it\tau})$ should collapse for different values of ${\it\tau}$ in the self-similar region when normalized by ${\it\tau}^{{\it\xi}_{\infty }}$ (Benzi et al. Reference Benzi, Biferale, Fisher, Kadanoff, Lamb and Toschi2008; Wang et al. Reference Wang, Yang, Shi, Xiao, He and Chen2013b ). We plot in figure 5 the normalized PDF ${\it\tau}^{-{\it\xi}_{\infty }}P$ for case S104 within ${\it\tau}/{\it\tau}_{{\it\eta}}\in (0.2,1.0)$ and for case C069 within ${\it\tau}/{\it\tau}_{{\it\eta}}\in (0.2,8.0)$ , i.e. the region where ${\it\xi}_{p}$ is evaluated. Indeed, the tails of the normalized PDFs collapse. It is remarkable that for case C069 the peak of the normalized PDFs at ${\it\delta}u=0$ decreases for several orders as ${\it\tau}/{\it\tau}_{{\it\eta}}$ increases, while the tails are still very close to one another for the entire range of ${\it\tau}/{\it\tau}_{{\it\eta}}\in (0.2,8.0)$ .

Figure 5. Probability density functions of ${\it\delta}u/u_{\mathit{rms}}^{c}$ normalized by ${\it\tau}^{{\it\xi}_{\infty }}$ , where $u_{\mathit{rms}}^{c}$ is the r.m.s. value of the velocity component: (a) ${\it\tau}/{\it\tau}_{{\it\eta}}\in (0.2,1.0)$ for case S104 and (b) ${\it\tau}/{\it\tau}_{{\it\eta}}\in (0.2,8)$ for case C069. The colour changes from blue to red as ${\it\tau}$ increases.

It has been revealed that the exponents of the Eulerian structure functions saturate due to the appearance of the front-like or shocklet structures (Benzi et al. Reference Benzi, Biferale, Fisher, Kadanoff, Lamb and Toschi2008; Wang et al. Reference Wang, Yang, Shi, Xiao, He and Chen2013b ). From the Lagrangian point of view, when a passive tracer moves through a shocklet it gains a velocity increment which is determined by the strength of the shocklet. In a homogeneous flow, a reasonable assumption is that the probability of a tracer encountering a shocklet is proportional to the time lag ${\it\tau}$ , which gives ${\it\delta}u\sim \text{constant}\cdot {\it\tau}$ . Therefore, the saturated value of the exponent should be very close to 1 if the shocklets in the flow have similar strength and tracers only penetrate shocklets once in the time interval considered. This argument is consistent with our numerical results. For a short time interval, it is likely that tracers encounter shocklets at most once, and thus the local relative exponent should be close to 1. Indeed, in figure 3 we observe that with solenoidal driving force the local minimum of ${\it\chi}_{p}$ within the shocklet dip approaches 1 as $M_{t}$ increases, while the ${\it\chi}_{p}$ of case C069 equals 1 at approximately ${\it\tau}={\it\tau}_{{\it\eta}}$ . For the three cases with solenoidal driving force, the major portion of the flow field is dominated by vortical structures, and the LSFs are similar to those of the incompressible cases at large time scales. However, the effects of compressibility do emerge at higher $M_{t}$ . With hybrid driving, however, shocks are more pronounced and dramatically change the LSFs even at large time lags.

Now it is clear that the shocklet dip in the local slope curves is related to the shocklet structures in flow. It is worth pointing out that the vortex dip and the shocklet dip have different shapes and respond differently to their relevant control parameters. For incompressible flow, the location and shape of the vortex dip are universal and have minor dependences on the Reynolds number (Arnèodo et al. Reference Arnèodo, Benzi, Berg, Biferale, Bodenschatz, Busse, Calzavarini, Castaing, Cencini, Chevillard, Fisher, Grauer, Homann, Lamb, Lanotte, Lévèque, Lüthi, Mann, Mordant, Müller, Ott, Ouellette, Pinton, Pope, Roux, Toschi, Xu and Yeung2008). The shocklet dip reported here, on the other hand, expands over time scales as $M_{t}$ increases or shocklets become dominant. This is caused by the different natures of the extreme velocity differences induced by vortex filaments and shocklets. For a tracer trapped by a vortex filament, the largest velocity difference occurs at the time lag during which the tracer travels half of the circle around the vortex filament, regardless of the strength of the vortex filament. However, when a tracer penetrates a shocklet, the velocity differences over a time lag are predominantly determined by the strength of the shocklet, provided that the tracer is located on different sides of the shocklet at the beginning and end of that time lag. As the shocklet becomes stronger, its influence lasts for a longer time period.