A three-dimensional gas-kinetic flux solver for simulation of viscous flows with explicit formulations of conservative variables and numerical flux

A truly three-dimensional (3D) gas-kinetic flux solver for simulation of incompressible and compressible viscous flows is presented in this work. By local reconstruction of continuous Boltzmann equation, the inviscid and viscous fluxes across the cell interface are evaluated simultaneously in the solver. Different from conventional gas-kinetic scheme, in the present work, the distribution function at cell interface is computed in a straightforward way. As an extension of our previous work (Sun et al., Journal of Computational Physics, 300 (2015) 492–519), the non-equilibrium distribution function is calculated by the difference of equilibrium distribution functions between the cell interface and its surrounding points. As a result, the distribution function at cell interface can be simply calculated and the formulations for computing the conservative flow variables and fluxes can be given explicitly. To validate the proposed flux solver, several incompressible and compressible viscous flows are simulated. Numerical results show that the current scheme can provide accurate numerical results for three-dimensional incompressible and compressible viscous flows.


Introduction
In the last few decades, the gas-kinetic scheme has been developed in both continuum [1][2][3][4][5][6][7][8] and rarefied flow regimes [9][10][11][12]. Unlike the traditional Riemann solver [13][14][15], the gas-kinetic scheme reconstructs the solution for the continuous Boltzmann equation at local cell interface. As the continuum assumption is avoided in the continuous Boltzmann equation, the gas-kinetic scheme can be applied in both continuum and rarefied flow problems, which is one of the advantages as compared with traditional Riemann solver. Another advantage is that the gas-kinetic scheme can compute the inviscid and viscous fluxes simultaneously from the solution of Boltzmann equation. In contrast, the traditional Riemann solver can only evaluate the inviscid flux and an additional step is required to compute the viscous flux by smooth function approximation.
Currently, there are two common types of gas-kinetic schemes: the kinetic flux vector scheme (KFVS) and gas-kinetic Bhatnagar-Gross-Krook (BGK) scheme. In KFVS, the Boltzmann equation without collision term, which is also called the collisionless Boltzmann equation, is solved. Basically, there are two stages in KFVS: free transport and collision. In the free transport stage, the collisionless Boltzmann equation is solved to calculate the flux at the interface. In the collision stage, the artificial collisions are added in the calculation of initial Maxwellian distribution at the beginning of next time step. The KFVS has been demonstrated to have good positivity property for simulation of flows with strong shock waves [2]. However, owing to the fact that the numerical dissipation in KFVS is proportional to the mesh size, the KFVS usually gives more diffusive results than the Godunov or flux difference splitting (FDS) scheme [16], and is not able to give accurate Navier-Stokes solutions except for cases in which the physical viscosity is much larger than the numerical viscosity. Some of representative researches on KFVS include Pullin [17], Deshpande [18], Perthame [19], Mandal and Deshpande [20], and Chou and Baganoff [21].
One of the significant developments in gas-kinetic schemes is the gas-kinetic BGK scheme, which was firstly proposed by Prendergast and Xu [22], and further developed by Chae et al. [23], Xu [24] and other researchers. In this method, the BGK collision model is adopted in the solution process to obtain the numerical fluxes across the interface. As a consequence, the dissipation in the transport can be controlled by a real collision time, which is proportional to the dynamic viscosity. The gas-kinetic BGK scheme enjoys some intrinsic advantages. Firstly, it has been shown that the gas-kinetic BGK scheme is able to generate a stable and crisp shock transition in the discontinuous region with a delicate dissipative mechanism [24]. At the same time, an accurate Navier-Stokes solution can be obtained in the smooth region. What is more, it is demonstrated that the entropy condition is always satisfied in the gas-kinetic BGK scheme and the "carbuncle phenomenon" is avoided for hypersonic flow simulations [25]. However, the gas-kinetic BGK scheme is not completely free from drawbacks. It is usually more complicated and inefficient than conventional computational fluid dynamics (CFD) schemes. This is because in the gas-kinetic BGK scheme, a number of coefficients related to the physical space should be calculated to evaluate the distribution function at each interface and each time step. Moreover, to the best of our knowledge, there is still no work of the gas-kinetic BGK scheme which can give explicit formulations for evaluating the conservative variables and numerical fluxes.
Recently, a straightforward way to evaluate the distribution function was proposed by Sun et al. [1], which is named as gas-kinetic flux solver (GKFS). Different from the gaskinetic BGK scheme [24], the non-equilibrium distribution function at cell interface is approximated by the difference of equilibrium distribution functions between the cell interface and its surrounding points in GKFS. To be specific, the equilibrium distribution functions at the surrounding points of the cell interface are firstly given by interpolation from the conservative variables at cell centers. Then, the equilibrium distribution function at cell interface can be evaluated by a streaming process from the surrounding points. After the above steps, the non-equilibrium distribution function at cell interface can be simply calculated and the explicit formulations for computing the conservative flow variables and fluxes can be derived. It has been proven that GKFS can give the same results and only requires about 60% of the computational time as compared with the conventional gas-kinetic BGK scheme [1]. Inspired by the previous work of GKFS [1], a 3D GKFS is developed in this work. In the scheme, the 3D Navier-Stokes equations are discretized by the finite volume method and the numerical flux across the interface is evaluated by the local solution of 3D Boltzmann equation. Therefore, the present scheme can be viewed as a truly 3D flux solver. At the same time, a coordinate transformation is made at the local cell interface to transform the velocities in the Cartesian coordinate system to the normal and tangential directions of interface. In this way, all the cell interfaces can be treated using the same way for evaluation of conservative variables and numerical fluxes. Like the two-dimensional (2D) case, the non-equilibrium distribution function is approximated by the difference of the equilibrium distribution functions between the cell interface and its surrounding points. It is indicated that the present work is the first time to give explicit formulations for evaluating the conservative variables and numerical fluxes for the 3D viscous flow problems. Like other gas-kinetic schemes, the present scheme can be applied to both incompressible and compressible viscous flow problems without any modification. To validate the developed scheme, both incompressible and compressible viscous test cases are solved, including 3D driven cavity flow, incompressible flow past a stationary sphere, flow around ONERA M6 wing and DLR-F6 wing-body configuration. Numerical results show that the present solver can provide accurate results for both incompressible and compressible flows.

