On the investigation of oblique shock‐wave/turbulent boundary‐layer interactions with a high‐order discontinuous Galerkin method

Shock‐wave/turbulent boundary‐layer interactions are still a challenge for numerical simulation. The shock capturing needs dissipation to avoid spurious oscillations while turbulence will be falsified by introducing dissipation. Especially, an accurate prediction of quantities such as the skin‐friction coefficient inside the interaction area of shock wave and turbulent flow is a critical point. In this article, we investigate a wall‐resolved large eddy simulation of oblique shock‐wave/turbulent boundary‐layer interactions by a high‐order discontinuous Galerkin scheme. The high‐order scheme handles the turbulent flow very well. The shock‐capturing is confined to the near shock region by switching locally to a finite volume second‐order TVD scheme on subcells. This strategy is completed with the application of a shock indicator to a filtered flow field. A global spanwise filter is applied to avoid switching on the shock‐capturing procedure in regions of under‐resolved turbulent structures. We validate our numerical results first at shock‐wave/laminar boundary‐layer interaction. The main simulation under consideration is a Mach 2$$ 2 $$ turbulent boundary‐layer with an inlet momentum‐thickness Reynolds number of 1628$$ 1628 $$ , interacting with an oblique shock that deflects the incoming flow by 8∘$$ {8}^{\circ } $$ . We employ a reformulated synthetic eddy method at the inlet to avoid the influence of recycling‐based turbulence generating schemes on the low‐frequency unsteadiness. The anisotropic linear forcing technique is adopted to further reduce the turbulence recovery length. Through the spectral analysis of wall pressure probes, a typical Strouhal number of around 0.03$$ 0.03 $$ is observed. We attribute the discrepancies between an experimental scaling law and our computation to the three‐dimensional sidewall effects in the experiment. With the assistance of numerical results from this article and other authors, a new scaling law for the spanwise‐periodic computations is suggested to quantify the difference between experimental and computed data.


INTRODUCTION
Shock-wave/boundary-layer interactions (SWBLIs) appear in fluid flow around vehicles in transonic, supersonic and hypersonic flow regimes. Prominent examples are transonic airfoils, supersonic air intakes and hypersonic air-breathing engines. In most cases, SWBLIs are accompanied by significant changes of pressure and temperature loadings. These phenomena may have a remarkable impact on the integrity and the fatigue of structures. 1,2 According to Dolling's famous review paper, 3 SWBLIs have been studied for more than 70 years. A wealth of review papers have been published after Dolling. [4][5][6][7] While substantial progress has been achieved in the last 70 years, there are still some open questions on SWB-LIs. For example, Gaitonde 7 argues that additional research is required in complex configurations, flow controls, heat transfer calculations, and skin friction predictions. One of our interests is the investigation of the accuracy in predicting the skin friction in SWBLIs with a high-order discontinuous Galerkin (DG) method. The subject of our study is one of the benchmarks in SWBLIs, that is, the oblique shock reflection over a flat plate.
Many successful experiments and computations for two-dimensional (2D) oblique shock-wave/laminar boundary-layer interactions (SWLBLIs) have been conducted in the last century, such as the experiments of Hakkinen et al. 8 and Degrez et al., 9 and the computational simulations of MacCormack 10 and Katzer. 11 In 1998, Knight and Degrez 12 stated that the aerodynamic and thermal loads of SWLBLIs could be accurately predicted with existing CFD technologies. However, the limitation of numerical accuracy at that time might had a negative influence on the results. For example, in order to ensure numerical stability, Katzer 11 introduced a numerical dissipation term into his oblique SWLBLI computation, and he found that the dissipation might impose an undesirable influence on the height of the separation bubble. Therefore, it is necessary to restudy the SWLBLIs with modern high-resolution, low-dissipation methods.
Three-dimensional (3D) oblique shock-wave/turbulent boundary-layer interactions (SWTBLIs) remain a challenging simulation task today. Oblique SWTBLIs can be categorized into two types according to shock reflection configurations. One is named regular reflection and the other is Mach reflection. 13 In general, the Mach reflection happens, where a large incident angle is imposed, small to moderate incident angles occur with regular reflection. Oblique SWTBLIs can also be classified as weakly or strongly interacting flows depending on the occurence of separation inside the turbulent boundary layer (TBL). 6 In this work, we focus on the regular reflection of the oblique SWTBLI with separation. For simplicity, we will refer to it as "oblique SWTBLI" in the rest of this article. The flow separation induced by strong shocks and the resultant flow unsteadiness are two dominant features of oblique SWTBLIs. As a result, without further improvements, time-averaged simulations based on the Reynolds-averaged Navier-Stokes (RANS) equations are not capable of accurately predicting even the mean-flow quantities inside the separated areas. 3 Methods with time-dependent modeling techniques like direct numerical simulation (DNS), large eddy simulation (LES) and hybrid LES/RANS are required for accurate oblique SWTBLIs simulations. Although wall-resolved LES are generally confined to small to moderate Reynolds numbers, it has been popular owning to its compromise between computational costs and accuracy. An early successful attempt to apply LES to the oblique SWTBLI has been conducted by Garnier and Sagaut. 14 More recent works on the LES of oblique SWTBLIs can be found in the work of Touber and Sandham 15 and Pasquariello et al. 16 Among the past LES computations of SWTBLIs, many adopted low-order spatial and temporal discretization schemes. These schemes are known to be dissipative without special treatments, and therefore might smear out important flow features, particularly in the near-wall interaction regions. This article seeks to apply a high-order DG method for the implicit large eddy simulation (ILES) to oblique SWTBLIs. The term ILES means no explicit subgrid scale dissipation model is employed in the governing Navier-Stokes equations to account for the influence of unresolved small turbulent scales, while we solely rely on the dissipation of interface Riemann solvers. This implicit approach has been widely tested and accepted within the DG community, and excellent turbulence results can be obtained. 17 Another deficiency of many previous LES for SWT-BLIs is that they adopted the recycling/rescaling method by Lund et al. 18 to generate the turbulence at inlets, see, for example, the simulation of Garnier and Sagaut. 14 It is well known that the recycling/rescaling technique introduces extra frequencies that can contaminate the separation-induced unsteadiness in SWTBLIs. 19 An ideal choice for the generation of the inflow turbulence are synthetic turbulent inflow approaches, for example, the digital filters method used by Touber and Sandham 15 and the synthetic eddy method (SEM) used by Agostini et al. 20 Compared with digital filters methods, SEMs are generally believed to have a shorter turbulence recovery length. 21,22 To further reduce the recovery length, we combine the SEM with the anisotropic linear forcing technique by de Laage de Meux et al. 23 As mentioned by Gaitonde,7 one of the future topics in SWTBLIs is the flow control. Understanding the contributing factors of SWTBLIs and quantifying their influences are the prerequisites for flow control. To reconcile the big variations in the interaction geometries and the aerodynamic conditions of past experiments, Souverein et al. 24 proposed a scaling law to quantify the relation between the interaction length and the separation criterion. Another significant purpose of the scaling law is to establish a practical criterion to assess the quality of numerical computations and help to identify possible numerical problems. However, as they are limited by the computational resources, DNS and LES usually adopt much lower Reynolds numbers than experiments. Moreover, it is difficult to accurately imitate the 3D sidewall effects of wind tunnels, whereas the 3D sidewall effects have been proven to have a marked impact on the interaction length scales. 25,26 Hence, an evaluation of the applicability of the experimental scaling law by Souverein et al. 24 into spanwise-periodic computations becomes necessary.
In recent years, high-order DG methods for the ILES of incompressible and compressible turbulent flows have aroused increasing interests. 27,28 DG methods can easily achieve high-order spatial and temporal accuracy on fully unstructured meshes. Moreover, high-order DG methods have been proven to have favorable dissipation and dispersion errors for turbulence simulations because numerous small-scale structures can be resolved. 17,27 The DG method in spectral collocation form combined with an explicit time integration is standout due to its high efficiency for massively parallel computing. Our code FLEXI, considered in this article, is an open-source CFD code based on the discontinuous Galerkin spectral element method (DGSEM) in space and a high-order explicit Runge-Kutta scheme in time, and it is designed to fully explore modern parallel architectures. 29 The main difficulty of applying the high-order DGSEM to the LES of under-resolved turbulence simulations is how to ensure computational robustness without excessive artificial damping. Gassner et al. 30 has addressed this problem by introducing the split forms of the advective fluxes. Subsequently, Flad and Gassner 17 introduced this method into FLEXI and demonstrated that this dissipation-free base line DGSEM could restore the nonlinear stability and high accuracy in under-resolved simulations by adding dissipation that mimics the physical subscale turbulent dissipation either by an explicit subgrid scale model or by a Riemann solver at element interfaces. However, if there exist strong discontinuities like shocks inside flow domains, the high-order split-form DGSEM still suffers the stability problem. To resolve this problem, Sonntag and Munz 31 implemented a hybrid DG/FV scheme into FLEXI to take advantage of the excellent shock-capturing property of FV schemes without losing the high order accuracy in smooth regions.
Low-order FV schemes introduce more numerical dissipation than the high-order DG method. In order to maintain a high-fidelity solution, the FV scheme is applied on subcells within a regular DG grid cell. Nevertheless, the total number of the FV grid cells in computational domains has to be as few as possible, especially in the near-wall regions. A difficulty is that the usual shock indicators cannot distinguish between the discontinuities of genuine shocks and the shocklets due to the under-resolution inside the TBL. The most popular shock indicator in SWTBLI is the Ducros indicator, 32 which separates shocks and turbulent vortices by evaluating whether the dilatation or the vorticity is the dominant factor. Recent improvements of the original Ducros indicator mainly focus on the decrease in its sensitivity to dilatation fluctuations outside the TBL, 33,34 and the increase in its capability to separate shocks and expansion regions. 35 However, due to the similar magnitude of the dilatation and the vorticity near the interaction area, it might be difficult for the Ducros indicator to distinguish the shocks inside the TBL. Therefore, adopting a Ducros-like indicator could lead to either nearly no detection of the shocks inside TBLs, or too many grid cells inside TBLs will be labeled as shocks. 36,37 The classic Jameson indicator 38 is seldom chosen due to its incapability of separating shocks and under-resolved turbulent structures. As a result, the excessive dissipation introduced for shock capturing will largely contaminate the solution in the interaction areas. 33 In this work, we propose to use spatial filters first to exclude the influence of turbulent structures while retaining the solution jumps due to shocks. In this way, any indicators, for example, the Jameson indicator, which can detect the solutions jumps, can easily be applied to obtain sharp shock profiles inside and outside TBLs. Nevertheless, a new stability problem might appear due to the frequent switching between the DG and FV grid cells in the interaction area. In order to stabilize the computation, the positivity preserving limiter (PPL) proposed by Zhang and Shu 39 is adopted. To the best of our knowledge, this is the first time investigating the high-order DGSEM for the ILES in oblique SWTBLIs. This article also aims to demonstrate that our high-order DG scheme is suitable for the complex simulation of the oblique SWTBLI. Finally, our validated results will be applied to assess the effectiveness of the scaling law by Souverein et al. 24 for quasi-2D configurations.
The format of this article is structured as follows. In Section 2, the governing equations are described. The numerical approaches, including the DGSEM framework and the hybrid DG/FV shock capturing scheme are presented in Section 3. The details of the modified Jameson shock indicator are provided in Section 3.2. Section 4 follows with the computational set-up and the results of the oblique SWLBLI. Subsequently, we perform a 3D oblique SWTBLI simulation in Section 5. Here, we validate the incoming TBL before presenting the instantaneous and mean flow of the oblique SWTBLI. Moreover, a spectral analysis technique is employed to examine the low-frequency unsteadiness of the oblique SWTBLI in Section 5.3.2. The evaluation of the experimental scaling law proposed by Souverein et al. 24 on the applicability into the spanwise-periodic computations is provided in Section 5.3.3. In the end, we summarize our results and conclusions in Section 6.

