# A Third-order Compact Gas-kinetic Scheme on Unstructured Meshes for Compressible Navier-Stokes Solutions

###### Abstract

In this paper, for the first time a compact third-order gas-kinetic scheme is proposed on unstructured meshes for the compressible viscous flow computations. The possibility to design such a third-order compact scheme is due to the high-order gas evolution model, where a time-dependent gas distribution function at a cell interface not only provides the fluxes across a cell interface, but also the time evolution of the flow variables at the cell interface as well. As a result, both cell averaged and cell interface flow variables can be used for the initial data reconstruction at the beginning of next time step. A weighted least-square reconstruction has been used for the construction of a third-order initial condition. Therefore, a compact third-order gas-kinetic scheme with the involvement of neighboring cells only can be developed on unstructured meshes. In comparison with other conventional high-order schemes, the current method avoids the use of Gaussian points for the flux integration along a cell interface and the multi-stage Runge-Kutta time stepping technique. The third-order compact scheme is numerically stable under CFL condition above . Due to the multidimensional gas-kinetic formulation and the coupling of inviscid and viscous terms, even with unstructured meshes the boundary layer solution and the vortex structure can be accurately captured in the current scheme. At the same time, the compact scheme can capture strong shocks as well.

###### keywords:

high-order scheme, gas-kinetic scheme, compact reconstruction, unstructured mesh, weighted least-square reconstruction.## 1 Introduction

In computational fluid dynamics, the second-order methods are generally robust and reliable, and they are routinely employed in the practical calculations. For the same computational cost, higher-order methods can provide more accurate solutions, but they are less robust and more complicated. In recent decades, there has been a continuous interesting and effort on the development of higher-order schemes. For engineering applications, the construction of higher-order numerical schemes on unstructured meshes becomes extremely demanding. Since a gigantic amount of publications have been devoted to the introduction and survey of higher-orders schemes, the current paper will mainly concentrate on the construction of the third-order compact gas-kinetic scheme on unstructured meshes.

The gas-kinetic scheme (GKS) has been developed systematically for the compressible flow computations GKS-Xu1 (); GKS-Xu2 (); GKS-Kumar (); GKS-Jiang (). An evolution process from kinetic to hydrodynamic scales has been constructed for the flux evaluation. The kinetic effect through particle free transport contributes to the capturing of the shock wave, and the hydrodynamic effect plays a dominant role for the resolved viscous and heat conducting solutions. In other words, the highly non-equilibrium of the gas distribution function in the discontinuous region provides a physically consistent mechanism for the construction of a numerical shock structure. In this sense, the GKS is close to the methodology of artificial dissipation approach, but with different dissipative mechanism. In smooth flow region, the hydrodynamic scale physics corresponding to the multi-dimensional central difference discretization captures the accurate viscous solutions. Due to the coupling of inviscid and viscous terms in the kinetic formulation, theoretically there is no difficulty for GKS to capture NS solutions in any structure or unstructured mesh. With the discretization of particle velocity space, a unified gas-kinetic scheme (UGKS) has been developed for the flow study in entire flow regime from rarefied to continuum ones UGKS-Xu (); UGKS-Luc (); UGKS-Guo ().

Recently, with the incorporation of high-order initial data reconstruction, a higher-order gas-kinetic schemes has been proposed in GKS-high1 (); GKS-high2 (); GKS-high3 (). The flux evaluation in the scheme is based on the time evolution of flow variables from an initial piece-wise discontinuous polynomials (parabola) around a cell interface, where higher-order spatial and temporal derivatives of a gas distribution function are coupled nonlinearly. The whole curves of discontinuous flow distributions around a cell interface interact through particle transport and collision in the determination of the flux function. Besides the evaluation of the time-dependent flux function across a cell interface, the higher-order gas evolution model also provides an accurate time-dependent solution of flow variables at a cell interface as well. Thus, it is feasible to develop a compact scheme with the consideration of time evolution of both cell averaged and cell interface flow variables. A compact third-order gas-kinetic scheme is proposed for the compressible Euler and Navier-Stokes equations on structure meshes with WENO-type reconstruction GKS-high4 (). However, this reconstruction technique is difficult to be used on unstructured meshes. Therefore, on the unstructured meshes, a weighted least-square reconstruction will be used in this paper. To the third-order accuracy, a quadratic distribution for the flow variables inside each cell needs to be determined. Based on the cell averaged and cell interface values of neighboring cells only, an over-determined linear system is formed. With the least-square solution for the system, the whole flow distribution can be fully determined. The shock detector can be also used as well to switch between higher-order (3rd) and lower order (2nd) reconstructions in different regions. In comparison with traditional schemes, the Gaussian points for the flux evaluation along the cell interface and the multi-stage Runge-Kutta technique are avoided in the current compact method. At the same time, the current third-order compact scheme is stable under the CFL condition .

