Parametric study of multiple shock-wave/turbulent-boundary-layer interactions with a Reynolds stress model

The flow of high-speed air in ducts may result in the occurrence of multiple shock-wave/boundary-layer interactions. Understanding the consequences of such interactions, which may include distortion of the velocity field, enhanced turbulence production, and flow separation, is of great importance in understanding the operating limits and performance of a number of systems, for example, the high-speed intake of an air-breathing missile. In this paper, the results of a computational study of multiple shock-wave/boundary-layer interactions occurring within a high-speed intake are presented. All of the results were obtained using the in-house computational fluid dynamics solver of Glasgow University, HMB3. First simulations of a Mach M=1.61\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M=1.61$$\end{document} multiple shock-wave/boundary-layer interaction in a rectangular duct were performed. The M=1.61\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M=1.61$$\end{document} case, for which experimental data is available, was used to establish a robust numerical approach, particularly with respect to initial and boundary conditions. A number of turbulence modelling strategies were also investigated. The results suggest that Reynolds-stress-based turbulence models are better suited than linear eddy-viscosity models. This is attributed to better handling of complex strain, in particular modelling of the corner separation. The corner separations affect the separation at the centre of the domain which in turn alters the structure of the initial shock and the subsequent interaction. Having established a robust numerical approach, the results of a parametric study investigating the effect of Mach number, Reynolds number, and confinement on the baseline solution are then presented. Performance metrics are defined to help characterize the effect of the interactions. The results suggest that reduced flow confinement is beneficial for higher-pressure recovery.


Introduction
The key to the design of supersonic intakes is compressing a large volume of air with minimal total pressure loss and flow distortion while avoiding the possibility of unstart over Communicated  a wide operating range. This design objective can be achieved in many different ways, for example, using an isentropic compression in a long duct, but the most reasonable, for a practical air vehicle, appears to be splitting the overall compression over a sequence of inclined shock waves. The overall efficiency of the intake then depends fundamentally on the nature of the shocks, their interaction with the intake boundary layer, and the resultant flow. Thus, for most flight conditions, formation of multiple shock-wave/boundary-layer interactions (SWBLIs) inside the intake is expected.
Multiple SWBLIs are often referred to as shock trains or pseudo-shocks. Although these terms are used, at times, interchangeably, the term shock train refers to the series of shocks and the pseudo-shock term refers to the entire region of pressure rise [1]. Referring to Fig. 1a, which shows a sketch of multiple normal SWBLI in a rectangular duct/intake, we can identify a number of distinct regions. Firstly, we have the pseudo-shock region, that part of the domain over which the shock train acts to produce the ultimate compression. Next, we identify a region downstream of the shock train termed the mixing region. The end of this region is often determined by the maximum of the wall pressure.
In addition to multiple SWBLI, oblique SWBLI and reflected shocks are often present within the intake. For the multiple normal SWBLI sketched in Fig. 1a, the first shock imparts an adverse pressure gradient on the incoming boundary layer developing on the walls of the rectangular duct. If the adverse pressure gradient and wall stresses are large enough, the boundary layer separates and a recirculating flow region forms. This region leads to the formation of oblique compression shocks (leading and trailing legs) which eventually join with the first shock to form a λ-shock structure. Higher Mach numbers downstream of the first shock are present near the corners of the duct, as the λ-shock structures on adjacent walls interact [2]. The result of their interaction is a weaker adverse pressure gradient at the corners of the duct compared to the centreline. From the bifurcation (or triple) point formed from the λ-shock structure, a secondary shear layer can develop in the form of a slip line. Downstream of the first shock, the flow is subsonic; however, the local streamline curvature re-accelerates the flow. Due to the curvature, the supersonic flow downstream of the trailing leg of the shock turns towards the wall. Depending on how sudden the turn is (governed by the ratio of boundary layer thickness δ r to the duct half-height h), the flow may again become supersonic and promote the formation of multiple SWBLIs (MSWBLIs). For high Mach numbers (M r ), the effect of flow confinement (δ r /h) on the MSWBLI is less pronounced [1]. At high Mach numbers, MSWBLIs form even under low flow confinement. As the Mach number or flow confinement increases, the SWBLI undergoes a transition from a single normal shock, to multiple normal shocks, to multiple oblique shocks. Multiple oblique shocks are almost exclusively observed at high Mach numbers. The transition from single SWBLI to multiple SWBLI is accompanied by an increase in the interaction length.
The importance of these phenomena, their prediction, and their control are paramount to effective intake design. For this reason, a robust numerical method capable of predicting multiple SWBLIs with an appropriate degree of accuracy is required. Not only must the method be able to adequately resolve shock-wave/boundary-layer interactions, but it should also be able to accurately resolve the corner flows present in ducts. Overestimation of the corner flows can lead to overestimation of pressure oscillations on the wall of the intake. Accurate prediction of the corner flow is also important at high levels of flow confinement. Such flows are known to be challenging for the linear eddy-viscosity models used in many industrial applications, as such models are known to be weak for flows containing complex strain, for example, the shock-wave/boundary-layer interaction itself or in the corner flow. However, higher-fidelity approaches, such as largeeddy simulation (LES), require reduction in the Reynolds number to be computationally tractable, and although they capture the initial pressure rise well, downstream pressure can be underpredicted, most likely due to incorrect boundary layer growth rate. Hybrid Reynolds-averaged Navier-Stokes (RANS)/LES methods such as the scale-adaptive simulation (SAS) and the detached-eddy simulation (DES) in different formulations rely on linear eddy-viscosity models in the RANS regions and suffer from many of the limitations of a pure RANS approach.
The objective of this paper is to explore the utility of nonlinear eddy-viscosity models for this flow. To do this, the effect of Mach number (M r ), Reynolds number (Re h ), and level of flow confinement (δ r /h) on an MSWBLI interaction is studied. Sections 2 and 3 discuss the numerical method and the problem set-up used in this work. The results in Sect. 4 are presented as follows: first, a validation of the approach is presented. This is followed by a parametric study that investigates the effect of lower pre-shock Mach number, lower Reynolds number, lower level of flow confinement, and reduced aspect ratio. Efficiency metrics are calculated for the shock train cases, each at different combinations of outlet pressure, inlet Mach number, and Reynolds number. The data are used to determine under which conditions the resulting shock train achieves the highest pressure recovery and lowest flow distortion, and conclusions are drawn as to which shock train structure is the most suitable for efficient flow deceleration.

