TwoStage Fourthorder Gaskinetic Scheme for Threedimensional Euler and NavierStokes Solutions
Abstract
For the onestage thirdorder gaskinetic scheme (GKS), success applications have been achieved for the threedimensional compressible flow computations GKShigh1 (). The highorder accuracy of the scheme is obtained directly by integrating a multidimensional timeaccurate gas distribution function over the cell interface within a time step without implementing Gaussian quadrature points and RungeKutta timestepping technique. However, for the further increasing the order of the scheme, such as the fourthorder one, the formulation becomes very complicated for the multidimensional flow. Recently, a twostage fourthorder GKS with high efficiency has been constructed for twodimensional inviscid and viscous flow computations GRPhigh (); GKShigh (), and the scheme uses the time accurate flux function and its time derivatives. In this paper, a fourthorder GKS is developed for the threedimensional flows under the twostage framework. Based on the threedimensional WENO reconstruction and flux evaluation at Gaussian quadrature points on a cell interface, the highorder accuracy in space is achieved first. Then, the twostage time stepping method provides the high accuracy in time. In comparison with the formal thirdorder GKS GKShigh1 (), the current fourthorder method not only improves the accuracy of the scheme, but also reduces the complexity of the gaskinetic solver greatly. More importantly, the fourthorder GKS has the same robustness as the secondorder shock capturing scheme. This scheme is applied to both inviscid and viscous, and low and high speed flow computations. Numerical results validate the outstanding reliability and applicability of the scheme for threedimensional flows, such as turbulent one. With the count of the degree of freedom in simulating a flow field, the current twostage fourthorder multidimensional GKS becomes one of the most reliable and efficient higherorder schemes among all shock capturing higherorder schemes for computational fluid dynamics.
keywords:
twostage fourthorder discretization, gaskinetic scheme, NavierStokes solutions, shock capturing scheme1 Introduction
Over the past half century, the computational fluid dynamics has been developed into a powerful tool for solving the fluid flow problems in industrial applications. Currently, there are a gigantic amount of numerical methods in literature. In comparison with welldeveloped secondorder schemes, the higherorder methods can provide more accurate solutions, but they are less robust and more complicated. For the highorder schemes, there are generally three parts, i.e. spatial discretization, temporal discretization, and flux solvers. For the spatial discretization, many highorder methods have been developed, including the discontinuous Galerkin (DG) DG1 (); DG2 (); DG3 (), essential nonoscillatory (ENO) ENO1 (); ENO2 (), weighted essential nonoscillatory (WENO) WENO (); WENOJS (); WENOZ (), and many others. For most of those methods, the exact and approximate Riemann solvers Riemannappro () are used for flux evaluation. Due to the firstorder evolution of Riemann solvers, the RungeKutta method is implemented to achieve higher order temporal accuracy TVDrk (), in which stage is needed for th order accuracy. Instead of Riemann solver, many schemes have been developed based on the timedependent flux function with highorder evolution, such as the generalized Riemann problem solvers GRP1 (); GRP2 (); GRP3 (); GRPhigh2 (), gaskinetic scheme GKSXu1 (); GKSXu2 (), and AEDR methods Riemannappro (); ADER (). Highorder temporal accuracy can be achieved in a onestage framework with the time integration of the timedependent flux function. Recently, in order to increase the accuracy and efficiency of these schemes a twostage fourthorder method has been developed for the timedependent flux solvers GRPhigh (); GKShigh (), where both the flux and temporal derivative of flux function are used in the construction of higherorder schemes. By combining the multistage multiderivative technique, a family of higherorder schemes has been constructed as well MSMDGKS ().
In the past decades, the gaskinetic scheme (GKS) based on the BhatnagarGrossKrook (BGK) model BGK1 (); BGK2 (); BGK3 () has been developed systematically for the compressible flow computations GKSXu1 (); GKSXu2 (); GKSXu3 (). Different from the traditional CFD methods based on the macroscopic governing equations, the main advantages of the gaskinetic scheme are the followings. (i) The inviscid and viscous coupling in the flux evolution GKSXu1 (); GKSXu2 (); (ii) Multidimensionality with the inclusion of both normal and tangential gradients of flow variables in the flux function across a cell interface GKShigh4 (); (iii) Compact stencils can be constructed with the use of the time accurate cell interface flow variables at the next time level GKShigh4 (); (iv) Extension to the whole flow regime from the rarefied to the continuum one UGKS1 (); UGKS2 (). Recently, with the highorder initial reconstruction, the thirdorder gaskinetic schemes have been constructed GKShigh2 (); GKShigh3 (); GKShigh4 (), in which the flux evaluation is based on the moments of spatialtemporal dependent gas distribution function. Highorder accuracy can be achieved in a onestage scheme without Gaussian point integration for spatial accuracy and RungeKutta technique for temporal accuracy. However, with the onestage gas evolution model, the formulation of gaskinetic scheme can become very complicated for the further development, such as the fourthorder method GKShigh5 (), especially for multidimensional computations. The twostage fourthorder temporal discretization for timeaccurate flux solvers in GRPhigh () provides a reliable framework to further develop the GKS into fourthorder and even higher accuracy with the implementation of the traditional secondorder or thirdorder flux functions GKShigh (); GKShigh6 (); MSMDGKS (). Most importantly, this scheme is robust, and works perfectly from the subsonic to the hypersonic viscous heat conducting flows. The robustness of the GKS in comparison with Riemann solver based CFD methods is solely due to the differences in the dynamical evolution model of the flux function. The GKS flux follows the dynamics from the particle free transport, to including collisional effect, and to the NS distribution function with the count of intensive particle collisions as the ratio of the time step over the particle collision time becoming large. The real physics used in the flux depends on the local xuliu (). In real NS computations for the compressible high Mach number flow, the ratio of is not too large as people think of for the validation of the NS modeling liukn (). However, for the Riemann solver based CFD methods, at the beginning of the step it is already assumed that there are infinity number of particle collision to generate distinguishable waves in the Riemann solution, and the collision needs to be reduced for the NS solutions. Theoretically, there is no such a physical process for the Riemann solverbased Godunovtype schemes for the NS equations. For the secondorder schemes, it is hard to distinguish the dynamical differences from the GKS and Riemann solver. However, for the higherorder schemes it seems that a reliable physical evolution model becomes more important due to the absence of large numerical dissipation in the secondorder schemes, and the delicate flow structures captured in higherorder schemes depend on the quality of the solvers greatly MSMDGKS ().
In this paper, with the twostage fourthorder discretization, a multidimensional fourthorder gaskinetic scheme is constructed for simulating threedimensional flows. Highorder accuracy in space is achieved by the threedimensional WENO method WENO (); WENOJS (); WENOZ () and Gaussian quadrature points for the numerical fluxes. In comparison with the formal threedimensional scheme GKShigh1 (), the current fourthorder scheme reduces the complexity of the gaskinetic flux solver greatly, and improves the robustness of scheme. Many numerical tests, including both inviscid and viscous, and low and high speed flow computations, will be used to validate the current fourthorder method. Numerical results show that the current scheme has the same reliability and applicability as the welldeveloped secondorder scheme, but is much more accurate and effective to capture the complicated flow structures. With the count of the degrees of freedom for the description of a flow field, the current scheme provides stateofart solutions from a higherorder scheme from the incompressible to the hypersonic flow simulations.
This paper is organized as follows. In Section 2, BGK equation and finite volume scheme are briefly reviewed. The general formulation for the twostage temporal discretization is introduced in Section 3, and the procedure of spatial reconstruction is given in Section 4. Section 5 includes numerical examples to validate the current algorithm. The last section is the conclusion.
2 BGK equation and finite volume scheme
The threedimensional BGK equation BGK1 (); BGK2 (); BGK3 () can be written as
(1) 
where is the gas distribution function, is the corresponding equilibrium state, and is the collision time. The collision term satisfies the compatibility condition
(2) 
where , the internal variables is equal to , , is the degrees of freedom, i.e. for threedimensional flows and is the specific heat ratio. In the continuum region, the gas distribution function can be expanded as
where . Based on the ChapmanEnskog expansion, the macroscopic equations can be derived GKSXu1 (); GKSXu2 (). For the Euler equations, the zerothorder truncation is taken, i.e. . For the NavierStokes equations, the firstorder truncation is
With the higher order truncations, the Burnett and superBurnett equations can be derived.
Taking moments of Eq.(1) and integrating over the control volume with , the semidiscretized form of finite volume scheme can be written as
(3) 
where are the conservative flow variables, is the cell averaged value over control volume and . For the threedimensional computation, the Gaussian quadrature for the numerical fluxes is used to achieve the accuracy in space, and the numerical fluxes in direction is given as an example
(4) 
where , are the Gauss quadrature points for and are corresponding weights, and is used in this paper. Based on the spatial reconstruction, which will be presented in the following section, the reconstructed point value and the spatial derivatives at each Gauss quadrature points can be obtained and the numerical fluxes can be provided by the flow solvers. Similarly, the numerical fluxes in the and directions can be obtained as well
In the gaskinetic scheme, the numerical fluxes at the Gauss quadrature point can be obtained by taking moments of the gas distribution function
(5) 
where is provided by the integral solution of BGK equation Eq.(1) at the cell interface
(6) 
and is the location of cell interface, is the particle velocity, is the trajectory of particles. For the secondorder scheme, the secondorder gaskinetic solver for the threedimensional flows GKSXu2 (); GKSXu1 () can be expressed as
(7) 
where the coefficients in Eq.(2) can be determined by the spatial derivatives of macroscopic flow variables and the compatibility condition as follows
where the superscripts or subscripts of these coefficients are omitted for simplicity, more details about the determination of coefficient can be found in GKSXu2 (). With the secondorder gaskinetic solver Eq.(2), the secondorder accuracy in time can be achieved by one step integration. In the onestage gas evolution model, the thirdorder gaskinetic solver has been developed as well GKShigh1 (); GKShigh2 (); GKShigh3 (); GKShigh4 ().
3 Fourthorder temporal discretization
In the further development of higherorder scheme, the formulation of onestage gaskinetic solver can become very complicated, such as the fourthorder method GKShigh5 (), especially for multidimensional computations. Recently, a twostage fourthorder timeaccurate discretization was developed for LaxWendroff flow solvers, particularly applied for hyperbolic equations with the generalized Riemann problem (GRP) solver GRPhigh () and gaskinetic scheme GKShigh (). Such method provides a reliable framework to develop a threedimensional fourthorder gaskinetic scheme with a secondorder flux function Eq.(2). Consider the following timedependent equation
with the initial condition at , i.e.,
where is an operator for spatial derivative of flux. The time derivatives are obtained using the CauchyKovalevskaya method,
Introducing an intermediate state at ,
(8) 
the corresponding time derivatives are obtained as well for the intermediate stage state,
Then, the state w can be updated with the following formula,
(9) 
It can be proved that for hyperbolic equations the above time stepping method Eq.(8) and Eq.(9) provides a fourthorder time accurate solution for at . The details of the analysis can be found in GRPhigh ().
The semidiscretized form of finite volume scheme Eq.(2) can be applied for the above twostage method . For the gaskinetic scheme, the gas evolution is a relaxation process from kinetic to hydrodynamic scale through the exponential function, and the corresponding flux is a complicated function of time. In order to obtain the time derivatives of the flux function at and with the correct physics, the flux function should be approximated as a linear function of time within a time interval. According to the numerical fluxes at the Gauss quadrature points Eq.(5), the following notation is introduced
In the time interval , the flux is expanded as the following linear form
The coefficients and can be determined as follows,
By solving the linear system, we have
Similarly, for the intermediate state can be constructed. The corresponding fluxes in the  and direction can be obtained as well. With these notations, the numerical scheme for the threedimensional flows can be constructed, and the detail procedure of twostage algorithm can be found in GKShigh ().
4 Spatial reconstruction
The above time evolution solution is based on the highorder initial reconstruction for macroscopic flow variables and WENO reconstruction WENO (); WENOJS (); WENOZ () is adopted for the spatial reconstruction. For the three dimensional computation, the reconstruction procedure for the cell interface is given as an example, and the stencil for reconstruction is given in Fig.1. The point value and and firstorder derivatives at the Gauss quadrature points , need to be constructed. The detailed procedure is given as follows

