2. Computational setup

The computational arrangement for the DNS is sketched in figure 1. Two turbulent boundary layers having the same thickness ( ${\it\delta}_{0}$ ) at the inflow of the computational domain, develop on the two sides of a splitter plate with zero thickness and merge past the trailing edge, forming a turbulent mixing layer. The mesh is uniform along the streamwise ( $x$ ) and spanwise ( $z$ ) direction, and stretched along the wall-normal direction to have ${\rm\Delta}y^{+}\lesssim 1$ at the wall. The computational box has a size of $L_{x}\times L_{y}\times L_{z}=120{\it\delta}_{0}\times 47{\it\delta}_{0}\times 12{\it\delta}_{0}$ , to be able to have fully developed upstream boundary layers, and follow the mixing layer evolution into the developed state. The computational domain has been selected based on previous experience (Foysi & Sarkar Reference Foysi and Sarkar2010; Sharma et al. Reference Sharma, Bhaskaran and Lele2011), to minimize numerical confinement effects, within the constraints of computational cost. The origin of the coordinate system is set at the trailing edge of the splitter plate, at $x_{0}\approx 40{\it\delta}_{0}$ downstream of the inflow. The physical parameters of the two incoming boundary layers are listed in table 1, where subscripts 1 and 2 are used to denote the upper and lower streams, respectively. While the free-stream conditions are made to fit those of case 2 of Goebel & Dutton (Reference Goebel and Dutton1991), the properties of the two boundary layers are somewhat different, mainly as we had to reduce the Reynolds numbers by about a factor of three (in the experiment of Goebel & Dutton Reference Goebel and Dutton1991, $\mathit{Re}_{{\it\delta}_{1}}\approx 51\,765$ , $\mathit{Re}_{{\it\delta}_{2}}\approx 56\,286$ ) to make the simulation feasible. Furthermore, in the experiment the two streams were brought together to form the mixing layer at an angle of $2.5^{\circ }$ across a thin splitter plate with a base height of approximately 0.5 mm (18 % of the averaged boundary layer thickness). As shown in the discussion below, these differences seem to have a very minor impact on the mixing layer development. The computational mesh includes $N_{x}\times N_{y}\times N_{z}=6144\times 301\times 768$ points, and the spacing in the wall-parallel directions (taken just upstream of the trailing edge) is ${\rm\Delta}x^{+}=4.42,{\rm\Delta}z^{+}=3.52$ for the upper stream, and ${\rm\Delta}x^{+}=9.83,{\rm\Delta}z^{+}=7.68$ for the lower stream, where the superscript $+$ indicates normalization by local wall units, based on the friction velocity $u_{{\it\tau}}$ and the viscous legthscale ${\it\delta}_{v}={\it\nu}_{w}/u_{{\it\tau}}$ , being ${\it\nu}_{w}$ the kinematic viscosity evaluated at the wall. In our experience (Pirozzoli & Bernardini Reference Pirozzoli and Bernardini2013), and based on the results of the next section, the mesh resolution is sufficient to have well-behaved boundary layer statistics.

Figure 1. Sketch of the computational domain; ${\it\delta}_{0}$ is the boundary layer thickness of the two streams at the inflow station. The subscript 1 is used to denote the upper (high-speed) stream, and 2 for the lower (low-speed) stream.

The numerical simulation has been carried out by means of a conservative finite-difference algorithm. The flow solver relies on central sixth-order discretization of the convective terms of the Navier–Stokes equations cast in fully split form (Pirozzoli Reference Pirozzoli2011), which guarantees discrete conservation of the total kinetic energy in the limit case of inviscid, incompressible flow, also in the presence of grid stretching in the coordinate directions. This approach allows robust spatial discretization of the convective terms without the addition of any spurious numerical dissipation in the form of upwinding or filtering. Being cast in locally conservative form, the method also guarantees excellent computational efficiency, facilitating hybridization with shock-capturing methods. A fifth-order weighted-essentially-non-oscillatory (WENO) scheme is used to replace central differencing in the presence of discontinuities, with switch controlled by the Ducros sensor (Bernardini, Pirozzoli & Grasso Reference Bernardini, Pirozzoli and Grasso2011). It is worth noting that in the present simulation WENO is only activated at about 0.36 % of the total number of mesh points, hence artificial viscosity has very little impact on the computed solution. The diffusive terms in the Navier–Stokes equations are expanded to Laplacian form for improved numerical stability, and also discretized with sixth-order central-difference formulae, to guarantee proper action of molecular diffusion at the smallest scales resolved on the computational mesh. The resulting semi-discrete system of equations is advanced in time by means of an optimized, fully explicit fourth-stage Runge–Kutta algorithm (Bernardini & Pirozzoli Reference Bernardini and Pirozzoli2009). Particular attention is devoted to the correct prescription of the inflow boundary conditions, which is a key ingredient in the simulation of spatially developing turbulent flows. In the present simulation the inflow conditions are prescribed through a recycling–rescaling procedure, suitably adapted to compressible flows (Bernardini et al. Reference Bernardini, Pirozzoli and Grasso2011), and the recycling station is placed sufficiently far from the inlet (about 30 ${\it\delta}_{0}$ ) to achieve full decorrelation from the inflow, as seen from inspection of the streamwise and spanwise two-point velocity correlations. In the absence of detailed information from experiments, the two walls are assumed to be adiabatic.