Numerical method
Numerical simulations of the steady flow fields have been performed using the Helicopter Multi-Block (HMB3) flow solver [3,4]. The flow solver is a three-dimensional, fully implicit, structured, multiblock unsteady Reynolds-averaged Navier-Stokes (URANS) code solving the URANS equations in integral form using the arbitrary Lagrangian-Eulerian (ALE) formulation for time-dependent domains, which may include moving boundaries. The equations are discretized using a cell-centred finite volume approach. The spatial discretization of these equations leads to a set of ordinary differential equations in time where i, j, and k represent the cell index, W and R are the vectors of conservative flow variables and flux residual, respectively, and V i, j,k is the volume of the cell i, j, and k. The Osher [5] approximate Riemann solver is used to evaluate the convective fluxes, while the viscous terms are discretized using a second-order central difference scheme. The Monotonic Upwind Scheme for Conservation Laws (MUSCL) approach developed by van Leer [6] is used to provide high-order accuracy in space. In regions where large gradients are encountered mainly due to shock waves, the alternative form of the van Albada limiter [7] is activated to avoid non-physical spurious oscillations. An implicit dual-time stepping method is employed to perform the temporal integration, where the solution is marching in pseudo-time iterations to achieve a fast convergence, which is solved using a first-order backward difference. The linearized system of equations is solved using the generalized conjugate gradient method with a block incomplete lower-upper (BILU) factorization as a pre-conditioner [8]. The solver offers several one-equation, two-equation, three-equation, and fourequation turbulence models. In addition, LES, DES, delayed DES (DDES), improved DDES (IDDES), and SAS methods are also available. The fully turbulent k-ω EARSM [9,10] turbulence model is used due to its good agreement with experimental data for this case.