According to one dimensional reconstruction, the cell averaged reconstructed value and cell averaged spatial derivatives can be constructed for the stencil shown in Fig.1.

With the onedimensional WENO reconstruction in the horizontal direction, the averaged value and the averaged spatial derivatives and over the interval with can be given.

With onedimensional WENO reconstruction in the vertical direction, the point value and spatial derivatives , and can be fully given at the Gaussian quadrature points .
In the computation, without special statement, the fifthorder WENOJS reconstruction WENOJS () is adopted for the flow with discontinuities and the linear scheme is used for the smooth flows to reduce the dissipation.
5 Numerical experiments
In this section, numerical tests for both inviscid and viscous flows will be presented to validate our numerical scheme. For the inviscid flow, the collision time takes
where and denote the pressures on the left and right sides of the cell interface. In the computation, and . For the viscous flow, we have
where is the viscous coefficient and is the pressure at the cell interface, and it will reduce to in the smooth flow regions. is the time step which is determined according to the CFL number, which takes in the computation. The reason for including artificial dissipation through the additional term in the particle collision time is to enlarge the kinetic scale physics in the discontinuous region for the construction of a numerical shock structure through the particle free transport and inadequate particle collision in order to keep the nonequilibrium property.
For the smooth flow, the WENO reconstruction can be used directly on the conservative flow variables. For the flow with strong discontinuity, the characteristic variables can be used in the reconstruction. Based on , where are the conservative variables, are the corresponding fluxes, and , the cell averaged and point conservative values can be projected into the characteristic field by , where is the matrix corresponding to right eigenvectors of . The reconstruction scheme is applied on the characteristic variables . With the reconstructed polynomials for characteristic variables, the conservative flow variables can be recovered by the inverse projection.
In the following, we present many numerical examples. Based on the same WENOJS reconstruction, the results presented below have very high accuracy in comparison with other solutions from same reconstruction.
5.1 Threedimensional Sod problem
This problem is a fully threedimensional extension of the Sod problem. The computational domain is , and a sphere of radius separates two different constant states
The exact solution of this problem can be given by the onedimensional system with geometric source terms
with
where the radial direction is denoted by , while is the radial velocity, is the number of space dimensions and . The twostage fourthorder gaskinetic scheme is used to solve this equation and the numerical results with cells are given as the reference solutions. In the computation, the uniform mesh with is used. The symmetric boundary condition is imposed on the plane with , and , and the nonreflection boundary condition is imposed on the plane with , and . The threedimensional density and pressure distributions are given in Fig.3, and the density and pressure distribution along the line are given in Fig.3. The numerical results agree well with the reference solutions.
5.2 Flow impinging on sphere
In this case, the inviscid hypersonic flows impinging on a unit sphere are tested to validate robustness of the current scheme with different Mach numbers with . In the computation, a mesh shown in Fig.4 is used, which representing the domain in the spherical coordinate . The density distributions in the whole domain, Mach number and pressure distributions at the plane with , and pressure profiles along the surface of sphere with for Mach number and are shown in Fig.4, Fig.6 and Fig.6, where the shock is well captured by the current scheme and the carbuncle phenomenon does not appear CasePandolfi (). As comparison, the inviscid hypersonic flows impinging on a twodimensional unit cylinder is also tested with mesh. The Mach number and pressure distributions for Mach number are shown in Fig.7.
5.3 RayleighTaylor instability
The RayleighTaylor instability occurs when an interface between two different fluids with different densities experiences a pressure gradient opposing the density gradient Casert (). This interface is unstable and any disturbance tends to grow, leading to the penetration of both fluids into each other. If the initial interface displacement is random, the RayleighTaylor instability usually evolves into complicated turbulent mixing. In this case, the threedimensional RayleighTaylor instability is tested in a rectangular box with domain with a gravity pointed downward . Boundary conditions are applied at the four sides, while symmetric boundary conditions are applied at the top and bottom walls. The instability develops from the imposed single mode initial perturbation
where is the height of the interface and is the box width. The dimensionless parameters is the Atwood number
where are densities of heavy and light fluids, respectively. In the computation, the Atwood number is used. The the fluid interface from a singlemode perturbation at are given in Fig.8, and the density distribution at , and with are given in Fig.9. As expected, the heavy and light fluids penetrate into each other as time increases. The light fluid rises to form a bubble and the heavy fluid falls to generate a spike.
5.4 Liddriven cavity flow
The liddriven cavity problem is one of the most important benchmarks for numerical NavierStokes solvers. As shown in Fig.10, the fluid is bounded by a unit cubic and driven by a uniform translation of the top boundary. The monatomic gas with is used, such that there is no bulk viscosity involved. Early threedimensional cavityflow calculations were carried out by De Vahl Davis and Mallinson CaseDavis () and Goda CaseGoda (). In this case, the flow is simulated with Mach number and all the boundaries are isothermal and nonslip. Numerical simulations are conducted with three Reynolds numbers of and using meshes for the domain . The velocity profiles along the vertical centerline line, velocity profiles along the horizontal centerline in the plane with and the benchmark data CaseAlbensoeder (); CaseShu () are shown in Fig.11. The simulation results match well with the benchmark data.
5.5 Compressible homogeneous turbulence
The highorder gaskinetic scheme is applied for the direct numerical simulations (DNS) of the compressible decaying homogeneous isotropic turbulence. The flow is computed within a square box defined as , and the periodic boundary conditions are used in all directions for all the flow variables DNS1 (); DNS2 (); DNS3 (). In the computation, the domain is discretized with a uniform Cartesian mesh cells . A divergencefree random initial velocity field is generated for a given spectrum with a specified root mean square
where is a volume average over the whole computational domain. The specified spectrum for velocity is given by
where is the wave number, is the wave number at spectrum peaks, is a constant chosen to get a specified initial kinetic energy. The initial volume averaged turbulent kinetic energy and the initial largeeddyturnover time is given by
The Taylor microscale Reynolds number and turbulence Mach number are given as
where is Taylor microscale
and can be determined from and with initialized and , and the dynamic viscosity is determined by
The time history of the kinetic energy, rootmeansquare of density fluctuation and skewness factor for velocity slope are defined as
In the computation, , and , and the uniform meshes with , and cells are used. The isosurfaces of criterion colored by velocity magnitude and the pressure distribution with at time are given in Fig.12. The time history of normalized kinetic energy , normalized rootmeansquare of density fluctuation and skewness factor with respect to are given in Fig.13. The numerical results agree well with the reference data. With fixed initial , the cases with are tested, and the time histories of normalized kinetic energy are given in Fig.13 as well. With the increase of , the dynamic viscosity increases, and the kinetic energy gets dissipated more rapidly.
5.6 TaylorGreen Vortex
This problem is aimed at testing the performance of highorder methods on the direct numerical simulation of a threedimensional periodic and transitional flow defined by a simple initial condition, i.e. the TaylorGreen vortex CaseBull (); CaseDebonis (). With a uniform temperature field, the initial flow field is given by
The fluid is then a perfect gas with and the Prandtl number is . Numerical simulations are conducted with two Reynolds numbers and . The flow is computed within a periodic square box defined as . The characteristic convective time . In the computation, , and the Mach number takes , where is the sound speed.
The volumeaveraged kinetic energy can be computed from the flow as it evolves in time, which is expressed as
where is the volume of the computational domain, and the dissipation rate of the kinetic energy is given by
The numerical results of the current scheme with mesh points for the normalized volumeaveraged kinetic energy and dissipation rate with Reynolds numbers and are presented in Fig.15 and Fig.17, which agree well with the data in CaseDebonis (); Casewang (). The isosurfaces of criterions colored by velocity magnitude at and are shown in Fig.15 for and in Fig.17 for . The evolution of flow structure is evident, starting from large vortices and decaying into more complex structures. Different from many other higherorder methods, the current scheme has no internal degrees of freedom to be updated within each cell.
6 Conclusion
In this paper, based on the twostage time stepping method, a fourthorder gaskinetic scheme is proposed for the threedimensional inviscid and viscous flow computations. With the threedimensional WENOJS reconstruction, a gaskinetic scheme with higherorder spatial and temporal accuracy is developed. In comparison with the classical methods based on the firstorder Riemann solver, for the same fourthorder accuracy in time the current scheme only uses two stages instead of four stages with the RungeKutta timestepping technique. As a result, the twostage GKS can be more efficient than the fourstage methods with the absence of two time consuming reconstructions. For the NavierStokes solutions, the current scheme doesn’t separate inviscid and viscous terms and they are treated uniformly from the same initial WENOtype reconstruction. The GKS can present very accurate viscous flow solutions due to its multidimensional flux function at a cell interface, where the gradients in both normal and tangential directions of flow variables participate in the gas evolution. The fourthorder GKS not only has the expected order of accuracy for the smooth flow, but also has favorable shock capturing property for the discontinuous solutions. Most importantly, the numerical tests clearly demonstrate that the current fourthorder scheme has the same robustness as the secondorder one. The scheme has been tested from the smooth flows to the flows with discontinuities, and from the low speed to the hypersonic ones. For the three dimensional NavierStokes solutions, the current scheme is one of the stateofart methods in the capturing of complicated flow structures.
Acknowledgements
The work of Pan is supported by the grants from NSFC (11701038) and China Postdoctoral Science Foundation (2016M600065), and the work of Xu is supported by Hong Kong research grant council (16206617, 16207715, 16211014) and the grants from NSFC (11772281, 91530319).
References
 (1) S. Albensoeder, H.C. Kuhlmann, Accurate threedimensional liddriven cavity flow, J. Comput. Phys. 206 (2005) 536558.
 (2) M. BenArtzi, J. Falcovitz, A secondorder Godunovtype scheme for compressible uid dynamics, J. Comput. Phys. 55 (1984), 132.
 (3) M. BenArtzi, J. Li, G. Warnecke, A direct Eulerian GRP scheme for compressible fluid flows, J. Comput. Phys. 218 (2006), 1943.
 (4) M. BenArtzi, J. Li, Hyperbolic conservation laws: Riemann invariants and the generalized Riemann problem, Numerische Mathematik, 106 (2007), 369425.
 (5) P.L. Bhatnagar, E.P. Gross, M. Krook, A Model for Collision Processes in Gases I: Small Amplitude Processes in Charged and Neutral OneComponent Systems, Phys. Rev. 94 (1954) 511525.
 (6) R. Borges, M. Carmona, B. Costa, W. S. Don, An improved weighted essentially nonoscillatory scheme for hyperbolic conservation laws, J. Comput. Phys. 227 (2008) 31913211.
 (7) J. R. Bull, A. Jameson, Simulation of the compressible TaylorGreen vortex using highorder flux reconstruction schemes, AIAA 20143210.
 (8) C. Cercignani, The Boltzmann Equation and its Applications, SpringerVerlag, (1988).
 (9) S. Chapman, T.G. Cowling, The Mathematical theory of NonUniform Gases, third edition, Cambridge University Press, (1990).
 (10) B. Cockburn, C. W. Shu, TVB RungeKutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework, Mathematics of Computation, 52 (1989), 411435.
 (11) B. Cockburn, C. W. Shu, The RungeKutta discontinuous Galerkin method for conservation laws V: multidimensional systems, J. Comput. Phys. 141 (1998), 199224.
 (12) J.Debonis, Solutions of the TaylorGreen vortex problem using highresolution explicit finite difference methods, AIAA Paper (2013) 20130382.
 (13) G. De Vahl Davis, G.D. Mallinson, An evaluation of upwind and central difference approximations by a study of recirculating flow, Comp. Fluids. 4 (1976) 2943.
 (14) Z. Du, J. Li, A Hermite WENO reconstruction for fourth order temporal accurate schemes based on the GRP solver for hyperbolic conservation laws J. Comput. Phys. 355 (2018) 385396
 (15) S. Gottlieb, C. W. Shu, Total variation diminishing rungekutta schemes, Mathematics of computation, 67 (1998) 7385.
 (16) K. Goda, A multistep technique with implicit difference schemes for calculating two or threedimensional cavity flows, J. Comput. Phys. 30 (1979) 7695.
 (17) A. Harten, B. Engquist, S. Osher and S. R. Chakravarthy, Uniformly high order accurate essentially nonoscillatory schemes, III. J. Comput. Phys. 71 (1987) 231303.
 (18) X.Y. He, R.Y Zhang, S.Y. Chen, G. D. Doolen, On the threedimensional RayleighTaylor instability, Phys Fluids 11 (1999) 1143.
 (19) J. Huang, K. Xu, P. Yu, A unified gaskinetic scheme for continuum and rarefied flows II: multidimensional cases, Commun. Comput. Phys. 12 (2012) 662690.
 (20) X. Ji, F.X. Zhao, W. Shyy, K. Xu, A Family of Highorder Gaskinetic Schemes and Its Comparison with Riemann Solver Based Highorder Methods, J. Comput. Phys. 356 (2018), 150¨C173.
 (21) G.S. Jiang, C. W. Shu, Efficient implementation of Weighted ENO schemes, J. Comput. Phys. 126 (1996) 202228.
 (22) J. Li, Z. Du, A twostage fourth order timeaccurate discretization for LaxWendroff type flow solvers I. hyperbolic conservation laws, SIAM J. Sci. Computing, 38 (2016) 30463069.
 (23) Q.B. Li, S. Fu, Highorder accurate gaskinetic scheme and turbulence simulation, SCIENTIA SINICA Physica, Mechanica Astronomica, 44 (2014) 278284.
 (24) Q. Li, K. Xu, S. Fu, A highorder gaskinetic NavierStokes flow solver, J. Comput. Phys. 229 (2010) 67156731.
 (25) W. Liao, Y. Peng, L.S. Luo, Gaskinetic schemes for direct numerical simulations of compressible homogeneous turbulence. Phys Rev E, 80 (2009) 046702.
 (26) C. Liu, K. Xu, G.Z. Zhou, Cell Size Effect on Computational Fluid Dynamics: The Limitation Principle for Flow Simulation, arXiv:1706.08200v4 [physics.compph].
 (27) J. Luo, K. Xu, A highorder multidimensional gaskinetic scheme for hydrodynamic equations, SCIENCE CHINA Technological Sciences, 56 (2013) 23702384.
 (28) N. Liu, H.Z. Tang, A highorder accurate gaskinetic scheme for one and twodimensional flow simulation, Commun. Comput. Phys. 15 (2014) 911943.
 (29) X.D. Liu, S. Osher, T. Chan, Weighted essentially nonoscillatory schemes, J. Comput. Phys. 115 (1994) 200212.
 (30) L. Pan, J. Li, K. Xu, A few benchmark test cases for higherorder Euler solvers, Numerical Mathematics: Theory, Methods and Applications, 10 (2017) 711736.
 (31) L. Pan, K. Xu, A thirdorder compact gaskinetic scheme on unstructured meshes for compressible NavierStokes solutions, Journal of Computational Physics, 318 (2016) 327348.
 (32) L. Pan, K. Xu, Q. Li, J. Li, An efficient and accurate twostage fourthorder gaskinetic scheme for the NavierStokes equations, J. Comput. Phys. 326 (2016), 197221.
 (33) L. Pan, K. Xu, A thirdorder gaskinetic scheme for threedimensional inviscid and viscous flow computations, Computers Fluids 119 (2015) 250260.
 (34) M. Pandolfi, and D. D’Ambrosio, Numerical Instabilities in Upwind Methods: Analysis and Cures for the ”Carbuncle” Phenomenon, J. Comput. Phys. 166 (2001) 271301.
 (35) W.H. Reed, T.R. Hill, Triangular mesh methods for the neutron transport equation, Technical Report LAUR73479, 1973, Los Alamos Scientific Laboratory, Los Alamos.
 (36) R. Samtaney, D.I. Pullin, B. Kosovic, Direct numerical simulation of decaying compressible turbulence and shocklet statistics. Phys Fluids, 13 (2001) 14151430.
 (37) C. Shu, L. Wang, Y.T. Chew, Numerical computation of threedimensional incompressible Navier.Stokes equations in primitive variable form by DQ method, Int. J. Numer. Meth. Fluids. 43 (2003) 345368.
 (38) C. W. Shu, S. Osher, Efficient implementation of essentially nonoscillatory shock capturing schemes, J. Comput. Phys. 77 (1988) 439471.
 (39) V.A. Titarev, E.F. Toro, ADER: arbitrary high order Godunov approach, J. Sci. Comput. 17 (2002) 609¨C618.
 (40) E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, (1997).
 (41) L. Wang, W. K. Anderson, Erwin, S. Kapadia, HighOrder Discontinuous Galerkin Method for Computation of Turbulent Flows, AIAA Journal, 53 (2015) 11571171.
 (42) K. Xu, Direct modeling for computational fluid dynamics: construction and application of unfied gas kinetic schemes, World Scientific (2015).
 (43) K. Xu, A gaskinetic BGK scheme for the NavierStokes equations and its connection with artificial dissipation and Godunov method, J. Comput. Phys. 171 (2001) 289335.
 (44) K. Xu, Gas kinetic schemes for unsteady compressible flow simulations, Lecure Note Ser. 199803, Von Karman Institute for Fluid Dynamics Lecture (1998).
 (45) K. Xu, J. Huang, A unified gaskinetic scheme for continuum and rarefied flows, J. Comput. Phys. 229 (2010) 77477764.
 (46) K. Xu, C. Liu, A paradigm for modeling and computation of gas dynamics, Phys. Fluids 29 026101 (2017).