In the following, unless otherwise specified, the velocity difference between the two streams, ${\rm\Delta}\unicode[STIX]{x1D650}=U_{1}-U_{2}$ , is used to normalize the velocity statistics, and lengths are made non-dimensional with respect to the arithmetic average of the 99 % thicknesses of the upper and lower boundary layers ( $\overline{{\it\delta}}=({\it\delta}_{1}+{\it\delta}_{2})/2$ ) evaluated just upstream of the trailing edge. The following scaled coordinates will then be used: ${\it\xi}=(x-x_{0})/\overline{{\it\delta}}$ , ${\it\eta}=y/\overline{{\it\delta}}$ , ${\it\zeta}=z/\overline{{\it\delta}}$ .

3. Velocity statistics

The statistical properties of the incoming boundary layers are provided in figure 2, where we show the mean velocity profile in van Driest scaling, $\overline{u}_{VD}=\int _{0}^{\overline{u}}(\overline{{\it\rho}}/\overline{{\it\rho}}_{w})^{1/2}\,\text{d}\overline{u}$ , and the density-scaled Reynolds stress components, $\unicode[STIX]{x1D619}_{ij}=(\overline{{\it\rho}}/\overline{{\it\rho}}_{w})\overline{u_{i}^{\prime }u_{j}^{\prime }}$ (here the overbar denotes Reynolds-averaged values, and the primes fluctuations thereof). In figure 2(a) the linear viscous sublayer law and the logarithmic overlap layer law are also shown for reference, with standard values of the log-law constants, $k=0.41$ , $C=5$ . The wake region is found to be stronger for the lower stream, which has higher friction Reynolds number, and consistently the streamwise velocity variance is also found to be somewhat higher than for the upper stream. Overall, the figure suggests that the incoming boundary layers are fully developed, thus ensuring proper initiation of the mixing process of the two streams.

Figure 2. Mean velocity profile (a) and Reynolds stresses (b) on the two sides of the splitter plate, upstream of the trailing edge ( ${\it\xi}=-1$ ). Solid lines indicate the upper stream, and dashed lines the lower stream. The dot-dashed lines in (a) are the viscous sublayer linear law ( $u^{+}=y^{+}$ ), and the logarithmic overlap layer ( $u^{+}=C+1/k\log y^{+}$ , with $C=5$ , $k=0.41$ ).

Figure 3. Mean spanwise vorticity in the $x{-}y$ plane (axes are not to scale). The sixteen contours are from $-1$ to 1, in colour scale from blue to red. The thick lines denote the lower ( $y_{L}$ ) and upper ( $y_{U}$ ) boundary of the mixing layer (in black), and the mixing layer mid-line ( $y_{0}$ , in white), according to the definition given by Goebel & Dutton (Reference Goebel and Dutton1991).

Mean vorticity contours are shown in a side view in figure 3. A weak compression wave is observed to form at the trailing edge on the high-speed (upper) side, whereas a weak expansion wave forms on the low-speed (lower) side, with subsequent upward deflection of the mixing layer. It is known in the incompressible case (Pope Reference Pope2000) that the deflection of the mixing layer axis occurs towards the side with less entrainment, which is the low-speed side, under ordinary circumstances. The slight upward deflection observed here is thus likely to be an effect of the non-unity density ratio. To visually appreciate the mixing layer evolution, in figure 3 we show the mixing layer mid-line ( $y_{0}$ , marked as a thick white line), and the position of the lower and upper boundary of the mixing layer ( $y_{L}$ , $y_{U}$ , thick black lines). These are defined by vertically scanning the mean velocity profiles as suggested by Goebel et al. (Reference Goebel, Dutton, Krier and Renie1990), as follows. For those profiles for which a location where $\overline{u}<U_{2}-0.1{\rm\Delta}\unicode[STIX]{x1D650}$ does not exist (i.e. the momentum deficit is low), the upper and lower boundaries of the mixing layer are defined as the transverse locations where the mean streamwise velocity equals $U_{1}-0.1{\rm\Delta}\unicode[STIX]{x1D650}$ and $U_{2}+0.1{\rm\Delta}\unicode[STIX]{x1D650}$ , respectively, and the mid-line is defined as the transverse position where the mean velocity equals $(U_{1}+U_{2})/2$ . For those profiles with strong momentum deficit, the lower edge of the mixing layer is defined to be the transverse location where the mean velocity equals $U_{2}-0.1{\rm\Delta}\unicode[STIX]{x1D650}$ , and the mixing layer centreline is defined as the average of $y_{U}$ and $y_{L}$ . This definition explains the jump in $y_{L}$ occurring at ${\it\xi}\approx 20$ .