Numerical domain, grids, and boundary conditions
A numerical domain of length L x /h = 46.32 was used for the shock train simulations of the experiment by Carroll [11]. The domain was made dimensionless by the duct half-height at x/h = 0 (h = 16.275 mm). Figure 1b shows a qualitative sketch of the numerical domain. At the inlet, either a uniform or 1/7th profile can be specified. Previous investigations showed that the length of the domain is sufficient to achieve the required boundary layer thickness at x r . If a thicker boundary layer is required at x r , a 1/7th profile can be specified or the domain can be extended. For all cases, a uniform profile was specified at the inlet. The Mach number at the inlet was higher than the Mach number before the start of the interaction to take into account the area reduction due to the boundary layer growth. The outlet pressure was also adjusted to match the experimental pressure at a specific station downstream of the onset of the interaction. The    adjustments of M u and p allowed to match the upstream and downstream conditions as closely as possible. At the outlet, a first-order extrapolation was performed except where the flow is subsonic. There, the outlet pressure p was specified. Adiabatic wall boundary conditions were used for all walls, and since previous investigations showed that the MSWBLI does not exhibit asymmetry, symmetry boundary conditions were applied at the x-z and x-y planes. This resulted in simulating only a quarter of the numerical domain. A previous study using grids A0-A3 [17] showed no significant differences in the wall pressure and the Reynolds stress profiles on finer grids. Grids A0-A2 featured refinement in the streamwise direction and grid A3 featured additional refinement in the wall-normal direction. Based on the results two additional grids B and C, having a similar level of refinement were created. Figure 2a shows grid B (mirrored across the z-x and y-x planes). Also shown in Fig. 2b, c are the cross sections of grid B and grid C at x/h = 0. As shown in Fig. 1b, the grids feature a refinement region in the streamwise (x/h) direction to resolve the pseudo-shock. All grids apart from grid C have an aspect ratio of h/w = 2.34. Grid C has an aspect ratio of h/w = 1.76. Grids B and C have identical streamwise and wall-normal spacings as the A1 grid and feature a longer refinement region in the streamwise direction to account for any shock train movement due to different inlet/outlet conditions and geometry. For all grids, the first grid point was located at y + 1. The average y + at the wall for grid A3 was ≈ 0.1. Table 1 lists the grid parameters. Table 2 lists the parameters for the shock train simulations of the experiment by Carroll [11]. The flow properties at the inlet and the onset of the interaction are reported in the table.

