Two-Stage Fourth-order Gas-kinetic Scheme for Three-dimensional Euler and Navier-Stokes Solutions
Liang Pan Kun Xu
math.NA physics.comp-ph




iapcm]Liang Pan
panliangjlu@sina.com

HKUST1,HKUST2]Kun Xucor1
makxu@ust.hk

[iapcm]Institute of Applied Physics and Computational Mathematics, Beijing, 100088, China
[HKUST1]Department of Mathematics, Department of Mechanical and Aerospace Engineering,
Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong
[HKUST2]HKUST Shenzhen Research institute, Shenzhen 518057, China
[cor1]Corresponding author




For the one-stage third-order gas-kinetic scheme (GKS), success applications have been achieved for the three-dimensional compressible
flow computations <cit.>.
The high-order accuracy of the scheme is obtained
directly by integrating a multidimensional time-accurate gas distribution function over the cell interface within a time step
without implementing Gaussian quadrature points and Runge-Kutta time-stepping technique.
However, for the further increasing the order of the scheme, such as the fourth-order one,
the formulation becomes very complicated for the multidimensional flow.
Recently, a two-stage fourth-order GKS with high efficiency has been constructed
for two-dimensional inviscid and viscous flow computations <cit.>,
and the scheme uses the time accurate flux function and its time derivatives.
In this paper, a fourth-order GKS is developed for the three-dimensional flows under the two-stage framework.
Based on the three-dimensional WENO reconstruction and flux evaluation at Gaussian quadrature points on a cell interface,
the high-order accuracy in space is achieved first. Then, the two-stage time stepping method provides the high accuracy in time.
In comparison with the formal third-order GKS <cit.>,
the current fourth-order method not only improves the accuracy of the scheme, but also reduces the complexity of the gas-kinetic solver greatly. More importantly, the fourth-order GKS has the same robustness as the second-order 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 three-dimensional flows, such as
turbulent one.
With the count of the degree of freedom in simulating a flow field, the current two-stage fourth-order multidimensional GKS becomes
one of the most reliable and efficient higher-order schemes among all shock capturing higher-order schemes for computational fluid dynamics.



two-stage fourth-order discretization, gas-kinetic scheme, Navier-Stokes solutions, shock capturing scheme