GOVERNING EQUATIONS
We solve the 3D compressible Navier-Stokes equations in conservative form where the solution ⃗ U denotes the vector of the conserved quantities ⃗ U = ( , v 1 , v 2 , v 3 , e) T , the subscript t represents the derivative with respect to time, and ∇ ⃗ x is the gradient operator in the physical coordinates ⃗ with i = 1, 2, 3. Here, represents the Kronecker delta. The viscous stress tensor ij is defined according to the Stokes hypothesis for a Newtonian fluid as where is the dynamic viscosity, obeying the Sutherland's law with the heat conductivity k The perfect gas law is used to close the system of equations. We follow the usual nomenclature for , (v 1 , v 2 , v 3 ) T , p, T, e, c v denoting density, velocity vector, pressure, temperature, specific total energy and specific heat at constant volume, respectively. In this work, the air is modeled with specific heat ratio = 1.4, the specific gas constant of R = 287 J (kg K) −1 , and a constant Prandtl number Pr = 0.72.

Discontinuous Galerkin spectral element method
To solve the system of governing equations (1), we subdivide the domain into nonoverlapping, unstructured, hexahedral elements with optionally curved surfaces to account for complex geometries, and then map each element onto a reference unit cube E(⃗) = [−1, 1] 3 . According to the mapping function ⃗ x(⃗), we get the corresponding Jacobian determinant , which will be used for transforming governing equations (1) into the reference space where ⃗  is the flux vector in the reference space. Afterwards, we multiply the transformed governing equations (8) by a test function Φ and integrate over the reference space E. Finally, we obtain the DG formulation of NSE in weak form using the divergence theorem, where the second and third integrals represent the surface and volume fluxes, respectively. In the second integral, ⃗  ⋅ ⃗ n denotes the outward pointing surface normal flux, which is composed of ⃗  ⋅ ⃗ n = ⃗  n − ⃗  n , with ⃗  n the element surface normal advective flux and ⃗  n the element surface normal viscous flux. For DGSEM, we employ a tensor-product basis N ijk of 1D Lagrange polynomials l N i of degree N, evaluated at Gauss or Gauss-Lobatto points, to approximate the solution and the three contravariant fluxes withÛ ijk (t) and m ijk (t) the time-dependent nodal solution and the flux values inside the reference element E. Furthermore, following the Galerkin ansatz that the test functions are chosen identically to the basis functions, we get the collocation of the integration and interpolation points. This leads to a significant reduction of operations for 2D and 3D cases, and finally results in a highly efficient numerical scheme. More detailed descriptions of DGSEM can be found in the book of Kopriva 40 and in the recent review paper of FLEXI by Krais et al. 29 In FLEXI, the element surface normal advective flux is approximated by Riemann solvers. Several well-known Riemann solvers have been implemented into FLEXI. Here, we choose the approximate Riemann solver with entropy fix by Harten and Hyman. 41 The element surface normal viscous flux is calculated with the lifting method proposed by Bassi and Rebay 42 commonly referred to as BR1 scheme. The five-state, fourth-order low-storage Runge-Kutta method of Carpenter and Kennedy 43 is adopted as the temporal integration scheme.
In order to stabilize the under-resolved turbulence simulation, the split form of DGSEM is utilized. This special DGSEM derives from the DG formulation of NSE in strong form on Legendre Gauss-Lobatto nodes. To obtain the strong formulation, a second employment of the divergence theorem is required in Equation (9). Afterwards, the two-point flux functions, for example, the formulation by Pirozzoli, 33 will replace the advective part of the standard physical flux function in the volume integral. More information about this split DG method can be found in the papers of Gassner et al. 30 and Flad and Gassner. 17

