Compressibility in suddenly expanded subsonic flows Open Access

Planar sudden expansion (PSE) flows are fascinating examples of asymmetric and unstable flow physics even when the geometry is fully symmetric.¹ This behavior could be ascribed to the Coanda effect,² but the exact mechanism is still not fully understood despite many years of research.

Suddenly expanded flows remain symmetric up to a certain Reynolds number, Re, with two recirculation regions of equal length formed on either side of the expanding channel. As the Re increases, the flow becomes asymmetric about its centerline and separation regions of unequal length arise and remain in the flow field even up to very high Reynolds numbers. Such flows are encountered in various engineering applications, such as hydraulic and fluidic devices, mixing equipment, combustion, and air ducts. Therefore, it is crucial to understand the mechanisms that dominate flows with the separation and reattachment of the free-shear layers.

Several PSE studies were performed in the laminar flow regime,^1,3,4 including both numerical and experimental investigations. For example, for an expansion ratio ER = 3 (the ratio of the main channel height to the inlet channel height), the flow is remarkably dependent on the Re and at high velocities three-dimensional (3D), even away from the channel corners. Furthermore, there is a strong relationship between the ER, the channel's aspect ratio AR, and the critical Reynolds number value Re_c, where the flow bifurcates to asymmetric separation.⁵ Several researchers, namely, Fearn et al.,⁴ Drikakis,¹ and Alleborn et al.,⁶ predicted that for ER = 3 the Re_c based on the half upstream height is $R e c ≃ 40.45 , 40$ and 40, respectively. More recently, Moallemi and Brinkerhoff⁷ showed using direct numerical simulation (DNS) that bifurcation initiates at higher Reynolds numbers as the expansion ratio decreases. Several recent numerical studies concerning sudden incompressible expansion flows at low Reynolds numbers have been published,^8–11 which primarily examine the symmetry breaking, bifurcation phenomena, and critical Reynolds number.

Three-dimensional numerical simulations using a finite volume method were performed by Tsui and Wang¹² to study the influence of different values of aspect ratio on the bifurcation of the flow through a symmetric, sudden expansion of ER = 3. For AR = 3, the flow remains symmetric for all Reynolds numbers under consideration. For AR = 1 the symmetry-breaking bifurcation occurs at $R e c ≃ 92$⁠, for AR = 4 at $R e c ≃ 61$⁠, and for AR = 8 its value reduces furthermore to about 58.5. Reducing the aspect ratio has a stabilizing effect on the flow. This observation agrees well with the numerical results of Schreck and Schafer.¹³ 

Transitional and turbulent flows through a PSE have been less investigated. The first study of turbulent flow over a PSE was undertaken by Abbott and Kline.¹⁴ In a symmetric plane channel, hot-film velocity measurements were performed with ER varying from 1.125 to 5 and AR from 2.5 to 5. They demonstrated that no changes regarding the reattachment length or the flow pattern for Reynolds numbers ranging from $4 × 10 4$ to 10⁵ occurred for a fully turbulent inlet profile.

The first numerical simulation of a turbulent flow over a two-dimensional (2D) PSE was carried out by Gagnon et al.¹⁵ The random vortex method was used to simulate a flow over a double backward-facing step. The geometric features of the channel and the Reynolds number (⁠ $R e = 10 5$⁠) were chosen to be the same as those used in the experimental study of Mehta.¹⁶ Mean velocity measurements were in good agreement with experiments. However, the turbulence intensities were not comparable, and this is mainly attributed to the three-dimensionality of the flow in the experiments conducted by Mehta.¹⁶ later, LDA measurements and numerical simulations of PSE were performed by De Zilwa et al.¹⁷ Spanwise measurements revealed that the mean velocity is uniform over more than 80% of the span.

Aloui and Souhar¹⁸ performed an experimental study of turbulent asymmetric flow in a flat channel with a sudden symmetric expansion of aspect ratio 0.18 and expansion ratio 2.27. Hot-film anemometer measurements were conducted at a Reynolds number of 32 000 based on the inlet bulk velocity and the upstream height of the duct. The flow is asymmetric downstream of the step, as in a 2D case, exhibiting two unequal size recirculation zones in both sides of the expansion. The regions of the recirculating flow are constituted of five vortices of different size lengths.

Studies concerning compressible flows through a sudden expansion are generally a few and far apart. Drewry¹⁹ performed an experimental study of compressible axisymmetric sudden expansion flow for an inlet Mach numbers below 1. Emmert et al.²⁰ on the other hand, performed a numerical study using compressible large eddy simulation of a sonic flow in a sudden planar expansion. Emmert et al.²¹ carried out a particle image velocimetry (PIV) study of a turbulent flow downstream of a PSE with ER = 3 and AR = 10. Two different Reynolds numbers were considered, $R e = 10 4$ and $R e = 4 × 10 4$⁠, to investigate their effect on the flow structure. It was shown that $R e = 10 4$ is not high enough to conclude that any further increase in Re does not affect the flow structure.

In this study, the compressibility effects on the structure of the turbulence are investigated for subsonic flows. We briefly present the governing equations in Sec. II, followed by a description of the numerical methodology, boundary conditions, and computational meshes employed in Sec. III. The results are discussed in Sec. IV. We present grid convergence studies and comparisons with experiments. Furthermore, we discuss the mean flow, separation bubbles and reattachment lengths, the pressure and friction coefficients along the channel's walls, and the Mach effects on the turbulence characteristics, including the Reynolds stress anisotropy and turbulence kinetic energy budgets. Finally, Sec. V provides a detailed summary of the key findings of the present study.

Figure 1 shows a schematic of the sudden expansion, which comprises of two channels, each of a different height. The flow domain consists of an inlet channel of height h and a downstream channel of height 3h. The characteristic length of the channel is the step height, h, with a value of 1. The total length of the domain is 84h, where the inlet channel has a length of 4h, and the downstream channel has a length of 80h. These particular geometrical features were chosen to ensure that the flow (a) is fully developed (turbulent) before reaching the expansion step, and (b) the buffer layer at the end of the domain has damped unsteadiness in the flow before exiting the domain. The expansion ratio, which is of great importance when simulating suddenly expanded flows, equals to 3, having the same value as that of Casarsa and Giannattasio.²¹ The aspect ratio (⁠ $A R = w / h$⁠) considered is 5. For the Ma = 0.5 case, the Reynolds number was also $R e = 10 4$ by scaling the dynamic viscosity accordingly; all other fluid flow properties remained the same for the two Mach numbers considered.