Two-Stage Fourth-order Gas-kinetic Scheme for Three-dimensional Euler and Navier-Stokes Solutions
    [
    Received: date / Accepted: date
=================================================================================================







§ 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 well-developed second-order schemes, the higher-order methods can provide more accurate solutions,
but they are less robust and more complicated. For the high-order schemes, there are generally three parts,
i.e. spatial discretization, temporal discretization, and flux solvers. For the spatial discretization, many high-order methods have been developed, including the discontinuous Galerkin (DG) <cit.>, essential non-oscillatory (ENO) <cit.>, weighted essential non-oscillatory (WENO) <cit.>, and many others.
For most of those methods, the exact and approximate Riemann solvers <cit.> are used for flux evaluation.
Due to the first-order evolution of Riemann solvers, the Runge-Kutta method is implemented to achieve higher order temporal accuracy <cit.>, in which n-stage is needed for n-th order accuracy. Instead of Riemann solver, many schemes have been developed based on the time-dependent flux function with high-order evolution, such as the generalized Riemann problem solvers <cit.>, gas-kinetic scheme <cit.>, and AEDR methods <cit.>. High-order temporal accuracy can be achieved in a one-stage framework
with the time integration of the time-dependent flux function.
Recently, in order to increase the accuracy and efficiency of these schemes
a two-stage fourth-order method has been developed for the time-dependent flux solvers <cit.>, where both the flux and temporal derivative of flux function are used in the construction of higher-order schemes.
By combining the multi-stage multi-derivative technique, a family of higher-order schemes has been constructed as well <cit.>.




In the past decades, the gas-kinetic scheme (GKS) based on the Bhatnagar-Gross-Krook (BGK) model <cit.> has been developed systematically for the compressible flow computations <cit.>.
Different from the traditional CFD methods based on the macroscopic governing equations, the main advantages of the gas-kinetic scheme are the followings. (i) The inviscid and viscous coupling in the flux evolution <cit.>; (ii) Multi-dimensionality with the inclusion of both normal and tangential gradients of flow variables in the flux function across a cell interface <cit.>; (iii) Compact stencils can be constructed with the use of the time accurate cell interface flow variables at the next time level <cit.>; (iv) Extension to the whole flow regime from the rarefied to the continuum one <cit.>. Recently, with the high-order initial reconstruction, the third-order gas-kinetic schemes have been constructed <cit.>, in which the flux evaluation is based on the moments of spatial-temporal dependent gas distribution function.
High-order accuracy can be achieved in a one-stage scheme without Gaussian point integration for spatial accuracy
and Runge-Kutta technique for temporal accuracy.
However, with the one-stage gas evolution model, the formulation of gas-kinetic scheme can become very complicated for the further development,
such as the fourth-order method <cit.>, especially for multidimensional computations.
The two-stage fourth-order temporal discretization for time-accurate flux solvers in <cit.> provides a reliable framework to further develop the GKS into fourth-order and even higher accuracy with the implementation of the traditional second-order or third-order flux functions <cit.>. 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 Δ t over the particle collision time τ becoming large. The real physics used in the flux depends on the local Δ t /τ <cit.>.
In real NS computations for the compressible high Mach number flow, the ratio of Δ t /τ is not too large as people think of for the
validation of the NS modeling <cit.>.
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 solver-based Godunov-type schemes for the NS equations.
For the second-order schemes, it is hard to distinguish the dynamical differences from the GKS and Riemann solver.
However,  for the higher-order schemes it seems that a reliable physical evolution model becomes more
important due to the absence of large numerical dissipation in the second-order schemes, and the delicate flow structures captured in higher-order
schemes depend on the quality of the solvers greatly <cit.>.



In this paper, with the two-stage fourth-order discretization, a multidimensional fourth-order gas-kinetic scheme is constructed for simulating  three-dimensional flows. High-order accuracy in space is achieved by the three-dimensional WENO method <cit.> and Gaussian quadrature points for the numerical fluxes. In comparison with the formal three-dimensional scheme <cit.>, the current fourth-order scheme reduces the complexity of the gas-kinetic 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 fourth-order
method.  Numerical results show that the current scheme has the same reliability and applicability as the well-developed second-order 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 state-of-art solutions from a higher-order 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 two-stage 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.



§ BGK EQUATION AND FINITE VOLUME SCHEME

The three-dimensional BGK equation <cit.> can be written as

    f_t+uf_x+vf_y+wf_z=g-f/τ,

where f is the gas distribution function, g is the corresponding equilibrium state, and τ is the collision time. The collision term satisfies the compatibility condition

    ∫g-f/τψdΞ=0,

where ψ=(1,u,v,w,1/2(u^2+v^2+w^2+ξ^2)), the internal variables ξ^2 is equal to ξ^2=ξ_1^2+...+ξ_K^2, dΞ=dudvddwdξ^1...dξ^K, K is the degrees of freedom, i.e. K=(5-3γ)/(γ-1) for three-dimensional flows and γ is the specific heat ratio. In the continuum region, the gas distribution function can be expanded as

    f=g-τ D_ug+τ D_u(τ D_u)g-τ D_u[τ D_u(τ D_u)g]+...,

where D_u=∂/∂ t+u·∇. Based on the Chapman-Enskog expansion, the macroscopic equations can be derived <cit.>. For the Euler equations, the zeroth-order truncation is taken, i.e. f=g. For the Navier-Stokes equations, the first-order truncation is

    f=g-τ (ug_x+vg_y+wg_z+g_t).

With the higher order truncations, the Burnett and super-Burnett equations can be derived.

Taking moments of Eq.(<ref>) and integrating over the control volume V_ijk=x_i×y_j×z_k with x_i=[x_i-Δ x/2,x_i+Δ x/2], y_j=[y_j-Δ y/2,y_j+Δ y/2], z_k=[z_k-Δ z/2,z_k+Δ z/2],  the
semi-discretized form of finite volume scheme can be written as

    dQ_ijk/dt=ℒ(Q_ijk)=1/|V_ijk|[
    ∫_y_j×z_k(F_i-1/2,j,k-F_i+1/2,j,k)dydz
       +∫_x_i×z_k(G_i,j-1/2,k-G_i,j+1/2,k)dxdz+∫_x_i×y_j(H_i,j,k-1/2-H_i,j,k+1/2)dx
    dy],

where Q=(ρ,ρ U,ρ V,ρ W,ρ E) are the conservative flow variables, Q_ijk is the cell averaged value over control volume V_ijk and |V_ijk|=Δ xΔ yΔ z.
For the three-dimensional computation, the Gaussian quadrature for the numerical fluxes is used to achieve the accuracy in space, and the numerical fluxes in x-direction is given as an example

    ∫_y_j×z_kF_i+1/2,j,kdydz=Δ yΔ z∑_m,n=1^Mω_mn F(x_i+1/2,j_m,k_n,t),

where x_i+1/2,,m,n=(x_i+1/2,y_j_m,z_k_n), (y_j_m,z_k_n), m, n= 1,...,M are the Gauss quadrature points for y_j×z_k and ω_mn are corresponding weights, and M=2 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 F(x_i+1/2,j_m,k_n,t) can be provided by the flow solvers. Similarly, the numerical fluxes  in the y and z-directions can be obtained as well

    ∫_x_i×z_kG_i,j+1/2,kdxdz=Δ xΔ z∑_m,n=1^Mω_mn G(x_i_m,j+1/2,k_n,t),
    ∫_x_i×y_jH_i,j,k+1/2dxdy=Δ xΔ y∑_m,n=1^Mω_mn H(x_i_m,j_n,k+1/2,t).

In the gas-kinetic scheme, the numerical fluxes at the Gauss quadrature point can be obtained by taking moments of the gas distribution function

    F(x_i+1/2,j_m,k_n,t)=∫ψ u f(x_i+1/2,j_m,k_n,t,u,ξ)dudvdwdξ,

where f(x_i+1/2,j_m,k_n,t,u,ξ) is provided by the integral solution of BGK equation Eq.(<ref>) at the cell interface

    f(x_i+1/2,j_m,k_n,t,u,ξ)=   1/τ∫_0^tg(x',t',u,ξ)e^-(t-t')/τdt'+e^-t/τf_0(-ut,ξ),

and x_i+1/2,j_m,k_n=0 is the location of cell interface, u=(u,v,w) is the particle velocity, x_i+1/2=x'+u(t-t'), y_j_m=y'+v(t-t'), z_k_n=z'+w(t-t') is the trajectory of particles.   For the second-order scheme, the second-order gas-kinetic solver for the three-dimensional flows <cit.> can be expressed as

    f(x_i+1/2,m,n,t,u,ξ)=   (1-e^-t/τ)g_0+((t+τ)e^-t/τ-τ)(a_1u+a_2v+a_3w)g_0
    
    +   (t-τ+τ e^-t/τ)A̅ g_0
    
    +   e^-t/τg_r[1-(τ+t)(a_1ru+a_2rv+a_3rw)-τ A_r)]H(u)
    
    +   e^-t/τg_l[1-(τ+t)(a_1lu+a_2lv+a_3lw)-τ A_l)](1-H(u)).