Parametric cases
To investigate the effects of the Reynolds number, Re h , and pre-shock Mach number, M r , on the shock train, the shock train was first simulated at inlet-outlet conditions corre-sponding to the experiment by Carroll [11]. Two additional simulations were then performed at the same inlet-outlet pressure ratio of p/ p u = 2.4776 and lower Mach number M u and Reynolds number Re h . For the lowest M u case, the pressure ratio had to be lowered, as the reference value resulted in an isolator unstart. The differences in the Mach and Reynolds numbers between the reference simulation and the simulations at lower M u and lower Re h were ΔM u ≈ 0.1 and ΔRe h ≈ 4.4 × 10 5 . A simulation at lower flow confinement (δ r /h) was also performed. To decrease the flow confinement (δ r /h) from the reference value (δ r /h = 0.3084), the outlet pressure was increased by ≈ 12 and ≈ 16%. The increase in the outlet pressure relocated the shock train upstream, where the boundary layer was thinner and hence the flow confinement was smaller. An additional case to investigate the effects of increased spanwise confinement (δ r /w) was also performed. The aspect ratio of the duct was reduced from w/h = 2.34 to w/h = 1.76 (25% reduction).
As changes of the inlet/outlet conditions or the geometry result in a different Mach number before the start of the interaction, M r , three additional cases were considered. The Mach number at the inlet, M u , for each case was adjusted until a Mach number before the start of the interaction of M r ≈ 1.61 was obtained. This allowed investigation of the isolated effect of Re h and δ r /h on the interaction. The pressure ratio p/ p u was kept constant at p/ p u = 2.4775 and p/ p u = 2.7876, respectively. Table 3 lists the parameters for the shock train parametric cases. The cases with matching pre-shock Mach number M r are marked in bold. In addition to the above cases, a Latin hypercube approach was used to generate combinations of inlet Mach number M u , outlet pressure percentage η, and Reynolds number Re h . The upper and lower limits for each parameter were set to 1.5 ≤ M u ≤ 3.0, 0.6 ≤ η ≤ 0.9, and 1.0 × 10 4 ≤ Re h ≤ 1.0 × 10 6 , where the outlet pressure percentage (η) was defined as the percentage of the pressure rise across an inviscid shock with a pre-shock Mach number of M u . The use of p and M u as parameters instead of δ r /h and M r was motivated by the fact that matching a specific δ r /h and M r requires iterative adjustments to M u and p which is not practical. Table 3 also lists the parameters for the generated shock train parametric cases.
All simulations were initialized with an inviscid shock having a strength of M u , placed at the end of the domain. At the inlet, the eddy-viscosity ratio for all cases was set to μ t μ = 10. Investigations showed that increasing the turbulence intensity at the inlet does not have a significant effect on the solution. Approximately 5 × 10 4 -1 × 10 5 implicit steps at Courant-Friedrichs-Lewy (CFL) number of 4 were required for the steady-state solutions to converge to at least five orders of magnitude in the flux residuals.

Validation cases
Previous numerical investigations of multiple shock wave boundary layer interactions showed that the fully turbulent k-ω EARSM turbulence model captures the wall pressure well. Figure 3 shows the wall pressure (a) and Mach number contours at the midplane for the A1 (b), A2 (c), and A3 (d) grids obtained with the k-ω EARSM turbulence model. The model slightly underpredicts the wall pressure at the beginning of the interaction which is attributed to the larger separation at the foot of the initial shock. Due to the ability of the model to account for the secondary flows, smaller corner separations are predicted by the model and the wall pressure, although slightly underpredicted, is in good agreement with the experiment. The linear k-ω [18], k-ω BSL [19], and k-ω SST [19] turbulence models were found to overpredict the size of the corner separations leading to oscillations of the wall pressure, not reported in the experiment [17]. The difference in the onset of the interaction x r between the A2 and A3 grids was ≈ 0.7% of the domain length. The obtained wall pressures on all grids agreed well when the axial coordinates were shifted such that x r = 0, where x r is the onset of the interaction. A slight reduction in the boundary layer thickness was present on finer grids. The decrease was turbulence-model-dependent and was more pronounced for the nonlinear k-ω EARSM model. Linear turbulence models showed no significant changes. Figure 4 shows the comparison between the density gradient from the A3 grid and the experiment by Carroll [11]. According to Edney's classification [20], the initial shock forms a type II shock pattern in the experiment. For detailed schematics of the six different shock patterns, refer to the work by Delery [21]. The A3 grid predicts a type I shock pattern, due to the larger separation at the foot of the shock predicted by the nonlinear turbulence model. In contrast, the linear turbulence models predict a type II shock pattern with a Mach stem. However, the wall pressures of the linear models feature oscillations not observed in experiments. Two-dimensional simulations of the problem with linear eddy viscosity turbulence models show good agree-  Table 3 ment in the shock spacing; however, they fail to account for the spanwise effects which are also important.
Several shock train models have been developed over the years including the shockless model by Crocco [22] followed by the diffusion model and modified diffusion models by Ikui et al. [23] and the mass averaging pseudo-shock model by Matsuo et al. [1]. These models can only predict the pressure rise across the shock train. An empirical quadratic equation for cylindrical ducts was developed by Waltrup and Billig [24]. The equation gives the pressure distribution and the distance over which the pressure rise is spread. Further modifications to the equation were made by Billig for square ducts. The equation for the pressure distribution is given by [25]: where θ r is the momentum thickness before the interaction, Re θ r is the momentum thickness Reynolds number, H = 2h is the duct height, and α is either 1/4 or 1/5 for circular or square ducts, respectively. Both circular and square duct equations have been applied. Figure 5 shows the wall pressure predictions. The square duct equation shows better agreement with experiments; however, the equation underpredicts the pressure gradient at the beginning of the interaction.