Shock-capturing method
In DG methods, discontinuities across element interfaces are allowed since the fluxes at grid element interfaces are calculated by Riemann solvers. Nevertheless, inside each element, the solution is approximated by high-order polynomials, which may result in instabilities at strong gradients. To handle the spurious oscillations of the high-order DGSEM solution inside the grid cells containing strong discontinuities, Sonntag and Munz 31 implemented into FLEXI an efficient shock-capturing scheme that is based on a finite volume approximation on subcells. For the integrity of this article, a brief summary of the hybrid DG/FV shock-capturing scheme is provided in the following paragraphs, where we mainly concentrate on the shock indicators. The general idea of the hybrid DG/FV shock-capturing is to replace each "troubled" DG element by a subcell element with equally-spaced subcells and to apply there an FV approximation. In order to maintain the same degrees of freedom (DoFs), the number of the FV cells inside each FV subcell element is determined by the polynomial degree of the original DG element. Take a 2D DG element with Gauss points and polynomial degree of N = 3 in the reference space E = [−1, 1] 2 as an example, see Figure 1. If the DG element is labeled as a "troubled" element, then it will be split into (N + 1) 2 equidistant FV cells. Each FV cell will then be resolved by a robust second-order FV scheme to stabilize the simulation. This switching to the subcells introduces a local grid-refinement and a better localization of shock waves. However in smooth parts of the flow, the second-order FV approach schemes introduce more numerical dissipation than the high-order DG method. Hence, it is essential to assure that the FV subcell elements appear only near the shocks. This is especially important for turbulent flows as the increased dissipation and the low accuracy of the second-order FV method may negatively influence the under-resolved turbulence. A shock indicator capable of accurately distinguishing between shocks and under-resolved turbulence is the key point for stable results with high resolution.
Several shock indicators to detect the "troubled" DG elements have been implemented into FLEXI. An indicator originates from the switching function of a dissipation term in the context of FV schemes by Jameson et al. 38 has been modified into our DG scheme. In this article, we call it the Jameson indicator. The indicator value for each node inside a DG element is of the form where p min,ijk is the minimal pressure of the neighboring nodal values of the node ijk and p max,ijk is the respective maximal value. By use of the indicator values at the nodes, the Jameson indicator value for a DG element reads as where V E denotes the volume of the DG element, and V ijk the volume of the ijkth FV sub cell. From Equations (12) and (13), it is apparent that the Jameson indicator detects shocks by checking the surrounding pressure difference of a node. The Jameson indicator should be suitable for laminar flows because no strong discontinuities exist in flow field except for the shocks. However, the Jameson indicator cannot distinguish the under-resolved turbulent structures from the shocks in the oblique SWTBLI, as shown in Figure 2, where red regions denote the FV subcell elements, and blue regions denote the DG elements. In Figure 2, except for the shock areas, a large amount of grid cells inside the TBL are also labeled as "troubled" ones. Therefore, without improvement, the Jameson indicator cannot be used in the computation of oblique SWTBLIs due to the excessive dissipation introduced by the FV subcell elements inside TBLs. According to the numerical test of Pirozzoli, 44 the Ducros indicator is capable of labeling real shocks. However, from the figure 4 in Pirozzoli, 44 we discover that the shock profile is totally outside the TBL, which for sure requires either a refined mesh or other special treatment to resolve the shocks existing inside the TBL. This problem becomes more obvious for strong SWTBLIs as the shock systems will penetrate deeply into the TBL. 16 The original Ducros indicator 32 is adapted for our DG framework. At a node inside a DG element, the Ducros indicator value is where ∇ ⋅ ⃗ v ijk is the dilatation, ∇ × ⃗ v ijk is the vorticity, and is a constant that prohibits division by zero. The  Ducros (i, j, k) in Equation (14) is dominated by the dilatation magnitude near shocks, and the vorticity magnitude near vorticies. The indicator value of each DG element is obtained from the weighted mean value inside the element, where V E and V ijk are the same as Equation (13). From Equation (15), we can find that the indicator value  Ducros is close to zero inside undisturbed TBLs and approaches one near inviscid shocks. From the instantaneous distributions of the dilatation in Figure 3 and the vorticity in Figure 4 of the oblique SWTBLI, we know there are roughly three ranges of the Ducros indicator values. One is outside the TBL, where the indicator value goes to one near shocks and ideally close to zero for the other areas. It must be noticed that due to the negligible vorticity outside the TBL, the Ducros indicator becomes very sensitive to any dilatation outside the TBL, and may mark the smooth regions as "troubled" areas. One easy solution is setting = (U ∞ ∕ 0 ) 2 as suggested by Pirozzoli,33 so that the influence of the dilatation with small amplitude is filtered out. The second range of the Ducros indicator value is inside the TBL but far from the interaction area containing shocks, where the indicator values are dominated by the vorticity, and close to zero. The third range of the indicator value is in the interaction area containing shocks: both the dilatation and the vorticity can be dominant due to their similar order of magnitudes. Furthermore, the separation-induced shock motion makes it more difficult to separate the shocks and the vortex structures in the interaction area. Generally there exist two strategies to choose a threshold value of the Ducros indicator for the switching between DG and FV subcell elements, either a large enough threshold value to make sure that the interaction area has no "troubled" cells, see for example the figure 15 in Fang et al., 36 or a small threshold value so that a large area close to shocks are marked as "troubled" fields, see for example the figure 1 in Jammalamadaka et al. 37 Most people have chosen the former strategy for better accuracy inside TBLs, but refined meshes are required to resolve the shocks inside the interaction area, especially while high-order methods are being adopted. 15,16,44,45 Another advantage of the former choice is that no DG/FV switching happens inside the interaction area, so that the stability problem due to the frequent switching will not appear. The latter strategy does not have too much advantages except that a coarse mesh inside the interaction area is allowed for better efficiency, but the substantial damping introduced in the  37 we can also discover that even the expansion fan regions are also marked as "troubled" area by adopting the modified Ducros indicator by Pirozzoli. 33 The erroneous detection of the expansion fan is difficult to skip if one simply uses the Ducros indicator or the modified Ducros indicator by Pirozzoli. 33 The reason is that the magnitude of dilatation square ( ∇ ⋅ ⃗ v ijk ) 2 of the expansion fan might be very close to the corresponding magnitude values of some weak shocks. One simple solution to exclude the wrong detection of the expansion fan is to convert the modified Ducors indicator by Pirozzoli 33 as which is organized in the context of our DG scheme but can be adapted to any other numerical method. The idea behind this modification is to take account of the physical fact that the dilatation of the expansion fan is positive, while the dilatation of shocks or compression waves are negative. This idea has already been successfully applied in the improvement of the original Ducros indicator by Larsson et al. 46 and Lusher and Sandham. 35 With the improved Ducros indicator in Equation (16), we find that it is still difficult to accurately capture shocks inside the interaction area in some cases. Some extra effort is required to further evaluate the magnitude of the dilatation and the vorticity in the interaction area for each specific case. Afterwards, a further small modification of the improved Ducros indicator in Equation (16) might be made as where the additional factor represents a scaling of the vorticity. The factor can used to further separate the effects of the dilatation and the vorticity in the interaction area. The modification on the Ducros indicator in Equation (17) can improve the shock-capturing accuracy, but it is not easy to be generalized for different cases.
In this article, we propose another simple way to separate the solution jumps from the vortex structures in the TBL and the shocks by filtering out the vortex structures firstly. Afterwards, we evaluate the Jameson indicator on the filtered flow field. Possible filters includes temporal filters and spatial filters. The biggest drawback of the temporal filters is that they do not allow a real-time evaluation of the indicator value. In the oblique SWTBLI simulation, an ideal shock indicator should be able to update the shock profiles along with time to account for the shock motions. One possible choice to use the temporal filter is to apply the Jameson indicator on a moving-averaged flow field, which is a kind of low-pass filter. Due to the fact that the shock profiles update with an apparent low frequency, whereas the free shear layers of the separated flow and the mixing layer in the TBL are oscillating with a medium to high frequency, a low-pass temporal filter is theoretically feasible to maintain sharp shock profiles along with time. However, according to our survey, it is difficult to obtain a sharp shock profile while keeping the computation stable. One possible solution is to labeling neighboring grid cells of a "troubled" element as FV subcell elements to maintain numerical stability, but obviously extra numerical dissipation will be introduced in this process. In comparison with the temporal filters, the spatial filters are more suitable for the oblique SWTBLI simulations. A simple spatial filter can be a globally spanwise average process, or for the high-order DG method, a more flexible choice can be the element-local spanwise average process. High-order DG methods generally have large element sizes, for which the element-local spanwise averaging can also significantly reduce the negative effects of the vortex structures. The average can be the mean value of all spanwise nodes, or the root-mean-square value of all spanwise nodes, and both work as low-pass filters. In our paper, we simply use the mean value of all spanwise nodes. Since our oblique SWTBLI simulation is a quasi-2D case, the global spanwise filter is more appropriate. We emphasize that the global and element-local spatial filters are not only limited to the Jameson indicator, but work for any other shock indicator, which can detect solution jumps, for example, the Persson indicator. 47 The idea behind the spatial filters is to exclude the influence of under-resolved turbulent structures while retaining the solution jumps nearby shocks so that we are able to locate all shocks accurately inside and outside the TBL. Here, we validate the effects of the global spanwise filter in the pressure field. In the instantaneous pressure field in Figure 5, we can see that many vortex structures exist inside the TBL. However, after the global spanwise-average procedure, we almost eliminate all the vorticies inside the TBL, see Figure 6, and thus successfully exclude the discontinuities due to the under-resolution in the TBL. Another favorable property of the global spanwise-average procedure is that even other fluctuations outside the TBL will be erased in this process. However, the capturing of unsteady shocks inside the TBL might lead to a stability problem. In order to stabilize the computation, we introduce the PPL of Zhang and Shu 39 into our DG framework. One of the biggest advantage of the PPL is that it can maintain the stability of DG methods without losing high-order accuracy. With the Jameson indicator, evaluated on the globally spanwise-averaged pressure field, we are able to capture all the shocks in the flow domain accurately, see the Figure 7, in which the red regions represent the FV subcell elements, and the blue regions represent the DG elements. In comparison with the distribution of the FV subcell elements and the DG elements based on the original Jameson indicator evaluated in the instantaneous pressure field in Figure 2, we can see a dramatic improvement on the shock capturing accuracy. The negligible amount of the FV subcell elements near the interaction area in Figure 7 ensures as less numerical dissipation as possible. In addition, comparing to the popular choice of the Ducros indicator with a big threshold value, our current shock indicator can easily capture the shocks inside the TBL. The limitation of the adopted global spanwise filter is that it only works in quasi-2D configurations such as the supersonic flows over compression ramps and the oblique shock reflection. The element-local spatial filter is more flexible and can be applied to some special 3D configurations, where a spanwise or approximately-spanwise direction can be defined for averaging process. For complex 3D configurations, considering the fact that the shocks generally have larger scales than the turbulence, the idea of averaging turbulence using a geometrical or temporal filter before calculating shock indicators should be feasible as well. However, much effort would be required to search a proper way of filtering the instantaneous solution. In this article, we focus on the quasi-2D oblique SWTBLI case. The extension of the idea of filtering to complex geometries is out of the scope of this article and is left for future considerations.