where the coefficients in Eq.(<ref>) can be determined by the spatial derivatives of macroscopic flow variables and the compatibility condition as follows

    ⟨ a_1⟩ =∂ W /∂ x,
    ⟨ a_2⟩ =∂ W /∂ y, ⟨
    a_3⟩ =∂ W /∂ z,
    ⟨ A+a_1u+a_2v+a_3w⟩=0,

where the superscripts or subscripts of these coefficients a_1, a_2, a_3, A are omitted for simplicity, more details about the
determination of coefficient can be found in <cit.>. With the second-order gas-kinetic solver Eq.(<ref>), the second-order accuracy in time can be achieved by one step integration.
In the one-stage gas evolution model, the third-order gas-kinetic solver has been developed as well <cit.>.



§ FOURTH-ORDER TEMPORAL DISCRETIZATION

In the further development of higher-order scheme, the formulation of one-stage gas-kinetic solver can become very complicated, such as the fourth-order method <cit.>, especially for multidimensional computations.
Recently, a two-stage fourth-order time-accurate discretization was developed for Lax-Wendroff flow solvers, particularly applied for hyperbolic equations with the generalized Riemann problem (GRP) solver <cit.> and gas-kinetic scheme <cit.>. Such method provides a reliable framework to develop a
three-dimensional fourth-order gas-kinetic scheme with a second-order flux function  Eq.(<ref>). Consider the following time-dependent equation

    ∂w/∂ t=ℒ(w),

with the initial condition at t_n, i.e.,

    w(t=t_n)=w^n,

where ℒ is an operator for spatial derivative of flux.
The time derivatives are obtained using the Cauchy-Kovalevskaya method,

    ∂w^n/∂ t=ℒ(w^n),   ∂/∂ tℒ(w^n)=∂/∂wℒ(w^n)ℒ(w^n).

Introducing an intermediate state at t^*=t_n+Δ t/2,

    w^*=w^n+1/2Δ tℒ(w^n)+1/8Δ t^2∂/∂
    tℒ(w^n),

the corresponding time derivatives are obtained as well for the intermediate stage state,

    ∂w^*/∂ t=ℒ(w^*), ∂/∂ tℒ(w^*)=∂/∂wℒ(w^*) ℒ(w^*).

Then, the state w can be updated with the following formula,

    w^n+1=w^n+Δ tℒ(w^n)+1/6Δ t^2(∂/∂
    tℒ(w^n)+2∂/∂
    tℒ(w^*)).

It can be proved that for hyperbolic equations  the above time
stepping method  Eq.(<ref>) and Eq.(<ref>)  provides a fourth-order time
accurate solution for w(t) at t=t_n +Δ t. The
details of the analysis can be found in <cit.>.




