An improved IP scheme for compressible Navier-Stokes equation

Based on the internal penalty discontinuous Galerkin method (IP-DG), an improved discretization method is proposed to solve the viscosity term of Navier-Stokes equation. The improved method can eliminate the calculation process of intermediate variables and variable gradients, and effectively improve the calculation efficiency. Because the calculation process of convection flux, viscous flux, and lifting operators can be performed at the current integration point, there is no need to calculate the derivative of the viscosity term with respect to the variable gradient, the improved method has very good local characteristics, and is easy to extend to other PDEs. The Couette flow verifies the design accuracy of the improved method, and the cavity flow problem proves the effectiveness of the method for viscous flow simulation. The example of the flow around a sphere shows that this method also has a good performance for three-dimensional problems.


Introduction
The high-order precision format has small numerical dispersion and dissipation, and is suitable for solving multi-scale flow problems such as turbulence and aeroacoustics problems. In recent years, highorder methods have been highly valued by CFD reachers. Among many high-order calculation methods, the discontinuous Galerkin (DG) method has high calculation accuracy and easy implementation of parallel algorithms. Besides, it has no need to expand the template, DG method has attracted much attention because of its good mathematical characteristics.
The DG method can only solve hyperbolic conservation equations with first-order partial derivatives such as Euler equations at the earliest, Subsequently, many scholars proposed different discrete schemes for high-order derivative terms. Bassi and Rebay reduced the original second-order equation to two first-order systems by introducing the gradient of the variable as an intermediate variable, and then discretized them according to the DG framework, thus proposed the BR1 scheme 1 . BR1 scheme increases the number of equations, and it turns out to be incompatible 2 . Then they further improved the BR1 scheme and proposed the BR2 scheme to improve the convergence and stability of the algorithm 3 . The BR2 scheme is widely used in the simulation of equations containing second-order partial derivatives. Bassi 4,5 used it to solve compressible equations, and Eyck 6 used it to solve nonlinear elasticity equations. Later, Parsani 7 apply the BR2 scheme to LES simulation; Another method of discretizing elliptic equations is called interior penalty (IP). The internal penalty scheme originated from the penalty boundary method of finite element, which was proposed by Lions when solving elliptic equations with Dirichlet boundary conditions. He converted the Dirichlet boundary to the Neumann boundary by introducing a penalty coefficient. The internal penalty method extends the idea of penalty boundary to the boundary of the internal unit, and establishes the information exchange  8 . Subsequently, the IP method was widely used in the DG format. Hartmann 9 is applied to the DG solution of the compressible Navier-Stokes equation. Rivière 10 used the IP-DG scheme to simulate the incompressible Navier-Stokes equation on a grid with dangling nodes.
However, when the IP scheme applies to multivariate equations, such as the Navier-Stokes equation, it is necessary to calculate the Jacobian matrix of the gradient of the viscosity term to the variable. The matrix is a fourth-order tensor with a very complicated form, which is not conducive to the generalization of other equations. Based on the work of Harmann 11 and Wang 12 , an improved highorder IP scheme is proposed to solve the compressible Navier-Stokes equation, and through the twodimensional and three-dimensional typical calculation Examples to verify the effectiveness of the method.

The improved IP scheme
The compressible Navier-Stokes equation can be written as: F respectively represent the convection term and the viscosity term. The IP scheme decomposes the NS equation into two equations by defining the viscosity term as an intermediate variable In the above formula, : represents the contraction of the tensor. g is the derivative of the gradient of the three direction components of the five components of the viscosity term with respect to the 5 variables, which is a fourth-order tensor. For the NS equation and the RANS-SA equation, the tensor can be calculated analytically. However, the form of the tensor is complicated, which is not conducive to expanding to other equations, and programming is inconvenient. For this reason, an improved scheme is proposed, which eliminates the need to solve the fourth-order tensor by defining a new lifting operator and greatly improve the efficiency of viscosity item calculation.
According to the solution steps of the DG framework, multiply both sides of the first equation of formula (2) by the test function  .
(3) T g is the transposition of the tensor, represents the variable flux,  represents the dyadic of the tensor. For the first-order tensors u and v .
Since the variable is discontinuous at the cell boundary, the variable flux is defined as  (2), multiply both sides by  , the partial integral can be obtained ˆc F and ˆv F represent inviscid flux and viscous flux.
In order to eliminate the intermediate variable G , let     Refer to the first equation in formula (2) Define the lifting operator on the element surface as Incorporating the above formula into formula (8), the semi-discrete form of the improved IP scheme can be obtained as The Lax-Friedrichs scheme is used for the calculation of the flux, and the viscous flux is defined as follows: Penalty factor Refer to the BR format. Generally, the face number of the unit is taken as the algebraic order of the discrete polynomial. The definition of is as follows.  (11) is solved by the coupling of Newton and GMRES method.

Couette flow
Couette flow describes the flow between two plates, one of which moves in its own plane at a constant speed relative to the other, so that the fluid between the two plates moves under the action of friction. When time is sufficient, the fluid movement reaches a steady state. Couette flow has accurate solutions and is widely used to verify the accuracy of programs and methods to simulate laminar viscous flow. The exact solution of the Couette flow is as follows: . The number of quadrilateral grid cells are 10*5, 20*10, 40*20, and the calculation is performed with a uniform flow field as the initial condition, boundary conditions apply to accurate solutions. Table 1 and Table 2 show the density error and numerical accuracy of the improved IP scheme and BR2 scheme under different accuracy conditions. It can be seen that each scheme has reached the design accuracy.  Figure 1 shows the simulation grid and boundary of the cavity problem. The calculation conditions are Mach=0.1, Re=1000, the top velocity is equal to 1m/s, and the rest adopt isothermal solid wall boundary conditions. Figure 2 shows the comparison of velocity distribution in the x direction at the position x=0.5 and the velocity distribution in the y direction at the position y=0.5 from improved IP scheme, BR2  13 . The second-order DG method shows obvious excessive dissipation. With the improvement of order, the results of third-order and fourth-order DG method have good agreement with Ghia with the dense grid, which reflects the high-precision Advantage. It also shows the effectiveness of the improved internal penalty method. Figure 3 shows the streamline obtained by the fourth-order DG format, compared with the result of Ghia.

Supersonic flow over sphere
The calculation of the three-dimensional NS equation is verified by the calculation example of the supersonic flow over sphere. The incoming Mach number is Ma=2.0, Re=300, the grid type selects highorder TRI6, the total number of grids is 431922, the third-order (p=2) precision is selected, and the Riemann flux selects the HLLC flux. The simulation pressure distribution and the spherical surface pressure contours are shown in Figure 4. It can be seen that the method in this paper can effectively capture the shock wave and the flow separation phenomenon at the same time. Figure 5 shows the comparison between the friction coefficient of the spherical symmetry plane calculated in this paper and the reference [14], and the two are basically the same.

Conclusion
Under high-order DG framework, an improved internal penalty scheme for the discretization of the viscous term of the compressible Navier-Stokes equation is proposed in this paper. The simulation results show that the new scheme can achieve the design accuracy and is also effective for threedimensional problems. In the next step, the scheme will be extended to the simulation of RANS equations in the next step.