This paper is organized as follows. In Section 2, the finite volume scheme on the unstructured mesh and third-order GKS are introduced. In section 3, the compact reconstruction on the triangular mesh is presented, and the techniques can be applied to rectangular mesh as well. Section 4 includes numerical examples to validate the current algorithm. The last section is the conclusion.

## 2 Finite volume gas-kinetic scheme

### 2.1 Finite volume scheme

The two-dimensional gas-kinetic BGK equation can be written as BGK-1 (),

(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 , , is the number of internal freedom, i.e. for two-dimensional flows, and is the specific heat ratio.

Based on the Chapman-Enskog expansion of the BGK model, the Euler and Navier-Stokes, Burnett, and Super-Burnett equations can be derived BGK-3 (); GKS-Xu1 (). In the smooth region, the gas distribution function can be expanded as

where . By truncating different orders of , the corresponding macroscopic equations can be derived. For the Euler equations, the zeroth order truncation is taken, i.e. . For the Navier-Stokes equations, the first order truncation is

(3) |

Based on the higher order truncations, the Burnett and super-Burnett eqautions can be obtained.

In the computation, the computational volumes are simply triangles. For a control volume , its boundary is given by three line segments

Thus, taking moments of the kinetic equation Eq.(1) and integrating with respect to time and space, the finite volume scheme can be expressed as

(4) |

where are the conservative variables, are the fluxes across the cell interface in the global coordinate, which is defined as

(5) |

where is the outer normal direction of the cell interface , and the tangential direction is denoted as . Eq.(4) is valid in any scale if the interface flux is properly defined, which is beyond the validity of the Navier-Stokes equations.

According to the coordinate transformation, the local coordinate for the cell interface is expressed as , where and , and the velocities in the local coordinate are given by

(6) |

For the gas distribution function in the local coordinate, and , then the line integral for the gas distribution function over the cell interface can be transformed as

(7) |

Thus, in the computation, the numerical fluxes in the local coordinate are obtained first by taking moments of the gas distribution function in the local coordinate

(8) |

where . According to Eq.(6) and Eq.(7), the fluxes in the global coordinate can be expressed as a combination of the fluxes in the local coordinate

(9) |

With the above numerical fluxes at the cell interface, the flow variables inside each control volume can be updated according to Eq.(4).

### 2.2 Gas-kinetic flux solver

In this section, the numerical flux will be presented in the local coordinate. For simplicity, all notations with tilde will be omitted here after.

In order to simulate the NS solutions, we need to model the interface flux function. For the distribution function at a cell interface, the integral solution of BGK equation Eq.(1) at the cell interface in the local coordinate can be written as

(10) |

where is the location of the cell interface, and are the trajectory of particles, is the initial gas distribution function, and is the corresponding equilibrium state. The target equations to be solved depend on the modeling of the initial condition term.

To construct a multidimensional third-order gas-kinetic solver, the following notations are introduced firstly

where is an equilibrium state. The dependence of these coefficients on particle velocity can be expanded as the following form GKS-Xu2 ()

For the kinetic part of the integral solution Eq.(10), the gas distribution function can be constructed as

(11) |

where is the Heaviside function, and are the initial gas distribution functions on both sides of a cell interface, which have one to one correspondence with the initially reconstructed polynomials of macroscopic flow variables on both sides of the cell interface. To construct a third-order scheme, the Taylor expansion for the gas distribution function in space and time at is expressed as

where . For the Euler equations, and the kinetic part of Eq.(10) can be obtained. For the Navier-Stokes equations, according to Eq.(3) and the notations introduced above, the distribution function is

where are the equilibrium states corresponding to the macroscopic variables given by the reconstruction procedure at both sides of cell interface. Thus, the corresponding kinetic part of Eq.(10) can be written as

(12) |

where are the equilibrium states at both sides of the cell interface, and the coefficients are defined according to the expansion of .

After determining the kinetic part , the equilibrium state in the integral solution Eq.(10) can be constructed as follows

(13) |

where is the equilibrium state located at interface, which can be determined through the compatibility condition Eq.(2)

(14) |

Based on Taylor expansion for the equilibrium state Eq.(2.2), the hydrodynamic part in Eq.(10) can be written as

(15) |

where the coefficients are defined from the expansion of the equilibrium state . The coefficients in Eq.(2.2) and Eq.(12) are given by

Substituting Eq.(2.2) and Eq.(12) into the integral solution Eq.(10), the gas distribution function at the cell interface can be obtained. The superscripts or subscripts of the coefficients in Eq.(12) and Eq.(2.2) are omitted for simplicity and they are determined by the spatial derivatives of macroscopic flow variables and the compatibility condition GKS-high2 () as follows

(16) |

where are the moments of gas distribution function, and defined by

In the following section, with the reconstruction procedure, the conservative value and at the center of cell interface corresponding to the equilibrium and the derivatives in Eq.(16) will be presented.

## 3 Compact reconstruction

This paper focuses on the high-order compact finite volume scheme. In the finite volume type schemes, to achieve higher-order accuracy, a reconstruction for the flow variables with high-order polynomials inside each cell is needed as the initial condition at the beginning of each time step. For the higher-order reconstruction, a large number of stencils is usually needed to determine all degrees of freedom through the WENO or least square techniques un-ENO (); un-WENO2 (); un-WENO3 (); k-exact-1 (); k-exact-3 (). In this section, the reconstruction will be done for the unstructured mesh with a compact stencil, which is shown in Fig.1. For simplicity, the whole reconstruction procedure is performed in a local coordinate relative to a cell interface, such as AB in Fig.1, which is consistent with the evaluation of a time-dependent gas distribution function at the cell interface.

In the gas-kinetic scheme, besides the numerical fluxes, the macroscopic pointwise values at a cell interface in the local coordinate can be obtained by taking moments of the gas distribution function,

(17) |

As shown in the last section, the whole curve of the polynomial of the macroscopic variables will participate the flow evolution, and the spatial and temporal derivatives of the gas distribution function are coupled nonlinearly. This point-wise value at the cell interface Eq.(17) is a solution of the evolution model, which can be used in the reconstruction stage at the beginning of next time step. Thus, in the following subsections, a third-order compact reconstruction will be presented for the unstructured mesh, in which the pointwise values at the cell interface and the cell averaged values shown in Fig.1 are used in the reconstruction.

The macroscopic variables for reconstruction is denoted by . For the smooth flow, the conservative variables will be directly used for reconstruction, i.e. . For the flow with discontinuity, in order to eliminate the spurious oscillation and improve the stability of the scheme, the compact reconstruction is based on the characteristic variables. Denote in the local coordinate. The Jacobian matrix can be diagnoalized by the right eigenmatrix , and the characteristic variables is defined as . For a cell interface, is the right eigenmatrix for and is the averaged conservative value from both side of cell interface. To the third order accuracy, the expansion of the macroscopic variable inside the cell can be expressed as

(18) |

where is the barycenter of , is the cell averaged value for , and

The cell averaged value for the base function over the triangle is denoted as

(19) |

and the point-wise value for the base function at the point is denoted as

(20) |

### 3.1 Initial data reconstruction

In this subsection, the weighted least-square reconstruction will be presented for the initial data reconstruction. As shown in Fig.1, three cell averaged values (blue square) form the neighboring cells and nine point-wise values (red circle) from the cell interface will be used in the weighted least square reconstruction.

For the third order expansion, with the definition of the cell averaged and point-wise values for the base function Eq.(19) and Eq.(20), we have

(21) |

where is the cell averaged value for the neighboring triangle , . For the nine cell interface points , , we have

(22) |

where is the point-wise value of at the point .

To solve the corresponding derivatives for , Eq.(21) and Eq.(22) can be written into an over-determined linear system

(23) |

Denote , , the above linear system is expressed as the matrix form

where is the coefficient matrix corresponding to Eq.(23).

To deal with the discontinuity, a diagonal matrix is introduced as the simple weight functions

where , , and . The derivatives can be obtained by solving the linear system

Generally, for most cases with Mach number , the weight function is enough to deal with the discontinuity. However, for strong discontinuity, the shock detection Shock-detection () technique is used in the current scheme. Analogous to the analysis of KXRCF detector Shock-detection (), for the third-order scheme, it is easy to distinguish the smooth region from the region near discontinuities as follows

where the index refers and the index refers , is the interpolated value at the center of and is the value at the center of extrapolated from . In the computation, the ”trouble cell” is detected according the following criterion

where is the area of the triangle, is a problem dependent coefficient, and is used in the computation. In those detected ”trouble cell”, the second order scheme with limiters are used. The above choice of weight functions may not be optimal and further study is needed.

With the derivatives , the whole flow distribution in the cell in Fig.2 can be obtained. For the smooth flow, no special treatment is needed. With , the interpolated value and the derivatives can be fully obtained in the cell . Similarly, the interpolated value and the derivatives in the cell can be obtained as well.

For the flow with discontinuity, the characteristic variables are reconstructed in the cell . With the derivatives , the interpolated value at the points in Fig.2 can be obtained. By the inverse projection, the conservative variables , where is the right eigenmatrix. Based on these point-wise values and their central difference, and can be obtained. Similarly, the interpolated value and the derivatives in the cell can be also obtained.

### 3.2 Reconstruction for equilibrium part

In this subsection, the reconstruction for the equilibrium part will be presented. This reconstruction will be based on the conservative variables . To the third-order accuracy, the Taylor expansion corresponding to equilibrium part at the center point of a cell interface is expressed as

(24) |

where is the conservative variable at the center point of cell interface based on the compatibility condition Eq.(14), and are corresponding derivatives.

As shown in Fig.3, with the reconstructed polynomials in and , the point values at those points can be determined, which has been obtained in the last subsection. Especially, we can get the point values at the interface (red) points at both sides of . By the compatibility condition Eq.(2), the reconstructed conservative variables at the cell interface can be determined. The derivatives can be obtained by the central difference of these point-wise values.

#### 3.2.1 Extension to rectangular mesh

For the rectangular mesh, the stencils are given in Fig.4. To reconstruct the polynomial for the rectangular , the cell averaged values , and point-wise values , at the cell interfaces can be used. Similar to the triangular case, we have the following matrix form for the over-determined linear system

where , . is the coefficient matrix and expressed as

By introducing the weight diagonal matrix , the derivative can be also obtained by solving the following linear system

The limiting process is also used for the flow with large discontinuity. In some cases of the numerical tests, the solutions from the compact scheme with rectangular mesh will be presented as well.

## 4 Numerical tests

In this section, numerical tests for both inviscid flow and viscous flow will be presented to validate the compact scheme. For the inviscid flow, the collision time takes

where and . For the viscous flow, we have

where and denotes the pressure on the left and right sides of the cell interface, is the viscous coefficient, is the pressure at the cell interface and . In the smooth flow regions, it will reduce to . The ratio of specific heats takes . is the time step which is determined according to the CFL condition. In the numerical tests, the CFL number takes a value of , even though the scheme works as well with a large CFL number. The value of is already more than two times of the time step used for the compact third-order DG method.

### 4.1 Accuracy test

The numerical order of the compact gas-kinetic scheme is tested in comparison with the analytical solutions of the Euler equations. The isotropic vortex propagation problem is presented to validate the accuracy for the solution of inviscid flow. The computational domain is taken to be . The free upstream is , and a small vortex is obtained through a perturbation on the mean flow with the velocity , temperature , and entropy . The perturbation is expressed as

where , , , , , and . In the computation, the unstructured meshes with mesh size and are used, and the errors and orders at are presented in Table.1, which shows a third-order accuracy of the current compact scheme.

### 4.2 One dimensional Riemann problem

In this case, two one-dimensional Riemann problems are tested to verify the capability in capturing the wave configurations. The mesh is presented in Fig.6, where the computational domain is , and mesh size is around . The first one is Sod problem, and the initial condition is given by

The second one is the Lax problem, and the initial condition is given as follows

To compare with the exact solution, points were extracted at for the Sod problem at and, for the Lax problem at . The density, velocity, and pressure distributions for the exact solutions and numerical results are presented in Fig.7, where the numerical results agree well with the exact solutions. The three dimensional density distributions for the two cases are given in Fig.6. In this case, the weighted least square reconstruction can deal with the discontinuity well, and the shock detection technique is not needed.

### 4.3 Flow impinging on a blunt body

In this case, the inviscid hypersonic flows impinging on a unit cylinder are tested to validate robustness of the current scheme. This problem is initialized by the flow moving towards a cylinder with different Mach numbers. The Euler boundary condition is imposed on the surface of cylinder, and outflow boundary condition on the right boundary. As mentioned in the reconstruction part, the weighted least square reconstruction is able to deal with the discontinuities at a Mach number . In this case, the flow with is tested without the detection of ”trouble cell”. The mesh and the pressure distribution for this case are also given in Fig.9, with mesh size , where the flow structure can be captured nicely in front of the cylinder. However, with a high Mach number, the weighted least square reconstruction is no longer able to capture strong discontinuities, and the shock detection technique is used to identify the trouble cells, where a second-order reconstruction is used in these cells. For the flow with , the mesh and the pressure distribution are shown in Fig.9 with mesh size . This test shows that the current scheme can capture the flow structure nicely in front of the cylinder and the carbuncle phenomenon does not appear Case-Pandolfi ().

### 4.4 Shock vortex interaction

The interaction between a stationary shock and a vortex for the inviscid flow is presented WENO2 (). The computational domain is taken to be