2.1 Configuration of the computational grid

An illustration of the computational domain is shown in figure 1. In determining the extent of the horizontal and vertical channels upstream of the entry of the jet, we have followed the guideline of Muppidi & Mahesh (Reference Muppidi and Mahesh2007). According to this, the distance of an inflow from the T-junction should be at least equal to the width of the corresponding channel, so that the velocity statistics in the vicinity of the T-junction remain unaffected by the inflow. Based on this guideline, the distance between the T-junction and the inflow of the horizontal channel is set to $5\unicode[STIX]{x1D6FF}$ , whereas the extent of the vertical channel in the $y$ -direction is set to $2\unicode[STIX]{x1D6FF}$ . Further, the extent of the horizontal channel downstream of the entry of the jet is set equal to $19\unicode[STIX]{x1D6FF}$ so as to capture as far as possible the mixing process between the two streams.

In the spanwise $z$ -direction, the incoming jet extends over the entire computational domain, and the flow in this direction is homogeneous. The size of the computational domain in the $z$ -direction should be large enough so as to avoid autocorrelation of the flow quantities. On the other hand, it has to be reasonably limited so as to keep the computational cost of the simulations at acceptable levels. For example, in their LES of flows in channels with periodic hills, Fröhlich et al. (Reference Fröhlich, Mellen, Rodi, Temmerman and Leschziner2005) opted for a size of $3\unicode[STIX]{x1D6FF}$ . Also, Flores & Jimenez (Reference Flores and Jimenez2010) suggested $3\unicode[STIX]{x1D6FF}$ as the minimum spanwise size for channel flow. This suggests that a similar size may be sufficient for the purposes of our study. However, in our case the incoming jet is laminar, which could lead to more persistent two-dimensional structures. For this reason, and in order to gain further insight into the size requirements, we conducted LES with 3 different spanwise sizes, $\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ , $(3/2)\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ and $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ , and compared the numerical results. As it turned out, the differences between the first- and second-order statistics of the flow quantities were very small. On the basis of these considerations, the spanwise size of the computational domain in our study is set to $\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ , which is deemed as a satisfactory compromise between required accuracy and computational cost.

Figure 1.

Illustration of the computational domain for the T-junction simulation. For visualisation purposes, the drawing is not to scale.

Given these dimensions, the computational grid for the DNS consists of $44.16\times 10^{6}$ computational nodes. Due to the significant flow inhomogeneities, in the streamwise and wall-normal directions we have employed non-uniform grid spacing. On the other hand, a uniform grid spacing is used in the homogeneous $z$ -direction. Henceforth, all lengths will be given in non-dimensional values, that is, normalised by $\unicode[STIX]{x1D6FF}$ . Finally, the origin of the axes is placed at the bottom wall, opposite to the downstream corner of the entry of the jet.

The computational grid can be divided in three sections. The first one covers the horizontal channel upstream of the entry of the jet, $-6<x\leqslant -1$ . In this section we have placed 100 computational nodes in the streamwise direction. More specifically, in the zone $-6<x\leqslant -2.78$ there are 36 nodes with constant spacing $\unicode[STIX]{x0394}x=0.091$ . Then, the grid is progressively refined, with a mild compression ratio of 0.95, up to $\unicode[STIX]{x0394}x=0.004$ . The second section is defined by $-1<x\leqslant 0$ and includes the vertical channel. Therein we have placed 100 nodes in the streamwise direction with a symmetric distribution with respect to the mid-plane of the vertical channel, $x=-0.50$ . For $-1<x\leqslant -0.5$ the grid spacing increases smoothly with a stretching factor of 1.03, up to $\unicode[STIX]{x0394}x=0.02$ . The third section covers the computational domain downstream of the entry of the jet and is discretised by 900 nodes in the streamwise direction. More specifically, in the interval $0<x\leqslant 0.9$ we have placed 85 nodes and the spacing increases with a constant stretching factor of 1.02 up to $\unicode[STIX]{x0394}x=0.022$ . Beyond the station $x=0.9$ , the streamwise meshing is uniform.

In the wall-normal direction, we have placed 300 cells across the horizontal channel, 140 nodes above and 160 nodes below the mid-plane. The choice of asymmetric meshing in the (wall-normal) $y$ -direction is motivated by the fact that in some regions of the bottom wall, the thickness of the boundary layer decreases considerably due to flow acceleration; see figure 7 and the relevant discussion below. Moreover, the fluid at the top wall is colder than the cross-flow and, therefore, has higher shear viscosity. Consequently, the computational requirements at the top wall are not as strict as those at the bottom wall. For these reasons, we opted for higher resolution at the bottom half of the horizontal channel. The centre of the first cell from the bottom wall, $y=0$ , is set at $y_{f.n.}=8.98\times 10^{-4}$ . In the bottom half of the horizontal channel, the spacing increases smoothly, with a constant stretching factor of 1.013, up to $\unicode[STIX]{x0394}y=0.015$ . The remaining 140 nodes are placed at the top half of the channel; therein the grid is refined smoothly as we approach the top wall with a constant stretching factor of 1.012. The first node off the top wall is located at $y_{f.n.}=1.33\times 10^{-3}$ . Along the vertical channel, $2<y<4$ , we have placed 150 cells. The grid spacing increases smoothly with $y$ , with a constant stretching factor of 1.018, up to $\unicode[STIX]{x0394}y=0.0195$ at the top inflow boundary.

Finally, in the spanwise direction, the domain is discretised via 128 uniformly distributed nodes, i.e.  $\unicode[STIX]{x0394}z=0.0245$ . An overview of the discretisation of the computational domain is provided in table 2.

Table 2.

Overview of the discretisation of the computational domain.