The semi-discretized form of finite volume scheme
Eq.(<ref>) can be applied for the above two-stage method
. For the gas-kinetic 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 t_n and t_*=t_n + Δ t/2 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.(<ref>),  the following notation is introduced

    𝔽_i+1/2,j,k(Q^n,δ)=∫_t_n^t_n+δF_i+1/2,j,k,(Q^n,t)dt
    =∑_m,n=1^Mω_mn∫_t_n^t_n+δF(x_i+1/2,j_m,k_n,t)dt.

In the time interval [t_n, t_n+Δ t], the flux is expanded as
the following linear form

    F_i+1/2,j,k(Q^n,t)=F_i+1/2,j,k^n+ ∂_t F_i+1/2,j,k^n(t-t_n).

The coefficients F_i+1/2,j,k(Q^n,t_n) and ∂_tF_i+1/2,j,k(Q^n,t_n) can be
determined as follows,

    F_i+1/2,j,k(Q^n,t_n)Δ t   +1/2∂_t
    F_i+1/2,j,k(Q^n,t_n)Δ t^2 =𝔽_i+1/2,j,k(Q^n,Δ t) , 
    1/2F_i+1/2,j,k(Q^n,t_n)Δ t   +1/8∂_t
    F_i+1/2,j,k(Q^n,t_n)Δ t^2 =𝔽_i+1/2,j,k(Q^n,Δ t/2).

By solving the linear system, we have

    F_i+1/2,j,k(Q^n,t_n)   =(4𝔽_i+1/2,j,k(Q^n,Δ t/2)-𝔽_i+1/2,j,k(Q^n,Δ t))/Δ t,
    ∂_t F_i+1/2,j,k(Q^n,t_n)   =4(𝔽_i+1/2,j,k(Q^n,Δ t)-2𝔽_i+1/2,j,k(Q^n,Δ t/2))/Δ
    t^2.

Similarly, F_i+1/2,j,k(Q^*,t_*), ∂_t
F_i+1/2,j,k(Q^*,t_*) for the intermediate state t^*=t_n+Δ t/2 can be constructed. The corresponding fluxes in the y- and z-direction can be
obtained as well. With these notations, the numerical scheme for the three-dimensional flows can be constructed, and the detail procedure of two-stage algorithm can be found in <cit.>.







§ SPATIAL RECONSTRUCTION

The above time evolution solution is based on the high-order initial reconstruction for macroscopic flow variables and  WENO reconstruction <cit.> is adopted for the spatial reconstruction.
For the three dimensional computation, the reconstruction procedure for the cell interface x_i+1/2,j,k is given as an example, and the stencil for reconstruction is given in Fig.<ref>.
The point value Q_l, Q_r and Q_0 and first-order derivatives at the Gauss quadrature points (x_i+1/2,y_j_m,z_k_n), m, n= 1,...,2 need to be constructed. The detailed procedure is given as follows


  * According to one dimensional reconstruction, the cell averaged reconstructed value (Q_l)_j-ℓ_1,k-ℓ_2, (Q_r)_j-ℓ_1,k-ℓ_2,(Q_0)_j-ℓ_1,k-ℓ_2, ℓ_1,ℓ_2=-2,...,2 and cell averaged spatial derivatives (∂_xQ_l)_j-ℓ_1,k-ℓ_2, (∂_xQ_r)_j-ℓ_1,k-ℓ_2,(∂_xQ_0)_j-ℓ_1,k-ℓ_2 can be constructed for the stencil shown in Fig.<ref>.

  * With the one-dimensional WENO reconstruction in the horizontal direction, the averaged value (Q_l)_j_m,k-ℓ_2, (Q_r)_j_m,k-ℓ_2,(Q_0)_j_m,k-ℓ_2 and the averaged spatial derivatives  (∂_xQ_l)_j_m,k-ℓ_2, (∂_xQ_r)_j_m,k-ℓ_2,(∂_xQ_0)_j_m,k-ℓ_2 and (∂_yQ_l)_j_m,k-ℓ_2, (∂_yQ_r)_j_m,k-ℓ_2,(∂_yQ_0)_j_m,k-ℓ_2 over the interval [z_k-ℓ_2-Δ z/2,z_k-ℓ_2+Δ z/2] with y=y_j_m can be given.

  * With one-dimensional WENO reconstruction in the vertical direction,  the point value (Q_l)_j_m,k_n, (Q_r)_j_m,k_n,(Q_0)_j_m,k_n and spatial derivatives (∂_xQ_l)_j_m,k_n, (∂_xQ_r)_j_m,k_n,(∂_xQ_0)_j_m,k_n,  (∂_yQ_l)_j_m,k_n, (∂_yQ_r)_j_m,k_n,(∂_yQ_0)_j_m,k_n
and (∂_zQ_l)_j_m,k_n, (∂_zQ_r)_j_m,k_n,(∂_zQ_0)_j_m,k_n can be fully given at the Gaussian quadrature points (x_i+1/2,y_j_m,z_k_n).