Boltzmann equation and conservative forms of moments for Maxwellian function
With Bhatnagar-Gross-Krook (BGK) collision model [26], the continuous Boltzmann equation in the three-dimensional Cartesian coordinate system can be written as where f is the real particle distribution function and g is the equilibrium particle distribution function. τ is the collision time, which is determined by the dynamic viscosity and the pressure. The right-hand side of the equation is the collision term which alters the distribution function from f to g within a collision time τ. Both f and g are functions of space (x, y, z), time (t) and particle velocity (u, v, w, ξ). The internal degree of freedom K in ξ is determined by the space dimension and the ratio of specific heat with the relation K + D = 2/(γ − 1), where D is the abbreviation of the dimension (D = 3 in three dimension) and γ is the specific heat ratio. The equilibrium state g of the Maxwellian distribution is where ρ is the density of the mean flow; U = (U, V, W) is the macroscopic velocity vector expressed in the x-, y-and z-directions; λ = m/(2kT) = 1/(2RT), where m is the molecular mass, k is the Boltzmann constant, R is the gas constant and T is the temperature. In the equilibrium state, ξ 2 is the abbreviation of With the Maxwellian distribution function in Eq. (2), the following 7 conservative forms of moments will be satisfied, which are used to recover Navier-Stokes equations by Eq. (1) through Chapman-Enskog expansion analysis: where u α , u β , u χ and U α , U β , U χ are the particle velocities and the macroscopic flow velocities in the α-, βand χdirections. p is the pressure and b = K + D represents the total degree of freedoms of molecules. dΞ = du α du β du χ dξ 1 dξ 2 ⋯dξ K is the volume element in the particle velocity space. The integral domain for u α , u β , u χ , ξ 1 , ξ 2 , …, ξ K is from −∞ to +∞. Eqs.