TWO-DIMENSIONAL OBLIQUE SHOCK-WAVE/LAMINAR BOUNDARY-LAYER INTERACTION (SWLBLI)
In this section, we employ the 2D simulation of an oblique shock wave with a laminar boundary layer performed by Katzer. 11 The numerical results of Katzer 11 have been extensively used for the validation of the numerical approaches in the last decades. The same simulation is computed with our DG method to present the desirable accuracy and efficiency.

Numerical set-up and boundary conditions
The incoming free-stream Mach number is Ma = 2.0. The free-stream temperature is T ∞ = 288.15 K. The Reynolds number at the inviscid shock impingement location is Re x imp = 300,000, with x imp = 169.1 ⋆ imp the distance from the leading edge of the laminar boundary layer to the inviscid shock impingement location. The ⋆ imp represents the boundary-layer displacement thickness at the x imp position. The incident shock angle is 32.58 • , leading to a shock strength of p 3 ∕p 1 = 1.4. A laminar boundary layer is prescribed in the inflow boundary, so that we can use a shorter streamwise distance of 83.9 ⋆ imp from the leading edge of the computational domain to the x imp position. The computational domain has the size of 297 ⋆ imp × 193.6 ⋆ imp in the streamwise and wall-normal directions. A diagram for the computational setup is presented in Figure 8 for clarity, with the length in the unit of ⋆ imp . The prescribed inflow profile is given by the solution of the 2D steady-state compressible laminar boundary layer equations with zero pressure gradient. We follow the classical way to transform the physical space to a reference space to seek a self-similar solution. Detailed descriptions can be found in the textbook of Schlichting. 48 The oblique shock is introduced at the upper part of the inflow boundary by imposing a jump condition according to the Rankine-Hugoniot relation. The flat-plate wall is treated as an adiabatic no-slip wall. The top boundary and the outlet are handled as supersonic outlets by adopting the "pressure outflow boundary condition" defined by Carlson. 49 The initial solution is set as the undisturbed laminar boundary layer flow for better solution convergence. The classic Jameson indicator in Equations (12) and (13) is employed for the DG/FV cells switching.
To ensure adequate grid resolution, grid convergence studies have been performed both for the DG scheme and the second-order FV scheme with the same amount of DoFs. Furthermore, according to the recent study of Sansica et al., 50 the computation should be run long enough to make sure the flow reach a fully-converged steady state. In this work, we run all computations until 50 flow through times, see Figure 9, where we evaluate the development of the total wall-friction force on the flat plate over time to guarantee the convergence. The "FV 45500DoFs" case in Figure 9 (left) corresponds to the same mesh resolution of Katzer, 11 that is, the same height of the first-layer grid cell Δy is ensured. The other cases represent different levels of refinement by a factor of 2, 3, 4, 5, 6 in both directions respectively. For example, the numbers of mesh grids in streamwise and wall-normal directions of the "FV 182000DoFs" case are two times of the corresponding "FV 45500DoFs" case. From Figure 9 (left), we can tell that the grid resolution adopted by Katzer 11 has not reached the necessary convergence level yet. In Figure 9 (right), we use the same mesh resolution as in Figure 9 (left). The difference is that we employ a polynomial degree of N = 4 inside each DG element. In Figure 9 (right), it is obvious that high-order DG method is more likely to obtain an accurate result given the same solution DoFs as the second-order FV method. For example, the "DG∕FV 45500DoFs" case in Figure 9 (right) has the same DoFs as the "FV 45500DoFs" case in Figure 9 (left), but the result of the "DG∕FV 45500DoFs" case is more close to the converged solution.
According to Wang et al., 51 a fair way to illustrate the efficiency of high-order methods is to look at the computational cost to achieve the same level of accuracy as the low-order method. In Figure 9, we discover that the "FV 409500DoFs" case in Figure 9 (left) and the "DG∕FV 182000DoFs" case in Figure 9 (right) have began to be converged, and therefore in the following Section 4.2, we focus on the comparison of the accuracy and efficiency between these two cases.