Parametric cases
The following section first considers the effects of Reynolds number, Mach number, confinement, and aspect ratio on the reference case. The pressure ratio p/ p u was kept constant, while M u , Re h , or w were varied. To investigate the effect of confinement, the p/ p u ratio was varied.

Effect of Reynolds and Mach numbers
The wall pressure and Mach number contours for the lower Mach and Reynolds results are shown in Fig. 6. The inletoutlet pressure ratio for both cases is constant and equal to the pressure ratio of the reference case. The shock train length increased considerably for the case at lower Re h . The Mach Fig. 7 Isosurfaces of M = 1.0 (green) and u/V u = 1 × 10 −3 (blue) for the reference, lower M u , lower Re h , and lower δ r /h cases given in Table 3 Fig. 8 Visualization of the wall shear stress using friction lines at the wall for the cases given in Table 3 number contours indicate that the supersonic core flow region extends further downstream. There is no distinct termination of this region as observed in the other two cases. The onset of the interaction begins at x r /h = 15.3. The difference with the reference case spans about 24% of the domain length. Although the location of the onset moves upstream for the lower Re h case, both the reference and the lower Re h cases have similar levels of flow confinement-δ r /h = 0.311 and δ r /h = 0.307 and the pre-shock Mach number difference between the two cases was about 3%. This shows that for a constant pressure ratio p/ p u and a reduced Reynolds number Re h the confinement is the dominant parameter in determining the onset of the interaction. Figures 7 and 8 show the sonic and separation isosurfaces and the shear stress visualized with friction lines on the wall. Both the reference and the lower Re h cases feature large separation on the top and bottom walls with less pronounced corner separations for the latter. Experimental studies of oblique and normal SWBLIs performed by Dupont et al. [26], Doerffer et al. [27], and Bruce et al. [28] report that the extent of the shock-induced separations at the centreline of the duct is strongly affected by the state of the flow at the corners of the duct. In the experiment by Bruce et al. [28], reduction in the corner separation by upstream suction of the boundary layer resulted in a separated region at the centreline in a previously attached flow field. The case at lower Mach number M u has larger corner separations and smaller separations on the upper and lower walls. The size of the centreline separation is reduced due to the increase in the corner separations which is in agreement with experiments. Both the reference and lower Re h cases feature an oblique initial shock structure. The case at lower M u shows an initial shock with a Mach stem. Such shock structure is usually observed at lower levels of confinement and lower pre-shock Mach numbers M r . Considering the reference case and the cases at lower M u and lower Re h , the pressure recovery p 0 / p 0,u at (x − x r )/h, 20 non-dimensional streamwise distances after the onset of the interaction, is highest for the case at lower Fig. 9 Wall pressure (top) and Mach number contours (bottom) for the reference and reduced δ r /h cases given in Table 3; dashed lines mark the location where the supersonic contours end Fig. 10 Shock train length L/h and onset x r /h with respect to confinement level; L/h is obtained from the midplane slices at the location where the supersonic contours end, dashed lines in Fig. 9 M u (0.825). The difference in pressure recovery between the reference and lower Re h cases is Δp 0 / p 0,u = 0.0052.