Macroscopic governing equations discretized by finite volume method
In this work, the 3D Navier-Stokes equations are solved using the finite volume discretization with the conservative variables defined at cell centers, which can be written as where W is the vector of conservative variables, Ω and N are the volume and number of interfaces of the control volume, respectively, F i and S i are the flux vector and the area of interface i. It should be noted that the numerical fluxes F i are reconstructed locally at cell interface from the conservative variables W at cell centers. In the gaskinetic scheme, the connection between the distribution function f and the conservative variables is where E ¼ 1 2 ðU 2 þ V 2 þ W 2 þ bRT Þ. φ α is the moment given by With the compatibility condition, Eq. (11) is equivalent to The above equation shows that the non-equilibrium distribution function has no contribution to the calculation of conservative variables.
After evaluation of conservative variables, the flux vector F can also be obtained from the distribution function It should be noted that Eq. (15) is the flux vector of x-direction in the Cartesian coordinate system. In the practical application such as curved boundary problems, we need to calculate the numerical flux in the normal direction of interface F n where u ′ is the particle velocity in the normal direction of interface. Suppose that n 1 = (n 1x , n 1y , n 1z ) is the unit vector in the normal direction of interface and n 2 = (n 2x , n 2y , n 2z ), n 3 = (n 3x , n 3y , n 3z ) are the unit vectors in the tangential directions. Then, the relationship between the particle velocities in the normal and tangential directions (u ′ , v ′ , w ′ ) and the particle velocities in the Cartesian coordinate system (u, v, w) are and similarly Substituting Eq. (18) into Eq. (12), we have φ α ¼ 1 0 0 0 0 0 n 1x n 2x n 3x 0 0 n 1y n 2y n 3y 0 0 n 1z n 2z n 3z 0 With the definition of a new moment and its corresponding flux vector the real flux vector F n can be obtained by substituting Eq. (19) into Eq. (16) and using Eq. (21) The above Eq. (22) shows that the calculation of F n is equivalent to the evaluation of F Ã n and the key issue is to obtain the gas distribution function f. In the next subsection, a 3D GKFS will be introduced to evaluate the gas distribution function f at cell interface.
3 Three-dimensional gas-kinetic flux solver As the flux vector F Ã n is evaluated at the local interface, a local coordinate system is applied in the derivation of distribution function f. It is known that the distribution function f can be separated into two parts, the equilibrium part f eq and the non-equilibrium part f neq with the relationship of Here, the equilibrium part f eq equals to With the Chapman-Enskog expansion analysis, the non-equilibrium distribution function can be approximated as Therefore, the gas distribution function truncated to the Navier-Stokes level becomes By applying the Taylor series expansion in time and physical space, the above equation can be simplified to where f(0, 0, 0, t + δt) is the gas distribution function at local interface; g(0, 0, 0, t + δt) and g(−u ′ δt, −v ′ δt, −w ′ δt, t) are the equilibrium distribution functions at local interface and its surrounding points, respectively. δt is the streaming time step. From Eq. (27), it can be seen that the non-equilibrium distribution f neq is calculated by the difference of equilibrium distribution functions between the interface and its surrounding points, which makes current GKFS be much more straightforward.
In the present work, the conservative variables in Eq. (10) are defined at cell centers. In order to solve Eq. (10) by marching in time, the numerical flux in the normal direction of each cell interface F Ã n should be evaluated first. Suppose that the conservative variables at cell centers and their first order derivatives are already known, the conservative variables at the left and the right sides of an interface can be easily given by interpolation. Then, the equilibrium distribution functions at these two sides of interface can be given via Eq. (2). After that, the second order approximation of g(−u ′ δt, −v ′ δt, −w ′ δt, t) at the time level t can be written as Where g l and g r are the equilibrium distribution functions at the left and the right sides of interface, respectively. Note that in Eq. (28), the equilibrium distribution functions at two sides of interface are not necessarily the same, which means that a possible discontinuity has been taken into account in the form. By substituting Eq. (28) into Eq. (27), we have where H(u ′ ) is the Heaviside function defined as & Equation (29) shows that the full information of distribution function at the interface can be decided once we have the equilibrium distribution function at cell interface and its surrounding points.

Evaluation of conservative variables W * at cell interface
It is known that the non-equilibrium distribution has no influence on the computation of conservative variables, and thus Eq. (14) can be adopted to calculate the conservative variables W * at local interface According to the compatibility condition (see Eq. (13)), by substituting Eq. (27) and Eq. (28) into Eq. (30), we have The above equation shows that the conservative variables at cell interface can be obtained by equilibrium distribution function of the surrounding points. By taking the limit δt → 0 [24], the conservative variables at cell interface can be calculated by The above equation means that the conservative variables at cell interface are simply computed by the reconstructed variables of left and right sides. With parameters defined in the Appendix, the conservative variables W * at cell interface are given by where "· l " and "· r " ("·" stands for any variable) denote the variables at the left and the right sides of interface, respectively.