Synthetic turbulent boundary conditions were implemented and used throughout the simulations based on the digital filter (DF) technique originally proposed by Klein et al.²² and further developed by Touber and Sandham.²³ The DF approach is used to produce a velocity signal in three directions by matching ad hoc first- and second-order statistical moments, length and time scales, and energy spectra. We selected the DF approach for generating artificial inflow data as the filtering operation applied only in 2D (rather than in 3D), making the whole process much faster and computationally efficient.

All DF input properties required were obtained from previously available DNS results of Iwamoto et al.²⁴ for a fully developed turbulent planar channel flow (database: https://thtlab.jp/index-orig.html). According to Pope,²⁵ the Reynolds number based on friction velocity, $R e τ$⁠, and mean velocity and channel height, Re_m, are related according to $R e τ = 0.09 R e m 1 / 0.88$⁠. Thus, for an $R e τ = 300$⁠, i.e., $R e m ≃ 10 4$⁠, as per the inflow conditions considered herein.

A buffer layer is employed at the outflow to avoid any numerical reflections, while in the spanwise direction, we prescribed periodic boundary conditions. Standard no-slip, viscous, and adiabatic wall boundary conditions were assigned to the surfaces of the domain normal in the y-direction, indicated by the translucent gray colored in Fig. 1.

We used the computational fluid dynamics (CFD) code $CNS3D$ (Compressible Navier–Stokes Solver in 3D). The code^26–29 uses the finite volume Godunov (upwind) method in conjunction with several approximate Riemann solvers. It includes high-resolution methods of up to 11th-order of accuracy in space and fourth-order of accuracy in time. In the present study, the Harten, Lax, and van Leer plus missing contact wave (HLLC) Riemann solver³⁰ is employed for the convective fluxes, in conjunction with a modified³¹ weighted essentially non-oscillatory (WENO) 11th-order spatial reconstruction scheme.³² For the viscous fluxes, a second-order central scheme was employed.

Very high-order WENO schemes (greater than seventh-order) are resilient to the low-Mach number dissipation inherent to compressible solvers, for Mach numbers as low as $M a ≃ 0.01$⁠.^33,34 A recent study³⁵ showed that both fifth-order monotone upstream-centered schemes for conservation laws (MUSCL)- and WENO-type schemes generally require a low-Mach number treatment³⁶ to remain accurate at very low-Mach numbers (⁠ $M a ≲ 0.1$⁠).

Moreover, high-order WENO schemes offer increased accuracy while improving parallel scalability performance.³⁷ Typically, a ninth-order WENO scheme is sufficient to resolve most turbulence scales—for a modest LES grid—and generally performs well in the context of ILES.³³ However, in cases for which grid stretching may lead to diminishing mesh resolution, a higher-order numerical scheme such as the 11th-order WENO scheme employed in this study will—to some extent—help lessen the impact.

The solution was advanced in time using a five-stage (fourth-order accurate) optimal strong-stability-preserving Runge–Kutta method.³⁸ Further details regarding $CNS3D$ can be found in past studies^31,33,34 and references therein.

The grid resolution at the channel inlet in the streamwise, wall-normal, and spanwise directions is $Δ x + = 26 , y w + = 1$ (corresponding to $y + = 0.5$ in the first cell-center off the wall), and $Δ z + = 15$⁠, respectively. $Δ x +$ gradually decreases to $y w +$ at the expansion corner for all grids examined. It then gradually coarsens again toward the exit boundary. For the fine grid in particular, it reaches approximately a value of $Δ x + ≃ 25$ at the expected $L 1 ≃ 14.4 h$ (Fig. 5). After that, it increases linearly toward the outlet, eventually reaching a value of $Δ x + ≃ 180$⁠. The grid resolution and numerical scheme employed were shown to yield accurate results³³ for wall-bounded turbulent flows at low-Mach numbers. Following typical resolution recommendations for LES and DNS simulations,^39–41 the present fine mesh corresponds to wall-resolved ILES.

For the fine mesh considered, the total simulation time was $300 h / U B$⁠, while statistics were obtained over the last $100 h / U B$⁠, where U_B is the bulk velocity of the channel inlet, and h is its height. The three mesh resolutions are presented in Table I.

Figure 2 shows the wall-normal profiles of the mean streamwise velocity component, U. The RMS turbulent fluctuation, $u ′$⁠, measured $x = 3 h$ downstream of the channel inlet (or six boundary layer heights $( h / 2 )$ downstream of the digital filter inflow) to ensure the turbulent flow has developed sufficiently. Both quantities are normalized by the bulk velocity in the (wall-normal) plane, U_B, and compared against the DNS results of Iwamoto et al.,²⁴ who performed a numerical simulation of a fully developed turbulent flow between two parallel flat plates at Re = 10 049. The ILES results agree reasonably well with the DNS data, particularly those related to the mean flow velocity component.

At the centerline of the inlet channel $y / h = 0.5$⁠, we found the streamwise turbulent fluctuations to be approximately 4.7% of U_B (⁠ $∼ 5 %$ deviation from the DNS results). The noticeable asymmetry of $u ″ / U B$ is associated with the inflow location relatively to the sudden expansion corner. Although the upstream effects deserve further investigation, we note that in any experimental or engineering device, the flow is not expected to be fully symmetric near to the entry point of the expansion corner due to upstream perturbations and wall roughness effects.

The maximum values of RMS velocity are observed near the lower and upper wall, at $y / h ≃ 0.1$ and $y / h ≃ 0.9$⁠, respectively. In these locations, the RMS velocity fluctuations are 15% of U_B (⁠ $∼ 9.3 %$ deviation from the DNS results with $u ′ = 16.4 %$ of U_B). Given the above, we consider the flow to be fully developed and turbulent in the xy-inlet plane, while overall the ILES results are in excellent agreement with the PIV measurements.²¹ 

Numerical measurements of the mean flow quantities (i.e., velocity components, fluctuating velocities, Reynolds stresses, and turbulence kinetic energy) were conducted at three different wall-normal locations in the x–y plane at the midplane of the channel (⁠ $z / h = 2.5$⁠), as shown in Fig. 3. The time-averaged streamwise and transverse velocity profiles are depicted in Fig. 4 using three different grid meshes. The velocity is normalized by the inlet bulk velocity, U_B, and the numerical results are compared against the experimental data.²¹ 

The Mach number based on the bulk velocity (U_B) at the channel inlet is 0.1, which after expansion reduces by a factor of one-third. Thus, compressibility effects are expected to be very weak so $ϕ ¯ ≡ ϕ ̃$⁠, and therefore the turbulent mass-flux becomes negligible, i.e., $a i ≡ − u i ″ ¯ ≡ ρ ′ u i ′ / ρ ¯ ¯ ≈ 0$⁠.

Mean $u ̃$ and $v ̃$ results along sample streamwise lines are plotted in Fig. 4. In the region of the shorter reattachment (⁠ $y = − 0.5 h$⁠), all grids are found to resolve the mean flow field well up to $x ≃ 5 h$⁠. Though the coarse grid resolution matches best with the experimental data, the two finest grid resolutions are grid converged and exhibit only a minor over-estimation of the streamwise velocities. This is evident for both (finer) grids up to $x / h = 16$⁠. The mean transverse velocity is satisfactorily predicted in the region of the shorter reattachment for all grid resolutions, but the finer grids perform overall better here.

In the core flow region (⁠ $y / h = 0.5$⁠), both velocities profiles exhibit a minor shift in the upstream direction, but the agreement is excellent, and the solution grid converged.

The Favre-averaged streamwise velocity is negative up to $x ≃ 12 h$⁠, in the region of the larger recirculation zone (⁠ $y = 1.5 h )$⁠. Similar to the $y = − 0.5 h$ line, the coarsest grid resolution best matches the experimental data, where the two finest overestimate the streamwise velocity magnitude compared to experimental data. Importantly, however, the two finest grid resolutions are grid converged. Besides, as the grid resolution increases, the mean velocity components agree better with the experimental data. The maximum absolute value of the streamwise velocity is observed at $x ≃ 6.3 h$ with $u / U B = − 0.21$⁠, which is considerably higher than the experimental value of $u / U B ≈ − 0.21$⁠.

The results obtained from the numerical simulations are compared against the experimental data of Casarsa and Giannattasio.²¹ The experimental study investigated a turbulent flow downstream of a planar sudden expansion using a 2D PIV technique. In Fig. 5, the time-averaged flow paths in the midplane of the channel are visualized by employing stream tracers for all three grid meshes. The 2D statistically symmetric flow encountered in the inlet channel emerges in the expansion, creating an asymmetric flow pattern deviated toward the lower sidewall. The high-velocity core is separated from two recirculation bubbles of considerably unequal length by two free-shear layers emanating from the expansion edges of the upstream channel. The shorter recirculation zone is encountered at the lower wall where the flow primarily impinges. The larger recirculation zone is placed at the upper wall. Two secondary vortices near the corners of the side walls directly after the expansion also form. These small-size vortices were also reported in Casarsa and Giannattasio.²¹ However, no further investigation on their origins and contribution to the primary large flow structure was carried out.

The primary reattachment and secondary separation lengths, determined by evaluating the locations at which the mean streamwise velocity component changes sign, are shown in Table II. The lengths of the shorter recirculation bubbles L₃ and L₄, and those of the larger recirculation bubbles L₁ and L₂, are normalized by the step height of the channel, h. The PIV data, which are also given in Table II, agree well with the numerical results. Notably, the reattachment length of the larger recirculation increases with increasing grid resolution.

The larger L₁ reattachment length obtained in the ILES could be attributed to the impact that the two free-shear layers (emanating from the expansion corner) have on the turbulent flow properties upstream of the expansion corner and into the smaller (inflow) channel section. As is evident by the asymmetry of the wall-normal profiles in Fig. 2, the ILES fluctuating streamwise velocity is weakened, particularly along the lower wall of the inflow channel section. It can be implied that a decrease in turbulence fluctuations is indicative of a lower Reynolds number to some extent. Since the PIV data of Casarsa and Giannattasio²¹ indicate that for a higher Reynolds number (⁠ $R e = 4 × 10 4$⁠), the large recirculation bubble reattachment length decreases (⁠ $L 1 = 14.17$⁠), reversely, a reduction in the turbulence intensity in the channel section can—justifiably—lead to an increase in L₁.

As far as the lengths of the secondary vortices are concerned, both are of the same order of the step height. The same feature was addressed by Spazzini et al.,⁴² where the unsteady behavior of a flow over a backward-facing step was investigated. They found that the length of the secondary vortex is of the order of one step height.

Abbott and Kline¹⁴ claimed that the reattachment length depends on the geometrical characteristics of the sudden expansion configuration, i.e., the expansion and aspect ratios. Therefore, a precise comparison of the reattachment lengths with previous experimental and numerical data is almost impossible because of the widely different geometries and inflow conditions used in these studies. Nonetheless, it is still worth addressing some of the differences with the previous results while at the same time analyzing the influence of the flow parameters on the reattachment lengths. In their pioneering work, Abbott and Kline¹⁴ stated that the reattachment lengths of the recirculation zones are strongly related to the expansion ratio. For the case of ER = 3, they obtained the following ranges: $L 1 / h = 11 − 15$ and $L 3 / h = 3.5 − 4$⁠, which agree reasonably well with the ILES results. They used a fully developed inlet velocity profile with the Reynolds number ranging between 20 000 and 50 000. Mehta¹⁶ for R = 3 and a fully developed inlet velocity profile at $R e = 1.25 × 10 5$ measured the reattachment lengths $L 1 = 15 h$ and $L 3 = 4.5 h$⁠, which are slightly lower and higher, respectively, than those found in the present study. This can be attributed to the low aspect ratio used for the channel configuration (AR = 0.25). The reattachment lengths found by Aloui and Souhar¹⁸ for a fully developed inlet velocity profile (⁠ $L 1 = 10.8$ and $L 3 = 5.4 )$ are not in good agreement with the other results presented in Table II. The relatively small expansion ratio (ER = 2.27) and the extremely low aspect ratio (AR = 0.18) are believed to be the main reasons for this discrepancy. De Zilwa, Khezzar, and Whitelaw¹⁷ obtained $L 3 = 3.4$ for the small recirculation area, which agrees reasonably well with the ILES findings herein. However, the reattachment length of the larger recirculation zone $L 1 ≃ 17$ was found to be much more prominent in their experiments (LDA technique). It is essential to underline that a uniform inlet velocity profile was used in their experiments along with ER = 2.86 and AR = 12.31, which can explain the difference between the experimental and numerical results. Finally, the reattachment lengths demonstrated in Escudier et al.⁴³ (⁠ $L 1 = 11.5$ and $L 3 = 3.13$⁠) are lower compared to the ILES results. However, care should be paid in the comparisons due to the different expansion ratios and initial conditions used in their LDA experiments (ER = 4 with a uniform inlet velocity profile at $R e = 5.55 × 10 4$⁠). The reattachment lengths (x_R) found in the literature of PSE flows along with the corresponding geometrical properties are all summarized and presented in Table III.

Along the upper wall, C_p exhibits a local minimum immediately after the step at $x ≃ 1.8 h$⁠. This is close to the position of the maximum backflow. Moreover, a local maximum and a second local minimum are observed at the streamwise locations of $x ≃ 3.8 h$ and $x ≃ 6.7 h$⁠, respectively. It should be noted that the local pressure maximum almost coincides with the reattachment point of the flow on the upper side of the wall. This last observation is in excellent agreement with the pressure coefficient results presented by El Khoury et al.⁴⁵ for a flow through a channel with a single thin-plate obstruction. De Zilwa et al.¹⁷ performed RANS calculations of a turbulent flow through a planar sudden expansion, calculating, among others, the wall pressure coefficient. The resulting shape of the wall-pressure curve along the upper wall is very similar to the variation of C_p along the wall with the smallest recirculation vortex shown in Fig. 6.

Measurements of the wall-pressure coefficient have also been reported by Escudier et al.⁴³ for a turbulent flow over a planar sudden expansion. Comparisons between ILES results and experimental data are presented in Fig. 6. The shape of the wall-pressure variation for the current study seems to be in excellent agreement with the experiment. However, its magnitude along the upper and lower wall of the channel differs from each other. The above is expected as both Re and ER are different between the two studies. In the experiment, the peaks of the wall-pressure coefficient are larger and shifted moderately upstream. This is consistent with the reattachment lengths reported by the same authors. The length of the shorter recirculation vortex was found to be larger in their experiments and equal to x_R = 4.7. In the recovery region, the pressure on the upper wall of the channel becomes identical to that on the opposite wall at $x ≃ 10 h$⁠, for both the numerical and experimental case, but with a higher C_p in the numerical study.

As previously stated, the secondary pressure minimum in the region of the shorter reattachment is located at $x ≃ 6.7 h$ downstream of the constriction. This local minimum is associated with the largest recirculation bubble, which is more pronounced in this cross-section area, forcing the flow to speed up along the upper wall region with an accompanying pressure drop. Moreover, this pressure drop can also be related to the impingement of the high-velocity core flow on the wall at that particular location.

The distributions of the normalized (Favre-averaged) Reynolds shear stress are shown in Fig. 8. The shear stress is consistent with the mean flow path of the streamwise velocity since the inversion of the stress takes place at points of maximum U velocity. In the region of the shorter reattachment (⁠ $y / h = − 0.5$⁠) at $x / h ≃ 2 − 5$⁠, the shear stress is negative with its maximum value, considerably larger than the experiment. However, further downstream, as the flow recovers from the initial deviation and its velocity decreases, the shear stress becomes approximately zero. For $x / h > 6$⁠, the numerical results are in excellent agreement with the PIV data. The turbulence kinetic energy increases rapidly up to $x / h ≃ 3.2$⁠, reaching its peak value, whereas, beyond that point, it decays gradually, as expected.

In the high core velocity region, the shear stresses are low up to the point where the measurement line crosses the free-shear layer of the upper wall, while it moves toward the opposite—lower wall. The numerical results agree well with the experiment with some small deviations downstream of the flow reattachment. The turbulence kinetic energy increases immediately after the step, reaching its peak value at the same location where the maximum value of the streamwise r.m.s. velocity is detected while maintaining its peak values up to $x ≃ 8 h$ before beginning to attenuate.

In the region of the larger recirculation zone (⁠ $y / h = 1.5$⁠), the ILES results are in excellent agreement with the PIV data, given the challenges of converging the solution data for this flow region, as the relatively larger fluctuations in the plotted numerical data imply. The magnitude of the shear stress $R ̃ 12$ is slightly lower in the present study, indicating lower turbulence levels at the core of the larger recirculation vortex than measured by Emmert et al.²¹ The same phenomenon was addressed by Sugawara et al.,⁴⁶ with the shear stress being approximately 0.15 or 15% of U_B at $x ≃ 8 h$⁠, which agrees well with the numerical results found at the same axial location. As far as the turbulence kinetic energy is concerned, its peak value is close to half that detected in the region of the shorter bubble, with dimensionless $k max ≃ 5$ at $x = 8 h$⁠, situated at the periphery of the core flow of the large recirculation area.

One method of examining the differences that arise in the turbulent flow properties as a result of compressibility effects is to examine the barycentric anisotropy invariant map of Banerjee et al.⁴⁷ To enhance the clarity between the different turbulence anisotropy states, we have instead implemented a modified form of the map first suggested by Emory and Iaccarino⁴⁸—a description of the procedure and underlying theory was detailed in Ref. 31.

In brief, the colored equilateral triangle legend in Fig. 9 is used to construct and visualize the barycentric anisotropy invariant map. The barycentric map componentality contours comprise of six distinct colors, three of which—RGB primary: red, green, blue—describe a limiting (triangle peaks) turbulence state, while the remaining three—RGB secondary: cyan, yellow, magenta—describe intermediate states resulting from the convex combination of the aforementioned limiting states.

Of the primary colors, the blue limiting state (⁠ $X 3 c$⁠) corresponds to isotropic turbulence, which is a three-component state in which turbulence fluctuations exist in all directions and are of equal magnitude. The green limiting state (⁠ $X 2 c$⁠) indicates the isotropic two-component axisymmetric limit in which turbulence fluctuations exist along two directions with equal magnitudes. The red limiting state (⁠ $X 1 c$⁠) is indicative of the one-component limit in which fluctuations only exist along one direction. Based on their visual description, the red, green, and blue limiting states designate turbulence that is shaped like a rod (or cigar-like), a circular disk, and a sphere, respectively.

The secondary colors are found at the mid-sides of the triangular barycentric anisotropy invariant map and define the so-called transition regions—otherwise known as intermediate states. The magenta intermediate state corresponds to the axisymmetric expansion limit in which one diagonal component of the Reynolds stress tensor is larger than the other two equal components. The cyan intermediate state represents the axisymmetric contraction limit in which one component is smaller than the other two equal components. And the yellow intermediate state is indicative of the two-component limit in which one component is zero and the two others are not equal. The latter (yellow), typically forms near solid walls since the wall-normal component of the fluctuations vanishes much faster than the other two components. Based on their visual description, the magenta, cyan, and yellow intermediate states designate turbulence shaped as a prolate spheroid, oblate spheroid, and an elliptical disk, respectively. As a result, both the intensity and type of turbulence anisotropy in the flow can be readily identified in Fig. 9.

Regarding the inlet channel section, either Mach-number simulations exhibit very similar turbulence anisotropy properties. It is evident that the turbulence anisotropy state gradually develops from the inflow. As expected, initially a two-component limit (yellow) forms near to the wall, while just below it, in closer proximity to the wall, a one-component limiting state (red) develops, indicative of the near-wall vortex streaks. According to Emory and Iaccarino⁴⁸ and the DNS data⁴⁹ for an incompressible turbulent channel flow, the upper region of the boundary layer should be predominantly in the axisymmetric expansion limit (magenta), slowly evolving toward the isotropic limit (blue, three-component state) at the channel center height. While the latter still appears to be the case here, for the former, a well-defined axisymmetric contraction state (cyan) can be seen instead. It is unclear whether this results from the proximity of the inflow to the expansion corner, or the flow establishing itself from the artificially generated turbulence generated by the digital filter at the inflow, or even a combination of both. Regardless, further investigation in the future is warranted to help establish the exact cause of this discrepancy.

Two very small one-component (red) turbulence regions can be distinguished within very close proximity to the lower and upper expansion corners. Both correspond to tiny (rod-like) recirculation bubbles that form at those locations along the spanwise direction.

In proximity to the vertical (yz) walls of the expansion corner, two long and thin regions can be distinguished. Close along the lower wall, and for both Mach numbers examined, the turbulence anisotropy componentality is in the axisymmetric expansion limit (magenta), or, in other words, it comprises predominantly of prolate spheroid-like structures. On the opposite side, along the upper vertical wall and near to the mixing layer, the turbulence anisotropy componentality is of the one-component limit (red) initially, which gradually transforms to axisymmetric contraction as the upper expansion corner is approached. It is thought that the flow acceleration in that vicinity—caused by the drag and downwards motion of the upper mixing layer—is the cause for the oblate spheroid structures becoming “stretched” and thus more tubular-like.

An interesting observation is the noticeable thinner structure of both of the aforementioned hook-like regions for the higher Mach-number case. On the other hand, the (upstream) inflow channel section seems to have been largely unaffected by compressibility effects, even at the expansion corner. However, as will be discussed in more detail later, this does not remain the case for the expanded channel, particularly the upper region associated with the L₁ and L₂ recirculation bubbles.

Regarding the lower mixing (free-shear) layer and the two separation bubbles corresponding to L₃ and L₄ (see Fig. 5), no significant compressibility effects on the turbulence anisotropy componentality are evident; it only differs at the L₄ bubble inner center region. The incompressible case (Ma = 0.1) appears to be nearer to the center region of the componentality contour legend triangle, with a slight incline toward axisymmetric expansion (magenta). On the other hand, for the compressible case (Ma = 0.5), the turbulence is almost isotropic (blue), with perhaps a slight incline toward axisymmetric contraction (cyan).

Comparing the outer region of the upper mixing layer and above it, several noticeable differences between the two Mach-number cases become immediately apparent:

1. The one-component limit (red) found in the vicinity of $x ≃ 1 h − 2 h$ and $y ≃ 1.5 h$⁠, part of the L₁ recirculation bubble region near to the L₂ bubble, transforms to a white-magenta “haze” colored region, indicating a shift toward, but not to, axisymmetric contraction. Similar to earlier, the one-component limit is thought to be caused by the flow acceleration (momentum magnitude increase) induced by the upper mixing layer, resulting in the pull/stretch of the vortices present and thus transforming them to tube-like structures. For the higher Mach-number case, the streamwise normal Reynolds stress component, $R ̃ u u$⁠, is found to be weaker relatively to the incompressible case in the same region. Careful examination of the streamwise pressure gradient reveals that the upstream flow acceleration in the counter-clockwise direction (along the upper wall of the L₁ recirculation bubble), is less for the compressible case. Consequently, the acceleration leading to the tube-like (one-component) turbulent structures is not as great and does not remain as an effective mechanism as for the incompressible case.

2. A thin and well-defined isotropic turbulence layer forms just above the upper mixing layer, extending into most of the L₁ separation bubble left half (⁠ $x = 0 − 5 h$⁠). The isotropic layer forms along a path at which the streamlines indicate the flow is changes direction and begins to move downstream (see, for example, Fig. 5), i.e., is at maximum shear and lowest velocity magnitude. This suggests that the isotropic state results from the mean flow deceleration and subsequent compression of the turbulence.

3. In the same general region and specifically the upper (left) half of the L₁ recirculation bubble, compressibility effects have also transformed what was predominantly an axisymmetric expansion (magenta) region to a mix of axisymmetric expansion (magenta) and isotropic turbulence (blue).

4. Similarly, at the lower (left) half of the L₁ recirculation bubble, compressibility effects have transformed a mix of axisymmetric contraction and expansion to a predominantly axisymmetric contraction (cyan).

5. Finally, the predominantly axisymmetric contraction (cyan) region at the right-half of the L₁ recirculation layer and further downstream, $x / h ≃ 10 − 24$⁠, has changed to a “mix” of axisymmetric expansion and contraction. This is typically a consequence of the time-averaging window size and indicates that the turbulence anisotropy componentality lies near the center of its contour triangle legend.

The above observations generally suggest that the turbulence in the expansion corner (left) side of the L₁ recirculation region undergoes a shift toward the left of the anisotropy componentality triangle contour legend. In other words, the turbulence in said region undergoes a contraction due to compressibility effects. Reversely, the far (right) side of the L₁ recirculation region undergoes a shift toward the right of the anisotropy componentality triangle contour legend. In other words, the turbulence in said region undergoes an expansion due to compressibility effects.

The thinner “hook”-like anisotropy state structures along the expansion corner walls, as well as the change of the anisotropy states at the center of the L₄ separation bubble, most of the L₂ separation bubble, and more extensively at the L₁ separation bubble half nearer to the expansion corner, are partly due to the decrease in the mean density and thus decrease in mean momentum. Further, examining the turbulent mass-flux components (⁠ $a i ≡ ρ ′ u i ′ ¯ / ρ ¯$⁠)—a measure of the cross correlation between the density and velocity fluctuations—shows that these are prevalent in the two mixing layers, while negligible elsewhere.

The turbulence kinetic energy budgets are examined next to analyze further the physics responsible for the differences in the turbulence properties due to compressibility effects.

We performed an analysis of the turbulence kinetic energy for both Mach-number cases considered, which sheds light on the dominant physical mechanisms that govern turbulence in PSE flows and provide some insight into turbulence modeling. To our knowledge, a similar study on the turbulence kinetic energy budget in PSE flows has not been presented in the past.

The terms in Table IV separately contribute toward the total rate-of-change of the turbulence kinetic energy per unit volume due to: (i) mean flow advection (⁠ $C$⁠), (ii) production due to mean flow deformation and rotation (⁠ $P$⁠), (iii) turbulent transport/diffusion (⁠ $T d$⁠), (iv) pressure diffusion (⁠ $T p$⁠), (v) turbulent viscous diffusion (⁠ $D v$⁠), (vi) turbulence dissipation into heat (ε), (vii) turbulent mass-flux arising from the coupling of density and velocity fluctuations (⁠ $M$⁠), and finally (viii) pressure-velocity correlations (Π_p).

Although the mean spanwise velocity component, $W ¯$⁠, and the corresponding velocity gradient, $∂ ϕ ¯ / ∂ z$⁠, are much lower in magnitude compared to their streamwise and lateral counterparts, all three velocity components and spatial directions are considered. A second-order central difference scheme approximates the derivatives applied throughout the sudden expansion domain apart from the boundaries. Along normal boundaries and overlapping regions (regions where two or more faces of different blocks overlap each other), a first-order difference scheme is used instead.

The actual value of the bulk inflow velocity, U_B, differs for each Mach number considered, i.e., its value is approximately 34 m/s and 170 m/s for the Ma = 0.1 and Ma = 0.5 cases, respectively. Hence, the turbulence kinetic energy budget terms considered are equivalent only in dimensionless space, i.e., are relative to the bulk inflow velocity and density used here to non-dimensionalize them. Consequently, the terms are scaled and compared relative to their ability to affect the mean flow, thus allowing for a more meaningful comparison between the two Mach numbers considered.

The results obtained for all eight budget terms are presented in Figs. 10–17, for each of which sub-figure (a) corresponds to the Ma = 0.1 case, while sub-figure (b) to Ma = 0.5. Starting with the advection budget term (Fig. 10), $C$⁠, it is immediately apparent that there are no evident differences as a result of compressibility effects; even though the mean density varies up to 8% in some regions, the advective flux remains essentially unchanged.

The sharp increase in $C$ at the two—upper and lower—mixing layers is a consequence of the rise in the production of turbulence kinetic energy (⁠ $P$⁠), as evident in Fig. 11. Since the streamlines in Fig. 5(c) show that the separated streamflow does not become narrower, the advection term is responsible for the gradual transfer of $k ̃$ from the mixing layers toward the center of the stream.

Another interesting feature is the two—upper and lower—positive regions that appear to separate outwards from the expansion corner and attenuate thereafter. The above two regions emanate from the interaction between the upper and lower mixing layers with the L₁ and L₂ recirculation bubbles, respectively. Finally, where the separated streamflow impacts the lower wall, $C$ results in a transfer of $k ̃$ from the near-wall region toward the center of the said stream.

The spanwise convection component is over one order of magnitude less than the streamwise and lateral convection terms. This is not surprising since the flow in the spanwise direction is homogeneous and hence mean flow properties will tend to zero as the sampling size becomes large. The same observation had also been addressed by Sideridis et al.⁵⁰ and Liu and Thomas.⁵¹ In both studies, the spanwise component of the convection term was neglected in the calculations due to its small magnitude.

Another important observation is the fact that the streamwise component of the advection term dominates the flow field in the region of the shorter reattachment, L₃ and L₄, mainly due to the sustaining streamwise velocity contribution. On the other hand, the lateral—or wall-normal—component becomes dominant in the larger recirculation zone region, L₁ and L₂, where the transverse velocity gradient is more pronounced.

Kasagi and Matsunaga⁵² carried out a 3D PIV measurement of the turbulence energy budget in a backward-facing step flow at Re = 5540 based on the channel step height and the upstream centerline velocity. They found that at $x / h = 1$ downstream of the expansion, the convection term played a significant role, contributing negatively to the kinetic energy budget. The general shape of the convection term is in excellent agreement with the ILES findings.

Regarding the turbulence kinetic energy production term, $P$⁠, compressibility effects again do not seem to have caused any notable differences. Unsurprisingly, the production of $k ̃$ is most significant at the mixing layers. The same feature has been addressed in several studies on turbulence kinetic energy budgets. Liu and Thomas⁵¹ performed an experimental study in which they measured the turbulence kinetic energy budget of a planar wake flow employing different pressure gradients. They found that the production term is the dominant term in the calculations. The longitudinal shear production $− u ′ υ ′ ¯ ( ∂ U ¯ / ∂ x + ∂ V ¯ / ∂ x )$ is the major contributor to the energy balance with positive values across the turbulent wake. The turbulence production budget term is the most essential, with the ILES results suggesting that it is over three times as large as the next largest, namely, the advection and turbulent diffusion budget terms.

Though $P$ increases rapidly within 1 h from the expansion corner, it begins to attenuate further downstream gradually. The only exception appears to be the apparent merging of the (large positive) production regions at the lower mixing layer and the reattachment near the lower wall. Regarding the latter region, compressibility effects seem to have enhanced the turbulence production.

Kasagi and Matsunaga⁵² also showed that turbulence production is almost three times larger than the turbulence advection and approximately two times larger than the dissipation term. They found that the maximum production rate decreased with increasing distance, and the major contributor over the entire flow region was the term $− u ′ υ ′ ¯ ∂ U ¯ / ∂ y$⁠. These last two observations are in excellent agreement with the results presented in Fig. 11. Regarding the turbulence dissipation, given in Fig. 15, its peak value along the free-shear layers is found to be five to seven times smaller than the turbulence production. This is in contrast to Kasagi and Matsunaga,⁵² who reported the value to be twice as small. The higher Reynolds number—almost twice as large—considered here results in weaker viscous forces that are responsible for the weaker dissipation. We note that past studies have also shed light on the physical vs numerical dissipation^26,53–55 for a variety of flows in the framework of ILES. Such a numerical analysis for the PSE flow case considered is beyond the scope of this paper.

The turbulent diffusion term, $T d$⁠, represents the transfer (advection) of the $u i ″ u i ″$⁠, which constitutes the turbulence kinetic energy, $k ̃$⁠, as a result of the velocity fluctuations $u j ″$⁠. Figure 12 shows that the turbulent diffusion budget term is of the same magnitude as the advection term, i.e., $T d ≈ P$⁠, and influences the same flow regions as the advective budget term. Thus, $T d$ is as influential to the development of the turbulence kinetic energy field as due to its advection by the mean flow. Also, compressibility effects cause no measurable differences again.

In Sec. IV B we mentioned that the transfer of mean momentum from the separated stream flow to the recirculation bubbles via the turbulent fluctuations in the free-shear layers could be—partly—responsible for the larger L₁ size. Here, $T d$⁠, shows just how of an effective mechanism turbulence fluctuations are in the—turbulent—diffusion of flow properties; in this case $k ̃$⁠.

Comparing the production and turbulent diffusion budget terms, given in Figs. 11 and 12, respectively, it is evident that $T d$ takes large negative values in regions at which $P$ takes—in contrast—large positive values. This indicates that the turbulent diffusion is responsible for a significant transfer of turbulence kinetic energy away from the free-shear layers and to the surrounding regions, particularly toward the interior of the separated channel stream.

Post-reattachment, turbulent diffusion remains relatively large, taking a negative value range in the upper-downstream quarter of the L₁ recirculation bubble, while a positive value range in the lower-downstream quarter. In the reattached channel stream, $T d$ becomes predominantly positive, with a much thicker negative region close to the lower wall than the upper wall, which is consistent with the channel section before the expansion corner. We believe that the turbulent boundary layer along the lower wall begins to reach equilibrium sooner than the upper wall, the latter being affected by the much larger L₁ recirculation bubble.

The lateral turbulent diffusion term, $− ∂ ( υ ″ u ″ 2 + υ ″ 3 + υ ″ w ″ 2 ) ¯ / ∂ y$⁠, is found to be the dominant diffusion mechanism, being approximately two times larger than the corresponding streamwise term, $− ∂ ( u ″ 3 + u ″ υ ″ 2 + u ″ w ″ 2 ) ¯ ∂ x$⁠, throughout the sudden expansion domain. The spanwise transport component appeared to be over one order of magnitude less than the lateral diffusion term.

The pressure diffusion, or velocity-pressure gradient term, is represented in Fig. 13. In most experimental studies found in the literature, the pressure diffusion term has either been neglected or extracted by difference. This implies that it is not a directly measurable quantity, and it cannot be obtained from hot-wiring or PIV data. Thus, comparisons of the ILES results with other experimental data cannot be generally made as the error encountered in such calculations is considerably large. Kasagi and Matsunaga⁵² did not consider the pressure diffusion term in their calculations. Instead, it was assumed to be small relative to the other turbulence kinetic energy budget terms and, finally, estimated the dissipation term as the resulting residual. The pressure diffusion term was also neglected in the estimation of the dissipation term in the experimental study of Panchapakesan and Lumley.⁵⁶ 

In the present study, the pressure diffusion (⁠ $T p$⁠) is found to be comparable in magnitude to the other turbulence kinetic energy terms, as shown in Fig. 13. More specifically, it is an order of magnitude smaller than the most crucial budget term, namely, the production term (⁠ $∼ P / 10$⁠), about three times smaller than the advection term (⁠ $∼ C / 3$⁠), and about two and a half times smaller than the turbulent diffusion term (⁠ $∼ T d / 2.5$⁠). This observation agrees well with the results of Sideridis et al.,⁵⁰ who argued that although the pressure-velocity correlation term might be smaller than the other terms in the energy budget, its magnitude nonetheless cannot be considered as negligible. Besides, Le et al.⁵⁷ found that the pressure diffusion is very significant in the near-wall region, where it balances the turbulence transport and dissipation terms. In their experimental measurements, though Liu and Thomas⁵¹ did not calculate the pressure diffusion term directly (it was inferred from the forced balance of the turbulence kinetic energy equation). Nevertheless, they stated that it is not negligible as it facilitates the transfer of turbulence energy from the free-shear layers of the wake to the outer region.

Apart from the free-shear layers, the pressure-velocity correlation term also contributes significantly to the energy budget between $x / h ≈ 2 − 4$⁠, particularly around the reattachment region of the lower free-shear layer corresponding to the shorter recirculation bubble. Thus, a significant amount of energy is transferred from the free-shear layer—where turbulence production is high—to regions where the production rate is relatively low.

Interestingly, the pressure diffusion term takes on values (is non-zero) at the same regions as the turbulent diffusion term, but with one crucial difference; their sign is opposite. In other words, the pressure diffusion term in PSE is to lessen the effect of the turbulent diffusion term. While this might not be the case at higher Reynolds numbers, the effect of the Mach number appears to be trivial, perhaps requiring high enough speeds reaching—or even surpassing—sonic conditions.

Viscous diffusion, $D v$⁠, is a dominant term in the near-wall regions, balancing mainly the turbulence dissipation term, which in turn contributes negatively to the turbulence kinetic energy budget. On the other hand, compared to other budget terms, it becomes considerably tiny in other regions of importance, particularly of high turbulence production. As Fig. 14 indicates, not surprisingly, the viscous diffusion magnitude is most significant in those regions at which the mean flow is under significant shear, i.e., the mixing/free-shear layers and near-wall region. The maximum (positive) value is attained where the lower free-shear layer impacts and reattaches onto the lower wall. Though not evident in Fig. 14 due to the proximity to the wall, $D v$ reaches a value of 0.05 (normalized by $ρ B U B 3 h − 1$⁠). On the contrary, just above the near-wall layer there exists a thicker—and thus more noticeable—negative layer where $D v$ reaches a value of 0.03. The above two values found at the lower wall reattachment region are approximately one order of magnitude larger than those found at the onset of the two free-shear layers at the expansion corner, and over two orders of magnitude larger than those in the remainder of the free-shear layer. Regarding the upper expansion region wall, the magnitude of $D v$ is found to be most prominent near the reattachment of the recirculation bubble (⁠ $x / h = 16 − 18$⁠), but one order of magnitude smaller than at the lower wall (similar to the onset of the free-shear layers at the expansion corner). The above indicates that the viscous diffusion term generally plays a minor role in redistributing the turbulence kinetic energy in PSE flows. The only exception is the near-wall region of the free-shear layer reattachment, particularly associated with the small recirculation bubble.

The dissipation term (ε) contributes negatively to the turbulence kinetic energy budget—it acts as a sink for $k ̃$—and typically becomes considerable in the near-wall region, where the viscous diffusion term balances it. It represents the decay of turbulence kinetic energy (into heat), and for that reason, it is preceded by a negative sign.

An important feature to consider is the ratio of the production to the dissipation term. In experiments carried out, the peak dissipation is approximately 64% and 62% of the production peak at $x / h = 1$ and $x / h = 2$⁠, respectively. For a backward-facing step flow, Lee et al.⁵⁷ found that in the recirculation region, the peak dissipation is about 60% of the production peak, reaching a maximum at the same point as in the free-shear layer. Kasagi and Matsunaga⁵² showed that the maximum dissipation value is approximately 61% of the production peak at $x / h = 1$ downstream of the expansion. Before the reattachment point, at $x / h = 4$⁠, they found that the ratio of the dissipation to the production term increases and approaches the value of $ε / P ≃ 0.70$⁠. In the ILES carried out herein, the dissipation is maximum at the reattachment location of the lower free-shear layer near the wall, i.e., near $x / h = 4$ where $ε ≃ 0.08 ρ B U B 2 / h$ just above $y + = 1$⁠. This suggests a dissipation over production ratio of one $ε / P ≃ 1.0$ at this location. The reason behind this discrepancy is twofold. The numerical calculation allows data to be obtained extremely near to the wall. At the same time, the Reynolds number considered in the ILES here is almost twice that in the experiments.⁵² The contribution of the dissipation term to the energy budget is significant, particularly near the walls and to a lesser extent at the free-shear layers.

Regarding compressibility effects, it is evident that the upper free-shear layer—associated with the large recirculation bubble—remains very similar, albeit slightly larger in extent. However, the dissipation along the lower free-shear layer sustains much longer. For the Ma = 0.1 case, the turbulence kinetic energy dissipation begins to attenuate downstream shortly after it peaks in both the upper and lower free-shear layers. For the Ma = 0.5 case, the dissipation in the lower free-shear layer is shown to persist and form a common “structure” with the increased dissipation located at the reattachment region of the lower free-shear layer.

As mentioned previously, the turbulent mass-flux budget term, $M$⁠, is strongly associated with compressibility effects. Indeed, Fig. 16 reveals that its magnitude is significantly smaller than most other terms, particularly for the Ma = 0.1 case. Specifically, $M$ is over two orders of magnitude smaller than the next weakest term, namely, the viscous diffusion term, while over four orders of magnitude relative to the largest term, namely, the production term. Increasing the bulk inflow Mach number to 0.5 increases the effect of the term by almost two orders of magnitude. Nonetheless, the term remains relatively small in comparison to the production term but becomes of the same order as the viscous diffusion term $D v$⁠. Despite the more substantial compressibility effects, the distribution (or pattern) of $M$ remains identical between the two Mach-number cases considered.

Interestingly, though all other budget terms investigated are asymmetric at the exit of the inflow channel in terms of their magnitude, mass-flux is the only term to be asymmetric in terms of its sign. Moreover, despite $M$ being predominantly adverse (dissipative) at both the upper and lower free-shear layers, it is positive at the upper expansion corner. The culprit is the difference in the sign of the mean pressure streamwise gradient between the upper and lower walls of the inlet channel exit. Above and along the lower wall at the inflow channel exit, the flow is “pulled” downwards and accelerates, thus resulting in a negative mean pressure streamwise gradient. On the other hand, at the upper wall, the flow is “pushed” against the upper free-shear layer and decelerates, thus resulting in a positive mean pressure streamwise gradient. In either case, the turbulent mass flux along both the upper and lower walls of the inlet channel is positive, i.e., $u ″ ¯ > 0$ for $x / h < 0$⁠.

It should also be noted that the mean pressure gradient component, $u i ″ ¯ ∂ τ ¯ i j / ∂ x j$⁠, is found to be several orders of magnitude larger than the mean stress component, $− u i ″ ¯ ∂ p ¯ / ∂ x i$⁠. The only exception being the inner boundary layer region that forms along the (no-slip) walls, particularly the lower wall (associated with the shorter recirculation bubble) and after $x / h = 5$ (post the lower free-shear layer reattachment).

Last but not least, the pressure strain (or dilation) turbulence kinetic energy term, $Π p$⁠, is considered. Figure 17 reveals that its role is predominantly to dissipate the turbulence kinetic energy, particularly at the free-shear layers and lower reattachment region. The pressure strain's largest negative value is $∼ − 0.039 ρ B U B 3 / h$⁠, which occurs at the reattachment region of the lower free-shear layer. A minimum value of $∼ − 0.035$ is reached along the lower free-shear layer itself. These values are about only two to two and a half times smaller than the maximum production value, which indicates $Π p$ plays a significant role in the development of PSE flows.

On the other hand, $Π p$ becomes weakly positive only at one—relatively small—region shortly before the reattachment of the lower free-shear layer, i.e., $Π p ≃ 0.001$ at $x / h ≃ 2.5$ and $d y ∼ 0.05 h$ from the lower wall.

Relative to the Ma = 0.1 case, compressibility effects in the Ma = 0.5 case are found to strengthen the pressure strain term only by approximately 10%. However, similar to the dissipation term, the pressure strain term is also found to attenuate slower along the lower free-shear layer, while intensifying at the lower reattachment region. Consequently, $Π p$ forms a more cohesive structure, having effectively merged the two otherwise distinct regions.

It is also important to point out that all turbulence kinetic energy budget terms decay with streamwise distance x, except those typically associated with the near-wall (viscous) region. Exceptions to the above are viscous diffusion, viscous dissipation, and, to a lesser extent, pressure diffusion. Though this downstream attenuating behavior is anticipated of the budget terms, it might be exaggerated in the present results due to the grid stretching employed, leading to high-aspect-ratio cells (large $Δ x$⁠) that artificially diffuse the turbulence field.

We performed a high-order ILES study of PSE flows and analyzed the turbulent flow characteristics and compressibility effects. We used an 11th-order WENO scheme for the spatial discretization in conjunction with a five-stage fourth-order of accuracy time integration scheme. The primary conclusions of this study are summarized below:

• ILES provides excellent agreement with the experiments.

• The reattachment lengths of the two smaller recirculation bubbles, L₂ and L₃, remain similar in size for the two Mach numbers considered here.

• However, the large recirculation bubble reattachment position, L₁, increases by $∼ 2 h$ (⁠ $∼ 13 %$⁠) for the higher Mach number.

• The PIV data of Emmert et al.²¹ indicate that for a smaller Reynolds number flow—and thus greater viscous dissipation—L₁ will become larger. The ILES data suggest that the compressible case exhibits increased dissipation of the turbulence kinetic energy due to the increased dissipative contributions by the viscous dissipation, turbulent mass-flux, and pressure strain budget terms, particularly at the developing free-shear layers.

• The turbulence kinetic energy streamwise profiles (⁠ $k ̃$⁠) show that compressibility effects lead to a redistribution of $k ̃$⁠.

• Above $y = 0.5 h$⁠, compressibility effects are shown to cause a reduction in $k ̃$ prior to $x / h ≃ 10 − 12$⁠, but an increase thereafter, of up to 35%.

• The L₄ recirculation bubble reattachment position increases with Mach number by $∼ 41 %$⁠. The pressure at the core of the larger neighboring L₃ recirculation bubble becomes significantly lower as the flow develops post the sudden expansion corner. This causes a “suction” effect; effectively “pulling” the smaller L₄ recirculation bubble toward it.

• The Reynolds stress anisotropy componentality is largely unaffected by compressibility at the smaller recirculation bubbles.

• However, above the free-shear layer associated with the larger recirculation bubbles, compressibility effects turn the turbulent structures from an axisymmetric expansion to an axisymmetric contraction in the left half of the L₁ recirculation bubble, while the exact opposite occurs in the right half. The latter does not noticeably influence the turbulence kinetic energy budgets in this region.

Results were obtained using the ARCHIE-WeSt High Performance Computer (www.archie-west.ac.uk) based at the University of Strathclyde.

The data that support the findings of this study are available from the corresponding author upon reasonable request.