Effect of confinement
As the outlet pressure is increased by 12% for the reference case, the onset of the interaction moves upstream, to The movement is equivalent to approximately 23% of the domain length. Figure 9 shows the wall pressure and Mach number contours for the reference and reduced δ r /h cases. The upstream movement of the shock is accompanied by a reduction in the shock train length and by changes of the initial shock structure. The differences in pre-shock Mach number and confinement between the cases amount to about 1.3% and 36.2%. As the confinement is reduced, the shock train becomes shorter and the spacing between subsequent shocks decreases (as observed by Carroll [11]). Further decrease in the confinement reduces the shock train length to approximately L/h = 8. Figure 10 shows the shock train length L/h versus the level of flow confinement and Fig. 11 the wall shear stress. The case at lower M u has similar confinement to the case at lower δ r /h-about 0.19. The difference in the pre-shock Mach numbers is about 7%. Nevertheless, both cases feature Fig. 11 Visualization of the wall shear stress using friction lines at the wall for the cases given in Table 3 Fig. 12 Wall pressure (a) and Mach number contours (b, c) for the reference and reduced w/h cases given in Table 3 shorter shock trains with an initial shock that has a Mach stem.
The pressure recovery for the lower δ r /h case is p 0 / p 0,u = 0.783, again taken at 20 dimensionless distance units after the start of the interaction. Figure 7 shows the sonic and separation isosurfaces for the reference and lower δ r /h cases. For the latter, the corner separations are larger. Since the size of the corner separation affects the centreline separation and the structure of the initial shock, for cases with larger levels of flow confinement (reference and lower Re h ), the corner separations are small, resulting in a larger separation at the centreline and no Mach stem.

Effect of aspect ratio
Reduction in the aspect ratio from w/h = 2.34 to w/h = 1.76 did not have a huge effect on the interaction. Figure 12 shows the wall pressure and Mach number contours for the reference and reduced w/h cases. A small Mach Fig. 13 Isosurfaces of M = 1.0 (shaded green) and u/V u = 1 × 10 −3 (shaded blue) for the reduced w/h case given in Table 3 stem is present for the latter, resulting in the appearance of supersonic tongues and subsonic core flow. The difference between the onset of the interaction for the cases is ≈ 2% of the domain length, and the difference in pre-shock Mach number is ≈ 0.5%. The spanwise confinement, defined as δ r /w where w is the half-width of the duct and δ r is the boundary layer thickness in the midplane perpendicular to the side wall, was δ r /w = 0.13 for the reference case and Fig. 14 Visualization of the wall shear stress using friction lines at the wall for the cases given in Table 3 δ r /w = 0.17 for the reduced w/h case. A 25% reduction of w/h resulted in 30.5% increase in δ r /w. Although the percentage increase is large the actual increase between the cases is not, Δδ r /h = 0.04. As stated earlier for the reduced w/h case, the pre-shock train Mach number and the onset of the interaction do not vary considerably. The main difference observed for the Mach number contours is the appearance of supersonic tongues and a subsonic core (shown in Fig. 12c). From the separation isosurface, shown in Fig. 13, it is observed that the corner separations are larger compared to the reference case (shown in Fig. 7a) and the centreline separation is smaller. Again, the smaller centreline separation results in an initial shock with a Mach stem. The wall shear stress, shown in Fig. 14, is qualitatively similar for both cases. Following the observations made from the separation isosurfaces, the centreline separation is reduced. The corner separations occupy approximately 4.06% and 9.65% of the duct width for the reference and reduced w/h cases, respectively. Figure 15 compares the wall pressure and Mach number contours of the lower δ r /h and Re h cases, both having M r of approximately 1.61. The inlet-outlet pressure ratio was kept constant, while the inlet Mach number varied. Note that for the lower δ r /h case the outlet pressure is 12% higher. As the pre-shock Mach number was approximately M r = 1.61 for all three cases, the sole effect of confinement on the interaction was investigated. For the lowest level of confinement, an initial shock with a Mach stem was observed. As the confinement was increased, the Mach stem reduced in size and disappeared (for both the reference and lower Re h cases). For the case at lower Re h , the confinement was δ r /h = 0.37 which resulted in weak crossing oblique shocks followed by a weak normal shock. From the wall pressure, it was observed that the increasing confinement smeared the pressure gradient mainly by altering the initial shock structure. The increase in the confinement level was accompanied by a reduction in the corner separations and an increase in the shock train length, as seen from Fig. 7a, d. The shock train is longest for the lower Re h case (shown in Fig. 7e). Figure 16b shows that similarly to the cases without matching Mach number, a reduction in the confinement δ r /h leads to an increase in the corner separations and a decrease in the centreline separation. Although slight, increase in the corner separations is also observed for the lower Re h case (shown in Fig. 16c). Figure 17 shows the Mach number contours for the remaining MSWBLI cases-10 to 14. For the shock train at M u = 1.736 (case 11), the initial shock features a Mach stem. The absence of a Mach stem is observed for the shock train at M u = 1.62 (case 10) although the pre-shock Mach number is lower. The difference in the pre-shock Mach number amounts to ≈ 8%. Comparing the confinement ratios shows that case 10 has higher confinement (δ r /h = 0.36) than case 11 (δ r /h = 0.31). The corner separations for case 10 are smaller, giving rise to a larger centreline separation which affects the initial shock structure. Cases 12, 13, and 14 have lower levels of flow confinement than cases 10 and 11. However, the shock train is formed by two crossing oblique shocks followed by a series of normal shocks. According to Matsuo et al. [1] for Mach numbers M r larger than 1.8-2.2, oblique shock trains are mostly observed, depending on the state of the boundary layers. The pre-shock Mach number for cases 12, 13, and 14 falls within this range. Case 12 features an initial shock with a very small Mach stem, whereas case 14 features two crossing oblique shocks terminated by a normal shock. The crossing of the oblique shocks for case 13 is significant, and as a result, there is no subsonic flow downstream of the crossing. In all cases, the downstream shocks in the shock train are concave facing upstream. All cases without a Mach stem do not have supersonic tongues (reference, lower Re h , case 10, case 13, and case 14). Cases featuring an initial shock with a Mach stem and a triple point have supersonic tongues (lower M u , lower δ r /h, case 11). The flow near the slip-line Mach number contours (bottom) for the reference and lower δ r /h and lower Re h cases given in Table 3 Fig. 16 Visualization of the wall shear stress using friction lines on the wall for the cases given in Table 3 emanating from the triple point of the initial shock remains supersonic for longer distances downstream. As these points move closer to the centreline, the core flow after the shock train remains supersonic. The efficiency metrics commonly used for high-speed intakes are the flow distortion (FD), and the total pressure recovery (TPR). The FD and TPR are given by:

