On the performance of WENO/TENO schemes to resolve turbulence in DNS/LES of high‐speed compressible flows

High‐speed compressible turbulent flows typically contain discontinuities and have been widely modeled using Weighted Essentially Non‐Oscillatory (WENO) schemes due to their high‐order accuracy and sharp shock capturing capability. However, such schemes may damp the small scales of turbulence and result in inaccurate solutions in the context of turbulence‐resolving simulations. In this connection, the recently developed Targeted Essentially Non‐Oscillatory (TENO) schemes, including adaptive variants, may offer significant improvements. The present study aims to quantify the potential of these new schemes for a fully turbulent supersonic flow. Specifically, DNS of a compressible turbulent channel flow with M = 1.5 and Reτ = 222 is conducted using OpenSBLI, a high‐order finite difference computational fluid dynamics framework. This flow configuration is chosen to decouple the effect of flow discontinuities and turbulence and focus on the capability of the aforementioned high‐order schemes to resolve turbulent structures. The effect of the spatial resolution in different directions and coarse grid implicit LES are also evaluated against the WALE LES model. The TENO schemes are found to exhibit significant performance improvements over the WENO schemes in terms of the accuracy of the statistics and the resolution of the three‐dimensional vortical structures. The sixth‐order adaptive TENO scheme is found to produce comparable results to those obtained with nondissipative fourth‐ and sixth‐order central schemes and reference data obtained with spectral methods. Although the most computationally expensive scheme, it is shown that this adaptive scheme can produce satisfactory results if used as an implicit LES model.