2D oblique SWLBLI results
Without further statement, the figures in this section show results from the "DG∕FV 182000DoFs" case. Figure 10 presents an instantaneous snapshot of the Schlieren, which is defined in this article as . In Figure 10, we can see that the laminar boundary layer separates due to the strong adverse pressure imposed by the incident shock. The transmitted incident shock impinges on the vertex of the separation bubble, and then reflects as an expansion fan leading to the reattachment on the wall. Two compression waves form upstream the separation and downstream the reattachment due to the change of the flow shape on the flat plate. Figure 11 presents the distribution of the DG grid elements and the FV subcell elements while the flow has reached a steady state. In the near-wall region of the laminar boundary layer, no FV subcell element exists and only the high-order DGSEM is active. Figure 12 shows the contour of the instantaneous pressure distribution in the flow field, where we can see the rapid rise of the pressure due to the upstream compression waves, the decrease of the pressure inside the expansion fan, and the gradual rise of the pressure owing to the weak compression waves downstream. The flat plateau of the pressure inside the recirculation region can be noticed in Figure 12 as well.
In Figure 13, we compare the scaled pressure p∕p ∞ and the skin-friction coefficient C f on the wall between the "DG∕FV 182000DoFs" and the "FV 409500DoFs" cases. Since the results of Sansica et al. 50 is overall close to the results of Katzer, 11 and the main difference appears at the separation bubble length and the minimum value of C f , we only include the main different results of Sansica et al. 50 (the red stars) in the skin-friction plot of Figure 13. The black dashed line provides the result of the "FV 45500DoFs" case to compare with the Katzer's solution (the blue stars). From Figure 13, we discover that the results of the second-order FV scheme in "FV 45500DoFs" case have an overall good agreement with Katzer, 11 which is expected due to their similarity in the grid resolution. According to the skin-friction coefficient distribution in the lower part of Figure 13, the "DG∕FV 182000DoFs" and "FV 409500DoFs" cases predict a bigger separation length than Katzer, 11 which is consistent with the latest discovery of Sansica et al. 50 (the red stars), where the DNS simulation predicts approximately 8.8% larger separation area length than Katzer. 11 In the current computation, both the "DG∕FV 182000DoFs" and "FV 409500DoFs" case predict a 5.87% larger separation length than Katzer, 11 and around 8.25% larger than the "FV 45500DoFs" case. Considering the small difference in the numerical schemes and computational setup, our solutions have already gained a good accuracy.
In Figure 13, the close results of "FV 409500DoFs" and "DG∕FV 182000DoFs" allow us to present the efficiency of the high-order DG method. According to the previous numerical experiments on the performance behavior of our shock-capturing method implementation, 29 different problem sizes (DoFs per core) have a big effect on the performance. Therefore, both computations are performed on the same test environment and the same problem size of approximately 355DoFs∕core. The "FV 409500DoFs" case is performed on 1152 cores with 9 min, that is, a total CPU time of 10,368 min. The "DG∕FV 182000DoFs" case is performed with 512 cores with 15 min, that is, a total CPU time of 7680 min, around 35% less than the "FV 409500DoFs" case.

THREE-DIMENSIONAL OBLIQUE SHOCK-WAVE/TURBULENT BOUNDARY-LAYER INTERACTION (SWTBLI)
In this section, we apply our DG framework to the three-dimensional turbulent boundary layer interacting with an oblique shock wave. We start with the simulation of the wall bounded flow without the oblique shock wave and demonstrate that the turbulent flow has been recovered before interacting with the shock by comparing with reference DNS data.

Numerical set-up and boundary conditions
The simulation set-up is based on the LES computation of Morgan et al. 52 and the experiment of Threadgill and Bruce. 53 The inflow Mach number is chosen to be Ma = 2.0. The free-stream temperature and pressure are T ∞ = 288.15 K and p ∞ = 101,306.34 Pa, respectively. The inflow Reynolds number is Re 0 = 1628, with 0 denoting the momentum boundary-layer thickness at the inlet. The flow deflection angle is = 8 • , which leads to an incident shock angle of = 37.2 • and the shock strength of p 3 ∕p 1 = 2.30. In order to be consistent with Morgan et al., 52 we specify a reference position x r located at 3 r upstream the inviscid shock impingement point x imp = 20 r , with r represents the boundary-layer thickness at this reference position. A schematic representation of the computational setup is as Figure 14, with the length in the unit of r for simplicity. The Reynolds number at x r is Re r = 2188, which is equivalent to the reference LES. However, due to the slightly bigger Mach number of 2.05 and the different boundary layer properties in the reference LES, we don't expect a fully quantitative agreement with the results of Morgan et al., 52 but the same trends of the wall pressure and skin-friction coefficient curves. The experiment of Threadgill and Bruce 53 has higher Reynolds numbers and the 3D sidewall effects of the wind tunnels, but the same Mach number and the same shock strength are utilized so that it is used for a supplementary validation by comparing with the dimensionless Strouhal number of the flow unsteadiness. Table 1 summarizes the characteristic properties of our flow conditions. The inflow boundary-layer thickness 0 and inflow displacement boundary-layer thickness * 0 are also showed in Table 1 for later comparison. Throughout this study we will use the boundary-layer parameters either from the inlet or the reference location.
The computational domain covers a rectangular box with dimensions of (0, 35 r ) × (0, 20 r ) × (0, 3 r ) and is discretized with 129 × 81 × 26 grid cells in the streamwise, wall-normal and spanwise directions, respectively. The mesh is designed to be composed of three parts in the streamwise directions, see Figure 15: the "recovery mesh" covering the streamwise distance of 0 ∼ 15 r used for the turbulence recovery, the "interaction mesh" covering the streamwise distance of 15 r ∼  25 r used for the study of the oblique SWTBLI, and a loosely stretched (from a grid resolution in wall units of 30 to 322 as showed in Table 2) "sponge mesh" covering the streamwise distance of 25 r ∼ 35 r used for a sponge zone to progressively damp the convected vortices to avoid reflections from the outlet. All three mesh components have equally distributed grid cells in the spanwise direction. Additionally, the same grid stretching law is applied in the wall-normal direction following where H j is the height of the jth element, l o is the height of the first element, and S fac is the stretching factor. From the wall to 1.2 r , a stretching factor of S fac = 1.05 is employed, and 38 cells are contained in this layer. The cell height is increased from Δy + min = 1.6 to Δy + = 9.8 in wall units. From 1.2 r to 5.0 r , a stretching factor of S fac = 1.15 is employed, leading to the maximum cell height of around Δy + = 79 in this layer. From 5.0 r to the top boundary, equidistant cell height is adopted. This discretization in combination with a polynomial degree of N = 6 leads to a grid resolution of Δy + min = 1.6, Δx + = 30 and Δz + = 16 at the reference position x r . The computational domain has a total number of around 93 million DoFs. The dimensions of the three parts of the mesh are summarized in Table 2. All the Δx + , Δy + and Δz + are computed based on the flow parameters at the reference position x r .
The synthetic eddy method (SEM) was firstly proposed by Jarrin et al. 54 to introduce turbulence into computational domains. The idea behind the SEM is that turbulence can be deemed as a superposition of coherent structures, and therefore one can generate randomly-distributed coherent structures with specific time and length scales as well as geometrical shapes at the inlet to model the turbulent fluctuations. Due to the random superposition of the vortices, SEMs do not generate spurious correlations of the inflow data, which fits well into the unsteadiness analysis of oblique SWTBLIs. Roidl et al. 22 improved the original SEM by defining three different wall-normal virtual volume boxes to take the different sizes of the vortices normal to the wall into account. The improved SEM also takes account of the compressibility effects of high-speed flows by enforcing the strong Reynolds analogy. We employ the improved SEM by Roidl et al. 22 at the lower part of the inflow boundary to generate turbulence. The anisotropic linear forcing (ALF) introduced by de Laage de Meux et al. 23 is limited to the forcing zone extending from 5.0 r to 10.0 r , indicated by the red region shown in Figure 15. The shock is introduced at the upper part of the inflow boundary, where the free-stream variables are prescribed upstream of the shock, and the variables downstream of the shock are calculated by the Rankine-Hugoniot relations. The flat plate is modeled as an adiabatic no-slip wall. Periodic boundary conditions are used in the spanwise directions. At the outlet boundary, a sponge layer is used and a supersonic exit condition is applied to all variables. The top boundary is also set as a supersonic exit condition. The reflected shocks and compression waves always exit the computational domain at the top boundary or the outlet boundary to eliminate the possible negative influence on the results, which is the same as the laminar case, see Figure 14.