A comparison of the evolution of the mixing layer apparent thickness, $b=y_{U}-y_{L}$ , with the experimental data of Goebel & Dutton (Reference Goebel and Dutton1991) is presented in figure 4(a). Very good agreement is recovered, even in the early development region. It is noteworthy that the position with a sudden drop of $b$ , which according to the previously given definition corresponds to a change from a strong velocity deficit state to a ‘canonical’ mixing layer state, is well captured. In figure 4 we also show the evolution of the vorticity thickness, defined as ${\it\Delta}_{{\it\omega}}={\rm\Delta}\unicode[STIX]{x1D650}/(\max _{y}\partial \overline{u}/\partial y)$ , which is the most frequently used indicator for mixing layer thickness in numerical simulations. As expected from its definition, ${\it\Delta}_{{\it\omega}}$ shows fast growth past the splitter plate trailing edge, associated with a sudden reduction of shear, and it attains a nearly linear behaviour at a distance of about $20\overline{{\it\delta}}$ past the trailing edge. In this region, which is typically associated with a self-similar state of the mixing layer (Pope Reference Pope2000), we find ${\it\Delta}_{{\it\omega}}\approx b$ , in agreement with the claims of Goebel & Dutton (Reference Goebel and Dutton1991).

Figure 4. Evolution of mixing layer thicknesses past the splitter plate, scaled by the reference upstream boundary layer thickness, $\overline{{\it\delta}}=({\it\delta}_{1}+{\it\delta}_{2})/2$ , in linear scale (a) and logarithmic scale (b). Open symbols: apparent thickness $b$ (squares), vorticity thickness ${\it\Delta}_{{\it\omega}}$ (circles). Solid squares: Goebel & Dutton (Reference Goebel and Dutton1991). Dashed line: ${\it\Delta}_{{\it\omega}}\sim {\it\xi}^{1.36}$ ; solid line: ${\it\Delta}_{{\it\omega}}\sim {\it\xi}^{0.587}$ .

Clearer scrutiny of the mixing layer thickness evolution is provided in figure 4(b), where the vorticity data are shown in a logarithmic scale, to highlight power-law behaviour. It is found that in the very early stages of mixing layer formation ( ${\it\xi}\lesssim 3$ ), the vorticity thickness increases faster than linearly, as

This behaviour can be tentatively explained as follows. Referring to figure 5, we assume that: (i) the merging boundary layers have the same thickness ${\it\delta}$ (see table 1), and have a power-law mean velocity profile, hence $\overline{u}_{i}(y)=U_{i}(y/{\it\delta})^{{\it\alpha}}$ for $y\leqslant {\it\delta}$ ; (ii) the mixing layer spreads symmetrically with respect to the geometric mid-plane, and ${\it\Delta}({\it\xi})$ is its local thickness; (iii) for $|y|\geqslant {\it\Delta}/2$ , the boundary layers are not affected by the entrainment process. From these hypotheses, taking as mixing length the local mixing layer thickness, and the local velocity difference across the mixing layer, $\overline{u}_{1}({\it\Delta}/2)-\overline{u}_{2}({\it\Delta}/2)$ as the typical eddy velocity scale, it follows that the relevant eddy viscosity at a given streamwise station is

Dimensional arguments easily yield the mixing layer thickness evolution equation

where $U_{0}=(U_{1}+U_{2})/2$ is an average convection velocity, whence it readily follows

Incidentally, after the mixing layer has grown to a thickness greater that the boundary layer thickness, the above analysis returns the classical linear growth law (Pope Reference Pope2000). However, until saturation is complete, (3.4) suggests that the growth rate should be faster than linear. At the low Reynolds number of the boundary layers under consideration, power-law exponents are expected to be in the range $0.182\lesssim {\it\alpha}\lesssim 0.196$ (Pope Reference Pope2000), which we have verified by fitting the DNS velocity profiles. The growth rate exponent predicted by (3.4) is then in the range between 1.22 and 1.24, not too far from the DNS data fit of (3.1).

Figure 5. Idealized sketch of mixing layer early development; ${\it\delta}$ is the initial thickness of the boundary layers (assumed to be the same), and ${\it\Delta}$ is the local mixing layer thickness.

It is also noteworthy that, although figure 4(a) seems to suggest that linear growth is achieved past ${\it\xi}\approx 20$ , figure 4(b) in fact favours power-law growth, with

In retrospect, a similar behaviour can be traced in the numerical simulation of Sandham & Sandberg (Reference Sandham and Sandberg2009), for the merging of a laminar and a turbulent stream. This non-canonical behaviour will be further commented on below.

Figure 6. Comparison of mean velocity and Reynolds stress components at various streamwise stations past the splitter plate: ${\it\xi}=3.65$ (a), ${\it\xi}=9.11$ (b), ${\it\xi}=19.69$ (c), ${\it\xi}=36.46$ (d), ${\it\xi}=54.68$ (e). Lines: present DNS; symbols: (Goebel & Dutton Reference Goebel and Dutton1991).