In the computation, without special statement, the fifth-order WENO-JS reconstruction <cit.> is adopted for the flow with discontinuities and the linear scheme is used for the smooth flows to reduce the dissipation.







§ 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

    τ=ϵΔ t+C|p_l-p_r/p_l+p_r|Δ t,

where p_l and p_r denote the pressures on the left and right sides of the cell interface. In the computation, ε=0.01 and C=1. For the viscous flow, we have

    τ=μ/p+|p_l-p_r/p_l+p_r|Δ t,

where μ is the viscous coefficient and p is the pressure at the cell interface, and it will reduce to τ=μ/p in the smooth flow regions. Δ t is the time step which is determined according to the CFL number, which takes 0.4 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 non-equilibrium 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 A_i+1/2,j=(∂ F/∂ Q)_Q=Q^*, where Q are the conservative variables, F(Q) are the corresponding fluxes, and Q^*=(Q_i,j,k+Q_i+1,j,k)/2, the cell averaged and point conservative values can be projected into the characteristic field by ω=RQ, where R is the matrix corresponding to right eigenvectors of A. 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 WENO-JS reconstruction, the results presented below
have very high accuracy in comparison with  other solutions from same reconstruction.




 §.§ Three-dimensional Sod problem

This problem is a fully three-dimensional extension of the Sod problem. The computational domain is (x,y,z)∈[0, 1]×[0, 1]×[0, 1], and a sphere of radius R =0.5 separates two different constant states

    (ρ, U, V, p)=
    
    (1, 0, 0, 1),           0<R<0.5,
    
    (0.125, 0, 0, 0.1),        0.5<R<1.

The exact solution of this problem can be given by the one-dimensional system with geometric source terms

    ∂ Q/∂ t+∂ F(Q)/∂
    r=S(Q),

with

    Q=[   ρ; ρ U; ρ E;     ],
    F(Q)=[      ρ U;  ρ U^2+p; U(ρ E+p);          ],
    S(Q)=-d-1/r[      ρ U;    ρ U^2; U(ρ E+p);          ],

where the radial direction is denoted by r, while U is the radial velocity, d is the number of space dimensions and γ=1.4. The two-stage fourth-order gas-kinetic scheme is used to solve this equation and the numerical results with 10000 cells are given as the reference solutions. In the computation, the uniform mesh with Δ x=Δ y=Δ z=1/100 is used. The symmetric boundary condition is imposed on the plane with x=0, y=0 and z=0, and the non-reflection boundary condition is imposed on the plane with x=1, y=1 and z=1. The three-dimensional density and pressure distributions are given in Fig.<ref>, and the density and pressure distribution along the line y=z=0 are given in Fig.<ref>. The numerical results agree well with the reference solutions.










 §.§ 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 γ=1.4. In the computation, a 40×40×40 mesh shown in Fig.<ref> is used, which representing the domain [-1.75,-1]×[-0.4π,0.4π]×[0.1π,0.9π] in the spherical coordinate (r,ϕ,θ).  The density distributions in the whole domain, Mach number and pressure distributions at the plane with ϕ=0, and pressure profiles along the surface of sphere with ϕ=0 for Mach number Ma=5 and 8 are shown in Fig.<ref>, Fig.<ref> and Fig.<ref>, where the shock is well captured by the current scheme and the carbuncle phenomenon does not appear <cit.>. As comparison, the inviscid hypersonic flows impinging on a two-dimensional unit cylinder is also tested with 40×80 mesh. The Mach number and  pressure distributions for Mach number Ma=5, 8 are shown in Fig.<ref>.












 §.§ Rayleigh-Taylor instability

The Rayleigh-Taylor instability occurs when an interface between two different fluids with different densities experiences a pressure gradient opposing the density gradient <cit.>. 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 Rayleigh-Taylor instability usually evolves into complicated turbulent mixing. In this case, the three-dimensional Rayleigh-Taylor instability is tested in a rectangular box with domain [0,0.25]×[0,0.25]×[0,1] with a gravity pointed downward g=(0, 0,-0.1).
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

    h(x,y)/W=0.05[cos(2π x/W)+cos(2π y/W)],

where h is the height of the interface and W is the box width. The dimensionless parameters is the Atwood number

    A=ρ_h-ρ_l/ρ_h+ρ_l,

where ρ_h, ρ_l are densities of heavy and light fluids, respectively. In the computation, the Atwood number A=1/3 is used. The the fluid interface from a single-mode perturbation at t=3, 4, 5, 6, 7, 8 are given in Fig.<ref>, and the density distribution at x=0, x=W/2 and x=y with t=3,4,5,6,7,8 are given in Fig.<ref>. 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.















 §.§ Lid-driven cavity flow