INTRODUCTION
High-speed compressible turbulent flows play a significant role in a variety of advanced fluid mechanics related applications and disciplines such as supersonic or hypersonic aerodynamics and combustion, aerospace propulsion, inertial confinement fusion, and astrophysics. [1][2][3][4][5] Those flows typically contain shockwaves and flow discontinuities and have been widely modeled using Weighted Essentially Non-Oscillatory (WENO) schemes 6,7 due to their high-order accuracy and sharp shock/discontinuity capturing capability. However, those schemes, which aim to stabilize the solution by adding numerical dissipation, may damp some scales of turbulence and as a consequence result in inaccurate solutions in the context of turbulence-resolving simulations, unless excessively fine spatial resolutions are used, which is not practical for direct numerical simulation (DNS) of complex flows. 8 Therefore, it is of great importance to critically evaluate and quantify the ability of such high-order schemes to accurately predict the wide range of flow structures in compressible turbulent flows.
There are not many studies available in the literature discussing and comparing the performance of modern high-order shock capturing schemes to resolve compressible turbulence, particularly in wall-bounded flows. Mossi and Sagaut 9 were among the first to investigate the performance of high-order shock capturing schemes for compressible turbulent flows. Specifically, they performed direct and large eddy simulations (LES) on quasi-incompressible (M = 0.5) and compressible (M = 1.5) Turbulent Channel Flows (TCFs) using second-order central and third-order Total Variation Diminishing (TVD) Roe schemes. They reported that the TVD scheme performed poorly in capturing turbulent structures mainly due to its significant numerical diffusion. Gerolymos et al 10 used high-order upwind, WENO, and mapped WENO 11 schemes to conduct DNS of subsonic (M = 0.35) and supersonic (M = 1.5) turbulent channel flows. It was noticed that very high-order upwind schemes could reproduce results obtained by pseudospectral schemes, given slightly higher spatial resolutions. However, noticeable numerical diffusion was reported for the WENO schemes. Gerolymos et al 10 also provided a brief review of some early DNS studies on TCFs using high-order methods.
The performance of the WENO scheme can be improved by using a dynamic balance between an upwind formulation for high-gradient and shock-containing regions and a central formulation in smooth regions such as the method proposed in a bandwidth-optimized WENO scheme, namely, WENO-SOL. 12 Taieb and Ribert 13 used the WENO-SOL scheme to perform DNS and LES of supersonic channel flows (M = 1.5). The fifth-order scheme performed well in the context of DNS. It was also used without and with Sub-Grid Scale (SGS) models to perform LES of TCFs. Without an SGS model, the WENO-SOL scheme was used as an Implicit LES (ILES) model and exhibited a satisfactory performance. The performance of the WENO-SOL scheme was improved when it was combined with a dynamic SGS LES model and produced a very close behavior to DNS. The worst performance of the WENO-SOL scheme was for when it was combined with a nondynamic non-dynamic SGS LES model mainly due to issues related to overdissipation. 13 Ribert et al 14 extended the WENO-SOL scheme and also the WENO scheme developed by Jiang and Shu (WENO-JS) 6 for real fluids in order to perform LES of turbulent channel flows under supercitical conditions. Under such conditions, where tiny variations in main flow parameters could significantly alter the fluid properties and impose considerable gradients, applying a shock capturing scheme could greatly improve the stability of the solution. A comparison of some fifth-and sixth-order shock capturing schemes of the WENO family to model low Mach (M = 0.01, 0.1 and 0.3) TCFs has also been reported by Matsuyama. 15 Various schemes including WENO-JS, WENO-Z, 16 Linear WENO (LWENO), 17 and WENO with a relative limiter (WENO-RL) 18 were investigated. 15 Again, the dissipative behavior of the WENO scheme was highlighted and the use of its modified versions, such as the WENO-RL scheme, which could noticeably reduce the dissipation, was emphasized.
In addition to channel flows, the performance of various high-order shock capturing schemes, mainly from the WENO family, to resolve compressible turbulent flows have been investigated in the context of turbulent boundary layers, isotropic decaying turbulence, and Shock-Boundary Layer Interactions (SBLIs). Similarly to the aforementioned TCF cases, different modifications were also applied to overcome the destructive dissipative behavior of the WENO scheme. Lagha et al 19 successfully performed DNS of turbulent boundary layers with Mach number values up to M = 20 using a fifth-order WENO scheme. Specifically, the authors used the WENO scheme at high-gradient regions and a fifth-order upwind scheme for smooth regions with global Lax-Friedrichs flux splitting. Supersonic turbulent boundary layers were also studied by Duan et al 20 using a seventh-order WENO scheme with limiters to reduce the numerical dissipation. Johnsen et al 21 studied the performance of various high-resolution shock capturing schemes on a variety of problems including the Taylor-Green Vortex (TGV) and shock-dominated problems such as the Shu-Osher problem, a shock-vorticity interaction, and the Noh problem all on underresolved grids. The authors applied and compared several methods including WENO, hybrid WENO/central, artificial diffusivity, an adaptive characteristic-based filter, and a shock fitting method. It was concluded that the WENO and the artificial diffusivity method (which used the magnitude of the strain-rate to activate the artificial viscosity) were not suitable for high-fidelity modeling of compressible turbulence. However, modified artificial diffusivity methods that use dilatation to activate the artificial viscosity were found to be able to resolve compressible turbulence. The hybrid WENO/central method was effective in reducing the numerical dissipation but quite challenging with respect to developing effective shock sensors. 21 Brehm et al 22 also compared the performance of various high-order shock capturing schemes on a variety of problems including the Shu-Osher problem, the isentropic vortex convection problem, double Mach reflection, and a TGV problem. The authors used central finite difference schemes with explicit artificial dissipation, a compact central finite difference scheme with localized artificial diffusivity, and various modified WENO schemes (a sixth-order WENO by Hu et al, 23 a fifth-order WENO by Ghosh and Baeder, 24 and WENO-Z) in both explicit and compact finite difference forms and on different grid resolutions. Based on their decaying isotropic turbulence case, the authors reported that the fifth-order compact WENO and the sixth-order WENO resolved the small-scale flow structures better compared to a fifth-order WENO-Z. It was also noted that the highest spectral resolution was obtained by the localized scheme followed by the fifth-order compact WENO and the sixth-order WENO. 22 Karami et al 25 used a hybrid fifth-order WENO/sixth-order central scheme to perform LES of various flow configurations including a compressible round jet and an underexpanded impinging jet. The authors proposed a new shock sensor based on a WENO smoothness indicator. A comprehensive study on the performance of various shock sensors for hybrid WENO-based schemes was recently conducted by Zhao et al. 26 The authors also proposed a new shock sensor for WENO-based hybrid schemes, which was not affected by the numerical dissipation caused by the WENO scheme. Hoffmann et al 27 used a hybrid Summation-By-Parts (SBP) central scheme and a fifth-order WENO scheme to perform LES of transitional flow over a flat plate with an impinging shock and also turbulent flow over a compression corner. The blending approach of the abovementioned upwind/central-WENO hybrid schemes requires a great care to ensure numerical conservation and stability at the interface of different schemes used. Non-linear filter schemes [28][29][30][31] are developed to avoid this interface deficiency of high-order hybrid schemes. Specifically, a hybrid scheme is applied to the inviscid fluxes and the nonlinear fluxes are typically evaluated by central discretization of an entropy conservative formulation.
In a prework to the present contribution, Lusher and Sandham 8 conducted a study of the performance of various high-order shock capturing schemes to model a compressible TGV and a transitional supersonic SBLI. The authors used fourth-order central, fifth-and seventh-order WENO-Z, fifth-order WENO-JS, and also various orders and formulations of the Targeted Essentially Non-Oscillatory (TENO) family scheme. 32-35 A significant improvement was obtained with the TENO schemes over the WENO schemes.
It is evident that there are limited works available in the literature that provide direct comparisons of the performance of various high-order shock capturing schemes to accurately resolve turbulence in compressible high-speed wall-bounded flows. Also, the majority of the earlier works are based on schemes from the WENO family, which have shown significant dissipative behaviors and failed to accurately resolve turbulent structures. The present work aims to compare and quantify the performance of the existing high-order shock capturing schemes, including the TENO family schemes, to model fully turbulent supersonic flows. The focus is on the capability of those schemes to resolve turbulent structures and predict statistical quantities. Therefore, a channel flow configuration is used to decouple the effect of flow discontinuities and turbulence. Specifically, it extends the earlier work of Lusher and Sandham 8 to supersonic TCFs. In addition, the present study evaluates the performance of such high-order schemes in the context of ILES by making direct comparisons against DNS and explicit LES approaches.