In determining the grid spacing in the near-wall regions we have followed the guidelines of Le, Moin & Kim (Reference Le, Moin and Kim1997) according to which at least one node should be placed in the interval $0<y^{+}<1$ , while the grid spacing in the streamwise and spanwise directions should be $\unicode[STIX]{x0394}x^{+}\approx 10$ and $\unicode[STIX]{x0394}z^{+}\approx 15$ , respectively. In figure 2 we show the position of the first node in the vertical direction and the grid spacing (in wall units) at the walls of the horizontal channel. More specifically, from figure 2(a) we can infer that at least two nodes are always placed in the interval $0<y^{+}<1$ . At the bottom boundary, $y_{f.n.}^{+}$ attains its maximum value of 0.44 at $x\approx 1$ and decreases significantly further downstream. In particular, for $x>3.0$ , $y_{f.n.}$ becomes smaller than 0.3. At the top boundary, $y_{f.n.}^{+}$ is always smaller than 0.3. Further, according to figure 2(b), the coarsest grid spacings at the bottom boundary are $\unicode[STIX]{x0394}x^{+}\approx 10.4$ and $\unicode[STIX]{x0394}z^{+}\approx 12.0$ , and they occur at the station $x=1$ . Further downstream, both grid spacings decrease and are consistently smaller than 8.0. At the top wall, the values of $\unicode[STIX]{x0394}x^{+}$ and $\unicode[STIX]{x0394}z^{+}$ are even lower and always smaller than 4.0. We therefore conclude that our computational mesh fulfils the above guidelines of Le et al. (Reference Le, Moin and Kim1997) for the near-wall regions.

Figure 2.

Grid spacing at the walls of the horizontal channel expressed in wall units. (a)  $y_{f.n.}^{+}$ , (b)  $\unicode[STIX]{x0394}x^{+}$ and $\unicode[STIX]{x0394}z^{+}$ . At the top wall and at $x\approx 0.5$ , $x\approx 1.4$ and $x\approx 3.9$ , the grid spacing drops to zero because the wall shear stress $\unicode[STIX]{x1D70F}_{w}$ vanishes due to reattachment of the flow.

For the assessment of the quality of the grid away from the walls, we compare the local grid spacing, $\unicode[STIX]{x0394}=(\unicode[STIX]{x0394}x\,\unicode[STIX]{x0394}y\,\unicode[STIX]{x0394}z)^{1/3}$ with the smallest length scales of the flow, namely, the Kolmogorov scale

(2.10) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}=\left(\frac{\unicode[STIX]{x1D708}^{3}}{\unicode[STIX]{x1D716}_{k}}\right)^{1/4}, & & \displaystyle\end{eqnarray}$$

and the temperature microscale $\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D703}}$ . For Prandtl numbers significantly higher than unity, $\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D703}}$ is given by the following relationship (Tennekes & Lumley Reference Tennekes and Lumley1972),

(2.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D702}\left(\frac{1}{Pr}\right)^{1/2}. & & \displaystyle\end{eqnarray}$$

However, in the flow under study, the reference Prandtl number is close to unity; $Pr=1.44$ , which corresponds to $\unicode[STIX]{x1D702}_{\unicode[STIX]{x1D703}}=0.83\,\unicode[STIX]{x1D702}$ . For this reason, we can examine the quality of the grid on the basis of the Kolmogorov scale $\unicode[STIX]{x1D702}$ . In figure 3 we present contour plots of the ratio $\unicode[STIX]{x0394}/\unicode[STIX]{x1D702}$ downstream of the entry of the jet for the DNS grid. According to this figure, the ratio $\unicode[STIX]{x0394}/\unicode[STIX]{x1D702}$ in our grid is close to 2, with a maximum value equal to 2.81. We therefore conclude that our grid resolution for the DNS is sufficient for the purposes of our study.

Finally, for the configuration of the computational mesh of our LES we have followed the guidelines of Zang (Reference Zang1991) so as to fully resolve the wall regions. The stretching factors above and below the entry of the jet are 1.1 and 1.07, respectively. The computational mesh of the LES consists of approximately $2\times 10^{6}$ nodes; in other words, it is 22 times coarser than the grid of the DNS. Although the LES mesh resolves the near-wall regions, it does not resolve all turbulent scales away from them; in fact, the ratio $\unicode[STIX]{x0394}/\unicode[STIX]{x1D702}$ for the LES mesh is an order of magnitude larger than for the mesh employed in our DNS.

Figure 3.

Contour plots of the ratio between the local grid spacing for the DNS and the Kolmogorov length scale, $\unicode[STIX]{x0394}/\unicode[STIX]{x1D702}$ , downstream of the entry of the jet.

2.2 Boundary conditions

As regards the implementation of boundary conditions, we have to prescribe a turbulent inflow at the inlet of the cross-flow. In our study, and for purposes of increased accuracy, we produce turbulent inflow data by performing an auxiliary simulation of channel flow simultaneously with the main simulation. In this auxiliary simulation, the streamwise dimension of the computational domain is $4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}$ , while the wall-normal and spanwise dimensions are set equal to the ones of the horizontal channel of the T-junction. To avoid spatial interpolation errors, the grid resolution is the same as in the corresponding inlet of the T-junction. The mesh of this auxiliary simulation consists of $5.76\times 10^{6}$ computational cells. Also, the same time step is used in the main and auxiliary simulations so as not to introduce time interpolation errors. An analogous procedure is followed for generating turbulent inflow data for the corresponding LES. In particular, we perform an auxiliary wall-resolved LES of channel flow on a grid consisting of $4.3\times 10^{5}$ cells. Also, since the incoming jet is laminar, at the top inflow of the vertical channel we impose a parabolic profile whose bulk velocity is equal to unity.

At the outflow boundary, we prescribe a convective boundary condition to all velocity components because it allows the exit of vortical structures from the computational domain with minimal distortion. When applied to a scalar field $\unicode[STIX]{x1D719}$ at an outflow plane normal to the $x$ -direction, this condition reads

(2.12) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+U_{conv}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}x}=0, & & \displaystyle\end{eqnarray}$$

where the convective velocity $U_{conv}$ is set equal to the bulk velocity at the outflow plane.

As regards the temperature field, we prescribe uniform profiles at the two inflows and the convective boundary condition (2.12) at the outflow plane. The walls of the T-junction are considered adiabatic and, therefore, a zero Neumann condition for the temperature is applied at them. Finally, since the flow is homogeneous in the $z$ -direction, we prescribe periodic boundary conditions in this direction.