The lid-driven cavity problem is one of the most important benchmarks for numerical Navier-Stokes solvers. As shown in Fig.<ref>, the fluid is bounded by a unit cubic and driven by a uniform translation of
the top boundary. The monatomic gas with γ = 5/3 is used, such that there is no bulk viscosity involved. Early three-dimensional cavity-flow calculations were carried out by De Vahl Davis and Mallinson <cit.> and Goda <cit.>. In this case, the flow is simulated with Mach number Ma=0.15 and all the boundaries are isothermal and nonslip. Numerical simulations are conducted with three Reynolds numbers of Re=1000, 400 and 100 using 65×65×65 meshes for the domain [-0.5, 0.5]×[-0.5, 0.5]×[-0.5, 0.5]. The u-velocity profiles along the vertical centerline line, v-velocity profiles along the horizontal centerline in the x-y plane with z=0 and the benchmark data <cit.> are shown in Fig.<ref>. The simulation results match well with the benchmark data.











 §.§ Compressible homogeneous turbulence

The high-order gas-kinetic 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 -π≤ x, y, z≤π, and the periodic boundary conditions are used in all directions for all the flow variables <cit.>. In the computation, the
domain is discretized with a uniform Cartesian mesh cells N^3. A divergence-free random initial velocity field u_0 is generated
for a given spectrum with a specified root mean square u'

    u'=<u·u/3>,

where <...> is a volume average over the whole computational domain. The specified spectrum for velocity is given by

    E(k)=A_0k^4exp(-2k^2/k_0^2),

where k is the wave number, k_0 is the wave number at spectrum peaks, A is a constant chosen to get a specified initial kinetic energy.  The initial volume averaged turbulent kinetic energy K_0  and the initial large-eddy-turnover time τ is given by

    K_0=3A_0/64√(2π)k_0^5,  τ=32/A_0(2π)^1/4k_0^-7/2.

The Taylor microscale Reynolds number Re_λ and turbulence Mach number Ma_t are given as

    Re_λ   =<ρ>u'λ/<μ>=(2π)^1/4/4ρ_0/μ_0√(2A_0)k_0^3/2,
    
    Ma_t   =√(3)u'/<c_s>=√(3)u'/√(γ T_0),

where λ is Taylor microscale

    λ^2=(u')^2/<(∂_1 u_1)^2>.

μ_0 and T_0 can be determined from Re_λ and Ma_t with initialized u' and ρ_0=1, and the dynamic viscosity is determined by

    μ=μ_0(T/T_0)^0.76.

The time history of the kinetic energy, root-mean-square of density fluctuation and skewness factor for velocity slope are defined as

    K(t)   =1/2<ρu·u>,
    ρ_rms(t)   =√(<(ρ-ρ)^2>),
    
    S_u(t)=   ∑_i<(∂_i u_i)^3>/<(∂_i u_i)^2>^3/2.

In the computation, A_0=1.3×10^-4, k_0=8, Re_λ=72 and Ma_t=0.5, and the uniform meshes with 64^3, 96^3 and 128^3 cells are used. The iso-surfaces of Q criterion colored by velocity magnitude and the pressure distribution with z=-π at time t/τ=1  are given in Fig.<ref>. The time history of normalized kinetic energy K(t)/K_0, normalized root-mean-square of density fluctuation ρ_rms(t)/Ma_t^2  and skewness factor S_u(t) with respect to t/τ are given in Fig.<ref>. The numerical results agree well with the reference data. With fixed initial Re_λ=72, the cases with Ma_t=0.2, 0.5, 0.6 are tested, and the time histories of normalized kinetic energy K(t)/K_0 are given in Fig.<ref> as well. With the increase of Ma_t, the dynamic viscosity increases, and the kinetic energy gets dissipated more rapidly.








 §.§ Taylor-Green Vortex

This problem is aimed at testing the performance of high-order
methods on the direct numerical simulation of a three-dimensional
periodic and transitional flow defined by a simple initial
condition, i.e. the Taylor-Green vortex <cit.>. With a uniform temperature field, the
initial flow field is given by

    u=   V_0sin(x/L)cos(y/L)cos(z/L),
    
    v=   -V_0cos(x/L)sin(y/L)cos(z/L),
    
    w=   0,
    
    p=   p_0+ρ_0V_0^2/16(cos(2x/L)+cos(2y/L))(cos(2z/L)+2).

The fluid is then a perfect gas with γ=1.4 and the Prandtl
number is Pr=0.71. Numerical simulations are conducted with two
Reynolds numbers Re=1600 and 280. The flow is computed within a
periodic square box defined as -π L≤ x, y, z≤π L. The
characteristic convective time t_c = L/V_0. In the computation,
L=1, V_0=1, ρ_0=1, and the Mach number takes M_0=V_0/c_0=0.1,
where c_0 is the sound speed.


The volume-averaged kinetic energy can be computed from the flow as it evolves in time, which is expressed as

    E_k=1/ρ_0Ω∫_Ω1/2ρu·udΩ,

where Ω is the volume of the computational domain, and the dissipation rate of the kinetic energy is given by

    ε_k=-dE_k/dt.

The numerical results of the current scheme with 128×128×128  mesh points for the normalized volume-averaged kinetic energy and dissipation rate with Reynolds numbers Re=1600 and 280 are presented in Fig.<ref> and Fig.<ref>, which agree well with the data in <cit.>. The iso-surfaces of Q criterions colored by velocity magnitude at t=5 and 10 are shown in Fig.<ref> for Re=1600 and in Fig.<ref> for Re=280. The evolution of flow structure is evident, starting from large vortices and decaying into more complex structures. Different from many other higher-order methods, the current scheme has no internal degrees of freedom to be updated within each cell.





§ CONCLUSION

In this paper, based on the two-stage time stepping method, a fourth-order gas-kinetic scheme is proposed for the three-dimensional inviscid and viscous flow computations.
With the three-dimensional WENO-JS reconstruction, a gas-kinetic scheme with higher-order spatial and temporal accuracy is developed.
In comparison with the classical methods based on the first-order Riemann solver,
for the same fourth-order accuracy in time the current scheme only uses two stages instead of four stages with the Runge-Kutta time-stepping technique. As a result, the two-stage GKS can be more efficient than the four-stage methods with the absence of two time consuming reconstructions.
For the Navier-Stokes solutions, the current scheme doesn't separate inviscid and viscous terms and they are treated uniformly from the same
initial WENO-type 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 fourth-order 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 fourth-order scheme has the same robustness as the second-order 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 Navier-Stokes solutions, the current scheme is one of the state-of-art 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).