A comparison of the velocity statistics with the experimental data of Goebel & Dutton (Reference Goebel and Dutton1991) is presented in figure 6 at various stations across the developing mixing layer, where the cross-stream coordinate is scaled with respect to $b$ , and the reference velocity is ${\rm\Delta}\unicode[STIX]{x1D650}$ . The scaled mean velocity profiles (left-most column) highlight the presence of significant momentum defect in the proximity of the plate trailing edge caused by flow retardation within the issuing boundary layers, which is gradually filled up as the flow evolves, although remnants are still visible (see also the analysis below) near the outflow of the computational domain. Overall, agreement with the experiments is quite remarkable, also given the previously noticed discrepancy in the Reynolds numbers and the assumption of a splitter plate with zero thickness. Similar considerations also apply to the streamwise turbulence intensity (second column from left), whose main peak is well predicted throughout, occurring in all cases near the mixing layer mid-line, where turbulence kinetic energy production is maximum. A shoulder in the streamwise velocity variance distribution is observed in the lower part of the mixing layer, approximately corresponding to the inflectional point of the momentum defect zone. Poorer agreement is found for the vertical Reynolds stress component (third column), which is overestimated with respect to experiments by up to 50 %, and which is likely to be due to difficulties in measuring the cross-stream velocity with laser Doppler velocimetry. This also implies, although to a lesser extent, discrepancies in the prediction of the shear stress peak (last column). The turbulence levels away from the mixing layer core are found to be substantially higher in experiments, presumably because of sound background wind tunnel turbulence. Similar results for the fully developed state are obtained from spatial simulations with synthetic inflow conditions (Salvadore, Bernardini & Botti Reference Salvadore, Bernardini and Botti2013).

The evolution of the maximum velocity deficit across the mixing layer, defined as $U_{2}-\min _{y}\overline{u}(y)$ , is shown in figure 7. Based on its definition, this property is expected to tend asymptotically to zero as $x\rightarrow \infty$ , where a canonical mixing layer behaviour should be recovered. Figure 7 shows a two-stage evolution whereby the momentum defect initially scales as ${\it\xi}^{-1/2}$ , and then decays exponentially past ${\it\xi}\approx 20$ . This behaviour can be explained by the prevalence of a wake-like behaviour in the trailing edge near field, recalling the scaling of momentum defect in a planar wake (Pope Reference Pope2000), which is then followed by rapid fill-up in the far field. It is noteworthy that linear stability analysis of mixing layers in the presence of momentum defect (Zhuang & Dimotakis Reference Zhuang and Dimotakis1995) suggested the existence of two modes of instability, namely a wake mode which is dominant in the presence of strong deficit, superposed onto a canonical mixing layer mode. The results of the stability analysis appear to be consistent with the observed two-stage evolution of the mixing layer properties.

Figure 7. Evolution of maximum velocity defect with trailing edge distance. The dashed line denotes the trend ${\it\xi}^{-1/2}$ .

Figure 8. Evolution of Reynolds stress peaks with trailing edge distance, in linear scale (a) and logarithmic scale (b). Symbols: $\overline{u^{\prime 2}}$ (circles), $\overline{v^{\prime 2}}$ (squares), $\overline{w^{\prime 2}}$ (triangles), $\overline{u^{\prime }v^{\prime }}$ (diamonds). The dashed line in (b) denotes the trend ${\it\xi}^{-0.33}$ .

The evolution of the Reynolds stress peaks with the distance from the trailing edge is depicted in figure 8. A sudden rise of all Reynolds stress components is observed past the trailing edge up to ${\it\xi}\approx 3$ , which is associated with flow rearrangement, and then followed by monotonic decay. This behaviour is overall consistent with the findings of Bell & Mehta (Reference Bell and Mehta1990) for a low-speed tripped mixing layer. However, those authors were unable to measure close enough to the trailing edge to clearly probe the near-field peaks. Most previous studies, including Bell & Mehta (Reference Bell and Mehta1990), Goebel & Dutton (Reference Goebel and Dutton1991), as well as the theory of self-similar mixing layers (Pope Reference Pope2000) suggest saturation of the Reynolds stresses with the trailing edge distance, when scaled by the mixing layer velocity difference. While tendency to saturation here may be inferred from a linear representation as in figure 8(a), the same data cast in a log-log representation (figure 8 b) do in fact suggest continuing decrease of the Reynolds stresses in the far-field region, according to a weak inverse power law of the streamwise coordinate, and no clear tendency to saturation. In particular, the inferred behaviour for the turbulent shear stress, which is the only component entering the mean streamwise momentum equation, is $-\overline{u^{\prime }v^{\prime }}\sim {\it\xi}^{-0.33}$ , for ${\it\xi}\gtrsim 10$ . This finding is consistent with the previously noticed nonlinear growth of the mixing layer thickness (recalling figure 4 b). In fact, analysis of the mean streamwise momentum equation (Townsend Reference Townsend1976; Pope Reference Pope2000) indicates that a necessary condition for self-similarity is

