Comparison of the performance of high-order schemes based on the gas-kinetic and HLLC fluxes
CComparison of the performance of high-order schemes based on thegas-kinetic and HLLC fluxes
Xiaojian Yang a , Xing Ji b , Wei Shyy a , Kun Xu a,b,c, ∗ a Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology,Clear Water Bay, Kowloon, Hong Kong b Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay,Kowloon, Hong Kong c Shenzhen Research Institute, Hong Kong University of Science and Technology, Shenzhen, China
Abstract
In this paper, a comparison of the performance of two high-order finite volume methods basedon the gas-kinetic scheme (GKS) and HLLC fluxes is carried out in structured rectangularmesh. For both schemes, the fifth-order WENO-AO reconstruction is adopted to achieve ahigh-order spatial accuracy. In terms of temporal discretization, a two-stage fourth-order(S2O4) time marching strategy is adopted for WENO5-AO-GKS scheme, and the fourth-order Runge-Kutta (RK4) method is employed for WENO5-AO-HLLC scheme. For theviscous flow computation, the GKS includes both inviscid and viscous fluxes in the evolu-tion of a single cell interface gas distribution function. While for the WENO5-AO-HLLCscheme, the inviscid flux is provided by HLLC Riemann solver, and the viscous flux is dis-cretized by a sixth-order central difference method. Based on the tests of forward Mach stepand viscous shock tube, both schemes show outstanding shock capturing property. Fromthe Titarev-Toro and double shear layer tests, WENO5-AO-GKS scheme seems to have abetter resolution than WENO5-AO-HLLC scheme. Both schemes show excellent robustnessin extreme cases, such as the Le Blanc problem. From the cases of the Noh problem and thecompressible isotropic turbulence, WENO5-AO-GKS scheme shows favorite robustness. Inthe compressible isotropic turbulence and three-dimensional Taylor-Green vortex problems,WENO-AO-GKS can use a CFL number up to 0 .
5, instead of 0 . Keywords:
WENO-AO reconstruction, gas-kinetic scheme (GKS), HLLC Riemann solver. ∗ Corresponding author
Email addresses: [email protected] (Xiaojian Yang), [email protected] (Xing Ji),
Preprint submitted to Elsevier September 9, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] S e p . Introduction The development of high-order schemes has been the main research direction in the cur-rent computational fluid dynamics. The targeting scheme should be accurate, robust, andefficient. The finite volume scheme is mainly composed of spatial reconstruction, flux eval-uation, and temporal discretization. The successful high-order reconstructions include theessentially non-oscillatory (ENO) and weighted essentially non-oscillatory (WENO) scheme[13, 18, 23]. There exists many modified versions of WENO, such as WENO-JS [18], WENO-Z [4], central WENO (CWENO) [19], WENO with adaptive order (WENO-AO) [1], multi-resolution WENO [45], etc.Besides the importance of initial reconstruction, the flux evaluation and temporal up-dating method also play important roles in the determination of the quality of the schemes.In the past decades, the gas-kinetic scheme (GKS) is mainly focusing on the time accurateflux function for capturing the Euler and Navier-Stokes solutions. The GKS is based onthe kinetic Bhatnagar-Gross-Krook (BGK) model and the Chapman-Enskog expansion isused for the flux evaluation [3, 8]. The scheme has been systematically developed for theflow computation from low-speed to hypersonic one [39, 40]. In GKS, a time-dependent gasdistribution function at the cell interface is obtained and covers a physical process from thekinetic free particle transport to the hydrodynamic NS wave propagation. In the smoothregion, GKS can accurately recover the Euler or Navier-Stokes solution. In the discontinuityregion, the particle free transport mechanism introduces the numerical dissipation withina shock layer and stabilize the numerical shock structure. Different from the traditionalCFD methods based on the macroscopic government equations directly, GKS has multi-scale property. Depending on the ratio of time step ∆ t over the particle collision time τ ,the flux function in GKS makes a smooth transition from the upwind flux vector splitting(kinetic scale) to the central difference (hydrodynamic scale). GKS has been adopted inmulticomponent flow [38, 26], acoustic computation [43], turbulence simulation [22, 7, 30],and hypersonic flow [21], etc. Furthermore, a unified GKS (UGKS) has been developedfor all flow regimes from rarefied to continuum one [41]. At the same time, in order todevelop high-order GKS, many techniques in CFD have been used in the kinetic schemes.The WENO reconstruction has been adopted to improve spatial accuracy [24]. Also, thehigh-order compact GKS on both structured and unstructured meshes have been developed[14, 17, 44]. Since the flux function in GKS is time-dependent, which provides not only thenumerical flux but also its time derivative. Therefore, multi-stage multi-derivative (MSMD)methods can be employed for time marching in GKS [16]. Particularly, a two-stage fourth-order (S2O4) temporal discretization for GKS has been developed with favorable numericalperformance [28, 27].In the CFD community, mostly the exact or approximate Riemann problems are usedin the flux construction [11]. One of the outstanding approximate Riemann solvers is theHLL flux [12]. In HLL, a configuration including two waves and three constant states isassumed. In order to improve the capacity of capturing contact surfaces in HLL solver, Toro [email protected] (Wei Shyy), [email protected] (Kun Xu)
2. WENO-AO-GKS and WENO-AO-HLLC
The WENO-AO reconstruction was proposed by Balsara et al. [1]. To meet the re-quirement for a fourth-order scheme in both space and time, the fifth-order WENO5-AOreconstruction is selected. Assume that Q is the cell-averaged variable, and Q is the recon-structed variable. Both Q and Q can be conservative or characteristic variables. To achievefifth-order spatial accuracy for Q , three sub-stencils S k , k = 0 , , Q li +1 / and right Q ri − / interface values at x i − / and x i +1 / . These three sub-stencils S k are, S = { I i − , I i − , I i } , S = { I i − , I i , I i +1 } , S = { I i , I i +1 , I i +2 } . S k , a unique quadratic polynomial p r k ( x ) is constructed by the require-ments, 1∆ x (cid:90) I i − j − k − p r k ( x )d x = Q i − j − k − , j = − , , , (1)and each p r k ( x ) achieves a third-order spatial accuracy in smooth flow region.On a large stencil S = { S , S , S } , a unique fourth-order polynomial p r ( x ) is obtainedby 1∆ x (cid:90) I i + j p r ( x )d x = Q i + j , j = − , − , , , . After determining the above reconstructions based on different stencils, p r ( x ) is definedagain as, p r ( x ) = γ [ 1 γ p r ( x ) − (cid:88) γ k γ p r k ( x )] + (cid:88) γ k p r k ( x ) , (2)where γ k , k = 0 , , , γ = γ Hi , γ = γ = (1 − γ Hi )(1 − γ Lo ) / , γ = (1 − γ Hi ) γ Lo , where γ Hi ∈ [0 . , .
95] and γ Lo ∈ [0 . , . (cid:80) γ k =1 and γ k > , k = 0 , , ,
3. In the current study, γ Hi = 0 .
85 and γ Lo = 0 .
85 are used.To deal with discontinuities, the WENO-Z type [4] non-linear weights are adopted, ω k = γ k [1 + δ ( β k + (cid:15) ) ] , where δ is the global smooth indicator, and it is defined as δ = 13 ( | β r − β r | + | β r − β r | + | β r − β r | ) = O (∆ h ) . More specifically, β k = β r k , k = 0 , ,
2, are the smooth indicator of sub-stencil S k , and β = β r is the smooth indicator of the large stencil S . The explicit formulas of β k can referto [1]. Besides, (cid:15) is a positive small number to avoid zero for denominator with (cid:15) = 10 − .Then, the normalized weights ω k can be defined as follows, ω k = ω k (cid:80) ω q . P AO (5 , ( x ) = ω [ 1 γ p r ( x ) − (cid:88) γ k γ p r k ( x )] + (cid:88) ω k p r k ( x ) . (3)The reconstructed left interface value Q li +1 / and the corresponding derivative become, Q li +1 / = P AO (5 , ( x i +1 / ) , ( Q lx ) i +1 / = P AO (5 , x ( x i +1 / ) . Similarly, the right interface value Q ri − / and its derivative can also be determined by, Q ri − / = P AO (5 , ( x i − / ) , ( Q rx ) i − / = P AO (5 , x ( x i − / ) . The reconstructed value and its normal derivative at the Gaussian quadrature points areobtained from the above procedure. Since GKS needs not only normal derivatives ( Q x ), butalso tangential derivative ( Q y , Q z ), the multi-dimensional reconstruction is performed. Thedetails are given in [15]. Here the GKS in 2D case is presented and the scheme in 3D can be obtained similarly.The two-dimensional BGK equation is written as [3], f t + u · ∇ f = g − fτ , (4)where u is the particle velocity, f is the gas distribution function, g is the correspondingequilibrium state, and τ is the collision time. The collision term satisfies the compatibilitycondition (cid:90) g − fτ ψ dΞ = 0 , (5)where ψ = (1 , u, v,
12 ( u + v + ξ )) T , the internal variables ξ = ξ + ... + ξ K , dΞ =d u d v d ξ ... d ξ K , K is the internal degree of freedom, i.e. K = (4 − γ ) / ( γ −
1) for two-dimensional flows, and γ is the specific heat ratio.In the continuum regime, the gas distribution function can be expanded as f = g − τ D u g + τ D u ( τ D u ) g − τ D u [ τ D u ( τ D u ) g ] + ..., where D u = ∂/∂t + u · ∇ . The corresponding macroscopic equations can be derived bytruncating on different order of τ . For example, when the zeroth-order truncation is taken,i.e. f = g , the Euler equations can be derived. When the first-order truncation is used, f = g − τ ( ug x + vg y + g t ) , (6)5he Navier-Stokes equations can be derived with τ = µ/p . The difficulties for the develop-ment of a reliable gas-kinetic scheme is the possible discontinuity of flow variables at thecell interface, where the above Chapman-Enskog expansion cannot be used directly for theflux evaluation, and the time evolution solution of the gas distribution function at the cellinterface has to be constructed properly from a piecewise discontinuous initial condition.Based on the conservation laws in a discretized space of control volume S ij = [ x i − ∆ x/ , x i + ∆ x/ × [ y j − ∆ y/ , y j + ∆ y/ W ij d t = − x ( F i +1 / ,j ( t ) − F i − / ,j ( t )) − y ( G i,j +1 / ( t ) − G i,j − / ( t )) , (7)where W ij = [ ρ, ρU, ρV, ρE ] T are the cell-averaged conservative variables. F i ± / ,j ( t ) and G i,j ± / ( t ) are the time-dependent numerical fluxes across the cell interfaces in x and y directions respectively. The fluxes can be obtained by a time-dependent gas distributionfunction f at the corresponding cell interface. To achieve the accuracy in space, the Gaussianquadrature is used. Taking the numerical fluxes in x directions F i +1 / ,j ( t ), for example, F i +1 / ,j ( t ) = 1∆ y (cid:90) y j +1 / y j − / F i +1 / ( y, t )d y = (cid:88) (cid:96) =1 ω (cid:96) F i +1 / ,j (cid:96) ( t ) , (8)two Gaussian quadrature points y j (cid:96) = y j + ( − (cid:96) − √ y , (cid:96) = 1 ,
2, and the correspondingweights ω = ω = 1 / F i +1 / ,j (cid:96) ( t ), (cid:96) = 1 ,
2, are numerical fluxes at the Gaussian quadrature points, F i +1 / ,j (cid:96) ( t ) = (cid:90) ψ uf ( x i +1 / , y (cid:96) , t, u, v, ξ )dΞ , (9)where f ( x i +1 / , y (cid:96) , t, u, v, ξ ), (cid:96) = 1 ,
2, are the gas distribution function at the Gaussianpoints. To obtain the numerical fluxes, the integral solution of BGK equation Eq.(4) atpoint ( x i +1 / , y (cid:96) ) and time t is used, f ( x i +1 / , y (cid:96) , t, u, v, ξ ) = 1 τ (cid:90) t g ( x (cid:48) , y (cid:48) , t (cid:48) , u, v, ξ ) e − ( t − t (cid:48) ) /τ d t (cid:48) + e − t/τ f ( − ut, − vt, u, v, ξ ) , (10)where ( x i +1 / , y (cid:96) ) = (0 ,
0) for the simplification of the notation, x = x (cid:48) + u ( t − t (cid:48) ) and y = y (cid:48) + v ( t − t (cid:48) ) are the trajectory of particles. f is the initial gas distribution function attime t = 0, and g is the corresponding equilibrium state.In the integral solution Eq.(10), the initial gas distribution function can be constructedas f = f l ( x, y, u, v ) H ( x ) + f r ( x, y, u, v )(1 − H ( x )) , (11)where H ( x ) is the Heaviside function, f l and f r are the initial gas distribution functions onthe left and right side of one cell interface, which can be determined by the corresponding6acroscopic variables. The initial gas distribution function f k , k = l, r , is constructed as f k = g k (cid:0) a k x + b k y − τ ( a k u + b k v + A k ) (cid:1) , where g l and g r are the Maxwellian distribution functions on the left and right hand sidesof a cell interface, and they can be determined by the corresponding conservative variables W l and W r . The coefficients a l , a r , b l , b r are related to the spatial derivatives in normaland tangential directions, which can be obtained from the corresponding derivatives of theinitial macroscopic variables, (cid:10) a l (cid:11) = ∂ W l /∂x, (cid:104) a r (cid:105) = ∂ W r /∂x, (cid:10) b l (cid:11) = ∂ W l /∂y, (cid:104) b r (cid:105) = ∂ W r /∂y, where (cid:104) ... (cid:105) means the moments of the Maxwellian distribution function, (cid:104) ... (cid:105) = (cid:90) ψ ( ... ) g dΞ . The non-equilibrium parts on the Chapman-Enskog expansion have no net contribution tothe conservative variables, (cid:10) a l u + b l v + A l (cid:11) = 0 , (cid:104) a r u + b r v + A r (cid:105) = 0 , and therefore the coefficients A l and A r , related to time derivatives, can be obtained. Afterthe determination of f , the equilibrium state g around the cell interface is modeled as, g = g (cid:0) ax + by + ¯ At (cid:1) , (12)where g is the local equilibrium at point ( x i +1 / , y (cid:96) ) and can be determined by the compat-ibility condition, (cid:90) ψ g dΞ = W = (cid:90) u> ψ g l dΞ + (cid:90) u< ψ g r dΞ , (cid:90) ψ ag dΞ = ∂ W /∂x = (cid:90) u> ψ a l g l dΞ + (cid:90) u< ψ a r g r dΞ , (13) (cid:90) ψ bg dΞ = ∂ W /∂y = (cid:90) u> ψ b l g l dΞ + (cid:90) u< ψ b r g r dΞ , and (cid:10) au + bv + ¯ A (cid:11) = 0 . After constructing the initial gas distribution function f and the equilibrium state g , andsubstituting Eq.(11) and Eq.(12) into Eq.(10), the time-dependent distribution function7 ( x i +1 / , y (cid:96) , t, u, v, ξ ) at a cell interface can be expressed as, f ( x i +1 / ,j (cid:96) , t, u, v, ξ ) =(1 − e − t/τ ) g + [( t + τ ) e − t/τ − τ ]( au + bv ) g +( t − τ + τ e − t/τ ) ¯ Ag + e − t/τ g r [1 − ( τ + t )( a r u + b r v ) − τ A r ] H ( u )+ e − t/τ g l [1 − ( τ + t )( a l u + b l v ) − τ A l ](1 − H ( u )) . (14)The collision time τ in Eq.(14) is defined by τ = µp + c | p l − p r p l + p r | ∆ t, for viscous flow computation, where p l and p r are the pressures on the left and right sidesof the cell interface, and p is the pressure at the interface from the equilibrium state. Here∆ t is the time step. For inviscid flow, the τ is given by τ = c ∆ t + c | p l − p r p l + p r | ∆ t, where c = 0 . c = 1 ∼ The two-stage fourth-order temporal discretization, originally developed for the gener-alized Riemann problem (GRP) solver [20], has been applied to GKS [28]. A fourth-ordertime-accurate GKS can be constructed using the second-order flux function Eq.(14). Forthe time-dependent equations, ∂ W ∂t = L ( W ) , (15)with the initial condition at t n , W ( t = t n ) = W n , (16)where L is the spatial operator of flux obtained in Eq.(7), a fourth-order temporal accuratesolution for W ( t ) at t = t n + ∆ t can be updated by, W ∗ = W n + 12 ∆ t L ( W n ) + 18 ∆ t ∂∂t L ( W n ) . (17) W n +1 = W n + ∆ t L ( W n ) + 16 ∆ t (cid:0) ∂∂t L ( W n ) + 2 ∂∂t L ( W ∗ ) (cid:1) . (18)The detailed proof is given in [20].The numerical fluxes and their time derivatives in the above equations, such as L ( W ni )8nd ∂∂t L ( W ni ), are determined by L ( W ni,j ) = − x [( F ) i +1 / ,j ( W n , t n ) − ( F ) i − / ,j ( W n , t n )] − y [( G ) i,j +1 / ( W n , t n ) − ( G ) i,j − / ( W n , t n )] , L t ( W ni,j ) = − x [ ∂ t ( F ) i +1 / ,j ( W n , t n ) − ∂ t ( F ) i − / ,j ( W n , t n )] − y [ ∂ t ( G ) i,j +1 / ( W n , t n ) − ∂ t ( G ) i,j − / ( W n , t n )] . (19)Similarly, the time derivatives for the intermediate state can be obtained, L t ( W ∗ i,j ) = − x [ ∂ t ( F ) i +1 / ,j ( W ∗ , t ∗ ) − ∂ t ( F ) i − / ,j ( W ∗ , t ∗ )] − y [ ∂ t ( G ) i,j +1 / ( W ∗ , t ∗ ) − ∂ t ( G ) i,j − / ( W ∗ , t ∗ )] . (20)In the gas-kinetic scheme, the flux Eq.(8) is a complicated function of time. To obtain thetime derivatives of the flux function used in the above two-stage fourth-order framework,the flux function is approximated as a linear function of time within a time interval. Thetime-dependent flux can be expanded as, F i +1 / ,j ( W n , t ) = F ni +1 / ,j + ∂ t F ni +1 / ,j ( t − t n ) , t ∈ [ t n , t n + ∆ t ] . (21)To get coefficients F ni +1 / ,j and ∂ t F ni +1 / ,j , the following notation of Eq.(8) is introduced, F i +1 / ,j ( W n , δ ) = (cid:90) t n + δt n F i +1 / ,j ( W n , t )d t = (cid:88) (cid:96) =1 ω (cid:96) (cid:90) t n + δt n (cid:90) u ψ f ( x i +1 / ,j (cid:96) , t, u, v, ξ )dΞd t. Take δ as ∆ t and ∆ t/
2, we have, F i +1 / ,j ( W n , t n )∆ t + 12 ∂ t F i +1 / ,j ( W n , t n )∆ t = F i +1 / ,j ( W n , ∆ t ) , F i +1 / ,j ( W n , t n )∆ t + 18 ∂ t F i +1 / ,j ( W n , t n )∆ t = F i +1 / ,j ( W n , ∆ t/ . Solving the above linear system, and we can obtain the expression of coefficients F i +1 / ,j ( W n , t n ) = (4 F i +1 / ,j ( W n , ∆ t/ − F i +1 / ,j ( W n , ∆ t )) / ∆ t,∂ t F i +1 / ,j ( W n , t n ) = 4( F i +1 / ,j ( W n , ∆ t ) − F i +1 / ,j ( W n , ∆ t/ / ∆ t . Similarly, the coefficients for the intermediate state F i +1 / ,j ( W ∗ , t ∗ ), ∂ t F i +1 / ,j ( W ∗ , t ∗ ) canbe determined as well. The fluxes in y-direction can be obtained through the same method.Then, the intermediate states W ∗ ij are updated by Eq.(17) and Eq.(19). The final states9 n +1 ij in Eq.(18) are determined through Eq.(19) and Eq.(20). The HLLC Riemann solver [36] is used to obtain the inviscid flux in the current WENO5-AO-HLLC scheme. Consider the following Riemann problem, W t + F x ( W ) = 0 , with the initial condition, W ( x,
0) = (cid:40) W L , x < , W R , x > , where W L and W R are the initial interface values. For the two-dimensional Euler equations,the conservative variables W and the corresponding fluxes F are, W = [ ρ, ρU, ρV, ρE ] T , F = (cid:2) ρU, ρU + p, ρU V, U ( ρE + p ) (cid:3) T . HLLC solver is an approximate Riemann solver, which consists of four constant states.Assume that the speeds of the slowest and fastest wave are S L and S R , and the speed of themiddle shear wave is S ∗ . Then, the HLLC solver can be written as follows, W ( x, t ) = W L , xt ≤ S L , W ∗ L , S L ≤ xt ≤ S ∗ , W ∗ R , S ∗ ≤ xt ≤ S R , W R , xt ≥ S R , (22)and the corresponding numerical flux can be defined as, F x +1 / = F L , ≤ S L , F ∗ L , S L ≤ ≤ S ∗ , F ∗ R , S ∗ ≤ ≤ S R , F R , ≥ S R , (23)where F ∗ K = F K + S L ( W ∗ K − W K ), K = L, R . The W ∗ K , K = L, R , is given by, W ∗ K = ρ K (cid:18) S K − U K S K − S ∗ (cid:19) S ∗ V KE K ρ K + ( S ∗ − U K ) (cid:104) S ∗ + p K ρ K ( S K − U K ) (cid:105) , (24)10here S ∗ is related to the speeds S L and S R , namely S ∗ = p R − p L + ρ L U L ( S L − U L ) − ρ R U R ( S R − U R ) ρ L ( S L − U L ) − ρ R ( S R − U R ) . There are many methods to estimate wave speeds S L and S R , and a pressure-based wavespeed estimate method proposed by Toro is adopted in the current work [34]. Firstly, weneed to estimate p ∗ , the pressure of the region x/t ∈ [ S L , S R ]. Based on the Two-RarefactionRiemann solver (TRRS), the estimated p ∗ is p ∗ = (cid:34) a L + a R − γ − ( U R − U L ) a L /p zL + a R /p zR (cid:35) /z where z = ( γ − / (2 γ ), and γ is the specific heat ratio. Then, the speeds S L and S R arecoming from the exact wave-speed relations in the exact Riemann solver, S L = U L − a L q L , S R = U R − a R q R , where a L , a R are the sound speeds of initial left and right state, and q K , K = L, R , are q K = , p ∗ ≤ p K , (cid:20) γ + 12 γ ( p ∗ /p K − (cid:21) / , p ∗ > p K . For viscous flow problems, the viscous fluxes in Navier-Stokes equations are needed. Tocalculate the viscous fluxes in the current WENO5-AO-HLLC scheme, both the conservativevariables Q i +1 / and the corresponding derivatives ( Q x ) i +1 / at the cell interface need to beconstructed by the cell averaged conservative variables Q . In this paper, a sixth-order centraldifference method is applied for the calculation of viscous fluxes. The conservative variablescan be written as follows, Q i +1 / = 160 ( Q i − − Q i − + 37 Q i + 37 Q i +1 − Q i +2 + Q i +3 ) , and the corresponding derivatives are,( Q x ) i +1 / = 1180∆ x ( − Q i − + 25 Q i − − Q i + 245 Q i +1 − Q i +2 + 2 Q i +3 ) . For two-dimensional problems, the dimension-by-dimension strategy is adopted [42, 15].The reconstructed value Q i +1 / ,j l at the Gaussian quadrature point j l , the correspondingnormal derivative ( Q x ) i +1 / ,j l , and tangential derivative ( Q y ) i +1 / ,j l can be obtained by thefourth-order polynomial p r ( y ) based on the above Q i +1 / and ( Q x ) i +1 / . Then, all terms in11he viscous fluxes can be fully determined. A similar procedure can be easily extended tothree-dimensional problems. To improve the robustness of WENO5-AO-HLLC scheme, theconservative variables at the cell interface Q i +1 / are obtained by simple averaging of theleft and right interface values of WENO5-AO reconstruction in some challenging cases. Considering the fourth-order temporal accuracy in WENO5-AO-GKS scheme, the clas-sical fourth-order Runge-Kutta method (RK4) is adopted for time integration in WENO5-AO-HLLC scheme for achieving the 4th-order temporal accuracy. The RK4 time marchingmethod reads, W = W n + 12 ∆ t L ( W n ) , W = W n + 12 ∆ t L ( W ) , W = W n + ∆ t L ( W ) , W n +1 = W n + 16 (cid:0) ∆ t L ( W n ) + 2∆ t L ( W ) + 2∆ t L ( W ) + ∆ t L ( W ) (cid:1) , with L defined in Eq.(19).
3. Numerical performance
In the following test cases, for the inviscid flow the time step is determined by,∆ t = CFL × ∆ x ( | U | + C ) Max , where C is sound speed. For viscous flow, the time step is given by,∆ t = CFL × Min (cid:20) ∆ x ( | U | + C ) Max , ρ ∆ x µ (cid:21) . Titarev-Toro problem is an inviscid flow problem with a shock wave impinging into ahigh-frequency density perturbation [32]. This problem consists of a main shock, a highgradient smooth post-shock region and multiple shocklets developed later. To representthese flow structures, a high-order scheme is needed. The initial condition is given by,( ρ, U, p ) = (cid:40) (1 . , . , . , − . ≤ x ≤ − . , (1 + 0 . πx ) , . , . , − . < x ≤ . . The computational domain is [ − ,
5] with a mesh of 1000 cells. Two CFL numbers, 0.5 and1.0, are employed for both WENO5-AO-GKS and WENO5-AO-HLLC, and the results at12he output time t = 5 . x d e n s i t y
4 2 0 2 40.811.21.41.61.8
ReferenceWENO5 AO GKSWENO5 AO HLLC x d e n s i t y ReferenceWENO5 AO GKSWENO5 AO HLLC
Figure 1: Titarev-Toro problem by WENO5-AO-GKS scheme and WENO5-AO-HLLC scheme. Densitydistribution with mesh number 1000 at t = 5 .
0. Left figure shows the whole domain; right figure shows theenlarged domain. The CFL number is 0.5. x d e n s i t y
4 2 0 2 40.811.21.41.61.8
ReferenceWENO5 AO GKSWENO5 AO HLLC x d e n s i t y ReferenceWENO5 AO GKSWENO5 AO HLLC
Figure 2: Titarev-Toro problem by WENO5-AO-GKS scheme and WENO5-AO-HLLC scheme. Densitydistribution with mesh number 1000 at t = 5 .
0. Left figure shows the whole domain; right figure shows theenlarged domain. The CFL number is 1.
Le Blanc problem is a class of 1-D Riemann problems with initially high ratios for densityand pressure [31]. Therefore, an extremely strong rarefaction wave is generated in the high-13ressure region. The initial condition here is chosen as,( ρ, U, p ) = (cid:40) (10 M , , M ) , ≤ x ≤ . , (1 , , , . < x ≤ . Here Le Blanc problem with initial pressure ratio 10 and 10 was calculated by WENO5-AO-GKS scheme and WENO5-AO-HLLC scheme, and the profiles of density, temperature,and pressure at t = 0 .
12 are presented in Figure 3. For this case, CFL number is 0.5. Forboth two schemes, there exist discrepancy in the vicinity of the shock wave, which has alsobeen observed in the previous research, especially in the coarse mesh case [31]. Both schemespresent a similar performance in this case.
Noh problem consists of two strong shocks moving from center to left and right siderespectively [25]. The initial condition is as follows,( ρ, U, p ) = (cid:40) (1 , , − ) , ≤ x ≤ . , (1 , − , − ) , . < x ≤ . The computational domain is [0 , γ = 5 /
3. The output time is t = 1 .
0. The results are presented inFigure 4. It is worth noting that WENO5-AO-HLLC scheme blows up for this problem whileWENO5-AO-GKS scheme can work well. Besides, as a comparison, the WENO5-AO-LFscheme is adopted for this problem. The WENO5-AO-LF scheme means that, only HLLCsolver in WENO5-AO-HLLC scheme is replaced by Lax-Friedrich solver. The results showthat both two schemes can resolve the shock very well, although the density profiles exista weak dip at the central region for both schemes. Besides, the result of WENO5-AO-GKSscheme shows a weaker dip.
The forward step problem proposed by Woodward and Colella [37] is an inviscid testcase. A uniform flow with
M a = 3 blows towards a wind tunnel containing a step. Thiswind tunnel size, is [0 , × [0 , ρ, U, V, p ) = (1 , , , /γ ) , where γ = 1 .
4. The supersonic inlet and outlet boundary condition is employed for theleft and right boundary respectively, while other boundaries are set as reflective boundaryconditions. It is worth remarking that, the ghost cells near the corner of the step [0 . , . U is given by the value obtained through applying the reflec-tive boundary condition for the upper flow region; velocity V is given by the value obtained14 d e n s i t y ReferenceWENO5 AO GKSWENO5 AO HLLC x d e n s i t y ReferenceWENO5 AO GKSWENO5 AO HLLC x t e m p e r a t u r e ReferenceWENO5 AO GKSWENO5 AO HLLC x t e m p e r a t u r e ReferenceWENO5 AO GKSWENO5 AO HLLC x ve l o c i t y ReferenceWENO5 AO GKSWENO5 AO HLLC x ve l o c i t y ReferenceWENO5 AO GKSWENO5 AO HLLC
Figure 3: Le Blanc problem with initial pressure ratio 10 (left three figures) and 10 (right three figures)by WENO5-AO-GKS scheme and WENO5-AO-HLLC scheme. For all figures, CFL number is 0.5, the meshnumber is 200 and the output time is t = 0 . d e n s i t y ReferenceWENO5 AO GKSWENO5 AO LF x d e n s i t y ReferenceWENO5 AO GKSWENO5 AO LF x p r ess u r e ReferenceWENO5 AO GKSWENO5 AO LF x ve l o c i t y ReferenceWENO5 AO GKSWENO5 AO LF
Figure 4: Noh problem by WENO5-AO-GKS scheme and WENO5-AO-LF scheme. The density, pressure,and velocity profiles respectively with CFL number 0.5 are shown. For all figures, the mesh number is 400and the output time is t = 1 .
0. The WENO-AO-HLLC fails for this test case. ρ and pres-sure p are given by algebraically averaging the corresponding values obtained through thereflective boundary condition for the upper and left flow regions. The CFL number 0 . x = ∆ y = 1 / , / , / t = 4 .
0. Theresults show that both WENO5-AO-GKS scheme and WENO5-AO-HLLC scheme performwell when adopting a fine mesh. In the top region, both the triple-point structure and thevortex sheet can be captured clearly. x y x y x y x y x y x y Figure 5: Mach 3 forward step problem by WENO5-AO-GKS scheme (left) and WENO5-AO-HLLC scheme(right). Density distribution with different mesh size at t = 4 .
0. The mesh size of the top figure, middlefigure, and bottom figure are 1/120, 1/240, and 1/360. The CFL number is 0.8 for both WENO-AO-GKSand WENO5-AO-HLLC. 30 equally spaced contours from 0.2 to 4.7 are plotted.
Laminar boundary layer is a standard test case for viscous flow [39]. A plane withthe characteristic length L = 100 is placed from 0 to 100. The computation domain is[ − , × [0 , ×
32. At the start point of plane, the minimal cell mesh ∆ x and ∆ y are 0.1and 0.12 separately. The inlet flow is described by,( ρ, U, V, p ) = (1 , . , , /γ ) , where γ = 1 .
4. In the case, kinematic viscosity coefficient is ν = 1 . × − , and thus Re = U ∞ L/ν = 1 . × and M a = 0 .
15. Besides, the adiabatic non-slip boundary17ondition is adopted on the plate, while the symmetric slip boundary condition is used forthe bottom boundary of [ − , ys = y √ Re/x ,and the non-dimensional velocity us = U/U ∞ , vs = V √ Re x /U ∞ , respectively. The valuesin the legend represent the location x/L . From the results, both WENO5-AO-GKS andWENO5-AO-HLLC are capable of capturing the velocity profile well in the boundary layerwith several mesh cells. Close to the leading edge, WENO5-AO-GKS gives a slightly better vs solution than WENO5-AO-HLLC at the location x/L = 0 . Figure 6: Mesh with 120 ×
32 cells for laminar boundary layer case.
Double shear layer is a viscous problem involving a pair of doubly-periodic shear layers[5]. When the numerical method is not enough to resolve the flow field, non-physical vortexeswill appear in the evolution stage. The “thin” shear layer problem is studied in [5], and theinitial U velocity is given by, U = (cid:40) tanh ( k ( y − . , ≤ y ≤ . , tanh ( k (0 . − y )) , . < y ≤ , and the initial V velocity, density, and pressure are given as follows, V = δ sin (2 πx ) , ρ = 1 , p = ρU M a γ , s u s BlasiusWENO5 AO GKS 0.050WENO5 AO GKS 0.105WENO5 AO GKS 0.165WENO5 AO GKS 0.265 ys vs BlasiusWENO5 AO GKS 0.050WENO5 AO GKS 0.105WENO5 AO GKS 0.165WENO5 AO GKS 0.265 ys u s BlasiusWENO5 AO HLLC 0.050WENO5 AO HLLC 0.105WENO5 AO HLLC 0.165WENO5 AO HLLC 0.265 ys vs BlasiusWENO5 AO HLLC 0.050WENO5 AO HLLC 0.105WENO5 AO HLLC 0.165WENO5 AO HLLC 0.265
Figure 7: Laminar boundary layer by WENO5-AO-GKS scheme (top two) and WENO5-AO-HLLC scheme(bottom two). For all figures, Re = 1 . × , M a = 0 .
15, CFL number is 0.5, and the mesh number is120 × k = 100, the perturbation size δ = 0 .
05, the Machnumber
M a = 0 .
15, and the specific heat ratio γ = 1 .
4. Besides, the kinetic viscosity is ν =5 . × − . Periodic boundary condition is employed for all boundaries. The computationaldomain is [0 , × [0 , ×
256 in this case. Linear reconstruction isemployed for both schemes in this test.The vorticity contours Ω = ( ∂V∂x − ∂U∂y ) at t = 0 . . , . x y x y Figure 8: Two-dimensional double shear flow by WENO5-AO-GKS scheme (left) and WENO5-AO-HLLCscheme (right) : vorticity. The CFL number is 0.5, the output time is t = 0 .
8, and the mesh number is256 × Viscous shock tube problem is a viscous flow problem with a strong shock [10]. Theinteraction of reflected shock from the right wall and the viscous boundary layer producesa series of complex flow structures, such as the typical λ − shape shock configuration. Theinitial condition is given by,( ρ, U, V, p ) = (cid:40) (120 , , , /γ ) , < x ≤ . , (1 . , , , . /γ ) , . < x < , where γ = 1 . P r = 0 .
73. The computational domain is [0 , × [0 , . Re = 200 is tested. The output time is t = 1 .
0. For the boundary condition, theupper boundary is asymmetric boundary, and others are the non-slip adiabatic wall. Thedensity contours at Re = 200 are shown in Figure 9, where both WENO5-AO-GKS and20ENO5-AO-HLLC can capture the main flow structures. The λ − shape structure, thevortices within the boundary layer, and the slip line in the lower right region are capturedclearly. The density profiles along the bottom wall are shown in Figure 10 with the localenlargement. The results show that both schemes have similar resolution. The results ofWENO5-AO-GKS on both coarse and fine meshes seem to be closer than those of WENO5-AO-HLLC. x y x y x y x y Figure 9: Viscous shock tube problem with Re = 200 by WENO5-AO-GKS scheme (left) and WENO5-AO-HLLC scheme (right): density distribution. For all cases, the CFL number is 0.2. For the top two figures,the mesh number is 500 × × The three-dimensional advection of density perturbation is adopted for accuracy test.The initial condition is, ρ ( x, y, z ) = 1 + 0 . π ( x + y + z )) ,U ( x, y, z ) = 1 , V ( x, y, z ) = 1 , W ( x, y, z ) = 1 , p ( x, y, z ) = 1 . The computational domain covers [0 , × [0 , × [0 , ρ ( x, y, z, t ) = 1 + 0 . π ( x + y + z − t )) ,U ( x, y, z, t ) = 1 , V ( x, y, z, t ) = 1 , W ( x, y, z, t ) = 1 , p ( x, y, z, t ) = 1 . The L error and convergence order of WENO5-AO-GKS and WENO5-AO-HLLC at t = 2 . d e n s i t y WENO5 AO GKS 500250WENO5 AO GKS 1000500WENO5 AO HLLC 500250WENO5 AO HLLC 1000500 x d e n s i t y WENO5 AO GKS 500250WENO5 AO GKS 1000500WENO5 AO HLLC 500250WENO5 AO HLLC 1000500
Figure 10: Viscous shock tube problem of Re = 200 by WENO5-AO-GKS scheme and WENO5-AO-HLLCscheme: density profile along the bottom wall ( y = 0). The right is the enlarged figure. For all cases, theCFL number is 0.2. CFL 0.20 0.60 1.00Mesh L Error Order L Error Order L Error Order5 × × × ×
10 2.234252e-03 4.36 1.665667e-03 4.48 3.741869e-03 3.8720 × ×
20 7.589204e-05 4.88 6.007786e-05 4.79 2.052167e-04 4.1940 × ×
40 2.470600e-06 4.94 2.640903e-06 4.51 1.220770e-05 4.0780 × ×
80 8.596794e-08 4.84 1.424737e-07 4.21 7.525234e-07 4.02
Table 1: 3-D accuracy test: L Error and convergence order by WENO5-AO-GKS scheme with differentCFL numbers.
CFL 0.20 0.60 1.00Mesh L Error Order L Error Order L Error Order5 × × × ×
10 3.390217e-03 4.14 3.658550e-03 4.08 5.859740e-03 3.7220 × ×
20 1.209088e-04 4.81 1.312929e-04 4.80 2.588044e-04 4.5040 × ×
40 3.900971e-06 4.95 4.563747e-06 4.85 1.318523e-05 4.2980 × ×
80 1.283506e-07 4.93 1.845443e-07 4.63 7.688684e-07 4.10
Table 2: 3-D accuracy test: L Error and convergence order by WENO5-AO-HLLC scheme with differentCFL numbers. .3.2. Compressible isotropic turbulence The decaying compressible isotropic turbulence is a case to evaluate the robustness ofdifferent schemes [29, 6]. The definitions of flow variables are introduced first. The turbulentfluctuating velocity U (cid:48) is, U (cid:48) = (cid:28) U + U + U (cid:29) / , where (cid:104)· · · (cid:105) means the space average over the whole computation domain. Then, turbulenceMach number M a t is given by, M a t = (cid:104) U + U + U (cid:105) / C = √ U (cid:48) C , where C is the local sound speed. Taylor microscale λ is defined by, λ = ( U (cid:48) ) (cid:10) ( ∂U /∂x ) (cid:11) , and the corresponding Taylor Reynolds number Re λ is Re λ = ρU (cid:48) λµ , where µ is the dynamic viscosity coefficient determined by µ = µ (cid:18) TT (cid:19) . . In this case, the velocity spectrum is given by, E ( k ) = A k e ( − k /k ) , where A is the initial kinetic energy, k is the wave number, and k is the peak value of k .The initial turbulent kinetic energy K and the initial large-eddy-turnover time τ can beobtained as follows, K = 3 A √ πk ,τ = (cid:114) A (2 π ) / k − / . K ( t ) and root-mean-square of density fluctuation ρ rms ( t ) are defined as K ( t ) = (cid:104) ρU + ρU + ρU (cid:105) ,ρ rms ( t ) = (cid:113)(cid:10) ( ρ − ρ ) (cid:11) . In this case, there are strong shocklets and shock-vortex interactions in the flow field, espe-cially at a high turbulence Mach number
M a t . Therefore, it is challenging for high-orderscheme to simulate high M a t flow. The simulations will cover a wide range of M a t to com-pare the robustness of WENO5-AO-GKS and WENO5-AO-HLLC. The mesh adopted inthis case is 128 . Other parameters take the values Re λ = 72, A = 1 . × − , and k = 8 . K ( t ) /K and root-mean-square of densityfluctuation ρ rms ( t ) /M a t are shown in Figure 11. Both WENO5-AO-GKS and WENO5-AO-HLLC perform well for a wide range of M a t from 0.5 to 1.4. The reference data of M a t = 0 . M a t = 1 .
4, iso-surface of the second invariant of velocity gradienttensor Q = 25 colored by the local Mach number at t/τ = 1 . . M a t up to 0.6, which is much smallerthan 1 .
4. In the simulations, to improve the robustness of WENO5-AO-HLLC scheme, theconservative variables at the cell interface Q i +1 / are obtained by simple averaging of theleft and right interface values of WENO5-AO reconstruction. The above results show thatWENO5-AO-GKS is more robust than WENO5-AO-HLLC in this case. The three-dimensional Taylor-Green vortex is studied by WENO5-AO-GKS and WENO5-AO-HLLC. The computational domain is [ − πL, πL ] × [ − πL, πL ] × [ − πL, πL ], and the initialcondition is U = U sin ( x/L ) cos ( y/L ) cos ( z/L ) ,V = − U cos ( x/L ) sin ( y/L ) cos ( z/L ) ,W = 0 ,p = p + ρ U (cos (2 x/L ) + cos (2 y/L )) (cos (2 z/L ) + 2) / . The simulation has L = 1, U = 1, ρ = 1, and the Reynolds number Re = U L/ν = 280.The Mach number is
M a = U /C =0.1 and the sound speed is C = √ γRT . The meshnumber is 64 , and periodic boundary condition is imposed at all boundaries. The volume-averaged kinetic energy is defined as, E k = 1 ρ Ω (cid:90) Ω ρ ( U + V + W )2 dΩ , / τ ρ r m s / M a t WENO5 AO GKS 0.5WENO5 AO GKS 0.8WENO5 AO GKS 1.0WENO5 AO GKS 1.2WENO5 AO GKS 1.4WENO5 AO HLLC 0.5WENO5 AO HLLC 0.8WENO5 AO HLLC 1.0WENO5 AO HLLC 1.2WENO5 AO HLLC 1.4Reference 0.5 t/ τ K (t) / K WENO5 AO GKS 0.5WENO5 AO GKS 0.8WENO5 AO GKS 1.0WENO5 AO GKS 1.2WENO5 AO GKS 1.4WENO5 AO HLLC 0.5WENO5 AO HLLC 0.8WENO5 AO HLLC 1.0WENO5 AO HLLC 1.2WENO5 AO HLLC 1.4Reference 0.5
Figure 11: Compressible isotropic turbulence at different
M a t by WENO5-AO-GKS scheme and WENO5-AO-HLLC scheme. The normalized kinetic energy (left) and normalized root-mean-square of density fluctu-ation (right). The CFL number is 0.3 for both WENO5-AO-GKS scheme and WENO5-AO-HLLC scheme.For all cases, the mesh number is 128 .Figure 12: Compressible isotropic turbulence with M a t = 1 .
4: iso-surface of the second invariant of velocitygradient tensor Q = 25 colored with local Mach number by WENO5-AO-GKS scheme (left) and WENO5-AO-HLLC scheme (right). The mesh number is 128 and output time is t/τ = 1 . (cid:15) k = − d E k d t . The linear reconstruction is taken for both schemes in this test case. The results arepresented in Figure 13, and are compared with the reference solution of [27]. The CFLnumber is 0 . . .
4, WENO-AO-HLLC will generate large oscillation. t k i n e t i c e n e r g y ReferenceWENO5 AO GKS(CFL0.5)WENO5 AO HLLC(CFL0.3)WENO5 AO HLLC(CFL0.4) t d k ReferenceWENO5 AO GKS(CFL0.5)WENO5 AO HLLC(CFL0.3)WENO5 AO HLLC(CFL0.4)
Figure 13: Taylor-Green vortex problem with Re = 280 by WENO5-AO-GKS scheme and WENO5-AO-HLLC scheme: the kinetic energy (left) and the dissipation rate (right). The CFL number is 0.5 forWENO5-AO-GKS scheme and 0.3 for WENO5-AO-HLLC scheme. For both cases, the mesh number is 64 .
4. Computational efficiency
The computational efficiency of WENO5-AO-GKS and WENO5-AO-HLLC is comparedin 2-D and 3-D cases. For both schemes, the main computational cost includes two parts,reconstruction, and evolution. For the reconstruction, WENO5-AO-HLLC needs only point-wise conservative variables, while the derivatives are also needed in WENO-AO-GKS. How-ever, additional reconstruction through central difference method for the viscous terms isrequired in WENO5-AO-HLLC. For the evolution stage, the GKS flux is more expensivethan HLLC, but GKS uses two stages instead of four stages in HLLC to achieve 4th-ordertime accuracy.The viscous shock tube is used to test the computational efficiency. The mesh pointsin the test are 1000 × CPU time ( s ) Time ratioWENO5-AO-GKS 154.91 1.00WENO5-AO-HLLC 196.47 1.27 Table 3: 2-D computational efficiency test of viscous shock tube problem. The mesh number is 1000 × CPU time ( s ) Time ratioWENO5-AO-GKS 476.04 1.00WENO5-AO-HLLC 403.03 0.85 Table 4: 3-D computational efficiency test of compressible isotropic turbulence problem with
M a t = 0 . .The shown CPU time is obtained for 10 time steps by a single Intelcore i7-9700 @ 3.00GHz.
5. Conclusion
A comparison of performance for two high-order schemes, namely WENO5-AO-GKS andWENO5-AO-HLLC, is presented. Both schemes use the fifth-order WENO-AO reconstruc-tion, the differences are mainly coming from the flux functions and the temporal updatingschemes. In GKS, due to the time accurate flux and its time derivative the multistageand multiderivative (MSMD) is used to update the solution. The two-stage fourth-ordertemporal discretization achieves a 4th-order temporal accuracy. For HLLC, four stagesRunge-Kutta method is used for the time accuracy. In WENO-AO-GKS, both inviscid andviscous flux terms can be evaluated from a single time-dependent gas distribution function.In WENO5-AO-HLLC, HLLC provides inviscid flux and a sixth-order central differencemethod is used to discretize the viscous flux. In the 3D accuracy test, both schemes canachieve the expected order of accuracy, and WENO5-AO-GKS shows a slightly smaller ab-solute L error. In terms of the shock and contact wave capturing, both schemes performwell and have similar robustness. With the same mesh and CFL number, WENO5-AO-GKSshows better accuracy in the double shear layer test. In the Noh problem, WENO5-AO-GKS presents favorable robustness. For the compressible isotropic turbulence and three-dimensional Taylor-Green vortex problem, WENO-AO-GKS can take a large time step27ith CFL number 0.5, instead of 0.3 for WENO5-AO-HLLC. For two-dimensional viscousshock tube problems, WENO5-AO-HLLC is (27%) more expensive than WENO-AO-GKS.While for the three-dimensional viscous test, WENO5-AO-HLLC is (15%) more efficientthan WENO5-AO-GKS. WENO-AO-GKS requires the reconstruction of flow variables inthe normal and two tangential directions on both sides of a cell interface in the 3D case.The multi-dimensional property and the coupling of inviscid and viscous fluxes in WENO5-AO-GKS have obvious advantages when the scheme is extended to the flow computationwith unstructured mesh. References [1] Dinshaw S Balsara, Sudip Garain, and Chi-Wang Shu. An efficient class of WENO schemes withadaptive order.
Journal of Computational Physics , 326:780–804, 2016.[2] P Batten, MA Leschziner, and UC Goldberg. Average-state Jacobians and implicit methods for com-pressible viscous and turbulent flows.
Journal of computational physics , 137(1):38–78, 1997.[3] Prabhu Lal Bhatnagar, Eugene P Gross, and Max Krook. A model for collision processes in gases I:Small amplitude processes in charged and neutral one-component systems.
Physical Review , 94(3):511–525, 1954.[4] Rafael Borges, Monique Carmona, Bruno Costa, and Wai Sun Don. An improved weighted essen-tially non-oscillatory scheme for hyperbolic conservation laws.
Journal of Computational Physics ,227(6):3191–3211, 2008.[5] David L Brown. Performance of under-resolved two-dimensional incompressible flow simulations.
Jour-nal of Computational Physics , 122(1):165–183, 1995.[6] Guiyu Cao, Liang Pan, and Kun Xu. Three dimensional high-order gas-kinetic scheme for supersonicisotropic turbulence I: criterion for direct numerical simulation.
Computers & Fluids , 192(104273),2019.[7] Guiyu Cao, Hongmin Su, Jinxiu Xu, and Kun Xu. Implicit high-order gas kinetic scheme for turbulencesimulation.
Aerospace Science and Technology , 92:958–971, 2019.[8] Sydney Chapman and Thomas George Cowling.
The mathematical theory of non-uniform gases: an ac-count of the kinetic theory of viscosity, thermal conduction and diffusion in gases . Cambridge universitypress, 1970.[9] Bing Chen, Yan Zhang, and Xu Xu. Numerical simulation of supersonic turbulent combustion flowsbased on flamelet model.
Journal of Propulsion Technology , 12:11, 2013.[10] Virginie Daru and Christian Tenaud. Evaluation of TVD high resolution schemes for unsteady viscousshocked flows.
Computers & Fluids , 30(1):89–113, 2000.[11] SK Godunov. A finite difference method for the computation of discontinuous solutions of the equationsof fluid dynamics.
Sbornik: Mathematics , 47(8-9):357–393, 1959.[12] Ami Harten, Peter D Lax, and Bram van Leer. On upstream differencing and godunov-type schemesfor hyperbolic conservation laws.
SIAM review , 25(1):35–61, 1983.[13] Ami Harten, Stanley Osher, Bj¨orn Engquist, and Sukumar R Chakravarthy. Some results on uniformlyhigh-order accurate essentially nonoscillatory schemes.
Applied Numerical Mathematics , 2(3-5):347–377,1986.[14] Xing Ji, Liang Pan, Wei Shyy, and Kun Xu. A compact fourth-order gas-kinetic scheme for the Eulerand Navier-Stokes equations.
Journal of Computational Physics , 372:446 – 472, 2018.[15] Xing Ji and Kun Xu. Performance enhancement for high-order gas-kinetic scheme based on weno-adaptive-order reconstruction.
Commun. Comput. Phys. , 28:539–590, 2020.[16] Xing Ji, Fengxiang Zhao, Wei Shyy, and Kun Xu. A family of high-order gas-kinetic schemes andits comparison with Riemann solver based high-order methods.
Journal of Computational Physics ,356:150–173, 2018.
17] Xing Ji, Fengxiang Zhao, Wei Shyy, and Kun Xu. A HWENO Reconstruction Based High-orderCompact Gas-kinetic Scheme on Unstructured Mesh.
Journal of Computational Physics , 109367, 2020.[18] Guang-Shan Jiang and Chi-Wang Shu. Efficient implementation of weighted ENO schemes.
Journal ofcomputational physics , 126(1):202–228, 1996.[19] Doron Levy, Gabriella Puppo, and Giovanni Russo. Central WENO schemes for hyperbolic systems ofconservation laws.
ESAIM: Mathematical Modelling and Numerical Analysis , 33(3):547–571, 1999.[20] Jiequan Li and Zhifang Du. A two-stage fourth order time-accurate discretization for Lax–Wendrofftype flow solvers I. hyperbolic conservation laws.
SIAM Journal on Scientific Computing , 38(5):A3046–A3069, 2016.[21] Qibing Li, Song Fu, and Kun Xu. Application of gas-kinetic scheme with kinetic boundary conditionsin hypersonic flow.
AIAA journal , 43(10):2170–2176, 2005.[22] Wei Liao, Yan Peng, and Li-Shi Luo. Gas-kinetic schemes for direct numerical simulations of compress-ible homogeneous turbulence.
Physical Review E , 80(4):046702, 2009.[23] Xu-Dong Liu, Stanley Osher, Tony Chan, et al. Weighted essentially non-oscillatory schemes.
Journalof computational physics , 115(1):200–212, 1994.[24] Jun Luo and Kun Xu. A high-order multidimensional gas-kinetic scheme for hydrodynamic equations.
Sci. China, Technol. Sci , 56(10):2370–2384, 2013.[25] William F Noh. Errors for calculations of strong shocks using an artificial viscosity and an artificialheat flux.
Journal of Computational Physics , 72(1):78–120, 1987.[26] Liang Pan, Junxia Cheng, Shuanghu Wang, and Kun Xu. A two-stage fourth-order gas-kinetic schemefor compressible multicomponent flows.
Communications in Computational Physics , 22(4):1123–1149,2017.[27] Liang Pan and Kun Xu. Two-stage fourth-order gas-kinetic scheme for three-dimensional Euler andNavier-Stokes solutions.
International Journal of Computational Fluid Dynamics , 32(10):395–411, 2018.[28] Liang Pan, Kun Xu, Qibing Li, and Jiequan Li. An efficient and accurate two-stage fourth-ordergas-kinetic scheme for the Euler and Navier–Stokes equations.
Journal of Computational Physics ,326:197–221, 2016.[29] Ravi Samtaney, Dale I Pullin, and Branko Kosovi´c. Direct numerical simulation of decaying compress-ible turbulence and shocklet statistics.
Physics of Fluids , 13(5):1415–1430, 2001.[30] Shuang Tan, Qibing Li, Zhixiang Xiao, and Song Fu. Gas kinetic scheme for turbulence simulation.
Aerospace Science and Technology , 78:214–227, 2018.[31] Huazhong Tang and Tiegang Liu. A note on the conservative schemes for the Euler equations.
Journalof Computational Physics , 218:451–459, 2006.[32] Vladimir A Titarev and Eleuterio F Toro. Finite-volume WENO schemes for three-dimensional con-servation laws.
Journal of Computational Physics , 201(1):238–260, 2004.[33] SA Tokareva and Eleuterio F Toro. HLLC-type Riemann solver for the Baer–Nunziato equations ofcompressible two-phase flow.
Journal of Computational Physics , 229(10):3573–3604, 2010.[34] Eleuterio F Toro.
Riemann solvers and numerical methods for fluid dynamics: a practical introduction .Springer Science & Business Media, 2013.[35] Eleuterio F Toro. The HLLC Riemann solver.
Shock Waves , 29:1065–1082, 2019.[36] Eleuterio F Toro, Michael Spruce, and William Speares. Restoration of the contact surface in theHLL-Riemann solver.
Shock waves , 4(1):25–34, 1994.[37] Paul Woodward and Phillip Colella. The numerical simulation of two-dimensional fluid flow with strongshocks.
Journal of computational physics , 54(1):115–173, 1984.[38] Kun Xu. BGK-based scheme for multicomponent flow calculations.
Journal of Computational Physics ,134(1):122–133, 1997.[39] Kun Xu. Gas-kinetic schemes for unsteady compressible flow simulations.
Lecture series-van KaremanInstitute for fluid dynamics , 3:C1–C202, 1998.[40] Kun Xu. A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificialdissipation and Godunov method.
Journal of Computational Physics , 171(1):289–335, 2001.[41] Kun Xu and Juan-Chen Huang. A unified gas-kinetic scheme for continuum and rarefied flows.
Journal f Computational Physics , 229(20):7747–7764, 2010.[42] Rui Zhang, Mengping Zhang, and Chi-Wang Shu. On the order of accuracy and numerical performanceof two classes of finite volume WENO schemes. Communications in Computational Physics , 9(3):807–827, 2011.[43] Fengxiang Zhao, Xing Ji, Wei Shyy, and Kun Xu. An acoustic and shock wave capturing compacthigh-order gas-kinetic scheme with spectral-like resolution. arXiv preprint arXiv:2001.01570 , 2019.[44] Fengxiang Zhao, Xing Ji, Wei Shyy, and Kun Xu. Compact higher-order gas-kinetic schemes withspectral-like resolution for compressible flow simulations.
Advances in Aerodynamics , 1(1):13, 2019.[45] Jun Zhu and Chi-Wang Shu. A new type of multi-resolution WENO schemes with increasingly higherorder of accuracy.
Journal of Computational Physics , 375:659–683, 2018., 375:659–683, 2018.