Evaluation of numerical flux F Ã n at cell interface
As soon as the conservative variables at local interface W * are obtained, the equilibrium distribution function g(0, 0, 0, t + δt) can be known by Eq. (2). Then the numerical flux across the cell interface can be calculated via Eq. (29) Note that g(0, 0, 0, t + δt) is the equilibrium distribution function at the interface and time level t + δt, and g l , g r are the distribution functions at the left and the right sides of interface and the time level t. By taking the limit δt → 0, we have According to the work of Xu [24], ∂g/∂t can be expanded by where A 1 , A 2 , A 3 , A 4 and A 5 are the derivatives of macroscopic variables with respect to physical space, which will be determined from the compatibility condition, Thus, the flux expression in Eq. (38) can be written as Note that δt → 0 has been applied in Eq. (41). In the above equation, the only undetermined variables are the coefficients A 1 , A 2 , A 3 , A 4 and A 5 .
Substituting Eq. (29) into Eq. (11) and adopting the compatibility condition, we have Using Eqs. (39)-(40), the above equation can be written as Defining The explicit formulations of G 1 to G 5 are given in the Appendix. After a similar derivation process to the work of Xu [24], the coefficients A 1 , A 2 , A 3 , A 4 and A 5 can be determined by where Once the above coefficients are obtained, the numerical flux F Ã n across the interface can be calculated via Eq. (41). Similar to the conservative variables W * , the explicit expressions for numerical flux F Ã n can also be given as In Eqs. (51)-(55), the definitions of 〈⋅〉 and all coefficients can be found in the Appendix. They are expressed explicitly as the functions of conservative variables and their derivatives. In addition, since the structured mesh is used in this work, the spatial derivatives in Eqs. (51)-(55) can be approximated directly by the finite difference scheme.
From the above derivations, it can be seen that there are two major differences between the present solver and the gas-kinetic BGK scheme [24]. The first difference is that the local asymptotic solution to the Boltzmann equation (see Eq. (26)) is used to calculate the distribution function on the cell interface for the GKFS, while the local integral solution to the Boltzmann equation is utilized for the gas-kinetic BGK scheme. Another difference is that the non-equilibrium distribution function is approximated by the difference of equilibrium distribution functions on the cell interface and its surrounding streaming nodes in GKFS (see Eq. (27)), while in the gas-kinetic BGK scheme, the non-equilibrium distribution function is included in the initial distribution function around the cell interface. These differences lead to the numerical flux reconstructed by the GKFS being time-independent (see Eq. (41)), while that of the gas-kinetic BGK scheme is time-dependent. Since δt → 0 is adopted, the GKFS actually reconstructs the numerical flux at the time level t, as shown in Eq. (41). In the gas-kinetic BGK scheme, the numerical flux can be viewed as the integral average in the time interval [t, t + Δt]. From this point of view, the temporal accuracy of the flux in GKFS is O(Δt), while that of the gas-kinetic BGK scheme is O(Δt 2 ). But in fact, most of the conventional CFD schemes, such as the Roe scheme, HLL scheme and AUSM, also calculate the numerical flux at the time level t, which indicates that the temporal accuracy of the flux may not be important for solving the Euler/Navier-Stokes equations for general cases. In terms of simplicity, fewer coefficients are involved in GKFS than the gas-kinetic BGK scheme.

Determination of collision time τ
Theoretically, the collision time τ in Eq. (1) is proportional to the physical viscosity where μ is the dynamic viscosity and p is the pressure. However, the numerical dissipation in Eq. (56) might not be sufficient to get a stable solution in cases such as strong shock wave. Therefore, the effective viscosity should be a combination of both physical and numerical ones. Xu [24] proposed a simple and effective treatment to incorporate the numerical viscosity into the gas-kinetic BGK scheme, which is also adopted in the present work: where Δt is the time step in the solution of Navier-Stokes equations, p l and p r are the pressure at the left and the right sides of interface, respectively. The additional term of the above equation corresponds to numerical viscosity, which is applied to take the pressure jump with a thickness in the order of cell size into account.