Governing equations
The dimensionless governing equations of a compressible Newtonian fluid flow that conserve mass, momentum, and energy are given as: where represents the density, u i (i = 0, 1, 2) denotes the velocity component in the ith direction, E is the total energy, and p and denote the pressure and the Kronecker delta, respectively. The equations are solved in three spatial dimensions x i (i = 0, 1, 2). c i is a constant forcing term, which drives the flow, with a value of unity in the x 0 direction and zero in other directions. Moreover, ij and q j are the viscous stress tensor and the heat flux (based on Fourier's law) respectively, defined as: where denotes the dynamic viscosity, T is the temperature, and Pr is the Prandtl number. Also, Re * and M * denote the Reynolds and Mach numbers based on a friction velocity defined as u * = √ ⟨ wall ⟩∕⟨ ⟩ , where wall denotes the wall shear stress and ⟨ ⟩ is the bulk-averaged density. Angle brackets ⟨⟩ denote averages over the homogeneous spatial directions and time throughout this article. The dynamic viscosity, if required, is calculated as = T 0.7 . Also, since the governing equations are solved for the conserved variables ( u i and E), then the velocity field is obtained by u i = u i and the tem- Here the pressure of an ideal Newtonian fluid is obtained using an equation of state as: All quantities are nondimensionalized by the bulk-averaged density, the friction velocity (u * ), and the wall temperature. It is worth mentioning that the governing equations with the WENO and the TENO schemes are solved as presented in Equations (1) to (3) and when central schemes are used a skew-symmetric formulation is applied to the convective terms. 36

Flux reconstruction and shock capturing schemes
In the present work the governing equations are discretized by applying finite differencing within OpenSBLI, a high-order Computational Fluid Dynamics (CFD) framework. 37 Comprehensive discussions on the WENO and TENO schemes implemented in OpenSBLI can be found in References [8,38]. However, a brief discussion on the flux reconstruction and the aforementioned shock capturing methods in OpenSBLI is provided here for completeness. To discretize the governing equations, their fluxes need to be reconstructed at each discrete grid point. With the finite difference approach, the flux reconstruction can be performed in each direction independently. For example, a 1D hyperbolic equation of the form U t + f (U) x = 0 can be discretized by reconstructing the flux term f (U) x over mid-point (half-node) locations with grid spacing Δx as follows: where i − 1 2 and i + 1 2 denote mid-point locations between grid point i and grid points i − 1 and i + 1 in the x direction, respectively. The full numerical stencil is reconstructed from a convex combination of lower order stencils as shown for WENO in Reference [38] and TENO in Reference [8]. Midpoint fluxes in Equation (7) are calculated as a weighted sum of Essentially Non-Oscillatory (ENO) interpolations over a set of r candidate stencils as: wheref (r) i± 1 2 are evaluated using the classical ENO scheme of Shu 7 and r are the nonlinear weights. The order of the WENO scheme is equal to 2K − 1 and for instance K = 3 for a fifth-order scheme, which results in r = 0, 1, 2 candidate stencils. 38 One difference between WENO and TENO schemes stems from the way that the non-linear weights of Equation (8) are calculated. In WENO schemes r is defined as: where r for the WENO-JS and the WENO-Z schemes are defined as r = d r (ϵ+ r ) 2 and r = d r , respectively. The polynomial smoothness indicator r ensures that stencils containing shocks are weighted close to zero. The optimal linear weights for the smooth (shock-free) regions are denoted by d r . Also, k is a global smoothness indicator over the full numerical stencil. It should be noted that is typically considered to be 1 × 10 −6 to avoid division by zero. 8 The TENO scheme also introduces a discrete cut-off function that removes nonsmooth candidate stencils from the flux reconstruction. Smooth candidate stencils are included in the reconstruction with their ideal linear weight to reduce the numerical dissipation compared with the WENO scheme. 8 For the TENO schemes r is defined as: where C T is an adjustable cut-off parameter and denotes the smoothness measure defined as = and r = 0, ..., K − 3. 8 The parameters C and q control the levels of dissipation invoked by the non-linear weights and typically have values of 1 and 6, respectively. 8,32 It is worth mentioning that reducing the value of C T reduces the numerical dissipation (which is desirable for compressible turbulence) at the cost of increased spurious oscillations in the vicinity of flow discontinuities. Values of 1 × 10 −5 and 1 × 10 −7 are suggested for the fifth-and sixth-order TENO schemes, respectively. 8,32 A shock sensor (Φ) can be added to the TENO scheme to adjust the threshold of the C T parameter in different regions of the domain and discard fewer stencils from the reconstruction in smooth regions. 8,34,35 In this work, a modified version of the Ducros sensor 39 is used to develop an adaptive TENO scheme (TENO-A) as discussed in Reference [8]. This sensor, which was designed to distinguish between dilatation rate (at shocks) and vorticity (in turbulence), was found to be more robust in the OpenSBLI framework. 8 Nevertheless, the sensor type would not have a noticeable effect on the results of the present study due to the absence of shocks and flow discontinuities. C T is dynamically adjusted by the shock sensor 4 (1 + 4Φ) and ⌊⌋ is the Gauss bracket. 8,35 The adjustable parameters 1 and 2 are taken as 1 = 10.5 and 2 = 4.5 in the present study.
The viscous terms in Equations (1) to (3) are computed with a fourth-order central scheme for the internal domain and over periodic boundaries, for nonperiodic boundaries, a fourth-order one-sided boundary closure is used. 8,40 The solution is advanced in time using a low-storage three-stage explicit Runge-Kutta scheme. 41

WALE LES model
In the present study, the performance of an explicit LES model, specifically the Wall-Adapting Local Eddy-viscosity (WALE) SGS model, is evaluated against those of DNS as well as ILES approaches. The original WALE formulation developed by Nicoud and Ducros 42 is implemented in OpenSBLI. The WALE model considers the effects of the strain and rotation rates of the resolved flow field. Specifically, in this model the turbulent (eddy) viscosity is evaluated as: where S ij is the deformation tensor of the resolved flow field defined as S ij = 1 . Also, S d ij is the traceless symmetric part of the square velocity gradient defined as where g ij is the velocity gradient tensor defined as g ij = u i x j . Moreover, Δ is the subgrid characteristic length scale, which in this study is defined as Δ = V 1/3 , where V is proportional to the local grid spacing as V = ΔxΔyΔz. C w is the SGS model constant, which has a value of C w = 0.325 in the present study. It should be noted that overbar here denotes the resolved flow quantities.
To conduct LES simulations using the WALE model in OpenSBLI, the same governing equations as presented in Equations (1) to (3) are solved for the resolved flow field. However, the effect of the turbulent viscosity is added to the viscous stress tensor and the heat flux (Equations (4) and (5)) as ij = ( . The turbulent Prandtl number has a value of Pr t = 1.0 in the present study. It is worth mentioning that it was found that the results reported here were almost insensitive to the value of Pr t , at least within the range of values found in the literature for similar studies.

OpenSBLI framework
OpenSBLI is a Python-based automatic source code generation and parallel computing framework for finite difference discretization. 37 It generates C codes for the Oxford Parallel library for Structured mesh solver (OPS), an embedded domain-specific language (DSL) with associated libraries and preprocessors to generate parallel executables for applications on multiblock structured meshes. 43,44 OPS can target various computational architectures based on CUDA or OpenCL for GPUs, MPI or OpenMP for CPUs, and CUDA-MPI for multiple GPUs. 44 A comprehensive discussion on the automatic derivation Python interface of OpenSBLI and the automatic procedure of the OPS C code generation can be found in References [37,38]. The results presented in this article are obtained using various backends including a single GPU (CUDA), multiple GPUs (MPI+CUDA), and distributed computing with CPUs (MPI).

PROBLEM SPECIFICATIONS
DNS and LES of subsonic and supersonic plane channel flows with a domain size of 4 H × 2H × 4 3 H are conducted. H denotes the channel half height and is set to H = 1 in the present study. The streamwise (x 0 ), normal (x 1 ), and spanwise (x 2 ) directions are denoted as x, y, and z as shown in Figure 1. Also, u, v, and w denote the velocity components in x, y, and z directions, respectively. Streamwise and spanwise boundaries are periodic and isothermal no-slip walls are assigned to the boundaries in the normal direction. The grid is stretched in the y direction in order to accurately resolve the near wall region. A hyperbolic tangent stretching function 45 is used as follows: where L y denotes the domain length in the y direction (here L y = 2), f s is the stretching factor and governs the distance ratio between two adjacent grid points, j is the grid point number (starts from 0), and n y is the total number of grid points in the y direction. In this study a stretching factor of f s = 1.7 is used for TCF simulations unless otherwise stated. It should be noted that the center of the channel is at y = 0 as shown in Figure 1. The test cases studied here are based on an incompressible channel flow with a friction velocity Reynolds number of Re = 180 46 and a compressible channel flow with a friction velocity Reynolds number of Re = 222 and a friction Mach number of M = 0.082. 47 As mentioned earlier the bulk-averaged density ⟨ ⟩ and its corresponding ) are used to normalize the governing equations. Therefore, the Reynolds number in Equations (4) and (5) is defined as All simulations are initialized by superimposing a random fluctuation over a mean velocity profile in the streamwise direction and imposing random fluctuations in other directions. Specifically, the initial mean streamwise velocity (u) is a piecewise function as: where b = 5.5, is the von Karman constant and has a value of 2.5. Amplitude of the fluctuations is defined as: and sine and cosine disturbances are applied to the velocity profile in all directions. For the x, y, and z directions, the sine disturbances are defined as s x = sin , s y = sin( y), and s z = sin . Similarly, the cosine disturbances are defined as c x = cos , c y = 1 + cos( y), and c z = cos . L x , L y , and L z are the domain dimensions in the x, y, and z directions, respectively, as shown in Figure 1. Finally, the initial values for the velocity components (specified with the subscript 1) are defined as: The initial density is 1 = 1 throughout the domain and the initial energy field is imposed as ( E) i = with the initial pressure field defined as p 1 = 1 M * 2 . The incompressible case is solved using the fourth-order central scheme on a grid with a 129 × 129 × 129 resolution. The compressible case is solved with the fourth-and sixth-order central schemes, fifth-order WENO-JS and WENO-Z schemes, fifth-and sixth-order TENO scheme, and the sixth-order TENO-A scheme. In addition to the grid with 129 × 129 × 129 resolution, the compressible case is also solved on various finer grids as shown in  47 are adjusted in this work to reflect the differences in the normalization approach used. Specifically, the velocity components (and their fluctuating parts) of the reference data are multiplied by a factor of Various time steps from Δt = 5 × 10 −5 to 2 × 10 −4 are used based on the stability of the solution as provided in Table 1. Statistics are collected after a fully turbulent flow is established (at t = 50) and for a duration of T = 1000 (around 1500 flow-through times).

RESULTS AND DISCUSSION
In the following subsections, first the performance of the OpenSBLI computational framework with the central scheme is compared against two incompressible and compressible test cases from the literature 46,47 on a baseline grid resolution (129 × 129 × 129). Also, a comparison is provided on the performance of various central, WENO and TENO schemes on the baseline grid as shown in Table 1. Then the effect of the grid resolution is explored for the central and the TENO schemes as also shown in Table 1. Finally, a comparative study is conducted on DNS and LES modeling (explicit and implicit) based on the test cases shown in Table 5.

Subsonic channel flow
DNS of incompressible turbulent channel flows with friction Reynolds number values around Re = 180 has been reported by several authors since it was first investigated by Kim, Moin and Moser. 48 A review of the earlier simulations of this test case is provided by Vreman and Kuerten. 46 The authors also conducted DNS on a similar channel flow case using two different codes, a finite difference code and a spectral code and provided a comprehensive discussion and comparison against the literature data. In the present work the channel flow case with Re = 180 is used to evaluate the performance of OpenSBLI as a preliminary validation. This is test case 1 of Table 1 and is compared against the spectral results of Vreman and Kuerten. 46 Figure 2 shows a direct comparison between OpenSBLI results obtained with the fourth-order central scheme and the spectral results of Reference [46] based on the mean streamwise velocity and the root mean square (rms) of the fluctuations of the velocity components. OpenSBLI is a compressible solver therefore a friction Mach number of M * = 0.01 is used to mimic incompressible conditions. Also, to maintain a stable solution the solver requires a smaller time step (Δt = 5 × 10 −5 ) compared to incompressible solvers which typically require a time step that ranges from 10 −4 to 10 −3 , depending on the numerical method and spatial resolution. 46 Figure 2 includes OpenSBLI results for statistics collection duration T = 500 and T = 1000. It is evident that statistical convergence is obtained even with T = 500. Nevertheless, in the present work all results are evaluated at T = 1000. With respect to the mean streamwise velocity a very good agreement is found between OpenSBLI and the spectral DNS of Reference [46]. The mean centerline velocity is calculated to be ⟨u c ⟩ = 18.295 which is around 0.18% higher than that of the reference data. 46 However, there are slight deviations in the rms of velocity fluctuations, particularly close to the channel centerline, as shown in Figure 2. Specifically, the rms of velocity fluctuations at the channel centerline are u rms = 0.738, v rms = 0.588, and w rms = 0.575, which are 7.28%, 3.90%, and 6.35% lower than those of the spectral DNS of Reference [46], respectively. Also, the peak of the rms of velocity fluctuations are u rms = 2.558, v rms = 0.841, and w rms = 1.084, which are 3.94%, 0.11%, and 0.36% lower than those of the spectral reference data, respectively. These slight differences are mainly attributed to the spatial resolution, particularly in the streamwise direction, which will be discussed again later in this article.

Scheme comparison
Various high-order shock capturing schemes, including central, WENO, and TENO schemes, are used to perform DNS of the compressible turbulent channel flow case with Re = 222 and M = 0.082 (M = 1.5). Figure 3 shows a comparison based on various mean flow parameters, including the streamwise velocity, Mach number, density, and temperature between the present DNS, using various high-order schemes and a spectral reference DNS conducted by Coleman et al. 47 Values of the obtained friction Reynolds numbers (Re * s and Re s ) and the mean streamwise velocity at the channel centerline (⟨u c ⟩) are tabulated in Table 1. The central schemes show the closest agreement to the reference data for all mean flow parameters presented in Figure 3. Also, the TENO schemes show a much better performance compared with the WENO schemes, with the sixth-order TENO-A scheme showing the closest agreement to the reference data along with the central schemes. For instance, the fourth-and sixth-order central schemes give 0.258% and 0.179% lower mean streamwise velocities at the channel centerline compared to the reference data (⟨u c ⟩ = 18.925), respectively. The sixth-order TENO scheme and the sixth-order TENO-A scheme give a 0.634% higher and 0.147% lower mean streamwise velocity at the channel centerline compared to the reference data, respectively. On the other hand, the fifth-order WENO-JS and after that the fifth-order WENO-Z exhibit the worst performance in predicting the mean flow quantities shown in Figure 3. The latter schemes give 7.98% and 4.97% higher mean streamwise centerline velocities compared to the reference data, respectively.
A direct comparison based on the mean Reynolds stresses between various high-order schemes studied here on the 129 × 129 × 129 grid and the reference solution of Coleman et al 47 is shown in Figure 4. Moreover, corresponding peak values of the stresses shown in Figure 4 are provided in Table 2. A similar trend to what observed for the mean flow parameters is also seen for the Reynolds stresses. Specifically, the central schemes exhibit the closest   Figures 3 and 4 is also reported by Lusher and Sandham. 8 Particularly, the sixth-order TENO-A scheme performed very well in predicting statistical quantities on relatively coarse grids. 8 There are various methods for vortex identification, which can be used to visualize vortical structures in turbulent flows. 49 Among them, Q criterion is one of the widely used methods to visualize complex vortical features. Q identifies vortex cores in a three-dimensional flow as connected regions, where the vorticity prevails over the strain rate S. [49][50][51] Specifically, based on its definition there exists a vortex core, where Q = 1 2 (∥ ∥ 2 − ∥ S∥ 2 ) > 0. However, this is only valid for incompressible flows where Q is directly related to the second invariant of the velocity gradient (II ∇u ). For compressible flows the divergence of the velocity field is non-zero, hence II ∇u ≠ Q. Therefore, as discussed in References [50,51], the extension of Q to compressible flows suffers from ambiguity and fails as a correct vortex identification criterion. Recently, Kolář and Šístek 51 suggested a modification to Q, which introduces compressibility and kinematic corotational aspects required for vortex identification in compressible flows. This approach is named Q M and defines as Q M = Q+ II S ∕2, where II S is the second invariant of S. 51 In the present study Q M is used to make a qualitative comparison of the accuracy of the aforementioned high-order schemes to predict vortical turbulent structures. Specifically, iso-surfaces of Q M with an isovalue of Q M iso = 30 are visualized for the DNS cases of Table 1 with the 129 × 129 × 129 grid resolution as shown in Figure 5. It should be noted that Figure 5 only shows the iso-surfaces that are located within the lower half of the channel to provide a more clear presentation. Large vortical structures that are mainly elongated in the streamwise direction are seen in the flow field of all schemes studied here. However, the schemes with a higher order of accuracy and/or an inherently lower dissipation capture a higher number of such elongated structures and noticeably more medium-and small-size structures. Specifically, the non-dissipative central schemes show a dense field with a wide range of vortex sizes. The TENO-A and TENO schemes show slightly less features compared with the central schemes, but still capture a wide range of structures. Nevertheless, the TENO-A schemes seem to better capture the medium-and small-size vortices compared to their baseline TENO counterparts with the same order of accuracy. On the other hand, the WENO schemes perform poorly by failing to capture a wide range of structures compared to the central and TENO family schemes. In fact, the WENO schemes generally only capture the large elongated vortical structures. The performance shown in Figure 5 supports the earlier discussion provided based on Figures 3 and 4 for the accuracy of the high-order schemes. The elongated structures are formed when near-wall vortices are stretched by shear forces. The near-wall structures are created when the vortices formed in the viscous layer are ejected into the outer layers. This vortex ejection mechanism can be seen in all snapshots of Figure 5 in locations where the vortical structures are almost attached to the wall. A comparable behavior and sensitivity to the numerical scheme and the discretization order, to what is shown in Figure 5, are reported in Reference [52] F I G U R E 5 Q M iso-surfaces on the 129 × 129 × 129 grid with an iso-value of Q M iso = 30 colored by the streamwise velocity for different schemes at t = 50 (the compressible case with Re = 222). The iso-surfaces are only visualized for the lower half of the channel [Colour figure can be viewed at wileyonlinelibrary.com] for a weakly compressible TCF with various schemes including WENO. Specifically, with the WENO scheme, fifth-order of accuracy was not enough to capture the vortical structures and the flow statistics with an acceptable level of accuracy and only a ninth-order WENO scheme exhibited a reasonable performance. 52

Computational cost
The computational cost associated with a high-order scheme is of great importance particularly for high fidelity DNS. all the high-order schemes studied here. The wall time is for one NVIDIA Tesla P100 GPU of Imperial College London's High Performance Computing facility. The central schemes are the most efficient ones in terms of the computational cost. The fourth-order central scheme, which has the lowest cost, is selected as the reference to calculate the relative costs provided in Table 3. The sixth-order central scheme is only 8.2% more expensive than its fourth-order counterpart. The high-order central schemes, considering their excellent performance (as shown in Figures 3 and 4) and noticeably lower computational costs compared to the other schemes, are the most attractive schemes for DNS of compressible turbulent flows, where shockwaves and flow discontinuities do not exist. However, when shock-capturing is required the WENO schemes are the most efficient ones in terms of the computational cost, with 17.2% and 17.8% relative costs for the WENO-JS and WENO-Z schemes, respectively. The WENO-Z scheme is preferred over the WENO-JS since it has only 0.52% higher computational cost but exhibits a much better performance as seen in Figures 3 and 4. The fifth-order TENO scheme has only around 5% higher relative cost compared with the fifth-order WENO-Z scheme but it shows a much better performance and should be considered if an average accuracy and a moderate computational cost are required. The adaptive schemes evaluate the sock sensor and therefore have relatively higher computational costs compared to their standard counterparts. Nevertheless, the sixth-order TENO-A scheme is recommended for when shock-capturing is required and the priority is to obtain a more precise solution with respect to the turbulent structures. A similar computational cost trend as discussed here for the WENO and TENO schemes is also reported for a TGV problem using the OpenSBLI framework in Reference [8].
It is worth mentioning that the computational efficiency of WENO/TENO-based schemes can be noticeably increased by using a hybrid approach in which a corresponding optimal linear scheme is applied in smooth regions. 53,54 For instance, Fu 54 reported a near four times speedup for a hybrid TENO scheme compared with a pure TENO scheme for a 2D Riemann problem. Nevertheless, as shown in Table 3, the most computationally expensive TENO scheme used here has only around 45% higher cost compared to the baseline central scheme, which has significantly fewer operations in comparison to the aforementioned hybrid TENO scheme. This indicates that obtaining a significant speedup by applying a similar hybrid approach is not plausible here and the TENO scheme used in this study is relatively well optimized (also, further modifications may cause load balancing issues).

Grid refinement study
In this section the effect of the spatial resolution is studied by increasing the grid resolution from 129 to 257 in each direction separately for the fourth-order central scheme and the sixth-order TENO and TENO-A schemes. Also, a central test case on a 257 × 257 × 257 grid is examined. Figures 6 and 7 show the mean streamwise velocity and the mean normal stresses for the aforementioned schemes with different spatial resolutions. Corresponding peak values of the normal stresses are also tabulated in Table 4. The mean streamwise velocity is calculated fairly well with the central and TENO schemes as discussed earlier and therefore increasing the grid resolution does not provide a noticeable improvement as seen in Figures 6 and 7 and no significant improvements are gained by refining the spatial resolution in all directions at once as can be seen in Figure 6 and Table 4.
With the sixth-order TENO and TENO-A schemes, increasing the grid resolution in the streamwise direction has a significant effect on the ⟨u ′ u ′ ⟩ stress and the ⟨v ′ v ′ ⟩ stress. Specifically, it results in around 4.5% and 6.0% improvements in the peak ⟨u ′ u ′ ⟩ and ⟨v ′ v ′ ⟩ stresses, respectively. The ⟨w ′ w ′ ⟩ stress also benefits from increasing the resolution in the streamwise direction, however mainly for y + s ≥ 90. It can be concluded that a grid resolution of 257 × 129 × 129 with the central and TENO schemes produces satisfactory results for this compressible TCF case within a relatively affordable computational time. The TENO scheme requires 81.5% higher total wall time (for 250 000 iterations) on the 257 × 129 × 129 grid compared to the 129 × 129 × 129 grid.

Implicit and explicit LES
The main concept of an implicit LES approach is to use the numerical dissipation, arising from discretization schemes (such as TENO), as a kind of turbulent viscosity comparable to that of an explicit SGS LES model. In this section, we evaluate the possibility of using the central, TENO, and TENO-A schemes of OpenSBLI to perform ILES of compressible turbulent channel flows. Also, a direct comparison is conducted between the ILES approach and the WALE LES model, which is implemented in OpenSBLI in the present work. The compressible turbulent channel flow test case with Re = 222 and M = 0.082 is modeled on relatively coarser grids compared to those used for the DNS studies discussed earlier in this article. Specifically, various combinations of grid resolutions including 32 × 129 × 32 and 64 × 64 × 64, high-order schemes and the WALE SGS model are used as presented in Table 5. Figures 8 and 9 show direct comparisons based on the mean streamwise velocity and the normal Reynolds stresses of the test cases presented in Table 5. It should be noted that with explicit LES modeling some of the turbulent energy is held in the SGS model. Therefore, when making comparison against DNS, it might be more appropriate to also filter the DNS solution and compare the respective rms fluctuations. However, since in this study ILES modeling is also investigated and compared against DNS and WALE, such filtering is avoided.
The 32 × 129 × 32 grid resolution is used to decouple the effect of the wall y + and as given in Table 5 Figure 8 it is evident that the WALE model improves the performance of the fourth-order central scheme to predict the mean streamwise velocity. However, the central scheme with the WALE model overpredicts the peak ⟨u ′ u ′ ⟩ stress and underpredicts the peak of the other two normal stresses in comparison to when it is used without the SGS model. This is a well-known behavior and has been reported in the literature. 55 Although the TENO-A scheme improves the prediction of the mean streamwise velocity as shown in Figure 8, it does not show a satisfactory performance in predicting the normal Reynolds stresses on the 32 × 129 × 32 grid resolution. Moreover, from Figure 8 it is clear that by reducing C T of the TENO scheme, from 10 −7 to 10 −9 , its ability to predict the mean statistical parameters noticeably improves. In fact, by reducing C T the numerical dissipation of the scheme reduces. However, this results in a reduced ability to produce a non-spurious stable solution at the location of shockwaves and flow discontinuities.
Similarly to the 32 × 129 × 32 grid resolution, the combination of the WALE model and the central scheme also improves the performance of the code to predict the mean streamwise velocity on the 64 × 64 × 64 resolution as shown in Figure 9. A similar trend also exists for the normal stresses as shown in Figures 8 and 9. However, since the central scheme without the WALE model underpredicts the ⟨u ′ u ′ ⟩ stress on the 64 × 64 × 64 grid, adding the WALE model results in a very close prediction of ⟨u ′ u ′ ⟩ to that of the reference DNS as shown in Figure 9. Comparing against the central scheme with the WALE model, it is evident that the TENO-A scheme exhibits a similar trend for the grid resolutions studied here. From Figures 8 and 9 it can be concluded that the combination of the TENO-A scheme and the WALE model is not beneficial and should be avoided. However, the results show that the TENO-A scheme alone may be used as an ILES model for compressible TCFs. From Table 5 it is evident that with the 64 × 64 × 64 resolution, increasing the stretching factor f s in Equation (12) from 1.7 to 1.9 reduces Δy + s by around 24.5%. However, this improvement in the wall y + value does not seem to noticeably improve the performance as seen in Figure 9.

CONCLUSIONS
The performance of a selection of low-dissipative high-order shock capturing schemes to model compressible turbulent flows was quantified and investigated. The focus was on the capability of those schemes to resolve turbulent structures. Therefore, a channel flow configuration was used to decouple the effect of flow discontinuities and turbulence. Specifically, DNS studies were conducted on subsonic (Re = 180) and supersonic (M = 1.5 and Re = 222) turbulent channel flows using OpenSBLI, a Python-based automatic source code generation and parallel computing framework for finite difference discretization. The schemes studied here were fourth-and sixth-order central, fifth-order WENO-Z and WENO-JS, and fifth-and sixth-order TENO, and adaptive TENO (TENO-A) schemes. Also, a comprehensive study was conducted to evaluate the effect of the spatial resolution on the performance of the central and TENO schemes. Moreover, the accuracy of the TENO family schemes as ILES models was evaluated against DNS and the WALE SGS model. Significant improvements were achieved using the TENO scheme, and particularly the TENO-A scheme, over the WENO schemes with respect to the mean and fluctuating flow quantities in DNS of the compressible TCF. Also, the TENO schemes captured the three-dimensional vortical structures very close to the non-dissipative central schemes and much better than the WENO schemes that failed to resolve a wide range of structures particularly the medium-and small-size vortices.
It was observed that using the WALE model in conjunction with the WENO and TENO schemes resulted in undesirable additional dissipation and should be avoided. However, the TENO family schemes were found to be appropriate for ILES. For instance, reducing the value of the C T constant in the TENO scheme, from 10 −7 to 10 −9 , improved its performance as an ILES model. However, it was still more dissipative than the TENO-A scheme on a similar coarse grid. A further reduction in the C T value would be required to be able to appropriately use the TENO scheme as an ILES model for the compressible TCF case. On the other hand, the TENO-A scheme exhibited a reasonable performance, particularly for the mean streamwise velocity (⟨u⟩) and the ⟨u ′ u ′ ⟩ normal stress. This scheme showed the best potential to be applied as an ILES model, with only modest increases in computational costs, for shock-containing turbulent flows.