Case-Albensoeder S. Albensoeder, H.C. Kuhlmann, Accurate three-dimensional lid-driven cavity flow, J. Comput. Phys. 206 (2005) 536-558.

GRP1 M. Ben-Artzi, J. Falcovitz, A second-order Godunov-type scheme for compressible uid dynamics, J. Comput. Phys. 55 (1984), 1-32.

GRP2 M. Ben-Artzi, J. Li, G. Warnecke,  A direct Eulerian GRP scheme for compressible fluid flows, J. Comput. Phys. 218 (2006), 19-43.

GRP3 M. Ben-Artzi, J. Li, Hyperbolic conservation laws: Riemann invariants and the generalized Riemann problem, Numerische Mathematik, 106 (2007), 369-425.

BGK-1 P.L. Bhatnagar, E.P. Gross, M. Krook, A Model for Collision Processes in Gases I: Small Amplitude Processes in Charged and Neutral One-Component Systems, Phys. Rev. 94 (1954) 511-525.

WENO-Z R. Borges, M. Carmona, B. Costa, W. S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, J. Comput. Phys. 227 (2008) 3191-3211.

Case-Bull J. R. Bull, A. Jameson, Simulation of the compressible Taylor-Green vortex using high-order flux reconstruction schemes, AIAA 2014-3210.

BGK-2 C. Cercignani, The Boltzmann Equation and its Applications, Springer-Verlag, (1988).

BGK-3 S. Chapman, T.G. Cowling, The Mathematical theory of Non-Uniform Gases, third edition, Cambridge University Press, (1990).

DG-2 B. Cockburn, C. W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: general framework, Mathematics of Computation, 52 (1989), 411-435.

DG-3 B. Cockburn, C. W. Shu, The Runge-Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems, J. Comput. Phys. 141 (1998), 199-224.

Case-Debonis J.Debonis, Solutions of the Taylor-Green vortex problem using high-resolution explicit finite difference methods, AIAA Paper (2013) 2013-0382.

Case-Davis 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) 29-43.

GRP-high-2 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) 385-396


TVD-rk S. Gottlieb, C. W. Shu, Total variation diminishing runge-kutta schemes, Mathematics of computation, 67 (1998) 73-85.

Case-Goda K. Goda, A multistep technique with implicit difference schemes for calculating two- or three-dimensional cavity flows, J. Comput. Phys. 30 (1979) 76-95.

ENO-1 A. Harten, B. Engquist, S. Osher and S. R. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes, III. J. Comput. Phys. 71 (1987) 231-303.

Case-rt X.Y. He, R.Y Zhang, S.Y. Chen, G. D. Doolen, On the three-dimensional Rayleigh-Taylor instability, Phys Fluids 11 (1999) 1143.

UGKS-2 J. Huang, K. Xu, P. Yu, A unified gas-kinetic scheme for continuum and rarefied flows II: multi-dimensional cases, Commun. Comput. Phys. 12 (2012) 662-690.

MSMD-GKS X. Ji, F.X. Zhao, W. Shyy, K. Xu,  A Family of High-order Gas-kinetic Schemes and Its Comparison with Riemann Solver Based High-order Methods, J. Comput. Phys. 356 (2018), 150C173.