where $u_{0}$ is the typical velocity scale of turbulent fluctuations. The power-law exponents found in DNS for the turbulent shear stress and for the mixing layer thickness (such as from (3.5), which yields $\text{d}{\it\Delta}_{{\it\omega}}/\text{d}{\it\xi}\sim {\it\xi}^{-0.413}$ ) are in satisfactory agreement with (3.6). However, full self-similarity also requires similarity of the Reynolds stress balance equation, which can only be achieved if ${\it\Delta}\sim {\it\xi}$ , $u_{0}\sim \text{const.}$ (Townsend Reference Townsend1976). The present data then point to the attainment of a state of incomplete similarity, whereby equilibrium of the mean velocity field is reached prior to equilibration of the velocity fluctuations field. While the existence of such a state does not seem to be very surprising, to our knowledge it has not been theoretically predicted nor numerically studied before, and probably deserves some further investigation. Full similarity, hence saturation of the Reynolds stresses, is not reached in the present simulation as it probably requires a much longer fetch, and indeed the experiments of Bell & Mehta (Reference Bell and Mehta1990) showed evidence for saturation only at ${\it\xi}\gtrsim 150$ . Analysis of the Favre-averaged stresses shows differences of no more than 1 % with respect to the Reynolds averages, hence density variations do not seem to be responsible for the observed scalings. In any case, it is noteworthy (see table 2) that the mixing layer state (in terms of the Reynolds stress peaks) at the right-most station available is comparable with the asymptotic structure quoted in experiments (Bell & Mehta Reference Bell and Mehta1990; Goebel & Dutton Reference Goebel and Dutton1991).

Although the mixing layer growth rate is found to be significantly different from linear, we have nevertheless attempted to compare with existing experimental data, by locally fitting the vorticity thickness distribution in figure 4 with a straight line towards the end of the computational domain. The results of this exercise are summarized in table 2. As probably to be expected, the table highlights significant differences in the alleged asymptotic growth rates from incompressible experiments, as well as from experiments with the same flow conditions. In this respect we note that the generally accepted correlation for incompressible mixing layers (Brown & Roshko Reference Brown and Roshko1974)

by taking $C_{{\it\delta}}=0.165$ (Goebel & Dutton Reference Goebel and Dutton1991), would yield $\text{d}{\it\Delta}/\text{d}x|_{i}=0.0465$ . Hence, the present DNS would yield a growth rate reduction factor with respect to the incompressible case $(\text{d}{\it\Delta}/\text{d}x)/(\text{d}{\it\Delta}/\text{d}x)_{i}\approx 0.601$ , whereas the data of Goebel & Dutton (Reference Goebel and Dutton1991) yield $(\text{d}{\it\Delta}/\text{d}x)/(\text{d}{\it\Delta}/\text{d}x)_{i}\approx 0.816$ . Several (more or less empirical) formulae have been proposed to relate compressible growth rates to incompressible ones, one of the most well established being by Slessor, Zhuang & Dimotakis (Reference Slessor, Zhuang and Dimotakis2000), who stipulated the correction factor to be

where ${\it\alpha}\approx 4$ , ${\it\beta}\approx 0.5$ , and ${\it\Pi}_{c}$ is an ‘effective’ convective Mach number

thus also covering the case that the two streams have different specific heat ratios, ${\it\gamma}_{i}$ . Applied to the flow conditions under study here, ${\it\Pi}_{c}=0.52$ , which translates into a growth rate reduction factor $(\text{d}{\it\Delta}/\text{d}x)/(\text{d}{\it\Delta}/\text{d}x)_{i}=0.69$ , which is not too far from the present DNS. It should also be noted that the growth reduction factor observed here is also entirely consistent with the findings of Clemens & Mungal (Reference Clemens and Mungal1995), who reported $(\text{d}{\it\Delta}/\text{d}x)/(\text{d}{\it\Delta}/\text{d}x)_{i}=0.63$ at $\mathit{Ma}_{c}=0.42$ , and $(\text{d}{\it\Delta}/\text{d}x)/(\text{d}{\it\Delta}/\text{d}x)_{i}=0.56$ at $\mathit{Ma}_{c}=0.50$ .

Figure 9. Instantaneous contours of (a) streamwise velocity, (b) temperature and (c) pressure in the $x{-}y$ plane. 32 contour levels are shown in the respective ranges: $0.45<u/U_{1}<1.1$ , $0.65<T/T_{1}<1.2$ , $0.75<p/p_{1}<1.07$ . The colour scale ranges from blue (lower end) to red (higher end).

4. Flow organization

Contours of velocity, temperature and pressure from instantaneous flow fields are shown in figure 9, in a side view. The figure highlights well the slow process of mixing of the two streams which leads to strongly convoluted interfaces with the surrounding irrotational fluid. However, no clear sign of flow organization into coherent spanwise vortices (rollers) is observed, as in classical visualizations (Brown & Roshko Reference Brown and Roshko1974). This is a likely effect of flow compressibility which acts to suppress the Kelvin–Helmholtz instability. Indeed, in the experiments of Clemens & Mungal (Reference Clemens and Mungal1995), the formation of well-organized rollers at a similar convective Mach number was only observed at ${\it\xi}\approx 200$ , well downstream of the present computational outflow. The pressure field is found to be organized into small-sized blobs with positive/negative deviation from the free-stream value in the proximity of the trailing edge, which is a clear imprint of the small-sized vortices issued from the incoming boundary layers. Large-scale organization of the pressure field into alternating, positive and negative, delta-scaled blobs is observed further downstream, which is probably the precursor of the formation of rollers.