The turbulent boundary layer (TBL) without shock
Before analyzing the oblique SWTBLI, we validate the accuracy of the reformulated SEM by comparing the LES results of a turbulent boundary layer with the DNS data provided by Wenzel et al. 55 All the boundary conditions of the undisturbed TBL stay the same as in the oblique 3D SWTBLI case, except that no shock is introduced at the upper part of the inflow boundary, and the top boundary is set as a Dirichlet boundary condition with free-stream properties. The mesh for this undisturbed TBL is the same as the mesh defined in Table 2. The simulation is performed with 16,384 CPU cores, which leads to a load of about 5687 DoFs per core. Due to an efficient parallel implementation of the reformulated SEM to take advantage of the high-performance computer architecture, we get a performance index (PID) of about 1.55 s for the present computation. The PID is defined as describing the time that it takes for one core to update one DoF for one time stage of the Runge-Kutta method. Comparing with the standard parallel performance of our code implementation conducted by Krais et al., 29 the PID of 1.55 s in current simulation is about 94.0% larger than the baseline data. There are several main reasons for this performance loss. First, although we have distributed the computation of the SEM in the inflow boundary into several CPUs, load imbalance problem still exists in the whole computational domain. Second, according to our previous experience, the use of the ALF by de Laage de Meux et al. 23 accounts for about 5 ∼ 10% performance loss. Third, the baseline performance of Krais et al. 29 is obtained using the free stream flow as the test case, while the current simulation involves the communication overhead due to the collection of the turbulent statistics for analysis. The flow reached a quasi-steady state after around 70 convective time units t * = r ∕U ∞ . Figure 15 shows the instantaneous isosurfaces of the Q-criterion colored by the velocity magnitude at the quasi-steady state. The red box represents the forcing zone where the ALF is activated. From here on, the turbulence statistics are collected for further 265t * to compare with the DNS data. In total, the current undisturbed TBL simulation costs around 46,173 CPU hours. Figure 16 presents the comparison of the temporally and spatially averaged local statistics at 11 r location, which is 1 r after the forcing zone. In Figure 16 (left), we can see that the van Driest transformed mean velocity profile of the LES is in good agreement with the DNS data. In Figure 16 (right), the density-scaled Reynolds stresses of the LES and the DNS match each other well. Given the large discrepancy in grid resolutions between the LES and the DNS, the agreement in results has already proven the capability of the adopted method. Furthermore, we present the developments of the skin-friction coefficient and the shape factor in Figure 17 to state that the TBL has reached an equilibrium state before interacting with the shock, which is the premise of an accurate oblique SWTBLI result in the next section. The van Driest II transformation is utilized to compare the skin-friction coefficient in the present compressible flow with the analytical correlation C fi = 0.024Re −0.25 i of incompressible flows, as Wenzel et al. 55 and Pirozzoli and Bernardini 56 did. Figure 17 (left) presents the development of the transformed skin-friction coefficient C f i over the transformed momentum thickness Reynolds number Re i . It is evident that the skin-friction coefficient of the present LES quickly approaches to the reference DNS and analytical results after the forcing zone (the red dashed line). Figure 17 (right) provides the comparison of the shape factor as a supplementary assessment of the reliance of the current LES computation.

Instantaneous and mean-flow organization
The oblique SWTBLI simulation is performed with 16,384 CPU cores, which leads to the same load of about 5687 DoFs per Core as the undisturbed TBL above. A PID of around 1.84 s is achieved in the simulation, which corresponds to an approximately 18.7% computational cost increase compared to the performance calculating the previous shock-less boundary layer. According to the numerical experiments of Krais et al., 29 a performance degradation of 15% − 20% due to the hybrid DG/FV scheme is normal. It proves that the parallel implementation of the filtering process for the Jameson shock indicator is well-designed.
The oblique SWTBLI calculation reaches a quasi-steady state after 270t * . The instantaneous separated flow and the numerical Schlieren of the studied oblique SWTBLI is shown in Figure 18. The flow separation (the blue areas in Figure 18) emerges due to the strong adverse pressure gradient imposed by the incident shock. The reattachment of the separated flow on the wall owes to the energy transferred from the outer high-speed flows. 6 The appearance of the separation deforms the flow streamlines near the flat plate, and accounts for the upstream compression waves, which coalesce into the reflected shock intersecting with the incident shock. The transmitted incident shock encounters the separation regions and then reflects as curved compression waves and expansion fan. The reflected curved compression waves finally merge into the transmitted reflected shock. The expansion fan deflects the detached flow toward the wall, and eventually accelerates the reattachment. The downstream compression waves after the reattachment region promote the equilibrium state of the flow. The appearance of the reflected curved compression waves in Figure 18 is not described by Delery and Dussauge. 6 In order to account for the influence of the reflected curved compression waves, Agostini et al. 57 establishes a new inviscid equivalent scheme for oblique SWTBLIs.
A 3D snapshot of the oblique SWTBLI is shown in Figure 19. The shock systems are identified by the dilatation in gray blocks. Turbulent structures are depicted by the isosurface of the Q-criterion colored by the streamwise velocity magnitude. In Figure 19, we can see the apparent growth of the TBL thickness due to the shock systems. Additionally, the phenomenon that the shocks penetrate into the TBL can also be clearly detected. From the instantaneous density field in Figure 20, the formation of the large vortical structures after the shock systems is obvious. Figure 21 shows the contour of the instantaneous streamwise velocity field, where we can also discover the TBL thickness increase and the flow separation.
The following mean-flow data are obtained by temporally and spatially averaging the flow contours after the quasi-steady state for further 270t * . Figure 22 shows the comparison of the wall pressure p∕p ∞ and the skin-friction coefficient C f between the current LES and the "fine mesh" results of the LES in Morgan et al. 52 The dimensionless streamwise length (x − x imp )∕ r is adopted in the abscissa axis. Both LES predict a similar trend of the pressure distribution along the flat plate in Figure 22 (left). The higher pressure amplitude of the reference LES adheres to the fact that it has a bigger shock strength p 3 ∕p 1 . From the pressure distribution of the present LES in Figure 22 (left), we can discover an early flow separation in the streamwise direction. The early separation of the current LES is more clear in Figure 22 (right), where we find that the present LES predicts a separation at x s = 16.7 r and that the reference LES predicts a separation at 17.6 r . The larger separation length of our LES computation than the reference LES is consistent with the experimental conclusion of Souverein et al., 24 who proves that, given the same incident angle and approximate Reynolds numbers, a smaller Mach number leads to a bigger separation length. In Figure 22 (right), we can see that both LES acquire similar reattachment locations of x r = 19.3 r in reference LES and x r = 19.5 r in our computation. The interesting thing is that, in the study of Morgan et al., 52 different mesh resolutions also lead to nearly the same the prediction of the reattachment position. This phenomenon might result in the conclusion that the prediction of the separation location of oblique SWTBLIs is highly sensitive to the computational resolution and flow conditions, whereas the reattachment location is not. The overall good agreement of the developments of the skin-friction coefficient of both LES shows that our present LES can compete the resolution of the reference LES, even though our present LES calculation has a much coarser grid resolution. Additionally, we can apparently see the skin-friction coefficient "plateau" of the present LES in Figure 22 (right). According to the study of Morgan et al., 52 a skin-friction coefficient "plateau" can be observed only if the computational resolution is high enough. This further demonstrates the high accuracy of our DG method. Except for the different Mach number and shock strength, another main reason for the discrepancy in the values of the skin-friction coefficient of the two LES lies in the difference in the TBL properties, which is the most obvious upstream and downstream the interaction area in Figure 22 (right).
We will further compare our computation with several DNS and LES computations according to the scaling law in the Section 5.3.3, for which the mean interaction length L int is required. According to the definition, the mean interaction length corresponds to the distance between the mean position of the reflected shock x 0 and the inviscid shock impingement location x imp . All the values are evaluated from the time-and spanwise-averaged pressure plot in Figure 23, where x 0 = 16.2 r , x imp = 20 r and thus L int = 3.8 r .