3.1 Instantaneous flow structures

In this subsection we examine the instantaneous flow structures, as extracted from the DNS data. In our analysis we focus on the primary and secondary shear layers and assess their contribution to the mixing process in the T-junction. We also investigate the effect of the incoming jet on the flow structures near the top and bottom walls. All the figures presented herein show results at the same time instance.

In figure 8 we provide contour plots of the amplitude of the instantaneous spanwise vorticity $|\unicode[STIX]{x1D714}_{z}|$ and normalised temperature $\unicode[STIX]{x1D703}$ ,

(3.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}=\frac{T-T_{jet}}{T_{ref}-T_{jet}}. & & \displaystyle\end{eqnarray}$$

The results shown in this figure are taken at the $xy$ -plane defined by $z=1.24$ .

Figure 8.

Contour plots of the amplitude of the instantaneous spanwise vorticity $|\unicode[STIX]{x1D714}_{z}|$ and normalised temperature $\unicode[STIX]{x1D703}$ at $z=1.24$ . (a)  $|\unicode[STIX]{x1D714}_{z}|$ , (b)  $\unicode[STIX]{x1D703}$ . In (a) we observe the formation of a large spanwise roller at the secondary shear layer at $x\approx 1.8$ .

Both shear layers roll up quickly, leading to the formation and shedding of spanwise rollers such as the one spotted at the secondary layer and at $x\approx 1.8$ , see figure 8(a). These rollers merge with adjacent ones in a pairing process, or they are torn apart and their vorticity is absorbed by the surrounding vortical structures. Most of the surviving rollers are convected downstream but some of them interact with the large recirculation bubble and eventually get entrapped therein. Also, as the two shear layers spread, they start to interact and, eventually, merge. Downstream of the reattachment point of the flow, the merged shear layer continues to spread until it extends across the entire channel. In fact, downstream of the station $x=6.0$ , many of its rollers interact with the bottom wall. Additionally, we can observe numerous thin and elongated vortical structures close to the bottom wall for $4<x<7$ . Their formation is attributed to the deceleration of the flow caused by the adverse pressure gradient in this region, see figure 7(b).