Top view visualizations are shown in figure 10, which indeed indicates inception of large-scale organization in the developing mixing layer. Specifically, the streamwise velocity organizes itself (as in any shear flow) into low- and high-speed streaks with spanwise spacing which monotonically grows from the trailing edge proximity (as inherited from the upstream near-wall streaks) to the computational outflow. On the other hand, the temperature and the pressure fields exhibit the signature of very wide flow structures, which eventually fill up the whole computational domain. Hints of spanwise organization are clear even at ${\it\xi}\approx 7{-}8$ , and large eddies inclined by about $70^{\circ }$ with respect to the free stream are observed until ${\it\xi}\approx 40$ . Similar diagonal patterns were also observed in the flow visualizations of Clemens & Mungal (Reference Clemens and Mungal1995), although at slightly higher Mach number ( $\mathit{Ma}_{c}\gtrsim 0.62$ ). At ${\it\xi}\approx 40$ , the eddies become more or less parallel to the spanwise direction. It is worth recalling that temporal stability of compressible laminar mixing layers (Sandham & Reynolds Reference Sandham and Reynolds1991) indicates the presence of both two-dimensional and three-dimensional unstable modes, the former expected to be dominant at $\mathit{Ma}_{c}\lesssim 0.6$ , and the latter at higher convective Mach numbers. The results of the present simulation (at $\mathit{Ma}_{c}\approx 0.45$ ) apparently support a similar scenario in which two-dimensional disturbances eventually occur, although strong amplification of oblique disturbances is observed in a transitory stage. Since a box sensitivity study is not feasible, these observations should be taken with some caution, being potentially affected by the spanwise size of the computational domain (here, $L_{z}/{\it\Delta}_{{\it\omega}}\approx 3.7$ at the end of the computational domain).

Figure 10. Instantaneous contours of (a) streamwise velocity, (b) temperature and (c) pressure in the $x{-}z$ geometric symmetry plane. 32 contour levels are shown in the respective ranges: $0.45<u/U_{1}<1.1$ , $0.65<T/T_{1}<1.2$ , $0.75<p/p_{1}<1.07$ . The colour scale ranges from blue (lower end) to red (higher end).

To quantify the growth of the size of the turbulent eddies, in figure 11 we provide the spanwise spectral densities of the streamwise velocity and pressure, taken at the mixing layer mid-plane. As already observed from the flow visualizations, the spectra show an overall shift to lower wavenumbers (i.e. larger scales) on progressing downstream, and the formation of an inertial layer in the $u$ spectrum, with exponent not far from $-5/3$ . On the other hand, the pressure spectra do not seem to exhibit the $k^{-7/3}$ inertial layer which would be expected in isotropic turbulence at sufficiently high Reynolds number (Batchelor Reference Batchelor1951). This is presumably due to the action of the mean shear which activates ‘rapid’ source terms related to mean shear/turbulence interaction in the Poisson equation for pressure (Blake Reference Blake1986), or to low-Reynolds-number effects, as the Taylor-scale Reynolds number here is no larger than 80, all the way to the end of the computational domain.

Figure 11. Spanwise spectral densities of streamwise velocity (a) and pressure (b) in the mixing layer symmetry plane, at ${\it\xi}=1.24$ , 6.37, 14.91, 25.16, 35.42, 45.67, 55.92 (the arrow indicating increasing $x$ ). The dot-dashed line in (a) denotes the trend $E(k)\sim k^{-5/3}$ , and in (b) the trend $E(k)\sim k^{-7/3}$ .

Figure 12. Pre-multiplied and normalized spanwise spectral densities of streamwise velocity ( $k_{z}\hat{E}_{u}$ , a,c) and pressure ( $k_{z}\hat{E}_{p}$ , b,d), as a function of the trailing edge distance. Wavenumbers are scaled with respect to $\overline{{\it\delta}}$ in (a,b), and with respect to the local vorticity thickness, ${\it\Delta}_{{\it\omega}}$ in (c,d). Contour levels are shown from 0.01 to 0.24, in intervals of 0.01, from light to dark grey shades.

Figure 13. Integral length scales of $u$ (solid line) and pressure (dashed line) as a function of trailing edge distance, scaled with respect to $\overline{{\it\delta}}$ (a) and with respect to the local vorticity thickness ${\it\Delta}_{{\it\omega}}$ (b).