Efficiency metrics
where p 0,average is the total averaged pressure at the outlet and p 0,u is the total pressure at the inlet. The effect of pseudo-shock length on the total pressure recovery was investigated by Mahoney [29] who showed that maximum total pressure recovery is achieved when the throat length of the intake equals or is slightly greater than the pseudo-shock length. For maximum total pressure recovery at off-design conditions, e.g., increase/decrease in the design Mach number or nonuniform flow due to change in the angle of attack or sideslip, it was concluded that the length of the throat section should be sufficient to account for changes in the pseudo-shock length. Table 3 lists the TPR and FD efficiency metrics for all shock train cases, and Figs. 19 and 18 show quadratic surface fits for the FD and TPR. For the first 9 cases in Table 3, the efficiency metrics are evaluated at 20 dimensionless distance units after the onset of the interaction (x r ). The TPR Fig. 17 Mach number contours for cases 10-14 given in Table 3 Fig. 18 Quadratic surface fit for FD; R 2 = 0.982 for the reference case is approximately 11.9% lower than the reference TPR, taken as the TPR across a normal shock with a pre-shock Mach number of M r . As the confinement is decreased (lower δ r /h case with matching M r ), the TPR increases and the flow distortion decreases). The lower δ r /h case with matching M r has an initial shock with a Mach stem due to the smaller centreline separation and larger corner separations. The extent of the separation in the streamwise direction is about L sep /h = 1.30. Both the reference and low Re h cases have larger separations in the streamwise direc-tion, of L sep /h = 1.55 and L sep /h = 1.61, respectively. The increase in the separation is accompanied by a decrease in the TPR and an increase in the FD. As the separation is affected by the corner separations and itself affects the initial shock structure, the cases with smaller centreline separations exhibit higher TPR. Such cases are the lower δ r /h with and without matching M r , and the lower M u case. For the lower δ r /h case with a matching M r , the TPR increases by 6.2% and the FD decreases by 17.8% compared to the reference case. Similarly, FD and TPR are evaluated at 10 dimensionless distance units after the onset of the interaction for cases 10-14 in Table 3. As the M u , Re h , and p vary for each case, it is more difficult to draw similar conclusions. However, the variation of Re h is moderate, and all cases exhibit reduced TPR and increased FD. On average, the TPR is approximately 10% lower than the reference TPR for cases 10-14. Case 10, being at a similar Mach number to the reference case, shows that larger separation resulting in a triple point closer to the centreline reduces TPR and increases FD. As the Mach number is increased for cases 11-14, the TPR continues to decrease and the FD to increase. Even though case 11 features a triple point below the centreline, due to the higher Mach number the TPR and FD are worse than for the lower M u or lower δ r /h cases. For cases 13 and 14, the FD is greater than one. The increase in the FD is due to the decrease in the minimum stagnation pressure   (13 and 14), the difference between the maximum and minimum stagnation pressures at x/h = 10 is greater than the average stagnation pressure which results in FD > 1.
Flow distortion was also evaluated for the reference, lower Re h , and lower M u cases at seven streamwise stations. The first station coincides with the onset of the pressure rise, x r /h = 0. Subsequent stations are spaced at equal distances downstream with the last station placed at x/h = 10. Figure 20 shows that the lowest M u case has the lowest FD among the three cases. For the shock trains featuring an initial shock with a Mach stem, the FD has the highest rate of decrease. No significant differences in the rate of FD decrease or the FD itself are observed between the reference and the lower Re h cases which both have larger centreline separations and an initial shock without a Mach stem. Shock trains of this type were systematically observed to have a reduced TPR and increased FD. To minimize FD and maximise TPR, one should aim for a shock train featuring an initial shock with a Mach stem, because of the smaller separation at its foot.
Shock trains were observed for all combinations of Mach numbers, confinement levels, and Reynolds numbers. The lowest M u (M u = 1.49) case featured the shortest shock train. Further increase in the outlet pressure by approximately 3% for the lowest M u case resulted in a single shock. This shows that even for low pre-shock Mach numbers in the range of 1.44-1.5 and confinement levels greater than δ r /h = 0.15 shock trains are still present.