Prandtl number fix
It is well known that the Prandtl number in the gas-kinetic BGK scheme corresponds to unity [24]. Several approaches are available to make the Prandtl number be consistent with the real problem. BGK-Shakhov model [27] is one of these attempts, which adjusts the heat flux in the relaxation term. In the Shakhov model, the Shakhov equilibrium distribution function is given by where g is the Maxwellian distribution function in Eq. (2), Pr is the Prandtl number, c = u − U is the peculiar velocity and q is the heat flux It can be seen from Eq. (58) that the Prandtl number can be changed to any realistic value easily. However, considerable work has to be devoted to extend the current GKFS to the above Shakhov model.
Another alternative approach is to make correction for heat flux, which has been presented in [24].
where F 5 is the energy flux and q is the heat flux defined in Eq. (59). Since all momentums in Eq. (60) have been obtained in the evaluation of energy flux F 5 , there will not be much additional work in the above Prandtl number fix. Therefore, Eq. (60) is employed to adjust the Prandtl number in the present work.

Computational sequence
In this section, the basic solution procedure of the current 3D GKFS is summarized as follows: (1) Firstly, we need to calculate the derivatives of conservative variables and reconstruct the initial conservative variables at two sides of cell interface.

Numerical results and discussion
To validate the proposed 3D GKFS for simulation of incompressible and compressible viscous flows, the 3D lid-driven cavity flow, incompressible flow past a stationary sphere, flow around ONERA M6 wing and DLR-F6 wing-body configuration are considered. For temporal discretization to Eq. (10), four-stages Runge-Kutta method is applied in cases of 3D lid-driven cavity flow and flow past a stationary sphere. In compressible cases, the Lower-upper symmetric-Gauss-Seidel (LU-SGS) scheme [28] is adopted to accelerate the convergence and the Venkatakrishnan's limiter [29] is used to calculate the conservative variables at two sides of interface W L and W R in the reconstruction process. Specifically, W L and W R are computed by where W L c and W R c are the conservative flow variables at centers of the left and the right cells, respectively; ∇W L c and ∇W R c are their corresponding first-order derivatives. x L c , x R c and x b are the coordinates of the left cell center, the right cell center and the midpoint of cell interface, respectively. Ψ L c and Ψ R c are the limiter functions utilized in the left and the right cells, respectively. In addition, all the simulations were done on a PC with 3.10GHz CPU.
Before applying the GKFS to various fluid flow problems, its accuracy is first validated by the advection of density perturbation for three-dimensional flows [30]. The initial condition of this problem is set as ; v x; y; z ð Þ¼1; w x; y; z ð Þ¼1; p x; y; z ð Þ¼1: ð62Þ The exact solutions under periodic boundary condition are ρ x; y; z; t ð Þ¼1 þ 0:2 sin π x þ y þ z−3t ð Þ ð Þ ; u x; y; z; t ð Þ¼1 ; v x; y; z; t ð Þ¼1; w x; y; z; t ð Þ¼1; p x; y; z; t ð Þ¼1: ð63Þ Since this test case belongs to the inviscid flow, the collision time τ takes

Case 1: 3D lid-driven cavity flow
The 3D lid-driven cavity flows in a cube are simulated to test the capability of the proposed explicit GKFS for simulating 3D incompressible viscous flows. The non-uniform mesh of 81 × 81 × 81 is used for the cases of Re = 100 and 400. The mesh point in the x-direction is generated by Where κ i = (i − 1)/((i max − 1)/2), i and imax are the mesh point index and total number of mesh points in the x direction; η is the parameter to control the mesh Fig. 1 L 1 error of the density field for the advection of density perturbation stretching and is selected as 1.1 in this study. Similarly, the mesh point in the y-and zdirections is generated in the same way.
In the current simulation, the fluid density is taken as ρ = 1.0 and the lid velocity is chosen as U ∞ = 0.1. Initially, the density inside the cavity is constant and the flow is static. The lid on the top boundary moves along the x-direction. The no-slip wall condition is imposed at all boundaries. To quantitatively examine the performance of 3D GKFS, the velocity profiles of x-direction component u along the vertical centerline and y-direction component v along the horizontal centerline for Re = 100 and 400 are plotted in Fig. 2. For comparison, the results of Shu et al. [31] and Wu and Shu [32] are also included in the figure. It can be found that all the velocity profiles by current 3D GKFS agree very well with those of Shu et al. [31] and Wu and Shu [32], which demonstrates the capability of present solver for the simulation of 3D incompressible flows on non-uniform grids. To further show the flow patterns of 3D lid-driven cavity flow, the streamlines for Re = 100 and 400 at three orthogonal mid-planes located at x = 0.5, y = 0.5 and z = 0.5 are displayed in Fig. 3. The flow patterns along the mid-plane of z = 0.5 in Fig. 3 demonstrate that the primary vortices gradually shift toward the center position and the second vortices gradually move to the lower bottom wall when the Reynolds number is increased. In this process, the strength of these vortices is also enhanced, which can also be proven by the flow patterns along other two mid-planes. All these observations match well with those in Shu et al. [31].