Spectral analysis
The unsteadiness feature of the oblique SWTBLI is studied in this section by means of spectral analysis. Inside the interaction zone, 100 equally spaced wall pressure probes have been placed along the streamwise center line. The probes are We calculate the weighted power spectral densities (PSD) of the signals by normalizing each premultiplied PSD contour by the integral of the PSD over the whole frequency range, that is, f ⋅ PSD(f )∕ ∫ PSD(f )df , which can highlight the frequencies that contribute most in the flow field according to Pasquariello et al. 2 The weighted PSD of all probes is reported in Figure 24, where the mean reflected shock foot x 0 (the red line), the mean separation position x s (the blue line), and the mean reattachment location x r (the green line) are indicated by horizontal dashed lines. As the authors 15 said, numerical computations can capture the energetically significant high-frequency turbulence oscillations on the wall, and consequently a different high-frequency end with experiments would appear in the spectrum. For example, the reference experiment of Threadgill and Bruce 53 reports that the upstream boundary layer is dominated by a high-frequency energy content of around 10 kHz. However, in the current computation we can discovery higher-frequency contents, for example, in Figure 24, the frequency higher than S t = 1 in upstream boundary layer is more obvious (corresponding to a frequency of around 353 kHz). Hence, it is appropriate to interpret the Figure 24 as the map of the most dominant wall-pressure fluctuations as one moves along the streamwise direction. We can see that: (1) no extra low-frequency unsteadiness is introduced into the incoming flow due to the adoption of the reformulated synthetic eddy method, (2) a low-frequency unsteadiness of around S t = 0.03 is detected near the mean reflected shock foot x 0 , (3) the flow goes back to higher frequencies downstream the interaction area. The phenomenon is consistent with the discovery in the reference experiment of Threadgill and Bruce. 53 For a better understanding of the low-frequency phenomenon, in Figure 25, we provide the raw wall pressure fluctuations p ′ w (the black line) at four streamwise positions corresponding to the vertical coordinates (x − x imp )∕ r of Figure 24, that is, at upstream flow of (x − x imp )∕ r = −5.0, mean reflect shock foot at (x − x imp )∕ r = −3.8, inside the separation flow at (x − x imp )∕ r = −2.0, and at downstream flow of (x − x imp )∕ r = 5.0. Additionally, a low-pass filter is applied to the raw pressure signals to highlight the low-frequency content, which is shown in red line in Figure 25. It is clear that in the shock foot x 0 , a low-frequency signal appears. The current oblique SWTBLI requires around 196,608 CPU hours altogether.