The spectral densities are plotted in pre-multiplied form as a function of the distance from the trailing edge in figure 12. This type of representation has the advantage of clearly highlighting the size of the energy-containing eddies, and it has the property that equal areas underneath the curves correspond to equal energy, when a semi-logarithmic representation is used. Further, in order to scale out the effect of variation of the amplitude of the disturbances with ${\it\xi}$ , the spectra are scaled with respect to the total energy, i.e. $\hat{E}_{u}(k)=E_{u}(k)/\overline{u^{\prime 2}}$ . The normalized velocity and pressure spectra are shown in figure 12(a,b) as a function of $k\overline{{\it\delta}}$ , and confirm the scale change effect previously noticed in figure 11, in that the spectral peak shifts to larger scales with the downstream coordinate. Apparently, collapse of the pre-multiplied spectra is found if the wavenumber is scaled with respect to the local vorticity thickness, as shown in figure 12(c,d). In this representation, the spectral ridge of $u$ is very nearly constant starting from ${\it\xi}\approx 10$ , whereas the peak of $p$ is slightly shifted to higher wavenumbers. This is further confirmed from inspection of the integral length scales, as shown in figure 13, which highlight well the proportionality of the velocity length scale to the vorticity thickness (with ${\it\Lambda}_{z}/{\it\Delta}_{{\it\omega}}\approx 0.2$ ), and to a lesser extent of the pressure scale (with ${\it\Lambda}_{z}/{\it\Delta}_{{\it\omega}}\approx 0.6$ , still decreasing). It is noteworthy that the observed integral length scale of the $u$ -containing eddies is comparable to that observed in the outer part of turbulent boundary layers, in which ${\it\Lambda}_{z}/{\it\delta}\approx 0.2{-}0.3$ (Pirozzoli Reference Pirozzoli2012), which probably points to a universal scaling law for the momentum streaks in all shear flows, either free or wall bounded. Finally, the enhanced spanwise coherence of the pressure field as compared to the velocity field is consistent with the flow visualizations of figure 10.

5. Dynamic mode decomposition

The DMD (Schmid Reference Schmid2010) is a technique of spectral analysis which operates on a sequence of snapshots $\boldsymbol{u}(t_{n})$ obtained from either experiments or, as in our case, simulations to extract information based on Koopman’s analysis of dynamical systems (Rowley et al. Reference Rowley, Mezic, Bagheri, Schlatter and Henningson2009). Introducing a linear time advancement operator $\unicode[STIX]{x1D63C}$ as $\boldsymbol{u}(t_{n+1})=\unicode[STIX]{x1D63C}\boldsymbol{u}(t_{n})$ , the DMD algorithm extracts approximate eigenvalues and eigenmodes of a reduced-order representation of $\unicode[STIX]{x1D63C}$ . While having some connection with proper orthogonal decomposition (POD), DMD produces modes associated with a single frequency, whereas POD modes may in general have broadband frequency content. Furthermore, for linear systems, the DMD modes coincide with the global modes, and for perfectly time-periodic data the DMD modes coincide with the Fourier modes. The essence of the DMD can be viewed as follows, with more complete and rigorous descriptions provided in Rowley et al. (Reference Rowley, Mezic, Bagheri, Schlatter and Henningson2009) and Schmid (Reference Schmid2010). Snapshots of data are organized as columns into a matrix $\unicode[STIX]{x1D61D}_{1}^{N}=\{v_{1},v_{2},\dots ,v_{N}\}$ (of size $M\times N$ , say), where two subsequent snapshots $v_{j}$ and $v_{j+1}$ are separated by a constant time interval ${\rm\Delta}t$ . Following the idea underlying Krylov subspace techniques, and in particular the Arnoldi method, it is possible to approximate the linear mapping $\unicode[STIX]{x1D63C}$ from one snapshot to the next using an $(N-1)$ -dimensional snapshot basis

or in matrix form

where the reduced-order matrix $\unicode[STIX]{x1D64E}$ is of companion type. Its eigenvalues (Ritz values) approximate some of the eigenvalues of $\unicode[STIX]{x1D63C}$ and the associated eigenvectors provide the coefficients of the linear combination that is needed to express the modal structure within the snapshot basis. $\unicode[STIX]{x1D64E}$ is easily estimated through a least-square procedure, by performing a $QR$ decomposition of $\unicode[STIX]{x1D61D}_{1}^{N-1}$ to solve (5.2), thus obtaining $\unicode[STIX]{x1D64E}=\unicode[STIX]{x1D619}^{-1}Q^{T}\unicode[STIX]{x1D61D}_{2}^{N}$ .

The DMD algorithm has been initially applied to temporal snapshots of the pressure field taken in an $x{-}y$ plane (at $z=L_{z}/2$ ), in a subdomain containing the final part of the splitter plate and the initial part of the mixing layer, extending from ${\it\xi}=-2$ to ${\it\xi}=11.6$ . A total of 300 snapshots, with time separation ${\rm\Delta}t=0.2\overline{{\it\delta}}/{\rm\Delta}\unicode[STIX]{x1D650}$ , have been processed. The Ritz values ( ${\it\lambda}$ ) associated with the resulting reduced-order system $\unicode[STIX]{x1D64E}$ are shown in figure 14 in the complex plane. Essentially all eigenvalues lie on the unit circle, which is to be expected as the system is statistically stationary. The eigenvalues of $\unicode[STIX]{x1D64E}$ are transformed via the mapping ${\it\omega}=\log ({\it\lambda}/{\rm\Delta}t)$ , and shown in figure 15, where the real and imaginary parts represent the exponential growth or decay rate and the temporal frequency, respectively. Since the processed data are real, all eigenvalues are real or in complex-conjugate pairs, and the associated DMD spectrum is symmetric about the zero frequency. The size and colouring of the symbols reported in the figure are indicative of the coherence of the associated modes, which is useful to distinguish large-scale energetic structures from smaller-scale, less-energetic ones. Coherence is gauged by projecting each dynamic mode onto the POD basis ( $\unicode[STIX]{x1D650}$ ) computed from the singular value decomposition of the data sequence $\unicode[STIX]{x1D61D}_{1}^{N-1}=\unicode[STIX]{x1D650}{\it\Sigma}\unicode[STIX]{x1D652}$ . This decomposition provides the matrices $\unicode[STIX]{x1D650}$ and $\unicode[STIX]{x1D652}$ containing the spatial and temporal structures respectively, whose energy ranking is given by the singular values ${\it\Sigma}$ (Schmid Reference Schmid2010).