In regard to the instantaneous temperature field, shown in figure 8(b), we remark that thermal mixing is initiated by the spanwise rollers of the primary shear layer close to the upstream corner of the entry of the jet. The mixing process is significantly enhanced once the two shear layers start to interact with one another. Also, downstream of the large recirculation bubble, we can observe the effect of the surviving rollers as they transport hot fluid originating from the cross-flow toward the top wall. A video animation of the evolution of the temperature field in the $(xy)$ plane is provided as supplementary movie 1 (available at https://doi.org/10.1017/jfm.2018.256).

In figure 9 we show isosurface plots of the fluctuating pressure at the value $p^{\prime }=-0.25$ . Such plots are frequently used for the identification of large coherent structures. In our case, the value $p^{\prime }=-0.25$ has been chosen because it provides a clear visualisation of these structures. The growth of the spanwise rollers from the secondary shear layer can be inferred from the three-dimensional view of figure 9(a). These rollers are convected downstream and eventually are torn apart into smaller structures in a highly turbulent environment. Also, from figure 9(b) which shows a projection of the isosurface plot of $p^{\prime }$ at the $xz$ -plane, we can observe the formation of streamwise vortex tubes, known as ‘rib’ vortices, that connect two adjacent spanwise rollers. On the other hand, the smaller streamwise vortex tubes close to the upstream corner of the entry of the jet, $x=0.0$ , are part of the primary shear layer and eventually break up into smaller ones.

Figure 9.

Isosurface plots of the fluctuating pressure at the value $p^{\prime }=-0.25$ . (a) Three-dimensional view, (b) view from above ( $xz$ -plane). The large spanwise rollers shown in these plots are shed by the secondary shear layer. Also, the streamwise vortex tubes connecting two consecutive rollers correspond to ‘rib’ vortices.

In order to elucidate the role of these streamwise vortex tubes to the mixing process, in figure 10 we present contour plots of $\unicode[STIX]{x1D703}$ at two streamwise stations, namely, at $x=0.0$ and $x=3.0$ . At the first station, these tubes are located at the primary shear layer and are part of mushroom-like structures that are formed right at the interface of the two streams. A mushroom moving upwards is formed by fluid coming from the cross-flow and wraps around itself cold fluid from the incoming jet, thus forming a spiral (Lin & Corcos Reference Lin and Corcos1984). Inversely, a mushroom moving downwards is formed by fluid coming from the jet and wraps around itself hot fluid from the cross-flow. Such mushroom-like structures are typical of turbulent mixing layers (Corcos & Lin Reference Corcos and Lin1984; Bernal & Roshko Reference Bernal and Roshko1986; Rogers & Moser Reference Rogers and Moser1992) and have also been observed experimentally in the shear layers of T-junctions (Hirota et al. Reference Hirota, Mohri, Asano and Goto2010). A video animation of the evolution of the temperature field at this streamwise station is provided as supplementary movie 2. At the second station, $x=3.0$ , the primary and secondary shear layers have started to interact and these mushroom-like structures have already broken up. Instead, the vorticity field therein is dominated by the large spanwise rollers of the secondary shear layer.

Figure 10.

Contour plots of the instantaneous normalised temperature $\unicode[STIX]{x1D703}$ at two streamwise locations. (a)  $x=0.0$ , (b)  $x=3.0$ . In (a) we can observe the formation of mushroom-like structures of streamwise vorticity at the primary shear layer. However at $x=3.0$ , (b) these structures have already been broken up.

Our analysis of the instantaneous flow structures concludes with the examination of the near-wall regions. In figure 11 we present contour plots of the instantaneous streamwise velocity $u$ close to the walls of the horizontal channel. In the vicinity of the reattachment point of the flow, vortical structures emanating from the shear layers interact with the top wall, leading to the formation of several irregularly shaped streaks. On the other hand, further downstream, the adverse pressure gradient results in the widening and weakening of the velocity streaks. This feature is typical of regions of adverse pressure gradient; see, for example the discussion in Skote & Henningson (Reference Skote and Henningson2002) and Lee et al. (Reference Lee, Hu, Saha and Kazimi2009). At the bottom wall and in the region $4<x<5$ , there are patches of negative streamwise velocity. This is also due to the adverse pressure gradient which causes the flow to decelerate. It should be noted, however, that this concerns instantaneous values, whereas the mean streamwise velocities remain positive close to the bottom wall. Also, further downstream, we identify zones with numerous velocity streaks. These streaks are much longer and thinner than the ones observed near the top wall and are caused by the interaction of vortices from the merged shear layer with the bottom wall.

Figure 11.

Contour plots of the instantaneous streamwise velocity $u$ close to the two walls of the horizontal channel. (a) Top wall, $y=1.98$ , (b) bottom wall, $y=0.02$ .

Finally, in figure 12 we show a representation of the instantaneous vorticity field at the walls, based on the quantity $Q$ . This quantity, introduced by Hunt, Wray & Moin (Reference Hunt, Wray and Moin1988), measures the difference between the amplitudes of the symmetric and skew-symmetric parts of the velocity-gradient tensor, and it aids the visualisation of the vorticity field. In particular, positive values of $Q$ correspond to non-zero vorticity amplitude. Also, the higher the value of $Q$ is, the higher the vorticity amplitude gets. In figure 12 we present plots of the isosurface of the value $Q=4.0$ . At the top wall, we can observe numerous irregularly oriented vortical structures in the region covered by the large recirculation bubble. Downstream of the reattachment of the flow, the density of these structures gradually decreases. The dense patches of interwoven vorticity that can be observed in this region coincide with the widened velocity streaks observed in figure 11(a).

According to figure 12(b), near the bottom wall there are very few vortical structures for $x<4$ . This is the region where the cross-flow accelerates due to the Venturi effect. Downstream of this region, we observe a large number of vortical structures. Some of them are formed due to the adverse pressure gradient caused by the expansion of the cross-flow, while others are transported by the merged shear layer which spreads and interacts with the bottom wall. We also remark that, for $x>4$ , the density of the vortical structures at the bottom wall is higher than at the top wall. This difference may be attributed to two factors. The first one is that the local velocity close to bottom top wall is more than twice the local velocity close to the top one. The second factor is the lower temperature and, therefore, higher viscosity of the fluid at the top which enhances dissipation of the vortical structures.

Figure 12.

Isosurface plots of the $Q$ -criterion ( $Q=4.0$ ) at the wall regions. (a) Top wall, (b) bottom wall.

3.2 Flow statistics

We first present the first- and second-order statistics of the temperature field and also compare the LES predictions with the DNS data. Figure 13 shows the profiles of the averaged normalised temperature $\langle \unicode[STIX]{x1D703}\rangle$ at different streamwise stations. The main feature of these profiles is the transition zone from the high temperature at the bottom to the low temperature at the top. Close to the entry of the jet, this transition takes place across the primary shear layer. Also, at the stations $x=1.0$ and $x=3.0$ , i.e. at the large recirculation bubble, the temperature profile is not monotonic but attains a minimum at the core of the secondary shear layer. In other words, the temperature in the recirculation bubble is higher than in the secondary layer. This implies that the bubble does not consist solely of fluid from the incoming jet but it also contains fluid from the hot cross-flow.

Further downstream, the two shear layers merge and a monotonic temperature profile is established due to thermal mixing. As the flow evolves and thermal mixing progresses, the temperature gradients become weaker. Also, the temperature profile becomes essentially linear, as can be asserted from figure 13(f). Since the solid walls are adiabatically isolated, a channel flow of uniform temperature will eventually emerge; however, thermal equilibrium is not established within the limits of our computational domain.

Figure 13.

Profiles of the mean normalised temperature $\langle \unicode[STIX]{x1D703}\rangle$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES.

We also remark that the LES predictions agree well with the DNS data and their differences are minor. For example, the LES slightly overpredicts the temperature inside the recirculation bubble, while further downstream at ( $x=12.0$ and $x=16.0$ ) there is a small discrepancy in the prediction of the temperature at the bottom wall.

The profiles of the r.m.s. normalised temperature $\unicode[STIX]{x1D703}_{rms}$ are shown in figure 14. Their most characteristic feature is the peak due to the primary shear layer. At $x=1.0$ this peak is very sharp. However, as the shear layer spreads out, the value of the peak drops and its location moves closer to the bottom wall. In turn, this implies that thermal mixing is most intense in the region between $x=0.0$ and $x=4.0$ . In this figure we also observe that the predictions of the LES agree very well with those of the DNS, albeit the LES mildly underpredicts the amplitudes of the peaks of $\unicode[STIX]{x1D703}_{rms}$ at $x=1.0$ and $x=3.0$ .

Figure 14.

Profiles of the r.m.s. of the normalised temperature $\unicode[STIX]{x1D703}_{rms}$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES.

We proceed with the presentation of the velocity statistics. The profiles of the mean streamwise velocity $\langle u\rangle$ are shown in figure 15. The negative values at the top wall for $x=1.0$ and $x=3.0$ occur inside the recirculation bubble. Also, the high wall-normal gradients of $\langle u\rangle$ at these streamwise stations correspond to the secondary shear layer. Downstream of the reattachment of the flow, the mean velocity converges to a channel-flow profile, albeit at a slow rate as can be attested by the asymmetry of the profiles shown in figure 15(d–f). We also remark that the LES predictions of $\langle u\rangle$ compare quite favourably with the DNS data at all streamwise stations.

Figure 15.

Profiles of the mean streamwise velocity $\langle u\rangle$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES.

In figure 16 we present plots of the r.m.s. profile of the streamwise velocity, $u_{rms}$ . The maximum at $x=-0.50$ occurs at the core of the primary shear layer. At $x=1.0$ , due to the bending of the two streams, this shear layer is located lower down. The peak r.m.s. inside the shear layer occurs right below the mid-plane of the channel and its value is decreased due to the spreading of the layer. On the other hand, the strong peak at $y\approx 1.5$ corresponds to the secondary shear layer. We also observe that at $x=3.0$ the values of $u_{rms}$ have increased, which is attributed to turbulence generation inside the recirculation bubble and near the bottom wall due to the adverse pressure gradient. Further downstream, $u_{rms}$ becomes almost constant away from the walls and its profile starts to resemble the one in channel flows. By contrast, in the corresponding gas–gas flow, the wall-normal variations of $u_{rms}$ decay very slowly (Georgiou & Papalexandris Reference Georgiou and Papalexandris2017b ). This difference is due to the fact that the flow under study involves liquids and, therefore, buoyancy effects are less significant.

Figure 16.

Profiles of the r.m.s. of the streamwise velocity $u_{rms}$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES.

Figure 17.

Profiles of the r.m.s. of the wall-normal velocity $v_{rms}$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES.

Figure 18.

Profiles of the r.m.s. of the spanwise velocity $w_{rms}$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES.

The r.m.s. profiles of the wall-normal and spanwise velocities are shown in figures 17 and 18, respectively. At $x=-0.5$ , the peak of these profiles corresponds to the primary shear layer. Also, at the streamwise station $x=1.0$ , both $v_{rms}$ and $w_{rms}$ profiles exhibit a double-peak structure. The upper peak corresponds to the secondary and the lower peak to the primary shear layer. This structure was also observed in our simulations of gas–gas flow in T-junctions (Georgiou & Papalexandris Reference Georgiou and Papalexandris2017b ) but not in isothermal flows (Georgiou & Papalexandris Reference Georgiou and Papalexandris2017a ). We may therefore infer that this double peak is due to the role of the temperature as an active scalar. More specifically, the primary shear layer enhances significantly the mixing between the cross-flow and the incoming jet when these two streams are not at the same temperature. On the other hand, the large r.m.s. values inside the recirculation bubble, figures 17(c) and 18(c), are a common feature of both isothermal and variable-temperature flows in T-junctions. Further, from figures 16–18 we infer that the LES captures fairly well the r.m.s. profiles of all velocity components. Nonetheless, there are some discrepancies; most notably, it underpredicts the r.m.s. values inside the recirculation bubble by up to 25 %.

Next, in figure 19 we present results for the Reynolds shear stresses $\langle u^{\prime }v^{\prime }\rangle$ . At the streamwise station $x=1.0$ , they assume an S-shaped profile across the shear layers. Similar profiles have been predicted for shear layers at interfaces between pure fluids and porous media (Antoniadis & Papalexandris Reference Antoniadis and Papalexandris2015, Reference Antoniadis and Papalexandris2016). Herein, the large negative peak occurs at the core of the primary shear layer, whereas the large positive peak occurs at the core of the secondary layer. Further, the secondary peak at $x=1.0$ close to the top wall ( $y\approx 1.75$ ) is due to the secondary recirculation zone located next to the downstream corner of the entry of the jet, cf. figure 4.

Figure 19.

Profiles of the Reynolds stress $\langle u^{\prime }v^{\prime }\rangle$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES. For purposes of clarity, the limits of the horizontal axis in (c) and (d) are different than in the other figures.

Our discussion of the flow statistics concludes with the results for the turbulent heat fluxes $\langle \unicode[STIX]{x1D703}^{\prime }u^{\prime }\rangle$ and $\langle \unicode[STIX]{x1D703}^{\prime }v^{\prime }\rangle$ . Their profiles at different streamwise stations are shown in figures 20 and 21, respectively. We observe that for $x<3$ , the peaks of these profiles occur at the core of the primary shear layer. In other words, the primary shear layer provides the main mechanism for thermal mixing. By comparison, the contribution of the secondary shear layer is very small. Further downstream, the amplitudes of both fluxes drop significantly as the merged shear layer spreads out. Another interesting point is that inside the recirculation bubble ( $x=3.0$ ), both turbulent fluxes are negative, i.e. they are directed upstream and downwards. This implies that, at this station, heat is transferred from the outer parts of the bubble toward its core. This is due to vortical structures emanating from the shear layers that get trapped inside the bubble.

Finally, we remark that the LES predictions agree fairly well with the DNS, especially as regards $\langle \unicode[STIX]{x1D703}^{\prime }v^{\prime }\rangle$ . The only noticeable disagreement is that the LES overpredicts the values of $\langle \unicode[STIX]{x1D703}^{\prime }u^{\prime }\rangle$ inside the recirculation bubble, and in the bottom half of the channel at $x=6.0$ . Based on the overall comparison of the first- and second-order statistics, we infer that the LES predictions for the flow under study are satisfactorily accurate, although there are some discrepancies in regions of intense turbulence production.

Figure 20.

Profiles of the streamwise turbulent heat flux $\langle \unicode[STIX]{x1D703}^{\prime }u^{\prime }\rangle$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES.

Figure 21.

Profiles of the wall-normal turbulent heat flux $\langle \unicode[STIX]{x1D703}^{\prime }v^{\prime }\rangle$ at different streamwise stations of the horizontal channel. ——, DNS; — ⋅ ⋅ —, LES.

3.3 Budget of the turbulent kinetic energy

In this subsection we analyse the budget of the turbulent kinetic energy (TKE) based on the DNS data. In order to fix ideas, let $K$ and $\bar{K}$ , denote the kinetic energy and the mean-flow kinetic energy, i.e. $K=(1/2)(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{u})$ , and $\bar{K}=(1/2)(\langle \boldsymbol{u}\rangle \boldsymbol{\cdot }\langle \boldsymbol{u}\rangle )$ . The transport equation for TKE can be obtained by subtracting the equation for $\bar{K}$ from the one for $K$ (Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002). However, when applied to variable-density flows, this procedure is quite cumbersome because it involves the density fluctuations. For this reason, it is customary to employ density-weighted averaging (Favre Reference Favre1983). For a generic quantity $\unicode[STIX]{x1D719}$ , the Favre average $\{\unicode[STIX]{x1D719}\}$ is defined as

(3.3) $$\begin{eqnarray}\displaystyle \{\unicode[STIX]{x1D719}\}=\frac{\langle \unicode[STIX]{x1D70C}\unicode[STIX]{x1D719}\rangle }{\langle \unicode[STIX]{x1D70C}\rangle }, & & \displaystyle\end{eqnarray}$$

and the Favre decomposition is then given by

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}=\{\unicode[STIX]{x1D719}\}+\unicode[STIX]{x1D719}^{\prime \prime }. & & \displaystyle\end{eqnarray}$$