Assessment of the scaling law
As mentioned in Section 1, Souverein et al. 24 summarized many experimental data for oblique SWTBLIs and SWTBLIs over compression ramps into a scaling law to quantify the relation between the interaction length and the separation criterion. The scaling law is based on the experimental data with lowest Reynolds number of Re 0 = 5000, which is almost unreachable for DNS and LES because of the huge computational costs. On account of the important role of an accurate scaling law in understanding scientific principles and validating numerical results, it is necessary to evaluate the effectiveness of the scaling law by Souverein et al. 24 with extensive computational data. Preliminary work of Morgan et al. 58 shows that the Reynolds number has an effect on the interaction length L int . Recent works of Bermejo-Moreno et al. 25 and Wang et al. 26 reveal that the 3D sidewall effects of wind tunnels lead to a much bigger interaction length L int in comparison with the spanwise-periodic computations, which might totally destroy the applicability of the scaling law into numerical simulations. In this section, we analyze the applicability of the scaling law with our LES result together with some DNS and LES data from literature. Souverein et al. 24 proposes a power law between the scaled interaction length L * and the separation criterion S * e , that is, L * = 1.3 ⋅ (S * e ) 3.0 . They determine the scaled interaction length L * based on the mass conservation in a control volume covered by the incoming and the outgoing boundaries where L int , * 0 , , are the interaction length, the inflow displacement boundary-layer thickness, the incident shock angle, and the deflection angle as defined before. The separation criterion S * e is introduced as TA B L E 3 The DNS and LES data considered in Figure 26 Method with Δp s = p 3 − p 1 the shock intensity. p 3 is the pressure after the shock reflection system. q 0 = 1 2 ∞ U 2 ∞ is the dynamic pressure in the free stream. k is a factor to take the influence of the Reynolds number into account. Based on Souverein et al., 24 k = 3 is recommended for Re 0 > 10 4 and k = 2.5 for Re 0 ≤ 10 4 . In current LES computation, the scaled interaction length is L * = 2.92 and the separation criterion is Se * = 1.39. The computed L * is underestimated by 16.4% in comparison with the scaling law. The corresponding values in the experiment of Threadgill and Bruce 53 are L * = 3.22 and Se * = 1.39. The L * in the experiment is close to the scaling law, with a small underestimation of 7.8%. We believe that the underestimation of L * in current computation can be contributed to the 3D side-wall effects in the experiment of Threadgill and Bruce. 53 Several computations have already proven the underestimation of the related L int . 25,26,58,59 In order to provide an overall evaluation on the discrepancy, we look into more results of recent DNS and LES computations.
In Table 3, the subscript r refers to the defined upstream reference locations in the respective papers. In Pirozzol et al., 45 Pirozzoli and Bernardini, 44 and Nichols et al., 60 the same Mach number, approximate Reynolds numbers, and similar inflow conditions are adopted, and thus we suppose they have the same ratios of * r ∕ r = 0.276 and r ∕ in = 0.27, with in the inflow boundary-layer thickness. Wang et al. 26 does not provide the necessary mean reflected shock foot position for the spanwise-periodic case, we approximate it with the mean separation position for our comparison. In the DNS simulation of Priebe et al., 61 we obtain the mean reflected shock foot location 6.7 in and the shock impingement location 13.3 in from the provided time-and spawise-averaged numerical Schlieren. Hadjadj et al. 62 adopts a similar set-up with the experiment of Piponniau et al. 63 so that the inflow displacement boundary-layer thickness * 0 = 3.4 mm in the experiment is assumed to be the same as in the computation. From the collected data in Table 3, we can draw two preliminary conclusions. First, according to the results of Pirozzol et al. 45 and Nichols et al. 60 with deflection angles of = 8 • , 9.5 • , we find that smaller spanwise length scale L z might lead to a bigger L * . Second, comparing the DNS data in Pirozzoli and Bernardini 44 with those results with the same deflection angle = 8 • in Nichols et al. 60 and Morgan et al., 58 we can conclude that bigger Reynolds numbers can lead to bigger L * . However, the overall influences from different spanwise length scales and Reynolds numbers are too small to be counted in the evaluation of the scaling law. We present the comparison between all the computational data in Table 3 and the scaling law in Figure 26 (left), where we can see that the DNS and LES data fails to match the scaling law, especially in the range of Se * ∈ (1.3, 1.8). Additionally, the position of our result in Figure 26 (left) is consistent with the distribution of other DNS and LES data, which indirectly proves that our result is in a reasonable range. In view of the negligible influences of the Reynolds number and the spanwise length scales, it is justified to attribute the underestimation of L * to the difference in the numerical and experimental set-up. That is, the missing of the modeling of the 3D sidewall effects result in the smaller predicted L * (and thus smaller L int ), which is consistent with the conclusion of Bermejo-Moreno et al. 25 and Wang et al. 26 Hence, we can conclude that modeling sidewall effects is compulsory if comparable interaction dimensions and Strouhal numbers with experiments are expected.
However, modeling the 3D sidewall effects will largely increase the complexity of simulations and requires much more computational costs. Following the way Souverein et al. 24 did, we use the least square technique to fit a new power law based on the data in Table 3. The new power law is L * = 0.95 ⋅ (S * e ) 3.15 , see Figure 26 (right). The scaling law of Souverein et al. 24 is concluded from the experimental data with the separation criterion Se * < 1.8, so that further experimental studies on the extension of the scaling law are necessary. Our fitted power law for spanwise-periodic numerical computations is based on a small amount of data, and it mainly aims for an overall impression on the discrepancy between the experimental scaling law. In order to generate a reliable scaling law for spanwise-periodic numerical calculations, a subsequent extensive numerical verification is still required. Moreover, Equation (20) shows that different inflow positions have a remarkable influence on the calculation of the L * , which is especially significant in numerical simulations. However, Souverein et al. 24 does not point out how to choose the inflow position for the calculation of the * 0 in Equation (20). From our experience, an inflow position upstream the interaction area within a distance of 10 r ∼ 20 r should produce a good agreement with the new scaling law. We emphasize that although the new scaling law is proposed based on the numerical data of oblique SWTBLIs, it should be applicable to the SWTBLIs over compressible ramps owing to their high similarity.

SUMMARY AND DISCUSSION
We investigated the accuracy and efficiency of our high-order discontinuous Galerkin method in the canonical oblique shock-wave/boundary-layer interactions in laminar and turbulent flows.
In order to accurately capture the shocks inside and outside the turbulent boundary layer, we proposed to filter the instantaneous flow field before evaluating the shock indicator. Specifically, in this article, we combined a global spanwise filter with the Jameson shock indicator. In comparison with the popular Ducros indicator, the modified Jameson shock indicator allowed us to easily but accurately capture shocks inside the interaction area. Due to the unsteady shock motion inside the turbulent boundary layer, a positivity preserving limiter was implemented to avoid the numerical instabilities due to frequent switching between the DG and FV subcell elements. We utilized a reformulated synthetic eddy method at the inlet to generate the unsteady turbulent inflow boundary condition. The anisotropic linear forcing method of de Laage de Meux et al. 23 has been employed to further reduce the turbulence recovery length. The accuracy and efficiency of current LES have been exhibited by comparing with the reference DNS. The usage of the reformulated synthetic eddy method in combination with a long integration time guaranteed an accurate analysis of the separation-induced low-frequency unsteadiness. Through comparing the wall pressure p∕p ∞ and the skin-friction coefficient C f with the LES of Morgan et al., 52 we have shown that our high-order DG framework is comparable to the accuracy of the established LES scheme. We also verified our computation with the reference experiment of Threadgill and Bruce 53 on the unsteady behavior.
Based on the accurate prediction of our DG method for the interaction dimensions in the oblique SWTBLI, we assessed the scaling law proposed by Souverein et al. 24 Some DNS and LES data from recent papers utilizing various computational schemes in different flow conditions are also contained for validation. It is shown that although different spanwise length scales L z and Reynolds numbers Re will affect the computed scaled interaction length L * , the missing of the 3D sidewall effects in experiments is the main reason for the disagreement of the numerical results with the scaling law. We suggested a new scaling law for the numerical simulations with spanwise periodic boundary conditions.
In our computation, we have also demonstrated the reliability of our implementation in dealing with a massively parallel computation of approximately 93 million DoFs on modern high performance computers using up to 16,384 cores. With the efficient schemes and parallel architectures, we show the high efficiency of our DG methods to resolve the complex flow of the oblique SWTBLI on the order of hours.
Future work include the application of our DG method to more complicated three-dimensional cases such as the SWTBLIs of a single fin or double fins mounted on a flat plate. The flow speeds will also be increased to hypersonic Mach numbers to include the effects of heat transfer on the flow field. The study of the transition from regular reflection into Mach reflection due to the increase of the incident shock angles is still in an early stage, and thus more efforts are necessary in order to reveal more physics behind the transition. Deeper research in the controlling of the SWTBLIs is still required. Considering the maturity of the high-fidelity LES and the large database on SWTBLIs, the main emphasis of the efficiency improvement of SWTBLIs simulations should be committed to the application of some advanced schemes like the wall-modeled LES. The significant changes of the pressure loading and the appearance of the low-frequency unsteadiness in SWBLIs might damage the elastic structures on the surface of high-speed vehicles, or at least limit the performance of the structures. In the future, we will develop a high-order structural solver to study the fluid-structure interactions of the SWBLIs on flexible structures.

ACKNOWLEDGMENTS
The first author would like to acknowledge the support of the China Scholarship Council (CSC) for supporting his stay in Stuttgart-201906020172. The work was supported by the German Research Foundation (DFG) through Germany's Excellence Strategy-EXC 2075-390740016. Computer time is provided by the High Performance Computing Center Stuttgart (HLRS) on the HPE APOLLO (HAWK) system under the grant hpcdg. FLEXI has been developed, tested and applied by a number of people, both from the Numerics Research Group at the University of Stuttgart and elsewhere. We acknowledge their contributions and work. A possibly incomplete list of contributors can be found in the FLEXI Doxygen: https://www.flexi-project.org/doc/doxygen/html/. Open access funding enabled and organized by Projekt DEAL.