Figure 14. Ritz values of the reduced system detected by DMD for temporal snapshots of the pressure field in an $x{-}y$ plane between ${\it\xi}=-2$ and ${\it\xi}=11.6$ .

Figure 15. Spectrum of eigenvalues from the DMD algorithm applied to the pressure field in an $x{-}y$ plane between ${\it\xi}=-2$ and ${\it\xi}=11.6$ . The colour and symbol size of the eigenvalues indicate the coherence of the associated dynamic modes (more coherent modes are rendered with larger symbols and in shades of red).

The pressure fields associated with three representative dynamic modes (labelled in figure 15) are displayed in figure 16. Besides the one stationary mode corresponding to the mean field (not shown), the most coherent mode (a) is characterized by large blobs with alternate positive and negative pressure, covering approximately three oscillatory periods over the streamwise extent of the interrogation window, and growing in the streamwise direction. Interestingly, these organized modes arise immediately past the end of the splitter plate, and they are apparently not inherited from any large-scale structure embedded in the upstream boundary layers. A similar pattern is also retrieved in mode (b), which is characterized by structures with comparable scale, but leaning in the streamwise direction. Mode (c) in figure 16 has a similar shape as mode (a); however its associated eigenvalue has lower frequency, and the associated spatial scale is correspondingly larger. Furthermore, it is characterized by a spatial extent in the ${\it\eta}$ direction which is less than for mode (a).

Figure 16. Leading dynamic modes in an $x{-}y$ plane, with snapshots taken in the interval $-2\leqslant {\it\xi}\leqslant 11.6$ , corresponding to the labelled eigenvalues in figure 15. Contours of modal pressure in colour scale from blue (negative values) to red (positive values).

Figure 17. Leading dynamic modes in the $y=0$ plane, with snapshots taken in the interval $-2\leqslant {\it\xi}\leqslant 11.6$ . Contours of modal pressure in colour scale from blue (negative values) to red (positive values).

The same analysis has been repeated for snapshots taken in the geometric mid-plane in the same streamwise interval, and the shape of the leading modes is shown in figure 17. It should be noted that, since the DMD algorithm has been applied plane-wise, the modes in figure 17 do not necessarily correspond to those of figure 16. Here the most salient feature, clearly noticeable in the leading modes, is the appearance of substantial spanwise coherence of the pressure field emerging without apparent connection with the upstream turbulent boundary layers. At this stage, it appears that the leading pressure-containing eddies are strongly undulated in the spanwise direction. These eddies may be interpreted as the precursors of the formation of mixing layer rollers, which have been observed in mixing layers at low/moderate values of the convective Mach number (Brown & Roshko Reference Brown and Roshko1974; Clemens & Mungal Reference Clemens and Mungal1995). It is noteworthy that, while their presence cannot immediately be inferred from the instantaneous visualizations reported in figures 9 and 10, they are clearly detected by the DMD algorithm, which identifies these structures as the most significant from the dynamical point of view. To analyse the downstream evolution of these large-scale structures, the DMD analysis in the mid-plane has been repeated for snapshots in the intermediate region $12\leqslant {\it\xi}\leqslant 32$ (figure 18) and in the far-field region $32\leqslant {\it\xi}\leqslant 56$ (figure 19). It becomes evident that in the intermediate region the pressure eddies gain further spanwise coherence, organizing themselves into diagonal bands filling the whole computational box. Further downstream, the bands become almost aligned along the spanwise direction, and the spanwise modulation is much reduced. Correspondingly, the modal pressure in the $x{-}y$ plane exhibits a pattern typical of mixing layer rollers. The DMD modes here appear to provide a condensed form of the information contained in the instantaneous flow visualization of figure 10.

Figure 18. Leading dynamic mode in the $y=0$ plane, with snapshots taken in the interval $12\leqslant {\it\xi}\leqslant 32$ . Contours of modal pressure in colour scale from blue (negative values) to red (positive values). The solid lines are inclined by $70^{\circ }$ with respect to the streamwise direction.

Figure 19. Leading dynamic mode in the $x{-}y$ plane (a) and in the $y=0$ plane (b), with snapshots taken in the interval $32\leqslant {\it\xi}\leqslant 56$ . Contours of modal pressure in colour scale from blue (negative values) to red (positive values).