Case 2: incompressible flow past a stationary sphere
In this section, the 3D GKFS is applied to a benchmark case of incompressible flow past a stationary sphere. In this case, the flow is characterized by the Reynolds number defined by Re = ρU ∞ D/μ, where ρ and μ are the fluid density and dynamic viscosity, respectively. U ∞ is the free stream velocity and D is the sphere diameter. To simulate this test case with a simple Cartesian mesh, the implicit boundary condition-enforced immersed boundary method [33,34] is coupled with the present 3D GKFS. The computational domain is selected as a rectangular box of 30D × 20D × 20D in the x-, y-and z-directions. The sphere is initially placed at (10D, 10D, 10D), which is discretized by triangular elements with 1195 vertices. As shown in Fig. 4, a non-uniform Cartesian mesh with mesh size of 137 × 122 × 122 is used, in which a uniform mesh spacing of 0.02D is applied around the sphere. The no-slip condition on the curved boundary is imposed by correcting the velocity on the Cartesian mesh through the immersed boundary method [33,34]. Here, laminar flows at low Reynolds numbers of 50, 100, 150, 200 and 250 are considered.
At first, the drag coefficients at Re = 100, 200 and 250 are computed and compared quantitatively in Table 1 to verify the accuracy of the present solver. The numerical results of Johnson and Patel [35], Wu and Shu [32], Kim et al. [36] and Wang et al. [37] are also included in the table for comparison. It can be clearly observed that the present results match well with those in the literature.
Then, for the steady axisymmetric flow, the streamlines of flow past a sphere at Re ≤ 200 are depicted in Fig. 5. Since the flow is axisymmetric, only the streamlines on the x-y plane of symmetry are given. From the figure, a recirculation region appears behind the sphere and its length L s increases with Reynolds number. Quantitative comparison between the present results of L s and those of Johnson and Patel [35] and Gilmanov et al. [38] is made in Fig. 6. Good agreement can be found in the figure. When the Reynolds number is increased to 250, the phenomenon of steady nonaxisymmetric pattern shows up, which can be seen in Fig. 7. In the figure, the streamlines on the x-z plane remain symmetric. However, there are two asymmetric vortices on the x-y plane, which implies that the symmetry is lost in this plane. These results are in good agreement with previous investigations [35,38].

Case 3: flow around ONERA M6 wing
The ONERA M6 test case is chosen to validate the present solver for the simulation of compressible viscous flows with complex geometry. For numerical simulation, the freestream Mach number is taken as M ∞ = 0.8395, the mean-chord based Reynolds number is chosen as Re = 11.72 × 10 6 and the angle of attack is α = 3.06 ∘ . The computational mesh in the NASA website [39] is adopted in this work, which has 4 blocks and 316, 932 grid points. The mesh spacing of the first mesh point adjacent to the wing surface is 4.5 × 10 −5 . To take turbulent effect into consideration, the Spalart-Allmaras turbulent model [40] is applied. Figure 8 shows the pressure contours at the wing surface obtained from the present solver, in which the "λ" shape shock wave on the upper surface is clearly presented. The above phenomenon matches well with the result from sphere function-based gas-kinetic scheme [41]. To further validate the present results, the pressure coefficient distributions at selected span-wise locations obtained from the present solver are displayed in Fig. 9. The numerical results of WIND scheme [39] and the experimental results [42] are also included for comparison. As can be seen from the figure, the present results are close to those of WIND scheme [39] and compare well with the experimental data [42]. What is more, the pressure coefficient distributions at 65% and 80% spans show that the present results are much closer to the experimental results [42] as compared with the results from WIND scheme [39]. It demonstrates that the present solver captures the shock wave more precisely and controls the numerical dissipation well.
To further investigate the performance of GKFS for simulation of high speed flows, in this test, we change the Mach number to M ∞ = 5 while keep other parameters the same as the above case. Figure 10 shows the pressure contours and the pressure

