Compact Higher-order Gas-kinetic Schemes with Spectral-like Resolution for Compressible Flow Simulations
CCompact Higher-order Gas-kinetic Schemes with Spectral-likeResolution for Compressible Flow Simulations
Fengxiang Zhao a , Xing Ji b , Wei Shyy a , Kun Xu a,b, ∗ a Department of Mechanical and Aerospace Engineering, Hong Kong University of Science and Technology,Clear Water Bay, Kowloon, HongKong b Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay,Kowloon, HongKong
Abstract
In this paper, a class of compact higher-order gas-kinetic schemes (GKS) with spectralresolution will be presented. Based on the high-order gas evolution model in GKS, both theinterface flux function and conservative flow variables can be evaluated explicitly from thetime-accurate gas distribution function. As a result, inside each control volume both thecell-averaged flow variables and their cell-averaged gradients can be updated within eachtime step. The flow variable update and slope update are coming from the same physicalsolution at the cell interface. Different from many other approaches, such as HWENO,there are no additional governing equations in GKS for the slopes or equivalent degrees offreedom independently inside each cell. Therefore, based on both cell averaged values andtheir slopes, compact 6th-order and 8th-order linear and nonlinear reconstructions can bedeveloped. As analyzed in this paper, the local linear compact reconstruction can achievea spectral-like resolution at large wavenumber than the well-established compact schemeof Lele with globally coupled flow variables and their derivatives. For nonlinear gas dy-namic evolution, in order to avoid spurious oscillation in discontinuous region, the abovecompact linear reconstruction from the symmetric stencil can be divided into sub-stencilsand apply a biased nonlinear WENO-Z reconstruction. Consequently discontinuous solu-tions can be captured through the 6th-order and 8th-order compact WENO-type nonlinearreconstruction. In GKS, the time evolution solution of the gas distribution function at acell interface is based on an integral solution of the kinetic model equation, which covers aphysical process from an initial non-equilibrium state to a final equilibrium one. Since theinitial non-equilibrium state is obtained based on the nonlinear WENO-Z reconstruction,and the equilibrium state is basically constructed from the linear symmetric reconstruction,the GKS evolution models unifies the nonlinear and linear reconstructions in gas evolutionprocess for the determination of a time-dependent gas distribution function, which givesgreat advantages in capturing both discontinuous shock wave and the linear aero-acousticwave in a single computation due to dynamical adaptation of non-equilibrium and equilib-rium state from GKS evolution model in different regions. This dynamically adaptive modelhelps to solve a long lasting problem in the development of high-order schemes about thechoices of the linear and nonlinear reconstructions. Compared with discontinuous Galerkin1 a r X i v : . [ phy s i c s . c o m p - ph ] J a n DG) scheme, the current compact GKS uses the same local and compact stencil, achievesthe 6th-order and 8th-order accuracy, uses a much larger time step with CFL number ≥ . Keywords: gas-kinetic scheme, WENO reconstruction, linear reconstruction, compactscheme
1. Introduction
High-order methods with spectral resolutions and shock capturing capability are neededin many engineering applications, such as turbulent flows, aeroacoustics, and various complexflows with shock and boundary layer interactions. With the great potentials of the high-order methods in solution accuracy and computational efficiency, extensive effort has beenpaid to the development of high-order schemes in the past decades. However, the spectralaccuracy and shock capturing seem to be contradictory in the current CFD algorithms, wherethe linear schemes for the smooth flow and nonlinear schemes for discontinuous flow playdifferent roles under different flow conditions [1]. It seems hard to possess both propertiesin a single high-order method.For smooth flows, the compact schemes [26, 34] and discontinuous Galerkin methods(DG) [43, 10, 11] are very attractive. The compact finite difference scheme constructsimplicit relation between derivatives and flow variables on compact stencils. Dispersion anddissipation properties of the scheme have been fully analyzed [26], and the scheme influencesgreatly on the development of linear schemes. The reason for designing compact scheme isits flexibility in complex geometries and effective parallelization. The DG has a second-order scheme stencil, but achieves high accuracy by using high-order piecewise polynomials ∗ Corresponding author
Email addresses: [email protected] (Fengxiang Zhao), [email protected] (Xing Ji), [email protected] (Wei Shyy), [email protected] (Kun Xu) ≥ .
5) in all test cases in this paper. The current 6th-order and8th-order compact GKS take advantages of both nonlinear WENO-Z reconstruction for theintroduction of numerical dissipation in discontinuous region and linear reconstruction withmuch reduced dissipation in resolving small amplitude waves. The current compact GKSare ideal schemes for the study of multiscale and complex flow interactions, such as the flowtransition and turbulence.The idea of using flow variables and their gradients has been investigated in finite volumeor finite difference methods. The finite volume Hermite WENO (HWENO) schemes updateboth flow variables and their first derivative [40, 41, 59]. The HWENO approach is also usedin the hybrid schemes [2], and a monotonicity preserving strategy for detecting troubledzones is proposed. The HWENO schemes on unstructured grids have been extensivelyinvestigated [30, 31, 32]. The HWENO approach needs additional limiters as both DGand reconstructed DG methods. The major advantage of these HWENO approaches is thecompactness of stencils in the reconstructions. The current compact GKS is different fromthe HWENO approach is mainly on the update of gradients. In GKS, there is no explicitevolution equation for the gradients and the gradients are obtained rigorously from theevolution solutions at cell interfaces. Therefore, the gradients in GKS are not introducedas additional degrees of freedom with their own evolution equations. The cell-averagedgradients in GKS like cell averaged conservative variables, which are obtained from thesame gas distribution functions at cell interfaces. The robustness and improved resolution4n GKS is due to the physically reliable gradients from the evolution solution of the originalgoverning equation. In other words, the updated gradients in GKS are not coming from the”equations” by taking additional derivatives to the original ones, such as the HWENO andmulti-layer compact (MLC) schemes [1]. In the HWENO and MCL, the gradients have theirown evolution equations which are only reliable for the smooth flow, where new governingequations can be continuously created by taking spatial derivatives to the original equation.This paper is organized as follows. The GKS and MSMD method will be introduced inSection 2. In Section 3, the compact linear and nonlinear 6th-order and 8th-order recon-structions will be presented. The dissipation and dispersion will be analyzed as well. InSection 4, the compact GKS will be tested in wide range of flow problems from the strongshock interaction to the linear acoustic wave propagation. The last section is the conclusion.
2. Gas-kinetic scheme and multi-stage multi-derivative method
In the past years, the gas-kinetic scheme (GKS) has been developed systematically [54,53]. The gas evolution model used in this paper for the flux and interface flow variableevaluation is almost identical to that of the fourth-order compact scheme [21]. Therefore,only a brief introduction about GKS will be presented in this section.
The gas kinetic evolution model in GKS is based on the BGK equation [6], f t + u · ∇ f = g − fτ , (1)where f is the gas distribution function, g is the corresponding equilibrium state that f approaches, and τ is defined as the collision time. In two-dimensional case, the equilibriumstate is the Maxwellian distribution g = ρ ( λπ ) K +22 e λ (( u − U ) +( v − V ) + ξ ) , where λ = m/ kT , and m, k, T represent the molecular mass, the Boltzmann constant, andtemperature, K is the number of internal degrees of freedom, i.e. K = (4 − γ ) / ( γ − γ is the specific heat ratio. ξ is the internal variable with ξ = ξ + ξ + ... + ξ K . The collision term satisfies the following compatibility condition (cid:90) g − fτ ψ dΞ = 0 , (2)where ψ = ( ψ , ψ , ψ , ψ ) T = (1 , u, v,
12 ( u + v + ξ )), dΞ = d u d v d ξ ... d ξ K .The macroscopic mass ρ , momentum ( ρU, ρV ), and energy ρE at a cell interface can be5valuated from the gas distribution function, W = ρρUρVρE = (cid:90) ψf d Ξ , (3)and the flux in the x direction can be evaluated as well, F = (cid:90) uf ψd Ξ . (4)Therefore, the central point in GKS is to evaluate the time-dependent gas distributionfunction f at a cell interface.By direct modeling on the mesh size scale [54], the conservations of mass, momentumand energy in a control volume get to written asd W ij d t = − x ( F i +1 / ,j ( t ) − F i − / ,j ( t )) − y ( G i,j +1 / ( t ) − G i,j − / ( t )) , where W ij is the cell averaged conservative variables, F i ± / ,j ( t ) and G i,j ± / ( t ) are the timedependent fluxes at cell interfaces in x and y directions. Like the traditional high-order finitevolume scheme, the Gaussian quadrature points are used in order to get a high accuracy. Inthis paper, three-points Gaussian quadrature will be used on each cell interface in the 2Dcase.To model the interface gas distribution function, the integral solution of BGK equation(1) is used f ( x i +1 / , y j (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) ) /τ dt (cid:48) + e − t/τ f ( − ut, − vt, u, v, ξ ) , (5)where ( x i +1 / , y j (cid:96) ) = (0 ,
0) is the location of cell interface, and x i +1 / = x (cid:48) + u ( t − t (cid:48) ) and y j (cid:96) = y (cid:48) + v ( t − t (cid:48) ) are the trajectory of particles. Here f is the initial gas distributionfunction and g is the equilibrium state in space and time. The integral solution basicallystates a physical process from the particle free transport in f in the kinetic scale to thehydrodynamic flow evolution in the integral term of g . The contributions from f and g in thedetermination of f at the cell interface depend on the ratio of time step to the local particlecollision time, i.e., exp( − ∆ t/τ ), and the modeling of f and g in space and time is neededto evaluate the solution f from Eq.(5). For the continuum flow simulation, such as the NSsolutions, the determination of f and g depend only on the macroscopic flow variables andtheir initial reconstructions. In this paper, the WENO-Z method will be used as a nonlinearreconstruction in the determination of f , and the linear reconstruction is adopted in thedetermination of g . Therefore, the above integral solution not only incorporates a physicalevolution process from non-equilibrium to an equilibrium state, but alo from a nonlinearreconstruction to a linear one. This fact is critically important for the current scheme to6apture both nonlinear shock and linear acoustic wave accurately in a single computation.Based on the integral solution, a simplified third-order gas distribution function can beobtained [58], f ( x i +1 / , y j (cid:96) , t, u, v, ξ ) = g + Ag t + 12 a tt g t − τ [( a u + a v + A ) g + ( a xt u + a yt v + a tt ) g t ] − e − t/τ g [1 − ( a u + a v ) t ]+ e − t/τ g l [1 − ( a l u + a l v )) t ] H ( u )+ e − t/τ g r [1 − ( a r u + a r v )) t ](1 − H ( u )) , (6)where the terms related to g are from the integral of the equilibrium state and the termsrelated to g l and g r are from the initial term f in the Eq.(5). All other terms, such as A, a tt , a, ... , are coming from the derivatives of a Maxwellian distribution. All these coefficientin the above equation can be determined from the initially reconstructed macroscopic flowvariables. With the above time accurate gas distribution function at a cell interface, in order toimprove the temporal accuracy of the scheme, the multi-stages multi-derivative (MSMD)technique will be used for the evaluation of the flux function and the interface values [21].The implementation of MSMD needs the flux function and its time derivative. For conser-vation laws, the semi-discrete finite volume scheme is written asd W ij d t = − x ( F i +1 / ,j ( t ) − F i − / ,j ( t )) − y ( G i,j +1 / ( t ) − G i,j − / ( t )) := L ( W ij ) , where L ij ( W ) is the numerical operator for spatial derivative of flux, F and G are obtainedfluxes at the Gaussian quadrature points.A fourth-order temporal accurate solution for W ( t ) at t = t n + ∆ t can be obtained by W ∗ = W n + 12 ∆ t L ( W n ) + 18 ∆ t ∂∂t L ( W n ) , (7) W n +1 = W n + ∆ t L ( W n ) + 16 ∆ t (cid:0) ∂∂t L ( W n ) + 2 ∂∂t L ( W ∗ ) (cid:1) , (8)where L and ∂∂t L are related to the fluxes and the time derivatives of the fluxes evaluatedfrom the time-dependent gas distribution function f ( t ) at cell interfaces. And the middlestate W ∗ is obtained at time t ∗ = t n + ∆ t/
2. Again, with the time accurate gas distributionfunction f ( t ), along the same line with MSMD technique the gas distribution function f ata cell interface at t n +1 becomes f n +1 = f n + ∆ tf nt + 16 ∆ t ( f ntt + 2 f ∗ tt ) , f ∗ is for the middle state at time t ∗ = t n + ∆ t/ f ∗ = f n + 12 ∆ tf nt + 18 (∆ t ) f ntt . Therefore, based on the cell interface f n +1 j +1 / , the flow variables W n +1 j +1 / can be explicitlyobtained, i.e., W n +1 j +1 / = (cid:82) ψf n +1 j +1 / d Ξ, from which the slope inside each cell, such as in the1D case, can be updated as ( W x ) n +1 j = ( W n +1 j +1 / − W n +1 j − / ) / ∆ x. More detailed formulation for the above GKS part can be found in [21].
3. Compact WENO reconstruction
In the classical high-order finite volume method [17, 29], the pointwise value is recon-structed based on cell averages, and the Riemann solver is used for the flux evaluation atthe interface, and the Runge-Kutta time stepping method is used for the time accuracy.However, for the gas kinetic scheme, as presented in the last section the cell-averaged con-servative flow variables and cell-averaged slopes inside each cell can be updated. Based onthese updated data, a 6th-order and a 8th-order compact linear and nonlinear reconstruc-tion will be designed consistently in this section, which can be subsequently used in thedetermination of equilibrium state and initial non-equilibrium one in the gas kinetic scheme.
In the section 2, the value W i +1 / is provided by taking moments of the gas distributionfunction at the cell interface. The averaged slope in the cell I i can be obtained by theGaussian formula W (cid:48) i = 1∆ x (cid:90) I i ∂W∂x d x = 1∆ x ( W i +1 / − W i − / ) , (9)where W (cid:48) i is the cell averaged slope in the cell I i .By taking advantage of the cell averaged variables and their slopes, the compact stencil S i +1 / = { I i − , I i , I i +1 , I i +2 } shown in Fig.1 will be used for the reconstruction of the valueat the cell interface x i +1 / . On each cell I i , the cell averages and their slopes are denotedas W i and W (cid:48) i . Based on the large stencil S i +1 / , a series of high-order polynomials can beconstructed, P n ( x ) ≡ n (cid:88) k =0 a k x k , (10)where n is the order of the polynomial and n = 5 , i I i +1 I i +2 I i − Figure 1: The large stencil for compact reconstruction of the value at the interface x i +1 / . The circlesrepresent cell averages, and the gradients represent averaged slopes. The polynomial P ( x ) can be uniquely determined with the condition1∆ x (cid:90) I k P ( x )d x = W k , x (cid:90) I k ( dP ( x ) /dx )d x = W (cid:48) k , I k ∈ S i +1 / , k = i − , · · · , i + 2 . (11)In order to determine the P ( x ), the above condition has to be modified. For the cells I i and I i +1 , the following equations are strictly enforced1∆ x (cid:90) I k P ( x )d x = W k , x (cid:90) I k ( dP ( x ) /dx )d x = W (cid:48) k , I k ∈ S i +1 / , k = i, i + 1 , (12)and for the cells I i − and I i +2 the following conditions are satisfied in the sense of the leastsquare 1∆ x (cid:90) I k P ( x )d x = W k , (cid:90) I k ( dP ( x ) /dx )d x = ∆ xW (cid:48) k , I k ∈ S i +1 / , k = i − , i + 2 . (13)Thus the linear reconstruction P n ( x i +1 / ) at the cell interface x = x i +1 / can be written as P ( x i +1 / ) = 1600 ( W i − + 299 W i + 299 W i +1 + W i +2 − xW (cid:48) i − + 111∆ xW (cid:48) i − xW (cid:48) i +1 + 3∆ xW (cid:48) i +2 ) , (14) P ( x i +1 / ) = 1420 (25 W i − + 185 W i + 185 W i +1 + 25 W i +2 +6∆ xW (cid:48) i − + 102∆ xW (cid:48) i − xW (cid:48) i +1 − xW (cid:48) i +2 ) , (15)where P n ( x i +1 / ) has the truncation error of O (∆ x n +1 ). The first-order and second-order9erivatives of the linear reconstruction at the cell interface can be obtained from P n ( x ), P x ( x i +1 / ) = 11560∆ x ( − W i − − W i + 3411 W i +1 + 3 W i +2 +11∆ xW (cid:48) i − − xW (cid:48) i − xW (cid:48) i +1 + 11∆ xW (cid:48) i +2 ) ,P xx ( x i +1 / ) = 140∆ x ( − W i − + W i + W i +1 − W i +2 +3∆ xW (cid:48) i − − xW (cid:48) i + 51∆ xW (cid:48) i +1 − xW (cid:48) i +2 ) , (16) P x ( x i +1 / ) = 1108∆ x ( − W i − − W i + 270 W i +1 + 14 W i +2 − xW (cid:48) i − − xW (cid:48) i − xW (cid:48) i +1 − xW (cid:48) i +2 ) ,P xx ( x i +1 / ) = 14∆ x ( − W i − + 4 W i + 4 W i +1 − W i +2 − ∆ xW (cid:48) i − − xW (cid:48) i + 9∆ xW (cid:48) i +1 + ∆ xW (cid:48) i +2 ) . (17)In order to eliminate the spurious oscillation and improve the stability, the above recon-struction is based on the characteristic variables. The dimension-by-dimension strategy isapplied for the reconstruction in the two-dimensional case. The above compact reconstruc-tion is mainly used to determine the equilibrium state in the normal direction in the gaskinetic scheme. Due to multi-dimensionality in the GKS evolution model, in the tangentialdirection a linear 5th-order reconstruction is used in both 6th-order and 8th-order schemes.The above linear reconstruction is used to determine the equilibrium state in GKS. In thereconstruction of the initial non-equilibrium state, the WENO-based compact nonlinearreconstruction will be designed in later subsection. In this section, the resolution analysis of the above compact spatial reconstruction willbe presented. The linear advection equation is used for the resolution analysis, W t = − U W x . (18)The integral form on the cell I i becomes( W i ) t = − U ( W i ) x , (19)where W i is the cell averaged W ( x ) in the cell I i . Remark 1.
In order to make Fourier expansion for the cell averages W i , a continuousfunction W ( x ) is introduced [17] as follows. W ( x ) = 1∆ x (cid:90) x +∆ x/ x − ∆ x/ W ( x + y ) dy. (20)10 hus the cell average W i becomes the value of the continuous function W ( x ) at cell center x i . In order to analyse the resolution of spatial discretization in the current compact schemes,we can compare the Fourier coefficients of the derivative given by the spatial discretizationwith the Fourier coefficients of exact derivative. Using the reconstructed value at the cellinterface, the spatial discretization of ( W i ) x in Eq.(19) becomes( W i ) thx = 1600∆ x (298( W i +1 − W i − ) + ( W i +2 − W i − ) − x ( W (cid:48) i +1 + W (cid:48) i − ) + 3∆ x ( W (cid:48) i +2 + W (cid:48) i − ) + 222∆ xW (cid:48) i ) , (21)( W i ) thx = 1420∆ x (160( W i +1 − W i − ) + 25( W i +2 − W i − ) − x ( W (cid:48) i +1 + W (cid:48) i − ) − x ( W (cid:48) i +2 + W (cid:48) i − ) + 204∆ xW (cid:48) i ) . (22)For the purpose of Fourier analysis, the function W ( x ) is assumed to be periodic overthe domain[0 , L ], and W ( x ) can be decomposed into its Fourier coefficients W ( x ) = k = N/ (cid:88) k =0 c k e ( πikxL ) . It is convenient to introduce new variables ω = 2 πk ∆ x/L = 2 πk/N and s = x/ ∆ x [26].Then the Fourier expansion can be written as W ( x ) = k = N/ (cid:88) k =0 c k e iω ( k ) s , and the first-order derivative of W ( x ) is W (cid:48) ( x ) = k = N/ (cid:88) k =0 c k iωh e iw ( k ) s = k = N/ (cid:88) k =0 c (cid:48) k e iω ( k ) s , where c (cid:48) k = c k iωh . Suppose the numerical solution given by the compact schemes can bewritten as W h ( x ) = k = N/ (cid:88) k =0 ˆ c k e iω ( k ) s , W (cid:48) h ( x ) = k = N/ (cid:88) k =0 ˆ c (cid:48) k e iω ( k ) s . Since the current spatial discretization is based on the symmetrical stencil, it can be shownthat ˆ c (cid:48) k = iω (cid:48) c k . Substituting the average value in Eq.(21) and Eq.(22) with the form of11ourier expansion, we get the modified wavenumber ω (cid:48) , th ( ω ) = 1300 (298 sin( ω ) + sin(2 ω ) − ω cos( ω ) + 3 ω cos(2 ω ) + 111 ω ) , (23) ω (cid:48) , th ( ω ) = 1210 (160 sin( ω ) + 25 sin(2 ω ) − ω cos( ω ) − ω cos(2 ω ) + 102 ω ) . (24)In addition to the modified wavenumber, the error in the phase speed can be alternativelyused to evaluate the dispersive error. For the linear advection equation, the phase speed fora wavenumber ω is given by the current spatial discretization (exact time advancement isassumed) as ( c p ) h = ω (cid:48) ( ω ) /ω. (25) wavenumber m od i f i e d w ave nu m b e r exactcompact GKS-6thcompact GKS-8thcs-6thcs-8th wavenumber ph ase s p ee d exactcompact GKS-6thcompa Figure 2: Plots of modified wavenumber and phase speed vs wavenumber for different schemes. cs-6th isthe compact 6th-order tridiagonal scheme studied by Lele [26], and the scheme has the best resolution inthe series of 6th-order scheme; cs-8th is the compact 8th-order pentadiagonal scheme [26], and the cs-8thuses the same formula in spatial discretization as the compact 8th-order GKS.
The comparison of the dispersion characteristics of the current compact GKS with thecompact schemes by Lele [26] is shown in Fig.2. The traditional compact schemes use onlythe node values [26] or cell averages [12] as independent values in spatial discretization, andthe curves of modified wavenumber and phase speed in the range ω ∈ ( π/ , π ) deviate greatlyfrom the curves corresponding to the analytical solution. Therefore, these schemes have poorresolution for the wave with ω ∈ ( π/ , π ), because the wavenumber range ( π > ω > π/ x, x ), and such wave cannot be determined by the cellaverages (or point values) alone with less than four points in a wavelength. However, thecompact GKS can retain a good resolution because of the use of slopes in each cell, andthese slopes are evaluated from the gas evolution solutions at the cell interfaces which areindependent from the cell averaged values. In the traditional compact schemes by Lele, the12 d e n s i t y exactlinear upwind-7thlinear compact GKS-6thlinear compact GKS-8th t=1 x d e n s i t y exactlinear upwind-7thlinear compact GKS-6thlinear compact GKS-8th t=1 Figure 3: Limit resolution for linear wave: Results of linear upwind 7th-order scheme, linear compact6th-order and 8th-order schemes on the meshes with ∆ x = 1 / x = 1 / t = 1. Since there areonly two points for each wavelength on the left figure, the real ”exact” solution should be averaged over halfwave length as well. node values and derivatives are not fully independent. The coupling increases the error sincethe node values cannot resolve high wavenumber solution at all. As an example, we test theschemes for the initial condition ( ρ, U, p ) = (1+0 . πx ) , ,
1) in a computational domain[0 , t = 1 with meshes ∆ x = 1 / x = 1 / In gas dynamics, shock and contact discontinuities can appear. To deal with discontinu-ities, WENO reconstruction [29, 23] can be used. In section 3 .
1, the linear reconstructiongives a unique reconstructed value and its slope at the cell interface. However, in order tocapture the possible discontinuity at a cell interface, the values at the left and right sides ofthe cell interface have to be valuated separately. Therefore, for the nonlinear reconstructionthe sub-stencils have to be defined first. For simplicity, the WENO reconstruction procedureis given in detail for the construction of the left side value of the cell interface x = x i +1 / .The procedure for the right side value can be obtained similarly according to the symmetricproperty, which will be omitted here. The left side value by WENO reconstruction is givenby the candidate polynomials as follows W ni +1 / = l n (cid:88) k =0 δ nk w nk,i +1 / , n = 6 , , (26)13here W ni +1 / is the nth − order WENO reconstruction for the left value of the cell interface x i +1 / , l n is the number of candidate polynomials, δ nk is the WENO weight, and w nk,i +1 / isthe point value of the candidate polynomial w nk ( x ) at x i +1 / . In the later presentation, thesets of candidate polynomials w nk ( x ) are the same for different nth − order scheme, such as n = 6 and 8, so the superscript is omitted for simplicity. I i I i +1 I i +2 I i − S S S S S S S Figure 4: Sub-stencils in the reconstruction for left side value of the cell interface x i +1 / in the compact6th-order and 8th-order reconstruction: the circles represent cell averages, and the gradients represent cellaveraged slopes. In the current compact WENO reconstruction, three principles are considered. Firstly,if a discontinuity is located at one of the cell interfaces, theoretically the point-wise valueat the interface cannot be reliable to determine the cell averaged slope on the left and rightside cells. Secondly, because the smooth sub-stencils can play a dominant role in dealingwith discontinuities appearing in the large stencil, the averaged slope used in the sub-stencilshould be possibly away from the interface for the data reconstruction in order to make thesub-stencil sufficiently smooth. Thirdly, the order of the candidate polynomial can be higherwithout increasing the spatial size of sub-stencils in the compact schemes, and such higher-order candidate polynomials can achieve better resolution to solve discontinuities withoutspurious oscillation. Based on the above principles, the following sub-stencils are designed,14here the 6th-order and 8th-order reconstructions are based on the same sub-stencils, S = { W i − , W i , W (cid:48) i − } ↔ w ( x ) ,S = { W i − , W i , W i +1 } ↔ w ( x ) ,S = { W i , W i +1 , W i +2 , W (cid:48) i +2 } ↔ w ( x ) ,S = { W i − , W i , W i +1 } ↔ w ( x ) ,S = { W i , W i +1 , W i +2 } ↔ w ( x ) ,S = { W i − , W i , W i +1 , W (cid:48) i − , W (cid:48) i } ↔ w ( x ) ,S = { W i , W i +1 , W i +2 , W (cid:48) i +1 , W (cid:48) i +2 } ↔ w ( x ) .w ( x ) is a cubic polynomial, w ( x ) and w ( x ) are fourth-order polynomials, and othersare quadratic polynomials. According to the similar reconstruction conditions in Eq.(11), w k,i +1 / can be obtained as w ,i +1 / = 16 ( − W i − + 13 W i − hW (cid:48) i − ) ,w ,i +1 / = 16 ( − W i − + 5 W i + 2 W i +1 ) ,w ,i +1 / = 124 (5 W i + 32 W i +1 − W i +2 + 6 hW (cid:48) i +2 ) ,w ,i +1 / = 16 ( − W i − + 5 W i + 2 W i +1 ) ,w ,i +1 / = 16 (2 W i + 5 W i +1 − W i +2 ) ,w ,i +1 / = 130 (10 W i − + 19 W i + W i +1 + 3 hW (cid:48) i − + 21 hW (cid:48) i ) ,w ,i +1 / = 130 ( W i + 19 W i +1 + 10 W i +2 − hW (cid:48) i +1 − hW (cid:48) i +2 ) . (27)In the smooth region, the convex combination with δ nk = d nk recovers the reconstruction inEq.(14) and Eq.(15), which can be the condition to get the linear weights [23]. The linearweight d k of the 6th-order reconstruction can be obtained as d = 33700 , d = 11140 , d = 22175 , d = 11100 , d = 11100 , d = 37140 , d = 37140 , and the d k of 8th-order reconstruction are d = 398 , d = 598 , d = 449 , d = 798 d = 798 , d = 1749 , d = 1749 . The WENO-Z nonlinear weights are used in the current compact schemes and they aredefined as [7] δ nk = α nk (cid:80) l K m =0 α nm , α nk = d nk (cid:104) (cid:0) Z nref β k + (cid:15) (cid:1)(cid:105) , k = 0 , ..., l n , (28)15here Z nref is the local higher-order reference value, which is related to the accuracy of thescheme and is given in detail later. Here β k is the smooth indicator and defined as [23] β k = q k (cid:88) q =1 ∆ x q − (cid:90) x i +1 / x i − / (cid:0) d q d x q w k ( x ) (cid:1) dx, (29)where q k is the order of w k ( x ).In smooth region, the first two candidate polynomial w ( x ) and w ( x ) can be combinedinto a cubic polynomial which is symmetric counterpart of w ( x ). Then, current sub-stencilscan become basically symmetric for interface x i +1 / . In order to maintain the symmetry forthe nonlinear schemes, the smooth indicator β of w ( x ) corresponding to S is replaced bythe indicator of the cubic polynomial (cid:101) w ( x ) on (cid:101) S = { W i − , W i , W i +1 , W (cid:48) i − } . Even withoutshowing in this paper, some tests demonstrate that the current choice β can present a betterresolution in the numerical results with excellent robustness, The detailed formulae for all β k , k = 0 , ..., l n are given in Appendix. In this subsection, the accuracy of the compact reconstruction is analysed. The nonlinearreconstruction in Eq.(26) can be rewritten as W ni +1 / = l n (cid:88) k =0 d nk w k,i +1 / + l n (cid:88) k =0 ( δ nk − d nk ) w k,i +1 / . (30)The reconstruction can be split into two terms [18], i.e., the linear term and nonlinear term.For the linear part W n,opti +1 / ≡ (cid:80) l n k =0 d nk w k,i +1 / , the error can be given as W n,opti +1 / − W ( x i +1 / ) = A n ( x i +1 / )∆ x n + O (∆ x n +1 ) , where W ( x i +1 / ) is the exact solution at x = x i +1 / . In the second term in Eq.(30), thepoint value w k,i +1 / approximates W ( x i +1 / ) as follows w k,i +1 / = W ( x i +1 / ) + B k ( x i +1 / )∆ x q k + O (∆ x q k +1 ) . Substituting w k,i +1 / into Eq.(30) and taking l n (cid:88) k =0 δ nk = l n (cid:88) k =0 d nk = 1 into account, we have W ni +1 / = W n,opti +1 / + l n (cid:88) k =0 ( δ nk − d nk ) (cid:0) W ( x i +1 / ) + B k ( x i +1 / )∆ x q k + O (∆ x q k +1 ) (cid:1) = W n,opti +1 / + l n (cid:88) k =0 B k ( x i +1 / )( δ nk − d nk )∆ x q k + l n (cid:88) k =0 ( δ nk − d nk ) O (∆ x q k +1 ) .
16o achieve the formal order of accuracy for the 6th-order and 8th-order reconstructions, thefollowing sufficient condition is proposed. δ nk − d nk = O (∆ x n − n (cid:48) ) , n (cid:48) ≥ . (31)For the current nonlinear weights with WENO-Z weighting functions, with the followingformulation of the local higher-order reference value Z nref , the sufficient condition Eq.(31)can be satisfied when Z nref = | β − β ) + ( β − β ) | . (32)In order to prove that the sufficient condition can be met, according to the Taylor series of β k , k = 0 , , β = (cid:0) ( W (1) i ) + 1312 ( W (2) i ∆ x ) (cid:1) ∆ x − W (1) i W (3) i ∆ x + J ( W ( l ) i , ∆ x )∆ x ,β = (cid:0) ( W (1) i ) + 1312 ( W (2) i ∆ x ) (cid:1) ∆ x + 13 W (1) i W (3) i ∆ x + J ( W ( l ) i , ∆ x )∆ x ,β = (cid:0) ( W (1) i ) + 1312 ( W (2) i ∆ x ) (cid:1) ∆ x − W (1) i W (3) i ∆ x + J ( W ( l ) i , ∆ x )∆ x , l = 1 , , · · · we have Z nref = | J − J − J | ∆ x , where J , J and J are the functions of the Taylor expanded terms of W ( x ) as W ( x ) isalways continuous. Suppose O ( J k ) ∼ O (1), the sufficient condition Eq.(31) is satisfied by α nk = d nk (cid:104) (cid:16) Z nref β k + (cid:15) (cid:17)(cid:105) = d nk (cid:2) C + O (∆ x ) (cid:3) , k = 0 , ..., l n ,δ nk = α nk (cid:80) l n m =0 α rm = d nk + O (∆ x ) , where C ≥
4. Numerical tests
In this section, we are going to test the 6th-order and 8th-order compact gas kineticschemes. The cases include linear acoustic waves, blast wave, shock-shock interactions,shock acoustic wave interaction, as well as viscous flow computations. The mesh used inthis paper is rectangular one and the time step is determined by the CFL condition with aCFL number ( ≥ .
3) in all test cases. In all tests, the same linear reconstruction is used forthe equilibrium state and the nonlinear reconstruction for the non-equilibrium state. Thereis no any additional ”trouble cell” detection or any limiter designed for specific test. The gas17inetic evolution model basically presents a dynamical process from the nonlinear to linearone, and the convergence rate exp( − ∆ t/τ ) depends on the flow condition. The collisiontime τ for inviscid flow at a cell interface is defined by τ = (cid:15) ∆ t + C | p l − p r p l + p r | ∆ t, where ε = 0 . C = 1, and p l and p r are the pressures at the left and right sides of a cellinterface. For the viscous flow, the collision time is related to the viscosity coefficient, τ = µp + C | p l − p r p l + p r | ∆ t, where µ is the dynamic viscosity coefficient and p is the pressure at the cell interface. Insmooth flow regions, it will reduce to τ = µ/p . The reason for including pressure jump termin the particle collision time is to add artificial dissipation in the discontinuous region tokeep the non-equilibrium dynamics in the shock layer. The one-dimensional advection of density perturbation is tested first, and the initialconditions are given as follows ρ ( x ) = 1 + 0 . πx ) , U ( x ) = 1 , p ( x ) = 1 , x ∈ [0 , . The periodic boundary condition is adopted, and the analytic solutions are ρ ( x, t ) = 1 + 0 . π ( x − t )) , U ( x, t ) = 1 , p ( x, t ) = 1 . With the r th-order spatial reconstruction and 2-stage 4th-order temporal discretization, theleading term of the truncation error is O (∆ x r + ∆ t ). To keep the r th-order accuracy inthe test, ∆ t = C ∆ x r/ needs to be used for the r th-order scheme. In the computation, theuniform meshes with N mesh points are used. The L , L and L ∞ errors and convergenceorders at t = 2 for the 6th-order and 8th-order compact GKS are presented in Table 1 toTable 4, respectively. For the 8th-order scheme, the errors reduce quickly, and the order ofaccuracy does not attain 8 with N = 80 due to the limited round-off error. mesh length L error Order L error Order L ∞ error Order1/5 1.761e-003 1.976e-003 2.721e-0031/10 3.171e-005 5.80 3.527e-005 5.81 4.899e-005 5.801/20 5.017e-007 5.98 5.554e-007 5.99 7.803e-007 5.971/40 7.818e-009 6.00 8.676e-009 6.00 1.224e-008 5.991/80 1.220e-010 6.00 1.355e-010 6.00 1.915e-010 6.00 Table 1: 1-D accuracy test: errors and convergence orders of 6th-order compact linear scheme with∆ t = 0 . x / . esh length L error Order L error Order L ∞ error Order1/5 1.196e-004 1.344e-004 1.849e-0041/10 4.715e-007 7.99 5.301e-007 7.99 7.286e-007 7.991/20 1.806e-009 8.03 1.998e-009 8.05 2.798e-009 8.021/40 6.979e-012 8.02 7.744e-012 8.01 1.092e-011 8.001/80 5.588e-014 6.96 6.207e-014 6.96 9.614e-014 6.83 Table 2: 1-D accuracy test: errors and convergence orders of 8th-order compact linear scheme with∆ t = 0 . x . mesh length L error Order L error Order L ∞ error Order1/5 7.495e-003 8.448e-003 1.248e-0021/10 3.624e-004 4.37 4.307e-004 4.29 7.357e-004 4.081/20 7.877e-006 5.52 8.677e-006 5.63 1.461e-005 5.651/40 4.216e-008 7.55 5.440e-008 7.32 1.235e-007 6.891/80 2.119e-010 7.64 2.647e-010 7.68 7.407e-010 7.381/160 2.100e-012 6.66 2.279e-012 6.86 3.540e-012 7.71 Table 3: 1-D accuracy test: errors and convergence orders of 6th-order compact nonlinear schemes with∆ t = 0 . x / . mesh length L error Order L error Order L ∞ error Order1/5 5.154e-003 5.969e-003 8.959e-0031/10 3.825e-004 3.75 4.392e-004 3.76 6.522e-004 3.781/20 7.963e-006 5.59 9.252e-006 5.57 1.588e-005 5.361/40 4.250e-008 7.55 6.042e-008 7.26 1.503e-007 6.721/80 1.339e-010 8.31 2.201e-010 8.10 7.984e-010 7.561/160 2.827e-013 8.89 3.801e-013 9.18 1.206e-012 9.37 Table 4: 1-D accuracy test: errors and convergence orders of 8th-order compact nonlinear scheme with∆ t = 0 . x . To test the performance of the schemes for capturing high frequency waves, the Shu-Osher problem for density-wave shock interaction is tested [47]. The initial condition isgiven by ( ρ, U, p ) = (cid:26) (3 . , . , . , x ≤ , (1 + 0 . x ) , , , < x ≤ . The computational domain is [0 ,
10] and 200 uniform mesh points are used. The non-reflecting boundary condition is used at both ends. The computed density profile and localenlargement at t = 1 . d e n s i t y referencecompact GKS-6thcompact GKS-8th x d e n s i t y referencecompact GKS-6thcompact GKS-8th Figure 5: Shu-Osher problem: the density distribution and local enlargement by 6th-order and 8th-ordercompact GKS at t = 1 . [49], and the initial condition in this case is the following( ρ, U, p ) = (cid:40) (1 . , . , . , − < x ≤ − . , (1 + 0 . πx ) , , , − . < x < . The computational domain is [ − , t = 5, local enlargement, and the exact solution for the Titarev-Toro problem areshown in Fig.6. In order to show the importance and accuracy of the 6th-order and 8th-ordercompact reconstructions, the results from the same GKS, but with the non-compact standard7th-order WENO-JS reconstruction using the cell-averaged flow variables only, are includedas well. As shown in Fig.6, even with 7th-order WENO reconstruction the dissipation anddispersion errors are much larger than those from the 6th-order compact GKS. Based on thisobservation, it clearly indicates that the use of high-order evolution model for the update ofslope is favorable in the design of high-order schemes. The compactness of the scheme is alsophysically sounded because the CFL condition constrains the signal propagation locally tothe neighboring cells only within a time step. To design a reliable compact scheme dependson the high-order gas evolution model at the cell interface. Any scheme based on the firstorder Riemann solver cannot achieve such a goal to have scheme with the properties ofcompactness, large CFL number, robustness, and high efficiency.The third test is the Woodward-Colella blast wave problem [51], and the initial conditions20 d e n s i t y -4 -2 0 2 411.21.41.6 referenceWENO-7JScompact GKS-6thcompact GKS-8th x d e n s i t y -2 -1.8 -1.6 -1.41.31.41.51.61.7 referenceWENO-7JScompact GKS-6thcompact GKS-8th x d e n s i t y referenceWENO-7JScompact GKS-6thcompact GKS-8th x d e n s i t y referenceWENO-7JScompact GKS-6thcompact GKS-8th Figure 6: Titarev-Toro problem: the density distribution and local enlargement by 6th-order and 8th-ordercompact GKS at t = 5 with 1000 mesh points. are given as follows ( ρ, U, p ) = (1 , , , ≤ x < , (1 , , . , ≤ x < , (1 , , , ≤ x ≤ . The computational domain is [0 , t = 3 . d e n s i t y referencecompact GKS-6thcompact GKS-8th x d e n s i t y
65 70 75 803456 re x d e n s i t y x d e n s i t y
65 70 75 803456 referencecompact GKS-6thcompact GKS-8th
Figure 7: Blast wave problem: the density distributions of 6th-order and 8th-order compact GKS for theblast-wave problem at t = 3 . The one-dimensional acoustic wave propagation in x -direction was proposed by Bai etal. in [1]. The test demonstrates the high order and high resolution of the compact GKSto compute acoustic wave propagating through a long-distance. The initial conditions are22 d e n s i t y pu r t e r b a t i on referencelinear compact GKS-6thcompact GKS-6thlinear compact GKS-8thcompact GKS-8th t=0.01 x d e n s i t y pu r t e r b a t i on referencelinear compact GKS-6thcompact GKS-6thlinear compact GKS-8thcompact GKS-8th t=1.0 Figure 8: One-dimensional acoustic problem: the distributions of density perturbation obtained by compactGKS at t = 0 .
01 (left, 1 period) and t = 1 . λ ρ /a ∞ = 4 . × − . x d e n s i t y pu r t e r b a t i on referencelinear compact GKS-6thcompact GKS-6thlinear compact GKS-8thcompact GKS-8th t=0.01 x d e n s i t y pu r t e r b a t i on referencelinear compact GKS-6thcompact GKS-6thlinear compact GKS-8thcompact GKS-8th t=1.0 Figure 9: One-dimensional acoustic problem: the distributions of density perturbation obtained by compactGKS at t = 0 .
01 (left, 1 period) and t = 1 . given as follows U = U ∞ + δU, δU = (cid:15)a ∞ cos( ωx ) , U ∞ = 0 ρ = ρ ∞ + δρ, δρ = (cid:15)ρ ∞ cos(2 ωx ) , ρ ∞ = 1 . pp ∞ = ( ρρ ∞ ) r , p ∞ = 101325 . a ∞ = (cid:114) γ p ∞ ρ ∞ , (cid:15) = 10 − is the magnitude of initial perturbation, and ω = 6 π is the wavenumber ofinitial perturbations in velocity. The specific heat ratio is γ = 1 .
4. The acoustic wave givenabove is approximately linear because of the very small (cid:15) , and an analytical solution [1] isgiven from the approximate acoustic wave equation, ρ ( x, t ) = ρ ∞ + 12 (cid:15)ρ ∞ [cos(2 ω ( x − a ∞ t )) + cos(2 ω ( x + a ∞ t ))+cos( ω ( x − a ∞ t )) − cos( ω ( x + a ∞ t ))] ,U ( x, t ) = 12 (cid:15)a ∞ [cos(2 ω ( x − a ∞ t )) − cos(2 ω ( x + a ∞ t ))+cos( ω ( x − a ∞ t )) + cos( ω ( x + a ∞ t ))] ,p ( x, t ) = p ∞ + 12 γ(cid:15)p ∞ [cos(2 ω ( x − a ∞ t )) + cos(2 ω ( x + a ∞ t ))+cos( ω ( x − a ∞ t )) − cos( ω ( x + a ∞ t ))] . (33)The computational domain is [0 , / ρ i = (cid:82) i +1 / i − / ρ ( x ) dx exactly, then convert them to the conservative variables. Similarly, we couldobtain the cell average of the derivatives, i.e. ˜ ρ x,i = (cid:82) i +1 / i − / ρ x ( x ) dx exactly, from whichthe derivatives of conservative variables can be obtained by the chain rules. Theoretically,this has a 2nd-order accuracy in space. However, it is enough for the comparison since theanalytical solution in Eq. (33) is also an approximate solution.Fig. 8 and Fig. 9 show the density distributions at t = 0 .
01 (left) and t = 1 . λ ρ /a ∞ = 4 . × − , where λ ρ is the initial wavelength of density perturbation. Thusthe distance of acoustic wave propagation is greatly larger than the initial wavelength ofdensity perturbation. The results from GKS with linear and nonlinear reconstructions arein good agreement. The results with 40 mesh points agree well with the analytical solution,where the small deviation near the extreme point is caused by the approximate analyticalsolution [1]. This case also shows that the 6th-order GKS cannot give the accurate solutionwith 20 mesh points after a long time wave propagation (100 periods) due to the relativelarge numerical error in comparison with the 8th-order scheme. The interaction between strong shock wave and turbulence is an important flow problemin gas dynamics, where the shock can cross the local wave packets which propagate fasteror slower than the mean velocity. To validate the current schemes for this kind of problems,we propose a test of shock with a Mach number
M a = 10 crossing a velocity perturbation.The initial condition is given as follows( ρ, U, p ) = (cid:26) (8 , . , . , x ≤ , (1 . , . π ( x − , , < x ≤ . d e n s i t y referenceWENO-7JScompact GKS-6thcompact GKS-8th x d e n s i t y referenceWENO-7JScompact GKS-6thcompact GKS-8th Figure 10: One-dimensional shock crossing velocity perturbation: the density distribution and local en-largement from WENO-7JS non-compact GKS, 6th-order and 8th-order compact GKS at t = 0 . x U referenceWENO-7JScompact GKS-6thcompact GKS-8th x U refer Figure 11: One-dimensional shock crossing velocity perturbation: the velocity distribution and local en-largement from WENO-7JS non-compact GKS, 6th-order and 8th-order compact GKS at t = 0 . The computational domain is [0 , t = 0 . Figure 12: Double Mach reflection: the density contours of 6th-order and 8th-order compact GKS with960 ×
240 mesh points.Figure 13: Double Mach reflection: the density contours of 6th-order and 8th-order compact GKS with1920 ×
480 mesh points. .5. Double Mach reflection problem In this subsection, the double Mach reflection problem is tested. The test was extensivelystudied by Woodward and Colella [51] for the inviscid flow. The computational domain is[0 , × [0 , x = 1 /
6. Initially a right-moving Mach 10 shock is positioned at ( x, y ) = (1 / , ◦ angle with the x-axis. The initial pre-shock and post-shock conditions are( ρ, U, V, p ) = (8 , . √ , − . , . , ( ρ, U, V, p ) = (1 . , , , . The reflecting boundary condition is used at the wall, and the exact post-shock conditionis imposed for the rest of bottom boundary. At the top boundary, the flow variables are setto describe the exact motion of the Mach 10 shock. In this case, the compact 6th-order and8th-order compact GKS are tested. The density distributions with 960 ×
240 and 1920 × t = 0 . In the following, the two-dimensional Riemann problems are considered. The first oneis the interaction of four shocks ←− S ←− S ←− S ←− S [24, 57], where the backward shock wavesconnecting the areas Ω l and Ω r are denoted by ←− S lr . The initial conditions are given as follows( ρ, U, V, p ) = (1 . , , , . , Ω : x > . , y > . , (0 . , . , , . , Ω : x < . , y > . , (0 . , . , . , . , Ω : x < . , y < . , (0 . , , . , . , Ω : x > . , y < . . This case is just the mathematical formation of the double Mach problem [51] and thesymmetric line x = y can be regarded as the rigid wall. The computational domain is[0 , × [0 , t = 0 . x = ∆ y = 1 / J − J − J − J − is tested [24, 57], where the backward contact discontinuities connectingthe areas Ω l and Ω r are denoted as J − lr . The initial conditions are given as follows( ρ, U, V, p ) = (1 , . , − . , , Ω : x > . , y > . , (2 , . , . , , Ω : x < . , y > . , (1 , − . , . , , Ω : x < . , y < . , (3 , − . , − . , , Ω : x > . , y < . . y compact GKS-6th x y compact GKS-6th x y compact GKS-8th x y compact GKS-8th Figure 14: Double Mach reflection: the local enlarged density distributions around the triple point of the6th-order (up) and 8th-order (down) compact GKS with 960 ×
240 (left) and 1920 ×
480 (right) mesh points.
The computational domain is [0 , × [0 , x = ∆ y = 1 /
500 are used. The numerical results from compact GKS are givenin Fig. 16 at t = 0 . Here we test the performance of the compact GKS for the capturing of viscous flow solu-tions. The lid-driven cavity problem is one of the most important benchmarks for validatingincompressible Navier-Stokes flow solvers. The fluid is bounded by a unit square and isdriven by a uniform translation of the top boundary. In this case, the flow is simulated with28 y compact GKS-6th x y Figure 15: 2D Riemann problem: the density distributions by 6th-order and 8th-order compact GKS forthe interaction of four shock waves at t = 0 . ×
500 mesh points.
Mach number
M a = 0 .
15 and all boundaries are isothermal and nonslip. The computa-tional domain [0 , × [0 ,
1] is covered with 65 ×
65 mesh points. Numerical simulations areconducted for two different Reynolds numbers, i.e., Re = 1000 and 3200. The streamlinescorresponding to Re = 1000 are shown in Fig. 17. The U -velocities along the center verticalline and V -velocities along the center horizontal line are shown in Fig. 18. The benchmarkdata [14] for Re = 1000 and 3200 are also presented, and the simulation results match wellwith these benchmark data. The results from 8th-order compact GKS are almost identi-cal to the reference solutions even for Re = 3200 case. The cavity case fully validates thehigher-order accuracy of compact GKS for viscous flow.29 igure 16: 2D Riemann problem: the density distributions from 6th-order and 8th-order compact GKS forthe interactions of four contact discontinuities at t = 0 . ×
500 mesh points.Figure 17: Lid-driven cavity flow: streamlines with 65 ×
65 uniform mesh points for Re = 1000 by 6th-orderand 8th-order compact GKS. The case is an interaction of a shock wave with a single vortex as studied by Inoue andHattori [20], where the fluid is viscous. The computational domain is [ − , × [ − , U θ ( r ) = M v re (1 − r ) / , U r = 0 , u compact GKS-6thcompact GKS-8thGhia’s data Re=1000 x v compact GKS-6thcompact GKS-8thGhia’s data Re=1000 y u compact GKS-6thcompact GKS-8thGhia’s data Re=3200 x v compact GKS-6thcompact GKS-8thGhia’s data Re=3200
Figure 18: Lid-driven cavity flow: U-velocity along vertical centerline and V-velocity along horizontalcenter-line with 65 ×
65 uniform mesh points at Re = 1000 and 3200. and the distribution of pressure and density upstream of shock are p ( r ) = 1 γ [1 − γ − M v e (1 − r ) ] γ/ ( γ − ,ρ ( r ) = [1 − γ − M v e (1 − r ) ] / ( γ − , where U θ and U r are the tangential and radial velocity respectively. Mach number M v of thevortex is M v = 0 .
25. The Mach number of shock wave is M s = 1 .
2. The Reynolds numberis Re = 800 defined by Re = ρ ∞ a ∞ /µ ∞ , where the subscript ∞ denotes the quantitydownstream of the shock wave. The initial location of vortex is ( x v , y v ) = (1 , x = 0. In the computation, the supersonic inflow boundary conditions31 y -15 -10 -5 0 5-10-50510 t=6compact GKS-6 x y -15 -10 -5 0 5-10-50510 t=6compact GKS-8th Figure 19: Two-dimensional acoustic problem: sound pressure contours obtained by 6th-order and 8th-order compact GKS at t = 6 with 700 ×
600 mesh points. Shock and vortex Mach numbers are M s = 1 . M v = 1. Here 150 equal-spaced sound pressure contours from − . .
05 are plotted. The dash linerepresents rarefaction region, and the solid line represents the compression region. r s ound p r ess u r e compact GKS-6threference, Inoue et al. r s ound p r ess u r e compact GKS-8threference, Inoue et al. Figure 20: Two-dimensional acoustic problem: radial distribution of the sound pressure (cid:52) p obtained by6th-order and 8th-order compact GKS at t = 6 with 700 ×
600 mesh points. The vortex center at t = 6 isapproximately located at ( − . , . at x = 8 as well as the periodic boundary conditions at y = ±
12 are imposed. The non-reflecting boundary conditions are adopted at x = − ×
600 cells. The sound pressure at t = 6 are shown in Fig. 19, where the soundpressure is defined as (cid:52) p = ( p − p ∞ ) /p ∞ . The pressure distribution in Fig. 19 clearly shows32hat the incident shock wave and two reflected shock waves are connected at the triple point,and the reflected shock waves emanate from the compression region of the incident shockwave. In addition, both the first and second sound waves have the clear quadruple naturewith opposite sign, which is similar with the result in [20]. Both the 6th-order and 8th-orderGKS perform well to resolve the shock wave and density wave, and the solution is sufficientlysmooth to resolve small amplitude waves. The radial distribution of the sound pressure (cid:52) p obtained by 6th-order and 8th-order compact GKS at t = 6 with 700 ×
600 mesh pointsare shown in Fig. 20. The radial distribution of the sound pressure is subtracted from thevortex center with the angle θ = 135 degree with the positive x direction. The vortex centerat t = 6 is approximately located at ( − . , . × /
10 cell size reduction at the shock region.
5. Conclusion
In this paper, a class of compact high-order gas-kinetic scheme with WENO reconstruc-tion has been proposed for the compressible Euler and Navier-Stokes solutions. Based on thegas-kinetic theory, the current scheme depends solely on the time-accurate solution of a gasdistribution function at a cell interface for the flux and conservative flow variable evaluation.With the update of pointwise values at cell interfaces, besides the cell-averaged conservativeflow variables, their cell-averaged gradients can be updated as well. Therefore, the com-pact linear and nonlinear reconstructions can be obtained from the local cell averaged flowvariables and their slopes. With symmetrical stencils, the 6th-order and 8th-order compactlinear reconstructions are given and used in the 1-D Fourier analysis. The analysis elucidatesspectral resolution of the current spatial discretization, even for very large wave number.The GKS flux makes the scheme stable even for the initial reconstruction from symmetricalstencils. In order to capture shock and other discontinuities, the nonlinear compact WENOreconstruction has been constructed as well. Compared with other compact schemes, thecurrent compact GKS is fully local and has explicit high-order discretization, without usingany formulation for the globally connected flow variable and their slopes. Equipped withboth linear and nonlinear compact reconstruction, the GKS gas evolution model can make asmooth transition in dynamics from the initial nonlinear reconstruction for non-equilibriumstate to the final linear reconstruction for equilibrium state. The transition rate depends onthe local flow conditions. As a result, the current compact high-order scheme can naturallycapture both shock discontinuities and small amplitude acoustic waves in a single computa-tion without using any trouble cell detection and additional limiters. Due to the existenceof the time-derivative of the flux function, the multi-stage and multi-derivative time step-ping technique has been applied here for the achieving high-order temporal accuracy of thescheme. More specifically, the two stage fourth-order time accuracy has been achieved. BothGKS and DG methods have the same compact stencils. But, the GKS uses a much largeCFL number ( ≥ .
3) for the 6th- and 8th-order schemes, and it is much more robust incapturing shock discontinuity and is much more efficient in simulating NS equations thanthe DG method. The extensive tests in this paper clearly demonstrate that the current 6th-33rder and 8th-order compact schemes advance the development of high-order CFD methodsto a new level of maturity. The success of the current compact schemes depends solely onthe high gas evolution model at a cell interface, which keeps and traces flow dynamics asclose as possible to a real physical time evolution process. Any scheme based on the firstorder Riemann solver, which gives a time-independent gas evolution at a cell interface fromtwo constant states whatever the high-order initial reconstruction is, may have difficultyto possess all these good properties of a high-order scheme, such as the compactness, largeCFL number, robustness, high efficiency, and the dynamic unification of linear and nonlinearreconstructions.
Acknowledgements
The current research is supported by Hong Kong research grant council (16206617) andNational Science Foundation of China (11772281, 91852114).
Appendix
For the WENO reconstruction, with the definition of smooth indicator in Eq.(29), thereconstructed polynomials from the sub-stencils for the determination of the left state of thecell interface x i +1 / are presented. Assume the candidate polynomial on a specific sub-stencilas w n ( x ) ≡ n (cid:88) k =0 a nk x k /k ! , (34)where n is the order of the polynomial and takes the values n = 2 , , w n ( x ) into the smooth indicator β in Eq.(29), the result can be expressed by the coefficient a nk as β =( a n · ∆ x ) + 1312 ( a n · ∆ x ) + 1043960 ( a n · ∆ x ) +3504532256 ( a n · ∆ x ) + 112 ( a n · ∆ x ) · ( a n · ∆ x ) + 780 ( a n · ∆ x ) · ( a n · ∆ x ) . For n = 2, a n = a n = 0, and for n = 3, a n = 0. Therefore, once the candidate polynomial isgiven, the smoothing indicator can be fully obtained.The polynomials of the sub-stencils for the 6th-order and 8th-order schemes are listedbelow. 34 ( x ) = 112 ( W i − + 11 W i + ∆ xW (cid:48) i − ) + 1∆ x ( − W i − + 2 W i − ∆ xW (cid:48) i − ) · x − x ( W i − − W i + ∆ xW (cid:48) i − ) · x , (cid:101) w ( x ) = 124 ( − W i − + 26 W i − W i +1 ) + 116∆ x ( − W i − + 20 W i + 3 W i +1 − xW (cid:48) i − ) · x +1∆ x ( W i − − W i + W i +1 ) · x x (3 W i − − W i + W i +1 + 2∆ xW (cid:48) i − ) · x ,w ( x ) = 148 (43 W i + 16 W i +1 − W i +2 + 6∆ xW (cid:48) i +2 ) + 116∆ x ( − W i + 60 W i +1 − W i +2 + 14∆ xW (cid:48) i +2 ) · x +12∆ x (5 W i − W i +1 + 11 W i +2 − xW (cid:48) i +2 ) · x x ( − W i + 4 W i +1 − W i +2 + 2∆ xW (cid:48) i +2 ) · x ,w ( x ) = 124 ( − W i − + 26 W i − W i +1 ) + 12∆ x ( − W i − + W i +1 ) · x + 1∆ x ( W i − − W i + W i +1 ) · x ,w ( x ) = 124 (23 W i + 2 W i +1 − W i +2 ) + 12∆ x ( − W i + 4 W i +1 − W i +2 ) · x + 1∆ x ( W i − W i +1 + W i +2 ) · x ,w ( x ) = 1960 ( − W i − + 1148 W i − W i +1 − xW (cid:48) i − − xW (cid:48) i ) + 18∆ x ( W i − − W i +1 + 10∆ xW (cid:48) i ) · x +14∆ x (19 W i − − W i + W i +1 + 6∆ xW (cid:48) i − + 12∆ xW (cid:48) i ) · x − x ( W i − − W i +1 + 2∆ xW (cid:48) i ) · x − x (5 W i − − W i − W i +1 + 2∆ xW (cid:48) i − + 4∆ xW (cid:48) i ) · x ,w ( x ) = 1960 (707 W i − W i +1 + 545 W i +2 − xW (cid:48) i +1 − xW (cid:48) i +2 )+18∆ x ( − W i + 8 W i +1 + 13 W i +2 − xW (cid:48) i +1 − xW (cid:48) i +2 ) · x +14∆ x (25 W i + 28 W i +1 − W i +2 + 60∆ xW (cid:48) i +1 + 18∆ xW (cid:48) i +2 ) · x − x (3 W i + 8 W i +1 − W i +2 + 10∆ xW (cid:48) i +1 + 4∆ xW (cid:48) i +2 ) · x x ( W i + 4 W i +1 − W i +2 + 4∆ xW (cid:48) i +1 + 2∆ xW (cid:48) i +2 ) · x , References [1] Z. Bai, X. Zhong, New very high-order upwind multi-layer compact (MLC) schemes with spectral-likeresolution for flow simulations, J. Comput. Phys. 378 (2019) 63-109.[2] D.S. Balsara, C. Altmann, C.-D. Munz, M. Dumbser, A sub-cell based indicator for troubled zones inRKDG schemes and a novel class of hybrid RKDG+HWENO schemes, J. Comput. Phys. 226 (2007)586-620.[3] M. Ben-Artzi, J. Falcovitz, A second-order Godunov-type scheme for compressible fluid dynamics, J.Comput. Phys. 55 (1984) 1-32.[4] M. Ben-Artzi, J. Li, Hyperbolic conservation laws: Riemann invariants and the generalized Riemannproblem, Numer. Math. 106 (2007) 369-425.[5] M. Ben-Artzi, J. Li, G. Warnecke, A direct Eulerian GRP scheme for compressible fluid flows, J. Comput.Phys. 218 (2006) 19-43.[6] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases I: small amplitudeprocesses in charged and neutral one-component systems, Phys. Rev. 94 (1954), 511-525.[7] R. Borges, M. Carmona, B. Costa, W. S. Don, An improved weighted essentially non-oscillatory schemefor hyperbolic conservation laws, J. Comput. Phys. 227 (2008) 3191-3211.
8] J.J. Choi, Hybrid spectral difference/embedded finite volume method for conservation laws, J. Comput.Phys. 295 (2015) 285-306.[9] A.J. Christlieb, S. Gottlieb, Z. Grant, D.C. Seal, Explicit strong stability preserving multistage two-derivative time-stepping scheme, J. Sci. Comput. 68 (2016) 914-942.[10] B. Cockburn, C. W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite elementmethod for conservation laws II: general framework, Mathematics of Computation, 52 (1989) 411-435.[11] B. Cockburn, C. W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws V:multidimensional systems, J. Comput. Phys. 141 (1998) 199-224.[12] X. G. Deng, H. X. Zhang, Developing high-order weighted compact nonlinear schemes, J. Comput.Phys. 165 (2000) 22-44.[13] M. Dumbser, D.S. Balsara, E.F. Toro, C.D. Munz, A unified framework for the construction of one-stepfinite volume and discontinuous Galerkin schemes on unstructured meshes, J. Comput. Phys. 227 (2008)8209-8253.[14] U. Ghia, K. N. Ghia, C.T Shin, High-Re solutions for incompressible flow using the Navier-Stokesequations and a multigrid method, J. Comput. Phys. 48 (1982) 387-411.[15] E. Hairer, G. Wanner, Multistep multistage multiderivative methods for ordinary differential equations,Computing 11 (1973) 287-303.[16] A. Harten, High resolution schemes for hyperbolic conservation laws, J. Comput. Phys. 49 (1983)357-393.[17] A. Harten, B. Engquist, S. Osher, S.R. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes, III, J. Comput. Phys. 71 (1987) 231-303.[18] 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.[19] D.J. Hill, D.I. Pullin, Hybrid tuned center-difference-WENO method for large eddy simulations in thepresence of strong shocks, J. Comput. Phys. 194 (2004) 435-450.[20] O. Inoue and Y. Hattori, Sound generation by shock-vortex interactions, J. Fluid Mech. 380 (1999)81-116.[21] X. Ji, L. Pan, W. Shyy, K. Xu, A compact fourth-order gas-kinetic scheme for the Euler and Navier-Stokes equations J. Comput. Phys. 372 (2018) 446-472.[22] X. Ji, F.X. Zhao, W. Shyy, K. Xu, A family of high-order gas-kinetic schemes and its comparison withRiemann solver based high-order methods, J. Comput. Phys. 356 (2018) 150-173.[23] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys. 126(1996) 202-228.[24] P. D. Lax, X.D. Liu, Solution of two-dimensional riemann problems of gas dynamics by positive schemes,SIAM J. Sci. Comput. 19 (1998) 319-340.[25] P.D. Lax, B. Wendroff, Systems of conservation laws, Commun. Pure Appl. Math. 13 (1960) 217-237.[26] S.K. Lele, Compact finite difference schemes with spectral-like resolution, J. Comput. Phys. 103 (1992)16-42.[27] J. Li, Z. Du, A two-stage fourth order time-accurate discretization for Lax-Wendroff type flow solversI. hyperbolic conservation laws, SIAM J. Sci. Comput. 38 (2016) 3046-3069.[28] J. Li, Q.B. Li, K. Xu, Comparison of the generalized Riemann solver and the gas-kinetic scheme forinviscid compressible flow simulations, J. Comput. Phys. 230 (2011) 5080-5099.[29] X.-D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115(1994) 200-212.[30] H. Luo, J.D. Baum, R. Lohner, A Hermite WENO-based limiter for discontinuous Galerkin method onunstructured grids, J. Comput. Phys. 225 (2007) 686-713.[31] H. Luo, Y. Xia, S. Li, R. Nourgaliev, C. Cai, A Hermite WENO reconstruction-based discontinuousGalerkin method for the Euler equations on tetrahedral grids, J. Comput. Phys. 231 (2012) 5489-5503.[32] H. Luo, Y. Xia, S. Spiegel, R. Nourgaliev, Z. Jiang, A reconstructed discontinuous Galerkin methodbased on a Hierarchical WENO reconstruction for compressible flows on tetrahedral grids, J. Comput.Phys. 236 (2013) 477-492.
33] J. Luo, K. Xu, A high-order multidimensional gas-kinetic scheme for hydrodynamic equations, SCI-ENCE CHINA Technological Sciences. 56 (2013) 2370-2384.[34] K. Mahesh, A family of high order finite difference schemes with good spectral resolution, J. Comput.Phys. 145 (1998) 332-358.[35] L. Pan, J.X. Cheng, S.H. Wang, K. Xu, A two-stage fourth-order gas-kinetic scheme for compressiblemulticomponent flows, Commun. Comput. Phys. 22 (2017) 1123-1149.[36] L. Pan, J.Q. Li, K. Xu, A few benchmark test cases for higher-order Euler solvers, Numer. Math. Theor.Meth. Appl. 10 (2017) 711-736.[37] L. Pan, K. Xu, A compact third-order gas-kinetic scheme for compressible Euler and Navier-Stokesequations, Commun. Comput. Phys. 18 (2015) 985-1011.[38] L. Pan, K. Xu, A third-order compact gas-kinetic scheme on unstructured meshes for compressibleNavier-Stokes solutions, J. Comput. Phys. 318 (2016) 327-348.[39] L. Pan, K. Xu, Q. Li, J. Li, An efficient and accurate two-stage fourth-order gas-kinetic scheme for theNavier-Stokes equations, J. Comput. Phys. 326 (2016) 197-221.[40] J. Qiu, C.-W. Shu, Hermite WENO schemes and their application as limiters for Runge-Kutta discon-tinuous Galerkin method: onedimensional case, J. Comput. Phys. 193 (2004) 115-135.[41] J. Qiu, C.-W. Shu, Hermite WENO schemes and their application as limiters for Runge-Kutta discon-tinuous Galerkin method II: Two dimensional case, Comput. Fluids. 34 (2005) 642-663.[42] D.C. Seal, Y. Guclu, A.J. Christlieb, High-order multiderivative time integrators for hyperbolic con-servation laws, J. Sci. Comput. 60 (2014) 101-140.[43] W.H. Reed, T.R. Hill, Triangular mesh methods for the neutron transport equation, Technical ReportLA-UR-73-479 (1973), Los Alamos Scientific Laboratory, Los Alamos.[44] Y.-X. Ren, M. Liu, H. Zhang, A characteristic-wise hybrid compact-WENO scheme for solving hyper-bolic conservation laws, J. Comput. Phys. 192 (2003) 365-386.[45] C.-W. Shu, High order weighted essentially nonoscillatory schemes for convection dominated problems,SIAM Rev. 51 (2009) 82-126.[46] C.-W. Shu, High order WENO and DG methods for time-dependent convection-dominated PDEs: Abrief survey of several recent developments, J. Comput. Phys. 316 (2016) 598-613.[47] C. W. Shu, S. Osher, Efficient implementation of essentially nonoscillatory shock-capturing schemes II,J. Comput. Phys. 83 (1989) 32-78.[48] E.M. Taylor, M. Wu, M.P. Martn, Optimization of nonlinear error for weighted essentially non-oscillatory methods in direct numerical simulations of compressible turbulence, J. Comput. Phys. 223(2007) 384-397.[49] V. A. Titarev and E. F. Toro, Finite volume WENO schemes for three-dimensional conservation laws,J. Comput. Phys. 201 (2014) 238-260.[50] E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, third edition, Springer (2009).[51] P. Woodward and P. Colella, Numerical simulations of two-dimensional fluid flow with strong shocks,J. Comput. Phys. 54 (1984) 115-173.[52] K. Xu, A slope-update scheme for compressible flow simulation, J. Comput. Phys. 178 (2002) 252-259.[53] K. Xu, A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificialdissipation and Godunov method, J. Comput. Phys. 171 (2001) 289-335.[54] K. Xu, Direct modeling for computational fluid dynamics: construction and application of unified gaskinetic schemes, World Scientific (2015).[55] K. Xu, J.C. Huang, A unified gas-kinetic scheme for continuum and rarefied flows, J. Comput. Phys.229 (2010) 7747-7764.[56] L. Zhang, W. Liu, L. He, X. Deng, H. Zhang, A class of hybrid DG/FV methods for conservation lawsI: Basic formulation and onedimensional systems, J. Comput. Phys. 231 (2012) 1081-1103.[57] T. Zhang, Y. Zheng, Conjecture on the structure of solutions of the Riemann problem for two-dimensional gas dynamics systems, SIAM J. Math. Anal. 21 (1990) 593-630.[58] G. Zhou, K. Xu, and F. Liu, Simplification of the flux function for a high-order gas-kinetic evolutionmodel, J. Comput. Phys. 339 (2017) 146-162.
59] J. Zhu, J. Qiu, Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuousGalerkin method, III: unstructured meshes, J. Sci. Comput. 39 (2009) 293-321.59] J. Zhu, J. Qiu, Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuousGalerkin method, III: unstructured meshes, J. Sci. Comput. 39 (2009) 293-321.