An important feature of this decomposition is that the fluctuating component $\unicode[STIX]{x1D719}^{\prime \prime }$ is not centred, i.e. $\langle \unicode[STIX]{x1D719}^{\prime \prime }\rangle \neq 0$ . Relationships between the standard and Favre averaging operations can be found in Huang, Coleman & Bradshaw (Reference Huang, Coleman and Bradshaw1995) and Chassaing et al. (Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002). Let $k$ denote the turbulence kinetic energy of the density-weighted velocity fluctuations,

(3.5) $$\begin{eqnarray}\displaystyle k={\textstyle \frac{1}{2}}(\boldsymbol{u}^{\prime \prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime \prime }). & & \displaystyle\end{eqnarray}$$

Then the Favre-averaged turbulent kinetic energy, $\{k\}=(1/2)\{\boldsymbol{u}^{\prime \prime }\boldsymbol{\cdot }\boldsymbol{u}^{\prime \prime }\}$ , satisfies the following transport equation (Chassaing et al. Reference Chassaing, Antonia, Anselmet, Joly and Sarkar2002),

(3.6) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}\langle \unicode[STIX]{x1D70C}\rangle \{k\}}{\unicode[STIX]{x2202}t} & = & \displaystyle -\underbrace{\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\langle \unicode[STIX]{x1D70C}\rangle \{k\}\{\boldsymbol{u}\})}_{C_{k}}-\underbrace{\langle \unicode[STIX]{x1D70C}\rangle \{\boldsymbol{u}^{\prime \prime }\boldsymbol{u}^{\prime \prime }\}\boldsymbol{ : }\unicode[STIX]{x1D735}\{\boldsymbol{u}\}}_{P_{k}}\nonumber\\ \displaystyle & & \displaystyle -\,\underbrace{\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\langle \unicode[STIX]{x1D70C}\rangle \{k\,\boldsymbol{u}^{\prime \prime }\})}_{T_{k}}-\underbrace{\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\langle p\rangle \langle \boldsymbol{u}^{\prime \prime }\rangle )}_{\unicode[STIX]{x1D6F7}_{k}}-\underbrace{\unicode[STIX]{x1D735}\boldsymbol{\cdot }\langle p^{\prime }\boldsymbol{u}^{\prime \prime }\rangle }_{\unicode[STIX]{x1D719}_{k}}\nonumber\\ \displaystyle & & \displaystyle +\,\underbrace{\langle p\,\unicode[STIX]{x1D735}\boldsymbol{\cdot }\langle \boldsymbol{u}^{\prime \prime }\rangle \rangle }_{\unicode[STIX]{x1D6F1}_{k}}+\underbrace{\frac{1}{Re}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\langle \unicode[STIX]{x1D749}\boldsymbol{\cdot }\boldsymbol{u}^{\prime \prime }\rangle }_{D_{k}}-\underbrace{\frac{1}{Re}\langle \unicode[STIX]{x1D749}\boldsymbol{ : }\unicode[STIX]{x1D735}\boldsymbol{u}^{\prime \prime }\rangle }_{\unicode[STIX]{x1D716}_{k}}.\end{eqnarray}$$