Case 4: DLR-F6 wing-body configuration
The DLR-F6 wing-body configuration is a generic transport aircraft model from the 3rd AIAA CFD drag prediction workshop (DPW III) [44]. At first, numerical simulations are conducted at a free-stream Mach number of M ∞ = 0.75, a mean-chord based Reynolds number of Re = 3 × 10 6 and an angle of attack α = 0.49 ∘ . The geometry and computational mesh from the NASA website [45] are utilized in the current work.
Owing to the limitation of the computer's memory, only the coarse mesh with 26 blocks and 2,298,880 cells is used. Figure 11 is the pressure contour of DLR-F6 wingbody obtained by present GKFS. The separation bubble at the intersection of wing and body is clearly recognized in Fig. 12, which is in line with the observations of Vassberg et al. [46]. To make a quantitative comparison, the pressure coefficient distributions at selected span-wise locations obtained by present 3D GKFS are compared with the experimental results [47] and numerical results of Vassberg et al. [46] and Yang et al. [48] in Fig. 13. It can be observed that the current results are close to those of Vassberg et al.  [46] and Yang et al. [48] and all of them basically agree well with the experimental measurement [47].
To further verify the force coefficients of current solver for the DLR-F6 wing-body, another test case is simulated with the free stream condition of Mach number M ∞ = 0.75, Reynolds number Re = 5 × 10 6 and angle of attach α = 0 ∘ . Table 2 shows the present results of force coefficients, including lift coefficient C l , pressure drag coefficient C d, p , friction drag coefficient C d, f , total drag coefficient C d and moment coefficient C M . The results of present solver are close to the results of LBFS [48] and can essentially match well with the reference data of Vassberg et al. [46]. This paper presents a three-dimensional GKFS for simulation of incompressible and compressible viscous flows. The present work is the extension of our previous work [1], where a new gas-kinetic scheme is presented to simulate two-dimensional viscous flows. In this work, the non-equilibrium distribution function is evaluated by the difference of equilibrium distribution functions at cell interface and its surrounding points. As a result, the distribution function at the interface can be simply derived and the formulations of the conservative variables and fluxes at cell interface can be explicitly given. Since the solution of 3D continuous Boltzmann equation is reconstructed locally at cell interface, the present scheme can be viewed as a truly 3D flux solver for viscous flows. To consider general 3D cases, a local coordinate transformation is made to transform the velocities in the global Cartesian coordinate system to the local normal and tangential directions at each cell interface. In this way, all the interfaces can be treated using the same way. Several numerical experiments are conducted to validate the proposed scheme, including 3D lid-driven cavity flow, incompressible flow past a stationary sphere, compressible flow around ONERA M6 wing and DLR-F6 wing-body configuration. Numerical results showed that the proposed flux solver can provide accurate numerical results for three-dimensional incompressible and compressible viscous flows.

Moments of Maxwellian Distribution Function
In the paper, some notations are taken to simplify the formulations. In this appendix, the notations for the moments of Maxwellian distribution function are introduced. At first, the Maxwellian distribution function for 3D flows is given as (Eq. (2)) Following the idea of [2], the notation for the moments of g is defined as It should be noted that in Section 3, the conservative variables and numerical fluxes are derived at local interface and the transformation of velocities is made at each cell  interface. Therefore, the moments of normal and tangential velocities (u ′ , v ′ , w ′ ) are presented to keep pace with the formulations in the paper. As the integration of three components of the particle velocity (u ′ , v ′ , w ′ ) are similar to each other, only the moments of u ′ are presented here. When the integral of velocity is from −∞ to +∞, the moments of u ′n and ξ n are ðA:8Þ