A new paradigm of dissipation-controllable, multi-scale resolving schemes for compressible flows
aa r X i v : . [ phy s i c s . c o m p - ph ] J u l A new paradigm of dissipation-controllable, multi-scale resolving schemes forcompressible flows
Xi Deng a , Zhen-hua Jiang b, ∗ , Peter Vincent a , Feng Xiao c , Chao Yan b a Department of Aeronautics,Imperial College London,SW7 2AZ, United Kingdom. b College of Aeronautics Science and Engineering, Beijing University of Aeronautics and Astronautics, Beijing 100191, PR China. c Department of Mechanical Engineering, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo, Japan.
Abstract
The scale-resolving simulation of high speed compressible flow through direct numerical simulation (DNS) or largeeddy simulation (LES) requires shock-capturing schemes to be more accurate for resolving broadband turbulence androbust for capturing strong shock waves. In this work, we develop a new paradigm of dissipation-controllable, shockcapturing scheme to resolve multi-scale flow structures in high speed compressible flow. The new scheme employsa polynomial of n-degree and non-polynomial THINC ( Tangent of Hyperbola for INterface Capturing) functions ofm-level steepness as reconstruction candidates. These reconstruction candidates are denoted as PnTm. From thesecandidates, the piecewise reconstruction function is selected through the boundary variation diminishing (BVD) al-gorithm. Unlike other shock-capturing techniques, the BVD algorithm e ff ectively suppresses numerical oscillationswithout introducing excess numerical dissipation. Then, a controllable dissipation (CD) algorithm is designed forscale-resolving simulations. This novel paradigm of shock-capturing scheme is named as PnTm-BVD-CD. The pro-posed PnTm-BVD-CD scheme has following desirable properties. First, it can capture large-scale discontinuousstructures such as strong shock waves without obvious non-physical oscillations while resolving sharp contact, ma-terial interface and shear layer. Secondly, the numerical dissipation property of PnTm-BVD-CD can be e ff ectivelycontrolled between n + + λ . Thirdly, with λ = . + Keywords: shock capturing, compressible turbulence, multi-scale resolving, dissipation controllable, implicit LES ∗ Corresponding author: Dr. Zhen-Hua Jiang (Email: [email protected])
Preprint submitted to Elsevier July 16, 2020 . Introduction
Numerical simulation of compressible flow involving multi-scale flow features remains one of unsolved issues thatare of great real life application importance. For example, in high speed compressible turbulent flow the interactionsbetween shock waves and turbulence raise challenges to the design of numerical methods because the contradic-tory properties of numerical schemes are required in dealing with discontinuous large-scale and smooth small-scalefeatures simultaneously. In nowadays community of computational fluid dynamics, the multi-scale flow features incompressible turbulence are resolved with so-called scale-resolving simulation through direct numerical simulation(DNS) or large eddy simulation (LES). The scale-resolving simulation requires high resolution numerical schemeswhich are low-dissipative to resolve small-scale structures in transitional and turbulent flows, and meanwhile are ro-bust to stabilize solutions around large-scale features such as strong shocks, contact discontinuities and shear layers.A great deal of e ff ort has been devoted to developing such numerical methods, and at present the shock-capturingschemes are the most popular methods among broad algorithms in the literature.The shock-capturing methods can be classified as two types of quite di ff erent methods: the central schemes andthe upwind schemes. The central schemes consider the flow solutions with smooth solution profiles and treat thediscontinuities as continuous solutions with large gradient. The non-dissipative central schemes are preferable forthe simulations of smooth broadband turbulence involving small-magnitude features. When simulating flows withsharp discontinuities, the artificial viscosity / di ff usivity is usually introduced in the classic central schemes to suppressspurious oscillation across various discontinuities. Nevertheless, the artificial viscosity should be carefully treated inorder to avoid overwhelming the small-scale turbulent structures [1–3]. Contrary to the central schemes, the upwindschemes consider the flow solutions with discontinuous solution profiles. Since the discontinuous solutions exist at theelement interface, the Riemann problems are usually solved to account for the physics of wave propagation. Popularupwind schemes, such as total variation diminishing (TVD) [4], essentially non-oscillatory (ENO) [5, 6] and weightedENO (WENO) [7] methods, are able to capture shocks without oscillations. Unfortunately, these traditional upwindschemes are found to be too dissipative to e ff ectively resolve the small flow structures which play essential role in theturbulent interaction process [8]. In previous studies of [9, 10], the authors demonstrated the importance of reducingnumerical dissipation for turbulence simulations by showing that the second-order central scheme performed betterthan the seventh-order upwind-biased WENO scheme in preserving turbulence energy and spectrum for the explicitLES (ELES) and DNS of compressible turbulence.Having realized the deficiency of traditional shock-capturing schemes which are based on dissipative upwindschemes, the researchers proposed several methods in aims of simulating compressible turbulence flows. In theso-called hybrid schemes for compressible turbulence simulations [11–15, 17, 18], the large-scale features such asdiscontinuities are usually captured with nonlinear upwind-biased shock-capturing schemes like TVD and WENO,while the central flux or upwind flux with low-dissipation / dispersion characteristics is introduced to reduce the nu-merical di ff usion on the smooth solution region. Nevertheless, the performance of the hybrid schemes relies on the2hock indicator that identifies the discontinuity in the solution. The authors of [19] investigated a wide range of suchindicators and concluded that there was no universally better performing method for every problem. Similar withhybrid schemes, several adaptive central-upwind WENO schemes have been designed in [15, 16] by recovering thenon-dissipative central scheme in smooth region. However, as shown in a very recent work [16], the adaptive central-upwind schemes may not be robust enough to suppress numerical oscillations across the discontinuity. Besides thesecentral-upwind schemes, extensive works [20–29] were also done in improving the classical shock-capturing WENOmethod [7] to minimize the numerical dissipation for bandwidth waves. For instance, early works include optimizingthe method of calculating the weights in the WENO scheme [20–22] or adding downwind candidate stencil in thereconstruction process [23, 24] to make schemes more centralized. Recent works include monotonicity-preserving(MP) WENO [25, 26], compact-reconstruction WENO (CRWENO) [27] and targeted ENO (TENO) [28, 29] et al.It is worth mentioning that many shock-capturing schemes described above have also been extended to other highorder discretization framework and implemented as limiting strategies for high speed compressible flow simulations[30–33].As indicated in [34], the numerical dissipation of current shock-capturing schemes mainly comes from the inherentdissipation errors in upwind-biased interpolations, and extra numerical errors introduced by limiting processes whichare designed to suppress numerical oscillations. While there have been researches in improving the spectral propertiesof the shock-capturing scheme by optimizing dispersion and dissipation relation of the scheme [24, 35], most of themethods are often problem dependent. Also as pointed out in [12], a small amount of dissipation is needed in order tosuppress numerical instabilities caused by under-resolved scale of high wavenumber structures. Therefore, it wouldbe beneficial to design high order numerical schemes with controllable dissipation especially for DNS or LES ofcompressible turbulence. In recognition of this point of view, the work in [36, 37] developed minimized dispersionand controllable dissipation (MDCD) schemes with the ability of adjusting dissipation without a ff ecting the optimizeddispersion properties of the method.Recently, a novel reconstruction approach named as boundary variation diminishing (BVD) algorithm [38] hasbeen proposed to construct shock capturing schemes in Godunov finite volume method [39–42]. In the BVD paradigm,the numerical dissipation is e ff ectively reduced by adaptively selecting the reconstruction functions to minimize thevariations at the cell boundary. Using proper BVD algorithms and BVD-admissible functions such as high orderpolynomials and jump-like THINC (Tangent of Hyperbola for INterface Capturing) function, a new class of high orderupwind-biased schemes named as PnTm-BVD (polynomial of n-degree and THINC function of m-level reconstructionbased on BVD algorithm) has been designed [40, 41]. Unlike conventional shock capturing strategies such as TVD andWENO schemes which introduce excessive numerical dissipations through limiting processes and thus undermine thespectral properties of underlying high order interpolations, the PnTm-BVD schemes significantly reduce the numericalerrors of limiting processes and are able to retrieve the dispersion and dissipation properties of high-order unlimitedpolynomials for all wave numbers. With this property, the PnTm-BVD has been applied in compressible turbulencesimulations in [34] through an approach known as implicit LES (ILES) [43–45] where numerical dissipation is used3o model under-resolved flow scale instead of explicit modelling. For instance, the recent work of [42] has shownthat the ILES of turbulent flows using higher order PnTm-BVD scheme, as high as 11th-order, is not only moreaccurate but also more e ffi cient by performing the method on the coarser meshes. However, as an upwind-biasedscheme, the inherent and non-tunable numerical dissipation in PnTm-BVD schemes still raises challenging issues toresolve small-scale structures in the DNS of broadband turbulence or to simulate the under-resolved flow with theILES approach.In order to solve multi-scale features in high speed compressible flows, current work has further explored theconcept of BVD by introducing the controllable dissipation (CD) algorithm. The BVD with controllable dissipation,denoted by BVD-CD, is then used as the guideline for designing a novel class of multi-scale resolving schemes. Thescheme is named as PnTm-BVD-CD (polynomial of n-degree and THINC function of m-level reconstruction based onBVD-CD algorithm) schemes. The PnTm-BVD-CD schemes employ n + / n + / central interpolationbased on a centered grid stencil and THINC functions with adaptive steepness as reconstruction candidates. Throughintroducing a tunable parameter, the final dissipation property of PnTm-BVD-CD schemes can be e ff ectively con-trolled between n + + ff ectively suppress spurious numericaloscillations across the strong shock; 2) The large-scale discontinuous features such as contact, material interface andshear layers can be resolved with less dissipation; 3) The scheme can recover to non-dissipative central interpolationfor smooth solution over all wave number, which is important to solve small scale in DNS as well as resolvable scalein ELES; 4) The e ff ectively controllable dissipation algorithm is believed to work for the ILES to solve the under-resolved flow scale. The rest of the paper is organized as follows. Section 2 briefly introduces the governing equationand the FV formulation applied in the work. Section 3 describes the details of BVD-CD algorithms to constructthe new dissipation controllable, multi-scale resolving shock capturing schemes. Numerical tests are carried out inSection 4 and some conclusions are drawn in Section 5.
2. Governing equation and finite volume method
To simulate high speed compressible viscous flow, the Navier-Stokes equations for a calorically perfect gas aresolved: ∂ρ∂ t + ∇ · m = , (1) ∂ m ∂ t + ∇ · ( m ⊗ u + p δ ) = ∇ · τ, (2) ∂ E ∂ t + ∇ · ( u E + p u ) = ∇ · ( u · τ − q ) , (3)where ρ is the density, m is the momentum, u is the velocity, p is the pressure, δ is the unit tensor, E is the total energy, τ is the viscous stress tensor and q is the heat flux. The equation of state for ideal gas is applied to closure the above4quation system as: p = ( γ − ρ e (4)in which γ is the ratio of specific heats and e is the specific internal energy. With Stokes’ hypothesis for a Newtonianfluid, the viscous stress tensor is calculated as τ = µ S − µ ( ∇ · u ) (5)where strain rate tensor S is defined as S =
12 ( ∇ u + ( ∇ u ) T ) and µ is the dynamic shear viscosity. Following Fourier’slaw, the heat flux is defined as q = − k ∇ T where k is the thermal conductivity and T is the temperature. For simulatingthe compressible turbulence, the accuracy of discretizing the convective terms has been considered critical [46]. Thus,the main challenge is how to design numerical scheme to solve the inviscid Navier-Stokes equations which is alsoknown as Euler equation systems. We use the one-dimensional scalar model equation to illustrate the finite volumemethod to solve Euler equation systems as ∂ q ∂ t + ∂ f ( q ) ∂ x = , (6)where q is the solution variables which can be in characteristic space and f ( q ) is the flux function.A standard finite volume method is applied to 1D scalar model equation of Eq. 6 to obtain the numerical solutions.We divide the computational domain into N non-overlapping cell elements, I i : x ∈ [ x i − / , x i + / ], i = , , . . . , N ,with a uniform grid spacing h = ∆ x = x i + / − x i − / . the volume-integrated average value ¯ q i ( t ) in cell I i is defined as¯ q i ( t ) ≈ ∆ x ˆ x i + / x i − / q ( x , t ) dx . (7)The semi-discrete version of Eq. (6) in the finite volume form can be expressed as an ordinary di ff erential equation(ODE) ∂ ¯ q ( t ) ∂ t = − ∆ x ( ˜ f i + / − ˜ f i − / ) , (8)where the numerical fluxes ˜ f at cell boundaries can be computed by a Riemann solver˜ f i + / = f Riemann i + / ( q Li + / , q Ri + / ) . (9)The left-side value q Li + / and right-side value q Ri + / at cell boundaries are obtained through reconstruction process.The Riemann flux can be generally written in a canonical form as f Riemann i + / ( q Li + / , q Ri + / ) = (cid:16) f ( q Li + / ) + f ( q Ri + / ) (cid:17) − | a i + / | (cid:16) q Ri + / − q Li + / ) (cid:17) , (10)where a i + / stands for the characteristic speed of the hyperbolic conservation law. Based on this canonical form, thecentral schemes lead to non-dissipative flux by considering the flow solutions as smooth solution profiles and regardingdiscontinuities as a continuous solution with large gradient. Moreover, by controlling the jump | q Ri + / − q Li + / | acrossthe cell boundary, the dissipation added to the flux can be controlled.5 . Dissipation controllable, multi-scale resolving schemes The solution quality is very dependent on how to approximate solution in each cell and reconstruct the q Li + / and q Ri + / which are used to calculate numerical flux. To solve smooth small-scale features, high order polynomialsare usually employed to approximate solutions. The chosen polynomials will determine the inherent dissipationand dispersion property of the scheme. In contrary, lower order or monotone functions are more preferable to presentdiscontinuous solutions. Following such idea, upwind-biased high order polynomials and jump-like monotone THINCfunction are employed as candidate interpolants. The final reconstruction function in each cell is selected from thesecandidate interpolants with the multi-step BVD algorithm which is in aims of minimizing numerical dissipation. Wewill review the candidate interpolants before the description of the multi-step BVD algorithm. P n : upwind-biased scheme with polynomial of n-degree A ( n + I i with a polynomial ˜ q Pni ( x ) of degree n . The n + ffi cients of the polynomial are determined by requiringthat ˜ q Pni ( x ) has the same cell average on each cell over an appropriately selected stencil S = { i − n − , . . . , i + n + } with n − + n + = n as 1 ∆ x ˆ x j + / x j − / ˜ q Pni ( x ) dx = ¯ q j , j = i − n − , i − n − + , . . . , i + n + . (11)As shown in [7, 47], 2 r − n − = n + = r −
1. The unknown coe ffi cients of polynomial of 2 r − q Pni ( x ), high order approximation for reconstructed values at the cell boundaries can be obtained by q L , Pni + = ˜ q Pni ( x i + ) and q R , Pni − = ˜ q Pni ( x i − ) . (12)In this work, we apply upwind schemes from fifth order ( r =
3) to ninth order ( r =
5) by using polynomials of 4th, 6th,and 8th degree as one of candidate interpolants. The explicit formulas of q L , Pni + / and q R , Pni − / for (n + • q L , P i + / =
130 ¯ q i − − q i − + q i +
920 ¯ q i + −
120 ¯ q i + , q R , P i − / =
130 ¯ q i + − q i + + q i +
920 ¯ q i − −
120 ¯ q i − . (13) • q L , P i + / = − q i − +
584 ¯ q i − − q i − + q i + q i + − q i + + q i + , q R , P i − / = − q i + +
584 ¯ q i + − q i + + q i + q i − − q i − + q i − . (14)6 q L , P i + / = q i − − q i − + q i − − q i − + q i + q i + − q i + + q i + − q i + , q R , P i − / = q i + − q i + + q i + − q i + + q i + q i − − q i − + q i − − q i − . (15) T m : THINC function with m-level steepness To present the solution in discontinuous regions, we introduce another candidate interpolation function is theTHINC interpolation [48, 49]. The piecewise THINC reconstruction function is written as˜ q Ti ( x ) = ¯ q min + ¯ q max " + θ tanh β x − x i − / x i + / − x i − / − ˜ x i !! , (16)where ¯ q min = min( ¯ q i − , ¯ q i + ), ¯ q max = max( ¯ q i − , ¯ q i + ) − ¯ q min and θ = sgn ( ¯ q i + − ¯ q i − ).The unknown ˜ x i is computed fromconstraint condition ¯ q i = ∆ x ˆ x i + / x i − / ˜ q Ti ( x ) dx . The jump thickness is controlled by the parameter β which determinesthe steepness. Given the reconstruction function ˜ q Ti ( x ), we calculate the boundary values q L , Ti + / and q R , Ti − / by q L , Ti + / = ˜ q Ti ( x i + / ) and q R , Ti − / = ˜ q Ti ( x i − / ) respectively.To present discontinuities with di ff erent steepness, we use THINC functions with β of m -level. Unlike the highorder polynomials, the THINC function can realize non-oscillatory as well as less-dissipative reconstructions fordiscontinuities. A THINC reconstruction function ˜ q Tki ( x ) with β k gives the reconstructed values q L , Tki + / and q R , Tki − / ,( k = , , . . . , m ). We will use m up to three in present study. In the P n T m − BVD − CD schemes, reconstruction values are determined from the candidate interpolants with k -stage ( k = m +
1) BVD-CD algorithm so as to minimize the jumps at cell boundaries as small as possible. Wedenote the reconstruction values at the cell boundary i + / k -th stage BVD-CD as q L ,< k > i + / and q R ,< k > i + / .The k -stage ( k = m +
1) BVD-CD algorithm is formulated as follows. (I): Initial stage ( k = : (I-I) As the first step, use the high-order upwind scheme as the base reconstruction scheme and initialize thereconstructed function as ˜ q < > i ( x ) = ˜ q Pn ( x ) i . (II) Limiting process at the intermediate BVD-CD stage ( k = , . . . , m − : (II-I) Set ˜ q < k > i ( x ) = ˜ q < k − > i ( x )(II-II) Calculate the TBV values for target cell I i from the reconstruction of ˜ q < k > i ( x ) T BV < k > i = (cid:12)(cid:12)(cid:12) q L ,< k > i − / − q R ,< k > i − / (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) q L ,< k > i + / − q R ,< k > i + / (cid:12)(cid:12)(cid:12) (17)7nd from the THINC function ˜ q Tki ( x ) with a steepness β k as T BV
Tki = (cid:12)(cid:12)(cid:12) q L , Tki − / − q R , Tki − / (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) q L , Tki + / − q R , Tki + / (cid:12)(cid:12)(cid:12) . (18)(II-III) Modify the reconstruction function for cells i − i and i + q < k > j ( x ) = ˜ q Tkj ( x ) , j = i − , i , i +
1; if
T BV
Tki < T BV < k > i . (19)(II-IV) Compute the reconstructed values on the left-side of x i + and the right-side of x i − respectively by q L ,< k > i + = ˜ q < k > i ( x i + ) and q R ,< k > i − = ˜ q < k > i ( x i − ) . (20) (III) The dissipation control process at BVD-CD stage ( k = m ) : Given the reconstruction values from the above stage m −
1, for each cell boundary i + modify the reconstructionvalues as q L ,< k > i + = λ q L , Pni + + (1 . − λ ) q R , Pni + ; q R ,< k > i + = λ q R , Pni + + (1 . − λ ) q L , Pni + ; if q L ,< m − > i + = q L , Pni + and q R ,< m − > i + = q R , Pni + (21) (IV) The large scale discontinuity sharpening process at the final BVD-CD stage ( k = m + : (III-I) Given the reconstruction values from the above stage m , compute the TBV using the reconstructed cellboundary values from previous stage by T BV < m > i = (cid:12)(cid:12)(cid:12) q L ,< m > i − / − q R ,< m > i − / (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) q L ,< m > i + / − q R ,< m > i + / (cid:12)(cid:12)(cid:12) , (22)and the TBV for THINC function of β m by T BV
Tmi = (cid:12)(cid:12)(cid:12) q L , Tmi − / − q R , Tmi − / (cid:12)(cid:12)(cid:12) + (cid:12)(cid:12)(cid:12) q L , Tmi + / − q R , Tmi + / (cid:12)(cid:12)(cid:12) . (23)(III-II) Determine the final reconstruction function for cell I i using the BVD algorithm as˜ q < m + > i ( x ) = ˜ q Tmi ; if
T BV
Tmi < T BV < m > i , ˜ q < m > i ; otherwise . (24)(III-III) Compute the reconstructed values on the left-side of x i + and the right-side of x i − respectively by q Li + = ˜ q < m + > i ( x i + ) and q Ri − = ˜ q < m + > i ( x i − ) . (25)Remark 1. With λ decreased from 1.0 to 0.5, the interpolation becomes more centralized. With this parameter,the dissipation property can be e ff ectively controlled between high order upwind scheme and non-dissipativecentral scheme. 8emark 2. The controllable dissipation is believed to work for solving under-resolved small-scale turbulence inILES, where controllable numerical dissipation can be used in replace of the explicit sub-grid scale (SGS)model. With λ = . n + λ = . T − BVD − CD, eighthorder scheme with P T − BVD − CD and tenth order scheme with P T − BVD − CD. According to previous study in[41], in all tests of the present study we use β = . β = . T − BVD − CD, and β = . β = . β = . n T − BVD − CD ( n = ,
8) schemes. As well known, the high speed compressible flow is challengingfor very low-dissipation schemes. Thus we mainly test the central scheme with λ = .
4. Numerical experiments
In this section, we verify the performance of proposed P n T m − BVD − CD schemes in simulating in both Euler andNS equations. We choose an extreme case with λ = . λ = . n + . The approximate dispersion relation (ADR) method described in [52] is applied to study the spectral property ofthe proposed scheme. The numerical dissipation property of the ninth order linear upwind schemes, WENOM andproposed P n T m − BVD − CD schemes are shown in Fig. 1. It is obvious that there is inherent numerical dissipation inupwind-biased schemes despite of very high order. Moreover,the spectral property of ninth order WENOM schemesis inferior to the same order linear scheme for high wavenumber regime. These discrepancies are partly caused byWENO non-linear adaptation which may regard high wavenumber waves as discontinuities. It is noteworthy thatdespite of having developed more accurate WENO schemes in recent year, to recover the upwind-biased scheme inwhole wavenumber regime hasn’t been fully solved through WENO methodology. On the contrary, through newlydevised BVD-CD processes the proposed schemes almost realize the non-dissipative property for all wavenumbers.Although it can be observed that there is still very small numerical dissipation in high wavenumber region for thesixth order scheme, the proposed eighth order and tenth order schemes retrieve corresponding non-dissipative high9rder central schemes. In following sections, we will show that the proposed schemes are also accurate and robust forshock capturing even with such non-dissipative central scheme for smooth regions.In Fig. 2, the dispersion property of proposed schemes are presented by plotting the real parts of modifiedwavenumbers. All proposed schemes show the almost same dispersion property as their corresponding high ordercentral schemes. Moreover, the dispersion properties of proposed schemes surpass those of very high order WENOMschemes especially in high wavenumber regions.In a summary, the proposed schemes are almost able to retrieve the non-dissipative central schemes for all re-solvable wavenumbers in smooth region, which are important for simulating turbulence involving broadband flowstructures through DNS or ELES. It can be expected that by adjusting parameter λ , dissipation property of proposedschemes can be exactly controlled between central scheme and their underlying upwind schemes. It is also noted thatthe dispersion property has been improved as order increased. To further improve the dispersion property which isalso important for turbulence flow, optimization processes of dissipation and dispersion [28, 29] or compact schemes[53, 54] will be considered to construct low-dispersion schemes in our future work. w I m a g Exact9th order UWWENOM9Present (a) 6th order scheme w I m a g Exact9th order UWWENOM9Present (b) 8th order scheme w I m a g Exact9th order UWWENOM9Present (c) 10th order schemeFigure 1: Approximate dissipation properties for di ff erent schemes are analyzed by imaginary parts of modified wavenumber. Comparisons amonglinear upwind schemes, WENOM and proposed P n T m − BVD − CD schemes.
The convergence rate of the proposed P n T m − BVD − CD schemes is studies with an advection test of smoothprofile on gradually refined grids. The initial smooth distribution was given by q ( x ) = sin (2 π x ) , x ∈ [ − , . (26)We ran the computation for one period (at t = .
0) and summarized the numerical errors and the convergence rates forP n T m − BVD − CD schemes in Table 1. As expected, the proposed P n T m − BVD − CD schemes achieved n + L and L ∞ errors have been significantly reduced through newly designed dissipation control processes and the convergencerate has been improved to one order higher. It can also be verified easily that the L and L ∞ errors from the proposed10 R ea l ExactWENOM96th CentralPresent (a) 6th order scheme w R ea l ExactWENOM98th CentralPresent (b) 8th order scheme w R ea l ExactWENOM910th CentralPresent (c) 10th order schemeFigure 2: Approximate dispersion property for di ff erent schemes are analyzed by real parts of modified wavenumber. Comparisons among centralschemes, WENOM and P n T m − BVD − CD schemes.Table 1: Numerical errors and convergence rate for linear advection test. Results are computed by the proposed P n T m − BVD − CD schemes.
Schemes Mesh L errors L order L ∞ errors L ∞ orderP T -BVD-CD 20 3 . × − . × −
40 5 . × − . × − . × − . × − . × − . × − T -BVD-CD 20 2 . × − . × −
40 1 . × − . × − . × − . × − . × − . × − T -BVD-CD 20 2 . × − . × −
40 2 . × − . × − . × − . × − . × − . × − n + Initial distributions containing critical points are more challenging for numerical schemes to distinguish smoothand truly discontinuous profiles. For example, at critical points where the high order derivatives do not simultaneouslyvanish, WENO-type schemes may not reach their formal order of accuracy. In this test, the initial condition is givenfollowing the work [50] by q ( x ) = sin( π x − sin( π x ) π ) , x ∈ [ − , . (27)The computation was conducted for four periods ( t = L and L ∞ of theproposed scheme in Table 2. As shown in results, the proposed algorithm keeps the expected convergence rate evenfor solution containing critical points. This test verifies again the smooth and discontinuous solutions can be e ff ectivelydistinguished by the BVD-CD algorithm. 11 able 2: Numerical errors and convergence rate for advection of a smooth profile containing critical points. Results are calculated by the proposedschemes. Schemes Mesh L errors L order L ∞ errors L ∞ orderP T -BVD-CD 20 3 . × − . × −
40 4 . × − . × − . × − . × − . × − . × − T -BVD-CD 20 4 . × − . × −
40 2 . × − . × − . × − . × − . × − . × − T -BVD-CD 20 9 . × − . × −
40 1 . × − . × − . × − . × − . × − . × − Shock capturing schemes should be able to solve profiles of di ff erent smoothness with high resolutions as well aswithout numerical oscillations. In this subsection, we simulated the propagation of a complex wave [7]. The initialprofile contains both discontinuities and smooth regions with di ff erent smoothness, which is given by q ( x ) = (cid:2) G ( x , β, z − δ ) + G ( x , β, z + δ ) + G ( x , β, z ) (cid:3) if | x + . | ≤ . , | x + . | ≤ . , − |
10 ( x − . | if | x − . | ≤ . , [ F ( x , α, a − δ ) + F ( x , α, a + δ ) + F ( x , α, a ))] if | x − . | ≤ . , , , (28)where functions F and G are defined as G ( x , β, z ) = exp h − β ( x − z ) i , F ( x , α, a ) = q max h − α ( x − a ) , i , (29)and the coe ffi cients are a = . , z = . , δ = . , α = . , β = log (cid:16) δ (cid:17) . (30)The computation was carried out for one period at t = . n T m − BVD − CD schemes with di ff erent order were presented in Fig. 3. It can be seen that all of schemes areessentially oscillation free and capture sharper discontinuities, which is important to solve shear layer, contact planand material interface with high fidelity. With polynomial degree increased, the extreme points of the initial profileare better resolved . Compared with upwind-biased scheme [41], the central scheme shows advantageous in solvingextremes. It is also noteworthy that P n T m − BVD − CD schemes produce almost same results for the discontinuitywhich is resolved by only four cells since the BVD-CD algorithm can properly choose THINC reconstruction functionacross discontinuities. Compared with other central-upwind shock-capturing scheme (see the Fig. 4 in [16]), thepresent solution is one of best with given cell elements. 12 y -0.5 0 0.500.20.40.60.81 Exactnumerical (a) 6th order scheme x y -0.5 0 0.500.20.40.60.81 Exactnumerical (b) 8th order scheme x y -0.5 0 0.500.20.40.60.81 Exactnumerical (c) 10th order schemeFigure 3: Numerical results for advection of complex waves. The numerical solutions at t = . x r ho exactnumerical (a) 6th order scheme x r ho exactnumerical (b) 8th order scheme x r ho exactnumerical (c) 10th order schemeFigure 4: Numerical results of Sod’s problem for density field at t = .
25 with 100 cells.
The Sod’s problem is employed here to test the performance of present schemes in solving Euler equation. Theinitial distribution on computational domain [0 ,
1] was specified as [55]( ρ , u , p ) = (1 , ,
1) 0 ≤ x ≤ . . , , .
1) otherwise . (31)The computation was carried out on a mesh of 100 uniform cells up to t = .
25. The numerical results calculatedfrom the proposed scheme were shown in Fig. 4 for density fields. We observe that P n T m − BVD − CD schemes cansolve the contact discontinuity and right moving shock waves without obvious numerical oscillations. Moreover, thecontact is resolved within only two cells. As expected, P n T m − BVD − CD schemes of di ff erent orders produce similarresults across the contact because the THINC function is selected by the BVD-CD algorithm.13 r ho exactnumerical (a) 6th order scheme x r ho exactnumerical (b) 8th order scheme x r ho exactnumerical (c) 10th order schemeFigure 5: Numerical results of Lax’s problem for density field at time t = .
16 with 100 cells.
We solve the Lax problem [6] which contains relatively strong shock in this subsection. The initial condition isgiven by ( ρ , u , p ) = (0 . , . , . ≤ x ≤ . . , . , . . (32)With the same number of cells as in the previous test case, we got the numerical results at t = .
16. The density fieldis plotted presented in Fig. 5. The proposed P n T m − BVD − CD schemes solve sharp contact and shock waves withoutnumerical oscillation. Compared with one of recently developed central-upwind schemes (Fig. 9 in [16]), the currentscheme produces less oscillatory and less di ff usive results. Dealing with high Mach number flow involving strong shock is very challenging for low-dissipative schemes.Here another Lax’s problem with strong discontinuities is used to test the robustness of di ff erent schemes on shockcapturing. The initial condition is prescribed as( ρ , u , p ) = (1 . , . , . , ≤ x ≤ . . , . , . , otherwise . (33)With initial high pressure ratio, a right-moving Mach 198 shock and a contact is generated. The computation lastsuntil time t = .
012 with 200 cell elements. The numerical solutions from di ff erent schemes for density field areplotted in Fig. 6. The zoomed regions between contact and strong shock are also presented in the figure. We compare n + n T m − BVD − CD with n + n T m − BVD − CD schemes give more faithful solutionwith reduced dissipation and suppressed overshooting. This test shows the proposed scheme is able to capture strongshock robustly even though it recovers to central scheme in smooth region.14 r ho exactnumericalWENOM5 (a) 6th order scheme x r ho exactnumericalWENOM7 (b) 8th order scheme x r ho exactnumericalWENOM9 (c) 10th order schemeFigure 6: Numerical results of strong Lax’s problem with Mach 198 for density field at time t = .
012 with 200 cells. Comparisons are madebetween WENOM and P n T m − BVD − CD at the same order. x r ho exactnumerical (a) 6th order scheme x r ho exactnumerical (b) 8th order scheme x r ho exactnumerical (c) 10th order schemeFigure 7: Numerical results of shock-density wave interaction problem at t = .
18 computed on 400 mesh elements.
A Mach 4 shock wave interacting with a density disturbance is simulating with proposed schemes. The initialcondition is set as ( ρ , u , p ) = (3 . , . , . , if 0 ≤ x ≤ . , (1 + . x − , , , otherwise . (34)The numerical solutions at t = .
18 computed on 400 mesh elements are shown in Fig. 7, where the reference solutionplotted by the solid line is computed by the classical 5th-order WENO scheme with 2000 mesh cells. The results showthe proposed schemes are able to solve problem containing shocks and complex smooth flow features.The advantage of introducing central flux is further demonstrated through another shock-density interaction prob-lem involving waves of higher wavenumber. Similar to [56], the initial condition is specified as( ρ , u , p ) = (1 . , . , . , if x ≤ − . , (1 + . π x ) , , , otherwise . (35)15 r ho -4 -2 0 2 411.21.41.6 exactnumerical (a) 6th order scheme x r ho -4 -2 0 2 411.21.41.6 exactnumerical (b) 8th order scheme x r ho -4 -2 0 2 411.21.41.6 exactnumerical (c) 10th order schemeFigure 8: Numerical results of another shock-density wave interaction problem involving high wavenumber waves. Numerical results at time t = . x r ho -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.51.351.41.451.51.551.61.651.7 exactnumericalWENOM9 (a) 6th order scheme x r ho -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.51.351.41.451.51.551.61.651.7 exactnumericalWENOM9 (b) 8th order scheme x r ho -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.51.351.41.451.51.551.61.651.7 exactnumericalWENOM9 (c) 10th order schemeFigure 9: A zoomed region of the Fig. 8. Comparisons were made between the P n T m − BVD − CD and the ninth order WENOM scheme.
The computation was carried out up to t = . − . , . n T m − BVD − CD schemes with λ = . .9. 2D viscous shock tube In this test, we employ 2D viscous shock tube to verify the performance of the proposed scheme in simulatingunsteady viscous flows. The initial condition is set in a computational domain of [0 , × [0 , .
5] as( ρ , u , v , p ) = (cid:16) , , , γ (cid:17) , if x ≤ . , (cid:16) . , , , . γ (cid:17) , otherwise . (36)The Prandtl number is Pr = .
72 and Reynolds number is Re = ∆ x = ∆ y = . n T m − BVD − CD schemes resolve the delicate structures produced by interaction between shock and boundary layer.The resolved structure of the primary vortex agrees well with the DNS result shown in the Fig. 5(a) of [64]. However,the primary vortex solved by the high order WENOM scheme is similar to the di ff usive result in the Fig. 5(d) of [64].Compared with upwind-biased high order WENOM scheme, the P n T m − BVD − CD schemes are more preferable forDNS of compressible turbulence flow.
In this subsection, we simulate the decaying compressible isotropic turbulence with eddy shocklets [8] throughthe ILES approach. In this problem, if a su ffi ciently high turbulent Mach number is provided weak shock waves willdevelop from the turbulent motions. Thus the coexist of shock waves and turbulence requires the shock-capturingscheme should be able to simultaneously handle shocklets and broadband turbulence motion. The computationaldomain is a cubic volume of side length 2 π L with periodic boundary condition. The initial condition consists of arandom isotropic velocity field velocity fluctuations satisfying a prescribed energy spectrum as E k = u rms , r π k k exp( − k k ) (37)where u rms , is the root mean square turbulence intensity and the wavenumber k =
4. The initial turbulent Machnumber and Taylor-scale Reynolds number are defined as M t , = √ u rms , / c , Re σ, = ρ u rms , σ /µ (38)where σ = k and µ is initial dynamic viscosity based on the initial temperature T . The power law is used tocalculate the viscosity. Initially, the constant density and pressure are prescribed with M t , = . Re σ, = t /τ = τ = σ/ u rms , is the turbulent time scale.17 y (a) P T -BVD-CD x y (b) P T -BVD-CD x y (c) P T -BVD-CD x y (d) 9th order WENOMFigure 10: Simulation results at t = . = n T m − BVD − CD and high order WENO scheme. The computation was conducted with a mesh spacing of ∆ x = ∆ y = . E ( k ) -4 -3 -2 -1 RefnumericalWENOM9 (a) P T -BVD-CD k E ( k ) -4 -3 -2 -1 RefnumericalWENOM9 (b) P T -BVD-CD k E ( k ) -4 -3 -2 -1 RefnumericalWENOM9 (c) P T -BVD-CDFigure 11: ILES of compressible isotropic turbulence. Velocity spectra at t /π = n T m − BVD − CDand WENOM schemes on 64 grids. The λ = .
75 is used for under-resolved ILES.
We first calculate the decaying compressible isotropic turbulence with a coarse mesh of 64 cells. With the coarsemesh, the turbulent model such as LES is required to simulate the under-resolved turbulent structure. In our simulation,the dissipation controllable P n T m − BVD − CD schemes with λ = .
75 are employed to conduct ILES. The additionaldissipation is used for ILES purpose. The numerical results for the velocity spectrum at t /τ = n T m − BVD − CD schemes of di ff erent order produce more accurate results than the ninth orderWENOM scheme. As order increasing, P T − BVD − CD obtains results close to the DNS result. Thus, the numericaldissipation introduced by a tunable λ can be employed to model the under-resolved scale through ILES. Comparedwith ninth order WENOM, the proposed P n T m − BVD − CD obtains more accurate results. More importantly, thenumerical dissipation of the new scheme can be controlled more e ff ectively, which allows investigation of adaptivecontrol to improve ILES results in our future work.Finally, we conduct the convergence test for the compressible isotropic turbulence with di ff erent parameter λ on the 128 mesh elements. The numerical results regarding to the velocity spectrum are presented in the Fig. 12where schemes with added dissipation through λ = .
75 and schemes with central flux through λ = . n T m − BVD − CD schemes provide faithful solutions for the LES andDNS of compressible turbulence flow.
5. Concluding Remarks
In this work, we propose a new paradigm of dissipation-controllable, multi-scale resolving scheme, named asPnTm-BVD-CD, for scale-resolving simulations. The desirable properties of the proposed schemes have been verifiedthrough benchmark tests involving multi-scale flow structures. Compared with high order upwind-biased WENOMschemes or other central-upwind WENO schemes, the PnTm-BVD-CD schemes resolve discontinuous large-scale19 E ( k ) -4 -3 -2 -1 Reflambda=0.75lambda=0.5 (a) P T -BVD-CD k E ( k ) -4 -3 -2 -1 Reflambda=0.75lambda=0.5 (b) P T -BVD-CD k E ( k ) -4 -3 -2 -1 Reflambda=0.75lambda=0.5 (c) P T -BVD-CDFigure 12: Convergence test for compressible isotropic turbulence problem with P n T m − BVD − CD using di ff erent λ . 128 grid elements areemployed. The DNS result marked by the circle is calculated with 256 from [8]. structures with less numerical dissipation as well as less oscillation. With the boundary variation diminishing algo-rithm and the newly designed dissipation-controllable algorithm, the numerical dissipation of PnTm-BVD-CD can bee ff ectively controlled between n + + λ . The n + λ = .
5, which is preferable for solving small-scale structures in DNS. Moreover, the under-resolved small-scale canbe solved with the proposed scheme through the implicit LES approach. Thus this work proposes an accurate androbust scheme for scale-resolving simulation of high speed compressible flow involving broadband turbulence andstrong shock waves.Since the dissipation can be e ff ectively controlled by the proposed scheme, we will explore the adaptive dissipationstrategy inspired by the work [65] to improve the ILES approach in our future work. The proposed methodology willalso be extended to finite di ff erence schemes, compact schemes as well as be equipped with positivity-preservingalgorithm such as [66, 67]. Acknowledgment
This work was supported in part by funds from Engineering and Physical Sciences Research Council (EP / R030340 / eferencesReferences [1] A.W. Cook, W.H. Cabot, A high-wavenumber viscosity for high-resolution numerical methods, J. Comput. Phys. 195 (2004) 594-601.[2] A.W. Cook, Artificial properties for large-eddy simulation of compressible turbulent mixing, Phys. Fluids 19 (2007) 055103 .[3] B. Fiorina , S.K. Lele, An artificial nonlinear di ff usivity method for supersonic reacting flows with shocks, J. Comput. Phys. 222 (2007)246-264.[4] A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (1983) 357-393.[5] A. Harten , B. Engquist , S. Osher , S. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes III, J. Comput. Phys.71 (1987) 231-323.[6] C.W. Shu, S. Osher, E ffi cient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys. 77 (1988) 439-471.[7] G. Jiang, C.W. Shu, E ffi cient implementation of weighted ENO schemes, J. Comput. Phys. 126(1) (1996) 202-228.[8] E. Johnsen, J. Larsson, A.V. Bhagatwala, et al. Assessment of high-resolution methods for numerical simulations of compressible turbulencewith shock waves, J. Comput. Phys. 229(4) (2010) 1213-1237.[9] Mittal, R., Moin, P., Suitability of upwind-biased finite di ff erence schemes for large-eddy simulation of turbulent flows. AIAA J. 35, 1415-1417 (1997)[10] Larsson, J., Lele, S.K., Moin, P., E ff ect of numerical dissipation on the predicted spectra for compressible turbulence. In: Annual ResearchBriefs, pp. 45-57. Center for Turbulence Research, Stanford (2007)[11] N.A. Adams, K. Shari ff , A high-resolution hybrid compact-ENO scheme for shock-turbulence interaction problems, J. Comput. Phys. 127(1996) 27-46.[12] S. Pirozzoli, Conservative hybrid compact-WENO schemes for shock-turbulence interaction, J. Comput. Phys. 178 (2002) 81-117.[13] Y. Ren, M. Liu, H. Zhang, A characteristic-wise hybrid compact-WENO scheme for solving hyperbolic conservation laws, J. Comput. Phys.192 (2003) 365-386.[14] D.J. Hill, D.I. Pullin, Hybrid tuned center-di ff erence-WENO method for large eddy simulations in the presence of strong shocks, J. Comput.Phys. 194 (2004) 435-450.[15] X.Y. Hu, Q. Wang, N.A. Adams, An adaptive central-upwind weighted essentially non-oscillatory scheme. J. Comput. Phys. 229 (2010),8952-8965.[16] H. Cong, L. L. Chen, A new adaptively central-upwind sixth-order WENO scheme. J. Comput. Phys. 357, (2018) 1-15.[17] Z. Jiang, C. Yan, J. Yu, Hybrid central-upwind finite volume schemes for solving the Euler and Navier-Stokes equations, Computers andMathematics with Applications 72(2016) 2241-2258.[18] Z. Jiang, C. Yan, J. Yu, E ffi cient methods with higher order interpolation and MOOD strategy for compressible turbulence simulations, J.Comput. Phys. 371 (2018) 528-550.[19] G. Li, J. Qiu, Hybrid weighted essentially non-oscillatory schemes with di ff erent indicators, J. Comput. Phys. 229 (2010), 8105-8129.[20] A. Henrick, T. Aslam, J. Powers, Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points, J.Comput. Phys. 207 (2005) 542-567.[21] R. Borges , M. Carmona , B. Costa, W.S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws,J. Comput. Phys. 227 (2008) 3191-3211.[22] Y. Shen, G. Zha, Improvement of weighted essentially non-oscillatory schemes near discontinuities, Comput. Fluids, 96 (2014) 1-9.[23] M.P. Mart´ın, E.M. Taylor, M. Wu, V.G. Weirs, A bandwidth-optimized WENO scheme for the e ff ective direct numerical simulation ofcompressible turbulence. J. Comput. Phys. 220 (2006) 270-289[24] E.M. Taylor, M. Wu, M.P. Mart´ın, Optimization of nonlinear error for weighted essentially nonoscillatory methods in direct numericalsimulations of compressible turbulence, J. Comput. Phys. 223 (2007) 384-397[25] A. Suresh , H.T. Huynh , Accurate monotonicity-preserving schemes with Runge-Kutta time stepping, J. Comput. Phys. 136 (1997) 83-99.[26] J. Fang , Z. Li , L. Lu, An optimized low-dissipation monotonicity-preserving scheme for numerical simulations of high-speed turbulentflows, J. Sci. Comput. 56 (2013) 67-95.[27] D. Ghosh, J.D. Baeder, Weighted Non-linear Compact Schemes for the Direct Numerical Simulation of Compressible, Turbulent Flows, J.Sci. Comput. 61 (2014) 61-89.[28] L. Fu , X.Y. Hu , N.A. Adams, A family of high-order targeted ENO schemes for compressible-fluid simulations, J. Comput. Phys. 305 (2016)333-359.[29] L. Fu, A low-dissipation finite-volume method based on a new TENO shock-capturing scheme, Comput. Phys. Commun. 235 (2019) 25-39.[30] C.W. Shu, On high order accurate weighted essentially non-oscillatory and discontinuous Galerkin schemes for compressible turbulencesimulations, Philosophical Transactions of the Royal Society A, 371 (2013) 20120172.[31] H. T. Huynh, Z. J. Wang, P. E. Vincent, High-Order Methods for Computational Fluid Dynamics: A Brief Review of Compact Di ff erentialFormulations on Unstructured Grids. Comput. Fluids, 98 (2014) 209-220.[32] B. Ahrabi, K. Anderson, J. Newman, An adjoint-based hp-adaptive stabilized finite-element method with shock capturing for turbulent flows,Comput. Methods Appl. Mech. Eng. 318 (2017) 1030-1065.[33] B. Xie, X. Deng, S.J. Liao, F. Xiao, High-order multi-moment finite volume method with smoothness adaptive fitting reconstruction forcompressible viscous flow, J. Comput. Phys. 394 (2019) 559-593.[34] X. Deng, Z. Jiang, F. Xiao, C. Yan, Implicit large eddy simulation of compressible turbulence flow with PnTm-BVD scheme, Appl. Math.Model. 77 (2020) 17-31.[35] Z.J. Wang, R.F. Chen, Optimized weighted essentially nonoscillatory schemes for linear waves with discontinuity, J. Comput. Phys. 174(2001) 381-404.[36] Z. Sun, Y. Ren, L. Cedric, et al. A class of finite di ff erence schemes with low dispersion and controllable dissipation for DNS of compressibleturbulence, J. Comput. Phys. 230 (2011) 4616-4635.
37] Q. Wang, Y. Ren, Z. Sun, Y. Sun, Low dispersion finite volume scheme based on reconstruction with minimized dispersion and controllabledissipation, SCI. CHINA. PHYS. MECH. 56 (2013) 423-431.[38] Z. Sun, S. Inaba, F. Xiao, Boundary variation diminishing (BVD) reconstruction: a new approach to improve Godunov schemes, J. Comput.Phys. 322 (2016) 309-325.[39] X. Deng, S. Inaba, B. Xie, K.M. Shyue , F. Xiao, High fidelity discontinuity-resolving reconstruction for compressible multiphase flows withmoving interfaces, J. Comput. Phys. 371 (2018) 945-966.[40] X. Deng , Y. Shimizu , F. Xiao , A fifth-order shock capturing scheme with two-stage boundary variation diminishing algorithm, J. Comput.Phys. 386 (2019) 323-349.[41] X. Deng , Y. Shimizu , B. Xie, F. Xiao, Constructing higher order discontinuity-capturing schemes with upwind-biased interpolations andboundary variation diminishing algorithm, Comput. Fluids, 200 (2020) 104433.[42] Z. Jiang, C. Yan, J. Yu, A higher order interpolation scheme of finite volume method for compressible flow on curvilinear grids, Communi-cations in Computational Physics (2020) (in press)[43] F.F. Grinstein, L.G. Margolin, W.J. Rider, Implicit Large Eddy Simulation, Computing Turbulent Fluid Dynamics, Cambridge UniversityPress, 2007.[44] C. Fureby, F.F. Grinstein, Large-eddy simulation of High-Reynolds-Number free and wall-bounded flows, J. Comput. Phys. 181 (2002) 68-97.[45] D. Drikakis, M. Hahn, A. Mosedale, B. Thornber, Large eddy simulation using high-resolution and high-order methods, Philos. Trans. R.Soc. 367 (2009) 2985-2997.[46] J.R. DeBonis, Solutions of the Taylor-Green Vortex Problem Using High Resolution Explicit Finite Di ff erence Methods, AIAA Paper (2014)2013-0382.[47] D.S. Balsara, C.W. Shu, Monotonicity prserving WENO schemes with increasingly high-order of accuracy, J. Comput. Phys. 160 (2000)405-452.[48] F. Xiao, S. Ii, C. Chen, Revisit to the THINC scheme: a simple algebraic VOF algorithm, J. Comput. Phys. 230 (2011) 7086-7092.[49] F. Xiao, Y. Honma, T. Kono, A simple algebraic interface capturing scheme using hyperbolic tangent function, Int. J. Numer. Methods Fluids.48 (2005) 1023-1040.[50] A.K. Henrick, T.D. Aslam, J.M. Powers, Mapped weighted essentially non-oscillatory schemes: achieving optimal order near critical points,J. Comput. Phys. 207 (2005) 542-567.[51] S. Gottlieb, L.A.J. Gottlieb, Strong stability preserving properties of Runge-Kutta time discretization methods for linear constant coe ffi cientoperators, J. Sci. Comput. 18 (1) (2003) 83-109.[52] S. Pirozzoli, On the spectral properties of shock-capturing schemes, J. Comput. Phys. 219 (2006) 489-497.[53] X. Deng, H. Zhang, Developing High-Order Weighted Compact Nonlinear Schemes, J. Comput. Phys. 165 (1) (2000) 22-44.[54] X. Liu, S. Zhang, H. Zhang, C. W. Shu, A new class of central compact schemes with spectral-like resolution II: Hybrid weighted nonlinearschemes, J. Comput. Phys. 284 (2015) 133-154.[55] G.A. Sod, A survey of several finite di ff erence methods for systems of nonlinear hyperbolic conservation laws, J. Comput. Phys. 27 (1978)1-31.[56] V.A. Titarev, E.F. Toro, Finite-volume WENO schemes for three-dimensional conservation laws, J. Comput. Phys. 201 (2004) 238-260.[57] C.W. Schulz-Rinne, Classification of the Riemann problem for two-dimensional gas dynamics, SIAM J. Math. Anal 24 (1993) 76-88.[58] A. Kurganov, E. Tadmor, Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers, Numer.Methods Partial Di ff erential Equations 18 (2002) 584-608.[59] M. Dumbser, O. Zanotti, R. Loub`ere, S. Diot, A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolicconservation laws, J. Comput. Phys. 278 (2014) 47-75.[60] C.Y. Jung, T.B. Nguyen, Fine structures for the solutions of the two-dimensional Riemann problems by high-order WENO schemes, Adv.Comput. Math. (2017) 1-28.[61] P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54 (1984) 115-173.[62] A. Rault, G. Chiavassa, R. Donat, Shock-vortex interactions at high Mach numbers, J. Sci. Comput. 19 (2003) 347-371[63] M. Dumbser, M. K¨aser, V.A. Titarev, E.F. Toro, Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinearhyperbolic systems, J. Comput. Phys. 226 (2007) 204-243.[64] V. Daru, C. Tenaud, Evaluation of TVD high resolution schemes for unsteady viscous shocked flows. Comput. Fluids, 30 (2001) 89-113.[65] X. Nogueira, L. Ramirez, et al. An a posteriori-implicit turbulent model with automatic dissipation adjustment for Large Eddy Simulation ofcompressible flows. Comput. Fluids, 197 (2020) 104371.[66] S. Clain, S. Diot, R. Loub`ere, A high-order finite volume method for systems of conservation laws-multi-dimensional optimal order detection(MOOD), J. Comput. Phys. 230 (2011) 4028-4050.[67] S. Diot, S. Clain, R. Loub`ere, Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured mesheswith very high-order polynomials, Comput. & Fluids 64 (2012) 43-63 (2012).erential Equations 18 (2002) 584-608.[59] M. Dumbser, O. Zanotti, R. Loub`ere, S. Diot, A posteriori subcell limiting of the discontinuous Galerkin finite element method for hyperbolicconservation laws, J. Comput. Phys. 278 (2014) 47-75.[60] C.Y. Jung, T.B. Nguyen, Fine structures for the solutions of the two-dimensional Riemann problems by high-order WENO schemes, Adv.Comput. Math. (2017) 1-28.[61] P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, J. Comput. Phys. 54 (1984) 115-173.[62] A. Rault, G. Chiavassa, R. Donat, Shock-vortex interactions at high Mach numbers, J. Sci. Comput. 19 (2003) 347-371[63] M. Dumbser, M. K¨aser, V.A. Titarev, E.F. Toro, Quadrature-free non-oscillatory finite volume schemes on unstructured meshes for nonlinearhyperbolic systems, J. Comput. Phys. 226 (2007) 204-243.[64] V. Daru, C. Tenaud, Evaluation of TVD high resolution schemes for unsteady viscous shocked flows. Comput. Fluids, 30 (2001) 89-113.[65] X. Nogueira, L. Ramirez, et al. An a posteriori-implicit turbulent model with automatic dissipation adjustment for Large Eddy Simulation ofcompressible flows. Comput. Fluids, 197 (2020) 104371.[66] S. Clain, S. Diot, R. Loub`ere, A high-order finite volume method for systems of conservation laws-multi-dimensional optimal order detection(MOOD), J. Comput. Phys. 230 (2011) 4028-4050.[67] S. Diot, S. Clain, R. Loub`ere, Improved detection criteria for the multi-dimensional optimal order detection (MOOD) on unstructured mesheswith very high-order polynomials, Comput. & Fluids 64 (2012) 43-63 (2012).