In the above equation, $C_{k}$ describes convection of TKE, while $P_{k}$ and $T_{k}$ describe production and transport of TKE, respectively. The terms $\unicode[STIX]{x1D6F7}_{k}$ and $\unicode[STIX]{x1D719}_{k}$ describe the work of the mean and fluctuating pressures through fluctuating motions, respectively. Also, $\unicode[STIX]{x1D6F1}_{k}$ represents the pressure correlation with the dilatation fluctuation. Finally, $D_{k}$ describes the viscous diffusion, which is the work of viscous stresses through fluctuating motions, whereas $\unicode[STIX]{x1D716}_{k}$ describes viscous dissipation. It should be mentioned that the work of the gravitational force does not enter the above equation.

It is also worth noting that since the fluid viscosity is temperature dependent, one could insert the expression (2.4) for $\unicode[STIX]{x1D749}$ in the last two terms of equation (3.6) and further expand these terms. This expansion results in terms involving triple moments and, as shown by Huang et al. (Reference Huang, Coleman and Bradshaw1995), allows for a more detailed description of the TKE budget. In particular, it provides the specific contribution of the fluctuating component of $\unicode[STIX]{x1D749}$ . Other authors such as Lessani & Zainali (Reference Lessani and Zainali2009) opted not to expand the last two terms of (3.6). In the present study we follow the approach of Lessani & Zainali (Reference Lessani and Zainali2009) for purposes of simplicity of the analysis. Conversely, Zonta et al. (Reference Zonta, Marchioli and Soldati2012) followed yet another approach and excluded altogether the contribution of viscosity fluctuations in the TKE budget.

In statistically stationary flows, the sum of all the terms at the right-hand side of equation (3.6) is zero. Indeed, the computed value of this sum from the DNS data is always very small and below 5 % of the amplitude of the largest term of (3.6), which is as an indication of the accuracy of our DNS. Further, since this computed sum is negligible, we have not included it in the figures for purposes of clarity of the plots.

In figure 22 we present the TKE budget across the channel at $x=-0.5$ . At this station, the main source of turbulence production is the primary shear layer, which can be confirmed by the high values of $P_{k}$ therein. Also, $T_{k}$ takes negative values inside this shear layer, which implies transport of turbulent kinetic energy from the core of the layer toward its edges. On the other hand, the turbulent convection term $C_{k}$ follows the opposite trend since it is positive at the core and negative at the edges of the primary shear layer. As concerns the turbulent dissipation $\unicode[STIX]{x1D716}_{k}$ , it is approximately 3 times smaller than $P_{k}$ . The high ratio between the production and dissipation terms, $P_{k}/\unicode[STIX]{x1D716}_{k}\approx 3$ , corroborates the fact that, with respect to turbulent kinetic energy, the flow at this region is far from equilibrium.