WENO-JS G.S. Jiang, C. W. Shu, Efficient implementation of Weighted ENO schemes, J. Comput. Phys. 126 (1996) 202-228.

GRP-high J. Li, Z. Du, A two-stage fourth order time-accurate discretization for Lax-Wendroff type flow solvers I. hyperbolic conservation laws, SIAM J. Sci. Computing, 38 (2016) 3046-3069.

DNS-3 Q.B. Li, S. Fu, High-order accurate gas-kinetic scheme and turbulence simulation, SCIENTIA SINICA Physica, Mechanica & Astronomica, 44 (2014) 278-284.

GKS-high-2 Q. Li, K. Xu, S. Fu,  A high-order gas-kinetic Navier-Stokes flow solver, J. Comput. Phys. 229 (2010) 6715-6731.

DNS-1 W. Liao, Y. Peng, L.S. Luo, Gas-kinetic schemes for direct numerical simulations of compressible homogeneous turbulence. Phys Rev E, 80 (2009) 046702.

liu-kn C. Liu, K. Xu, G.Z. Zhou, Cell Size Effect on Computational Fluid Dynamics: The Limitation Principle for Flow Simulation, arXiv:1706.08200v4 [physics.comp-ph].

GKS-high-3 J. Luo, K. Xu, A high-order multidimensional gas-kinetic scheme for hydrodynamic equations, SCIENCE CHINA Technological Sciences, 56 (2013) 2370-2384.

GKS-high-5 N. Liu, H.Z. Tang, A high-order accurate gas-kinetic scheme for one- and two-dimensional flow simulation,  Commun. Comput. Phys. 15 (2014) 911-943.

WENO X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys. 115 (1994) 200-212.

GKS-high-6 L. Pan, J. Li, K. Xu, A few benchmark test cases for higher-order Euler solvers, Numerical Mathematics: Theory, Methods and Applications, 10 (2017) 711-736.

GKS-high-4 L. Pan, K. Xu, A third-order compact gas-kinetic scheme on unstructured meshes for compressible Navier-Stokes solutions, Journal of Computational Physics, 318 (2016) 327-348.

GKS-high L. Pan, K. Xu, Q. Li, J. Li, An efficient and accurate two-stage fourth-order gas-kinetic scheme for the Navier-Stokes equations, J. Comput. Phys. 326 (2016), 197-221.

GKS-high-1 L. Pan, K. Xu, A third-order gas-kinetic scheme for three-dimensional inviscid and viscous flow computations, Computers & Fluids 119 (2015) 250-260.

Case-Pandolfi M. Pandolfi, and D. D'Ambrosio, Numerical Instabilities in Upwind Methods: Analysis and Cures for the "Carbuncle" Phenomenon, J. Comput. Phys.  166 (2001) 271-301.

DG-1 W.H. Reed, T.R. Hill, Triangular mesh methods for the neutron transport equation, Technical Report LA-UR-73-479, 1973, Los Alamos Scientific Laboratory, Los Alamos.

DNS-2 R. Samtaney, D.I. Pullin, B. Kosovic,  Direct numerical simulation of decaying compressible turbulence and shocklet statistics. Phys Fluids, 13 (2001) 1415-1430.

Case-Shu C. Shu, L. Wang, Y.T. Chew, Numerical computation of three-dimensional incompressible Navier.Stokes equations in primitive variable form by DQ method, Int. J. Numer. Meth. Fluids. 43 (2003) 345-368.

ENO-2 C. W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys. 77 (1988) 439-471.

ADER  V.A. Titarev, E.F. Toro, ADER: arbitrary high order Godunov approach, J. Sci. Comput. 17 (2002) 609C618.

Riemann-appro E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, (1997).

Case-wang L. Wang, W. K. Anderson, Erwin, S. Kapadia, High-Order Discontinuous Galerkin Method for Computation of Turbulent Flows, AIAA Journal, 53 (2015) 1157-1171.

GKS-Xu1 K. Xu, Direct modeling for computational fluid dynamics: construction and application of unfied gas kinetic schemes, World Scientific (2015).

GKS-Xu2 K. Xu, A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method, J. Comput. Phys. 171 (2001) 289-335.

GKS-Xu3 K. Xu, Gas kinetic schemes for unsteady compressible flow simulations, Lecure Note Ser. 1998-03, Von Karman Institute for Fluid Dynamics Lecture (1998).

UGKS-1 K. Xu, J. Huang, A unified gas-kinetic scheme for continuum and rarefied flows, J. Comput. Phys.  229 (2010) 7747-7764.

xu-liu K. Xu, C. Liu, A paradigm for modeling and computation of gas dynamics, Phys. Fluids 29 026101 (2017).