Conclusions and future work
The following conclusions have been drawn from this study: -The most suitable approach for simulating shock trains is to adjust both the Mach number at the inlet of the domain and the pressure at the outlet of the domain to match the Mach number M r and confinement δ r /h before the interaction.
-Linear eddy-viscosity models exhibit overprediction of the corner separations. Since corner separations affect the separation at the centreline, the overprediction often leads to the stronger flow re-acceleration resulting in pressure oscillations at the wall.
-Due to its capability of resolving the secondary (corner) flows, the k-ω explicit algebraic Reynolds stress model (EARSM) predicts smaller corner separation and results in a slight underprediction of the wall pressure with no wall pressure oscillations.
-For the reference case considered, reduction in confinement resulted in a shorter shock train and larger corner separations. The larger corner separations lead to a smaller separation at the centreline and an initial shock with a Mach stem.
-Reduction in the upstream Mach number has a similar effect, with the effect of Mach number on the length of the shock train being more pronounced.
-An opposite trend is observed for the reference case and the reduced Reynolds number cases. For both cases, the confinement is larger with smaller corner flows and larger centreline separation leading to an initial shock without a Mach stem.
-The core flow after the shock train, of the cases featuring an initial shock without a Mach stem, is supersonic. Cases featuring an initial shock with a Mach stem have supersonic tongues and subsonic core flow.
-The total pressure recovery (TPR) and flow distortion (FD) efficiency metrics show that shorter shock trains (at lower Mach numbers and lower confinement) featuring an initial shock with a Mach stem are more efficient, resulting in higher TPR and lower FD.
Future work includes investigations of shock trains in a more realistic intake geometry with a store forebody.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.