Figure 22.

Budget of the turbulent kinetic energy at the streamwise station $x=-0.5$ . (a) Bulk of the channel, (b) close up look at the bottom wall.

The TKE budget at $x=-0.5$ and close to the bottom wall is shown in figure 22(b). Therein, the boundary layer is subject to a favourable pressure gradient due to the aforementioned Venturi effect. The peak of the production term $P_{k}$ occurs at $y^{+}\approx 10$ ; this is somewhat lower than the location $y^{+}\approx 12$ of the corresponding peak in zero-pressure-gradient boundary layers (Kim, Moin & Moser Reference Kim, Moin and Moser1987; Spalart Reference Spalart1988). Further, the production term is mostly balanced by the work of the fluctuating viscous forces $D_{k}$ . Finally, as expected, closer to the wall all terms vanish except for $D_{k}$ which is balanced by the viscous dissipation.

The profiles of the TKE budget at $x=1.0$ are plotted in figure 23. The main features of the budget are the high amplitudes of the production and convection terms in the zone occupied by the secondary shear layer. In fact, the peak of the production term is an order of magnitude higher than the one at $x=-0.5$ , which confirms that the interaction of the incoming jet with the recirculation bubble is a major source of turbulence generation. Also, the negative values of the convection term $C_{k}$ imply that TKE is convected from the core of secondary shear layer toward its edges and the core of the bubble. At the bottom wall, the pressure gradient is still favourable. Therefore, the profiles of the terms of the TKE budget are similar to those at $x=-0.5$ .

Figure 23.

Budget of the turbulent kinetic energy at the streamwise station $x=1.0$ . (a) Bulk of the channel, (b) close up look at the top wall, (c) close up look at the bottom wall.

The TKE budget at $x=3.0$ is presented in figure 24. The spreading of the shear layers and their merging lead to even higher values of turbulence production above the mid-plane. Also, the peak value of $P_{k}$ is considerably higher; this is in accordance with the high r.m.s. velocities at this station that are plotted in figures 16–19. The production term is balanced by turbulence transport, convection and dissipation. Interestingly, the convection term $C_{k}$ takes negative values close to the mid-plane. However, inside the recirculation bubble it changes sign and becomes positive, which implies convective transport of turbulent kinetic energy from the reattachment region toward the centre of the bubble. This has also been observed in the simulations of flow in channel with constrictions by Fröhlich et al. (Reference Fröhlich, Mellen, Rodi, Temmerman and Leschziner2005). Also, near the bottom wall, the cross-flow is subject to the adverse pressure gradient which results to large turbulence generation; see figure 24(c). In fact, the near-wall peak of $P_{k}$ , located at $y^{+}\approx 12$ , is 3 times higher than the corresponding peak at the upstream station $x=1.0$ . Moreover, according to figure 24(c), the near-wall turbulence production is mostly balanced by $C_{k}$ , $D_{k}$ and $\unicode[STIX]{x1D716}_{k}$ .

Figure 24.

Budget of the turbulent kinetic energy at the streamwise station $x=3.0$ . (a) Bulk of the channel, (b) close up look at the top wall, (c) close up look at the bottom wall.

Figure 25.

Budget of the turbulent kinetic energy at the streamwise station $x=6.0$ .

Downstream of the reattachment point of the flow, the amplitudes of all the terms in the TKE equation (3.6) decrease with $x$ . In figure 25 we have plotted the TKE budget at $x=6.0$ . We observe that $T_{k}$ is positive above and negative below the mid-plane, which implies transport of turbulent kinetic energy away from the walls and toward the bulk of the channel. Additionally, the ratio $P_{k}/\unicode[STIX]{x1D716}_{k}$ is close but not equal to unity. At the same time, due to the spreading of the merged shear layer downstream of the recirculation bubble, the convective and transport terms are also significant and of the same order of magnitude as $P_{k}$ and $\unicode[STIX]{x1D716}_{k}$ . Moreover, nowhere in the flow field is $P_{k}$ counter-balanced by $\unicode[STIX]{x1D716}_{k}$ alone. This implies that, even when the ratio $P_{k}/\unicode[STIX]{x1D716}_{k}$ gets close to unity, the flow turbulence is not at equilibrium.

Finally, in figure 26 we provide the TKE budget far downstream, at the station $x=12.0$ . We can readily observe that as the flow evolves, the amplitudes of the various terms in the TKE budget decrease even further. Moreover, the profiles at the near-wall regions, which are not shown herein for purposes of economy of space, start to resemble those of zero-pressure-gradient boundary layers; see, for example, Kim et al. (Reference Kim, Moin and Moser1987) and Spalart (Reference Spalart1988).

Figure 26.

Budget of the turbulent kinetic energy at the streamwise station $x=12.0$ .

Figure 27.

Budget of the fluctuating temperature variance (FTV) at the streamwise station $x=-0.5$ .

3.4 Budget of the fluctuating temperature variance

In this subsection we analyse the budget of the fluctuating temperature variance (FTV) at different streamwise locations based on the DNS data. The transport equation for FTV, $(T^{\prime \prime }T^{\prime \prime })/2$ , can be derived by first multiplying every term of (2.3) with $T$ so as to obtain an equation for $T^{2}$ , then by Favre averaging the result in order to derive an equation for $\{T^{2}\}$ and, finally, by subtracting the equations for $T^{2}$ and $\{T^{2}\}$ . The resulting transport equation reads

(3.7) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2}\frac{\unicode[STIX]{x2202}(\langle \unicode[STIX]{x1D70C}\rangle \{T^{\prime \prime }T^{\prime \prime }\})}{\unicode[STIX]{x2202}t}=-\frac{1}{2}\underbrace{\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\langle \unicode[STIX]{x1D70C}\rangle \{T^{\prime \prime }T^{\prime \prime }\}\{\boldsymbol{u}\})}_{C_{T}}-\underbrace{\langle \unicode[STIX]{x1D70C}\rangle \{T^{\prime \prime }\boldsymbol{u}^{\prime \prime }\}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\{T\}}_{P_{T}}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\frac{1}{2}\underbrace{\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\langle \unicode[STIX]{x1D70C}\rangle \{T^{\prime \prime }T^{\prime \prime }\,\boldsymbol{u}^{\prime \prime }\})}_{T_{T}}+\underbrace{\frac{1}{Re\,Pr}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\langle T^{\prime \prime }\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}T\rangle }_{D_{T}}-\underbrace{\frac{1}{Re\,Pr}\langle \unicode[STIX]{x1D705}\unicode[STIX]{x1D735}T\boldsymbol{\cdot }\unicode[STIX]{x1D735}T^{\prime \prime }\rangle }_{\unicode[STIX]{x1D716}_{T}}.\end{eqnarray}$$

In the above equation, the terms $C_{T}$ , $P_{T}$ and $T_{T}$ describe convection, production and transport of fluctuating temperature variance, respectively. Further, $D_{T}$ describes viscous transport, whereas $\unicode[STIX]{x1D716}_{T}$ describes viscous dissipation of FTV. Similarly to the TKE equation (3.6) one could expand the last two terms in the above equation even further so as to obtain terms that involve the fluctuating component of the heat flux vector $\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}T$ . As in the case of the TKE equation, we opt not to follow this approach for purposes of simplicity of the analysis; in any case, in the flow under study the ratio between the conductivities of the two streams is small (approximately equal to 1.15).

Once the flow becomes statistically stationary, the sum of all the terms at the right-hand side of the above equation is zero. Indeed, the computed value of this sum based on the DNS data is small and below 20 % of the amplitude of the largest term.

In figure 27 we show plots of the terms at the right-hand side of (3.7) at the streamwise station $x=-0.5$ . We see that the various terms of the budget become significant only across the primary shear layer, which provides the main mechanism for thermal mixing. The peak of the production term occurs at the core of this shear layer and is counter-balanced by the transport term and, to a lesser extent, by the dissipation term. We also observe that near the edges of this shear layer, the convection and transport terms are of opposite sign and they essentially cancel each other out. At the same time, the values of the viscous transport term $D_{T}$ are very small compared to those of the other terms.

The same observations hold for the budget of the FTV at $x=1.0$ , which is shown in figure 28. Moreover, at this streamwise location the profiles of the various terms have ‘spread out’ and their values have decreased due to the spreading of the primary shear layer. In figure 29 we have plotted the FTV budget further downstream, at $x=6.0$ . We readily infer that the amplitudes of all terms, have decreased even further. On the other hand, the relative strength of the dissipation term $\unicode[STIX]{x1D716}_{T}$ has increased so that it almost counter-balances the production term $P_{T}$ . At the same time, the convective and transport terms have opposite signs and almost equal amplitudes so that they essentially cancel each other out. Finally, according to our simulation, as the thermal mixing evolves, all the terms in the FTV budget decrease substantially further downstream.

Figure 28.

Budget of the fluctuating temperature variance (FTV) at the streamwise station $x=1.0$ .

Figure 29.

Budget of the fluctuating temperature variance (FTV) at the streamwise station $x=6.0$ .

3.5 Spectral analysis of the velocity and temperature signals

We proceed with the spectral analysis of the velocity and temperature signals at two locations of interest. The first one is at $(x,y)=(-0.5,\,1.6)$ , i.e. beneath the entry of the jet and at the interface of the two streams. The second location is $(x,y)=(12.0,\,1.0)$ , i.e. far downstream of the recirculation bubble and at the mid-plane. Our procedure is the following. Once the flow has become statistically stationary, we record the time signals of the spanwise-averaged velocity components $\overline{u}$ and $\overline{v}$ and temperature $\overline{T}$ at the above locations for 40 non-dimensional time units. Then, we compute the power spectral densities, $S_{uu}(f)$ , $S_{vv}(f)$ and $S_{TT}(f)$ , of these signals so as to estimate the characteristic frequencies $f$ with which they oscillate. The frequency is then normalised with respect to the half-width of the horizontal channel and the bulk velocity of the cross-flow, which provides the corresponding Strouhal number $St=f\unicode[STIX]{x1D6FF}/U_{c}$ .

Figure 30.

Power spectral density of the spanwise-averaged signals $\overline{u}$ , $\overline{v}$ and $\overline{T}$ , at two locations in the horizontal channel. Each profile has been normalised with respect to its maximum value. (a)  $(x,\,y)=(-0.5,\,1.6)$ , (b)  $(x,\,y)=(12.0,\,1.0)$ .

Figure 30 shows the power spectral densities at these two locations. At the first one, $(x,y)=(-0.5,\,1.6)$ , they all exhibit three low- $St$ peaks, at $St=0.035,\,0.09$ and 0.15. Further, the peaks of $S_{uu}$ and $S_{vv}$ occur at the same frequencies. These peaks are thus attributed to the oscillatory motion of the primary shear layer and to its shedding of vortices. Also, the three low- $St$ peaks of $S_{TT}$ occur at slightly different Strouhal numbers than those of $S_{uu}$ and $S_{vv}$ ; this is due to the fact that thermal mixing at this location is still at an early stage and, therefore, the fluctuations of the temperature are not yet synchronised with those of the velocity. We further observe that the profile of $S_{TT}$ has a strong peak at a higher Strouhal number, $S_{t}\approx 0.5$ , in contrast to the profiles of $S_{uu}$ and $S_{vv}$ . Based on the video animations of the DNS (see the supplementary movies), we may attribute this peak to the wrapping of fluid of different temperature by the mushroom-like structures at the interface of the two streams such as the ones depicted in figure 10(a).

On the other hand, at the second location $(x,y)=(12.0,\,1.0)$ , the peaks of the power spectral densities occur at the same frequencies. In other words, the dominant frequencies of the oscillatory motion of $u$ , $v$ and $T$ coincide. The synchronisation of the power spectral densities is due to the extent of thermal mixing that has taken place. It also implies that once the two streams have been aligned sufficiently downstream of the recirculation bubble, the convective motions of spanwise vortices provide the main mechanism for thermal mixing.