Skip to main content

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

Abstract

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.

1 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 gas-kinetic 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.

2 Boltzmann equation, Maxwellian distribution function and Navier-Stokes equations

2.1 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

$$ \frac{\partial f}{\partial t}+u\frac{\partial f}{\partial x}+v\frac{\partial f}{\partial y}+w\frac{\partial f}{\partial z}=\frac{g-f}{\tau }, $$
(1)

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

$$ g=\rho {\left(\frac{\lambda }{\pi}\right)}^{\frac{K+3}{2}}{e}^{-\lambda \left({\left(u-U\right)}^2+{\left(v-V\right)}^2+{\left(w-W\right)}^2+{\xi}^2\right)}, $$
(2)

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 \( {\xi}^2={\xi}_1^2+{\xi}_2^2+\cdots +{\xi}_K^2. \)

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:

$$ \int gd\Xi =\rho, $$
(3)
$$ \int {gu}_{\alpha }d\Xi =\rho {U}_{\alpha }, $$
(4)
$$ \int g\left({u}_{\alpha }{u}_{\alpha }+{\xi}^2\right)d\Xi =\rho \left({U}_{\alpha }{U}_{\alpha }+ bRT\right), $$
(5)
$$ \int {gu}_{\alpha }{u}_{\beta }d\Xi =\rho {U}_{\alpha }{U}_{\beta }+p{\delta}_{\alpha \beta}, $$
(6)
$$ \int g\left({u}_{\alpha }{u}_{\alpha }+{\xi}^2\right){u}_{\beta }d\Xi =\rho \left[{U}_{\alpha }{U}_{\alpha }+\left(b+2\right) RT\right]{U}_{\beta }, $$
(7)
$$ \int {gu}_{\alpha }{u}_{\beta }{u}_{\chi }d\Xi =p\left({U}_{\alpha }{\delta}_{\beta \chi}+{U}_{\beta }{\delta}_{\chi \alpha}+{U}_{\chi }{\delta}_{\alpha \beta}\right)+\rho {U}_{\alpha }{U}_{\beta }{U}_{\chi }, $$
(8)
$$ {\displaystyle \begin{array}{l}\int g\left({u}_{\alpha }{u}_{\alpha }+{\xi}^2\right){u}_{\beta }{u}_{\chi }d\Xi \\ {}=\rho \left\{{U}_{\alpha }{U}_{\alpha }{U}_{\beta }{U}_{\chi }+\left[\left(b+4\right){U}_{\beta }{U}_{\chi }+{U}_{\alpha }{U}_{\alpha }{\delta}_{\beta \chi}\right] RT+\left(b+2\right){R}^2{T}^2{\delta}_{\beta \chi}\right\},\end{array}} $$
(9)

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χ12K is the volume element in the particle velocity space. The integral domain for uα, uβ, uχ, ξ1, ξ2, …, ξK is from −∞ to +∞. Eqs. (3)–(5) are applied to recover the fluid density, momentum and energy, respectively. Eqs. (6) and (7) are used to recover convective fluxes of the momentum equation and the energy equation. Eqs. (8) and (9) are to recover diffusive fluxes of the momentum equation and the energy equation.

2.2 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

$$ \frac{d\mathbf{W}}{dt}+\frac{1}{\Omega}\sum \limits_{i=1}^N{\mathbf{F}}_i{S}_i=0, $$
(10)

where W is the vector of conservative variables, Ω and N are the volume and number of interfaces of the control volume, respectively, Fi and Si are the flux vector and the area of interface i. It should be noted that the numerical fluxes Fi are reconstructed locally at cell interface from the conservative variables W at cell centers. In the gas-kinetic scheme, the connection between the distribution function f and the conservative variables is

$$ \mathbf{W}={\left(\rho, \kern0.5em \rho U,\kern0.5em \rho V,\kern0.5em \rho W,\kern0.5em \rho E\right)}^T=\int {\boldsymbol{\upvarphi}}_{\alpha } fd\Xi, $$
(11)

where \( E=\frac{1}{2}\left({U}^2+{V}^2+{W}^2+ bRT\right) \). φα is the moment given by

$$ {\boldsymbol{\upvarphi}}_{\alpha }={\left(1,\kern0.5em u,\kern0.5em v,\kern0.5em w,\kern0.5em \frac{1}{2}\left({u}^2+{v}^2+{w}^2+{\xi}^2\right)\right)}^T. $$
(12)

With the compatibility condition,

$$ \int {\boldsymbol{\upvarphi}}_{\alpha}\frac{g-f}{\tau }d\Xi =0, $$
(13)

Eq. (11) is equivalent to

$$ \mathbf{W}={\left(\rho, \kern0.5em \rho U,\kern0.5em \rho V,\kern0.5em \rho W,\kern0.5em \rho E\right)}^T=\int {\boldsymbol{\upvarphi}}_{\alpha } gd\Xi . $$
(14)

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

$$ \mathbf{F}=\int u{\boldsymbol{\upvarphi}}_{\alpha } fd\Xi . $$
(15)

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 Fn

$$ {\mathbf{F}}_n={\left({F}_1,\kern0.5em {F}_2,\kern0.5em {F}_3,\kern0.5em {F}_4,\kern0.5em {F}_5\right)}^T=\int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha } fd\Xi, $$
(16)

where u is the particle velocity in the normal direction of interface. Suppose that n1 = (n1x,  n1y,  n1z) is the unit vector in the normal direction of interface and n2 = (n2x,  n2y,  n2z), n3 = (n3x,  n3y,  n3z) 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

$$ {u}^{\prime }={un}_{1x}+{vn}_{1y}+{wn}_{1z},\kern1.5em {v}^{\prime }={un}_{2x}+{vn}_{2y}+{wn}_{2z},\kern1.5em {w}^{\prime }={un}_{3x}+{vn}_{3y}+{wn}_{3z}, $$
(17)

and similarly

$$ u={u}^{\prime }{n}_{1x}+{v}^{\prime }{n}_{2x}+{w}^{\prime }{n}_{3x},\kern1.5em v={u}^{\prime }{n}_{1y}+{v}^{\prime }{n}_{2y}+{w}^{\prime }{n}_{3y},\kern1.5em w={u}^{\prime }{n}_{1z}+{v}^{\prime }{n}_{2z}+{w}^{\prime }{n}_{3z}. $$
(18)

Substituting Eq. (18) into Eq. (12), we have

$$ {\boldsymbol{\upvarphi}}_{\alpha }=\left(\begin{array}{ccccc}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\\ {}0& 0& 0& 0& 1\end{array}\right){\left(1,\kern0.5em {u}^{\prime },\kern0.5em {v}^{\prime },\kern0.5em {w}^{\prime },\kern0.5em \frac{1}{2}\left({u^{\prime}}^2+{v^{\prime}}^2+{w^{\prime}}^2+{\xi}^2\right)\right)}^T. $$
(19)

With the definition of a new moment

$$ {\boldsymbol{\upvarphi}}_{\alpha}^{\ast }={\left(1,\kern0.5em {u}^{\prime },\kern0.5em {v}^{\prime },\kern0.5em {w}^{\prime },\kern0.5em \frac{1}{2}\left({u^{\prime}}^2+{v^{\prime}}^2+{w^{\prime}}^2+{\xi}^2\right)\right)}^T, $$
(20)

and its corresponding flux vector

$$ {\mathbf{F}}_n^{\ast }={\left({F}_1^{\ast },\kern0.5em {F}_2^{\ast },\kern0.5em {F}_3^{\ast },\kern0.5em {F}_4^{\ast },\kern0.5em {F}_5^{\ast}\right)}^T=\int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast } fd\Xi, $$
(21)

the real flux vector Fn can be obtained by substituting Eq. (19) into Eq. (16) and using Eq. (21)

$$ {\mathbf{F}}_n=\int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha } fd\Xi =\left(\begin{array}{ccccc}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\\ {}0& 0& 0& 0& 1\end{array}\right){\mathbf{F}}_n^{\ast }. $$
(22)

The above Eq. (22) shows that the calculation of Fn is equivalent to the evaluation of \( {\mathbf{F}}_n^{\ast } \) 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 \( {\mathbf{F}}_n^{\ast } \) 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 feq and the non-equilibrium part fneq with the relationship of

$$ f={f}^{eq}+{f}^{neq}. $$
(23)

Here, the equilibrium part feq equals to

$$ {f}^{eq}=g. $$
(24)

With the Chapman-Enskog expansion analysis, the non-equilibrium distribution function can be approximated as

$$ {f}^{neq}=-\tau \left(\frac{\partial }{\partial t}+\mathbf{u}\cdot \nabla \right){f}^{eq}=-\tau \left(\frac{\partial }{\partial t}+\mathbf{u}\cdot \nabla \right)g. $$
(25)

Therefore, the gas distribution function truncated to the Navier-Stokes level becomes

$$ f={f}^{eq}+{f}^{neq}=g-\tau \left(\frac{\partial g}{\partial t}+{u}^{\prime}\frac{\partial g}{\partial {n}_1}+{v}^{\prime}\frac{\partial g}{\partial {n}_2}+{w}^{\prime}\frac{\partial g}{\partial {n}_3}\right). $$
(26)

By applying the Taylor series expansion in time and physical space, the above equation can be simplified to

$$ {\displaystyle \begin{array}{l}f\left(0,0,0,t+\delta t\right)\\ {}=g\left(0,0,0,t+\delta t\right)-\frac{\tau }{\delta t}\left[g\left(0,0,0,t+\delta t\right)-g\left(-{u}^{\prime}\delta t,-{v}^{\prime}\delta t,-{w}^{\prime}\delta t,t\right)\right],\end{array}} $$
(27)

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 fneq 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 \( {\mathbf{F}}_n^{\ast } \) 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

$$ g\left(-{u}^{\prime}\delta t,-{v}^{\prime}\delta t,-{w}^{\prime}\delta t,t\right)=\left\{\begin{array}{c}{g}_l-\frac{\partial {g}_l}{\partial {n}_1}{u}^{\prime}\delta t-\frac{\partial {g}_l}{\partial {n}_2}{v}^{\prime}\delta t-\frac{\partial {g}_l}{\partial {n}_3}{w}^{\prime}\delta t,\kern2em {u}^{\prime}\ge 0,\\ {}{g}_r-\frac{\partial {g}_r}{\partial {n}_1}{u}^{\prime}\delta t-\frac{\partial {g}_r}{\partial {n}_2}{v}^{\prime}\delta t-\frac{\partial {g}_r}{\partial {n}_3}{w}^{\prime}\delta t,\kern2em {u}^{\prime }<0.\end{array}\right. $$
(28)

Where gl and gr 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

$$ {\displaystyle \begin{array}{l}f\left(0,0,0,t+\delta t\right)=g\left(0,0,0,t+\delta t\right)\\ {}\kern2em -\frac{\tau }{\delta t}\left[g\left(0,0,0,t+\delta t\right)-{g}_lH\left({u}^{\prime}\right)-{g}_r\left(1-H\left({u}^{\prime}\right)\right)\right]\\ {}\kern2em -\tau \left[\left(\frac{\partial {u}^{\prime }{g}_l}{\partial {n}_1}+\frac{\partial {v}^{\prime }{g}_l}{\partial {n}_2}+\frac{\partial {w}^{\prime }{g}_l}{\partial {n}_3}\right)H\left({u}^{\prime}\right)+\left(\frac{\partial {u}^{\prime }{g}_r}{\partial {n}_1}+\frac{\partial {v}^{\prime }{g}_r}{\partial {n}_2}+\frac{\partial {w}^{\prime }{g}_r}{\partial {n}_3}\right)\left(1-H\left({u}^{\prime}\right)\right)\right],\end{array}} $$
(29)

where H(u) is the Heaviside function defined as

$$ H\left({u}^{\prime}\right)=\left\{\begin{array}{c}0,\kern2.5em {u}^{\prime }<0,\\ {}1,\kern2.5em {u}^{\prime}\ge 0.\end{array}\right. $$

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.

3.1 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

$$ {\mathbf{W}}^{\ast }={\left(\rho, \kern0.5em \rho {U}^{\prime },\kern0.5em \rho {V}^{\prime },\kern0.5em \rho {W}^{\prime },\kern0.5em \rho E\right)}^T=\int g{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }d\Xi . $$
(30)

According to the compatibility condition (see Eq. (13)), by substituting Eq. (27) and Eq. (28) into Eq. (30), we have

$$ {\displaystyle \begin{array}{c}{\mathbf{W}}^{\ast }=\int {\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(0,0,0,t+\delta t\right)d\Xi =\int {\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(-{u}^{\prime}\delta t,-{v}^{\prime}\delta t,-{w}^{\prime}\delta t,t\right)d\Xi \\ {}\kern1em =\int {\int}_{u^{\prime }>0}{\boldsymbol{\upvarphi}}_{\alpha}^{\ast}\left({g}_l-\frac{\partial {g}_l}{\partial {n}_1}{u}^{\prime}\delta t-\frac{\partial {g}_l}{\partial {n}_2}{v}^{\prime}\delta t-\frac{\partial {g}_l}{\partial {n}_3}{w}^{\prime}\delta t\right)d\Xi \\ {}+\int {\int}_{u^{\prime }<0}{\boldsymbol{\upvarphi}}_{\alpha}^{\ast}\left({g}_r-\frac{\partial {g}_r}{\partial {n}_1}{u}^{\prime}\delta t-\frac{\partial {g}_r}{\partial {n}_2}{v}^{\prime}\delta t-\frac{\partial {g}_r}{\partial {n}_3}{w}^{\prime}\delta t\right)d\Xi .\end{array}} $$
(31)

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

$$ {\mathbf{W}}^{\ast }=\int {\int}_{u^{\prime }>0}{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_ld\Xi +\int {\int}_{u^{\prime }<0}{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_rd\Xi . $$
(32)

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

$$ \rho =\left({\rho}_l{a}_l+{\rho}_r{a}_r\right), $$
(33)
$$ \rho {U}^{\prime }=\left({\rho}_l{b}_l+{\rho}_r{b}_r\right), $$
(34)
$$ \rho {V}^{\prime }=\left({\rho}_l{V_l}^{\prime }{a}_l+{\rho}_r{V_r}^{\prime }{a}_r\right), $$
(35)
$$ \rho {W}^{\prime }=\left({\rho}_l{W_l}^{\prime }{a}_l+{\rho}_r{W_r}^{\prime }{a}_r\right), $$
(36)
$$ {\displaystyle \begin{array}{l}\rho E=\frac{1}{2}{\rho}_l\left[{c}_l+\left({V^{\prime}}_l^2+{W^{\prime}}_l^2+\left(b-1\right){RT}_l\right){a}_l\right]\\ {}\kern2.5em +\frac{1}{2}{\rho}_r\left[{c}_r+\left({V^{\prime}}_r^2+{W^{\prime}}_r^2+\left(b-1\right){RT}_r\right){a}_r\right],\end{array}} $$
(37)

where “·l” and “·r” (“·” stands for any variable) denote the variables at the left and the right sides of interface, respectively.

3.2 Evaluation of numerical flux \( {\mathbf{F}}_n^{\ast } \) 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)

$$ {\displaystyle \begin{array}{l}{\mathbf{F}}_n^{\ast }=\int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }f\left(0,0,0,t+\delta t\right)d\Xi \\ {}\kern2.5em =\int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(0,0,0,t+\delta t\right)d\Xi -\frac{\tau }{\delta t}\left[\int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(0,0,0,t+\delta t\right)d\Xi \right.\\ {}\kern2.5em -\int {\int}_{u^{\prime }>0}{u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_ld\Xi -\left.\int {\int}_{u^{\prime }<0}{u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_rd\Xi \right]\\ {}\kern2.5em -\tau \left[\frac{\partial }{\partial {n}_1}\int \left({\int}_{u^{\prime }>0}{u^{\prime}}^2{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u^{\prime }<0}{u^{\prime}}^2{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \right.\\ {}\kern4.5em +\frac{\partial }{\partial {n}_2}\int \left({\int}_{u^{\prime }>0}{u}^{\prime }{v}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u<0}{u}^{\prime }{v}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \\ {}\kern4.5em +\left.\frac{\partial }{\partial {n}_3}\int \left({\int}_{u^{\prime }>0}{u}^{\prime }{w}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u<0}{u}^{\prime }{w}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \right].\end{array}} $$
(38)

Note that g(0, 0, 0, t + δt) is the equilibrium distribution function at the interface and time level t + δt, and gl, gr 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

$$ {\displaystyle \begin{array}{l}-\frac{\tau }{\delta t}\left[\int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(0,0,0,t+\delta t\right)d\Xi \right.-\int {\int}_{u^{\prime }>0}{u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_ld\Xi -\left.\int {\int}_{u^{\prime }<0}{u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_rd\Xi \right]\\ {}=-\tau \int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast}\frac{\partial g\left(0,0,0,t\right)}{\partial t}d\Xi .\end{array}} $$
(39)

According to the work of Xu [24], ∂g/∂t can be expanded by

$$ \frac{\partial g\left(0,0,0,t\right)}{\partial t}=g\left(0,0,0,t\right)\left({A}_1+{A}_2{u}^{\prime }+{A}_3{v}^{\prime }+{A}_4{w}^{\prime }+{A}_5\varepsilon \right), $$
(40)

where A1, A2, A3, A4 and A5 are the derivatives of macroscopic variables with respect to physical space, which will be determined from the compatibility condition, \( \varepsilon =\frac{1}{2}\left({u^{\prime}}^2+{v^{\prime}}^2+{w^{\prime}}^2+{\xi}^2\right) \). Thus, the flux expression in Eq. (38) can be written as

$$ {\displaystyle \begin{array}{l}{\mathbf{F}}_n^{\ast }=\int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(0,0,0,t\right)d\Xi \\ {}\kern2.5em -\tau \int {u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(0,0,0,t\right)\left({A}_1+{A}_2{u}^{\prime }+{A}_3{v}^{\prime }+{A}_4{w}^{\prime }+{A}_5\varepsilon \right)d\Xi \\ {}\kern2.5em -\tau \left[\frac{\partial }{\partial {n}_1}\int \left({\int}_{u^{\prime }>0}{u^{\prime}}^2{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u^{\prime }<0}{u^{\prime}}^2{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \right.\\ {}\kern4.5em +\frac{\partial }{\partial {n}_2}\int \left({\int}_{u^{\prime }>0}{u}^{\prime }{v}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u<0}{u}^{\prime }{v}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \\ {}\kern4.5em +\left.\frac{\partial }{\partial {n}_3}\int \left({\int}_{u^{\prime }>0}{u}^{\prime }{w}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u<0}{u}^{\prime }{w}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \right].\end{array}} $$
(41)

Note that δt → 0 has been applied in Eq. (41). In the above equation, the only undetermined variables are the coefficients A1, A2, A3, A4 and A5.

Substituting Eq. (29) into Eq. (11) and adopting the compatibility condition, we have

$$ {\displaystyle \begin{array}{l}\frac{1}{\delta t}\left[\int {\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(0,0,0,t+\delta t\right)d\Xi -\int {\int}_{u^{\prime }>0}{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_ld\Xi -\int {\int}_{u^{\prime }<0}{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_rd\Xi \right]\\ {}=-\left[\int {\int}_{u^{\prime }>0}{\boldsymbol{\upvarphi}}_{\alpha}^{\ast}\left(\frac{\partial {g}_l}{\partial {n}_1}{u}^{\prime }+\frac{\partial {g}_l}{\partial {n}_2}{v}^{\prime }+\frac{\partial {g}_l}{\partial {n}_3}{w}^{\prime}\right)d\Xi \right.\\ {}\kern1.5em +\left.\int {\int}_{u^{\prime }<0}{\boldsymbol{\upvarphi}}_{\alpha}^{\ast}\left(\frac{\partial {g}_r}{\partial {n}_1}{u}^{\prime }+\frac{\partial {g}_r}{\partial {n}_2}{v}^{\prime }+\frac{\partial {g}_r}{\partial {n}_3}{w}^{\prime}\right)d\Xi \right].\end{array}} $$
(42)

Using Eqs. (39)–(40), the above equation can be written as

$$ {\displaystyle \begin{array}{l}\int {\boldsymbol{\upvarphi}}_{\alpha}^{\ast }g\left(0,0,0,t\right)\left({A}_1+{A}_2{u}^{\prime }+{A}_3{v}^{\prime }+{A}_4{w}^{\prime }+{A}_5\varepsilon \right)d\Xi \\ {}=-\left[\frac{\partial }{\partial {n}_1}\int \left({\int}_{u^{\prime }>0}{u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u^{\prime }<0}{u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \right.+\frac{\partial }{\partial {n}_2}\int \left({\int}_{u^{\prime }>0}{v}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u^{\prime }<0}{v}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \\ {}\kern2em +\left.\frac{\partial }{\partial {n}_3}\int \left({\int}_{u^{\prime }>0}{w}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u^{\prime }<0}{w}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \right].\end{array}} $$
(43)

Defining

$$ {\displaystyle \begin{array}{c}\frac{\partial }{\partial {n}_1}\int \left({\int}_{u^{\prime }>0}{u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u^{\prime }<0}{u}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi +\frac{\partial }{\partial {n}_2}\int \left({\int}_{u^{\prime }>0}{v}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u^{\prime }<0}{v}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \\ {}+\frac{\partial }{\partial {n}_3}\int \left({\int}_{u^{\prime }>0}{w}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_l+{\int}_{u^{\prime }<0}{w}^{\prime }{\boldsymbol{\upvarphi}}_{\alpha}^{\ast }{g}_r\right)d\Xi \end{array}}=\left(\begin{array}{c}{G}_1\\ {}{G}_2\\ {}{G}_3\\ {}{G}_4\\ {}{G}_5\end{array}\right). $$
(44)

The explicit formulations of G1 to G5 are given in the Appendix. After a similar derivation process to the work of Xu [24], the coefficients A1, A2, A3, A4 and A5 can be determined by

$$ {A}_5=-\frac{8{\lambda}^2}{\left(K+3\right)\rho}\left[{G}_5-{U}^{\prime }{G}_2-{V}^{\prime }{G}_3-{W}^{\prime }{G}_4-\left({\Re}_1-{U^{\prime}}^2-{V^{\prime}}^2-{W^{\prime}}^2\right){G}_1\right], $$
(45)
$$ {A}_4=-\frac{2\lambda }{\rho}\left({G}_4-{W}^{\prime }{G}_1\right)-{W}^{\prime }{A}_5, $$
(46)
$$ {A}_3=-\frac{2\lambda }{\rho}\left({G}_3-{V}^{\prime }{G}_1\right)-{V}^{\prime }{A}_5, $$
(47)
$$ {A}_2=-\frac{2\lambda }{\rho}\left({G}_2-{U}^{\prime }{G}_1\right)-{U}^{\prime }{A}_5, $$
(48)
$$ {A}_1=-\frac{1}{\rho }{G}_1-{U}^{\prime }{A}_2-{V}^{\prime }{A}_3-{W}^{\prime }{A}_4-{\Re}_1{A}_5, $$
(49)

where

$$ {\Re}_1=\frac{1}{2}\left({U^{\prime}}^2+{V^{\prime}}^2+{W^{\prime}}^2+\frac{K+3}{2\lambda}\right). $$
(50)

Once the above coefficients are obtained, the numerical flux \( {\mathbf{F}}_n^{\ast } \) across the interface can be calculated via Eq. (41). Similar to the conservative variables W, the explicit expressions for numerical flux \( {\mathbf{F}}_n^{\ast } \) can also be given as

$$ {F}_1^{\ast }=\rho {U}^{\prime }, $$
(51)
$$ {\displaystyle \begin{array}{l}{F}_2^{\ast }=\left(\rho {U^{\prime}}^2+p\right)-\tau \rho \left[{A}_1\left\langle {u^{\prime}}^2\right\rangle +{A}_2\left\langle {u^{\prime}}^3\right\rangle +{A}_3\left\langle {u^{\prime}}^2\right\rangle \left\langle {v^{\prime}}^1\right\rangle +{A}_4\left\langle {u^{\prime}}^2\right\rangle \left\langle {w^{\prime}}^1\right\rangle \right.\\ {}\kern3.5em +\left.\frac{1}{2}{A}_5\left(\left\langle {u^{\prime}}^4\right\rangle +\left\langle {u^{\prime}}^2\right\rangle \left\langle {v^{\prime}}^2\right\rangle +\left\langle {u^{\prime}}^2\right\rangle \left\langle {w^{\prime}}^2\right\rangle +\left\langle {u^{\prime}}^2\right\rangle \left\langle {\xi}^2\right\rangle \right)\right]\\ {}\kern3.5em -\tau \left[\frac{\partial \left({\rho}_l{d}_l+{\rho}_r{d}_r\right)}{\partial {n}_1}+\frac{\partial \left({\rho}_l{V}_l^{\prime }{c}_l+{\rho}_r{V}_r^{\prime }{c}_r\right)}{\partial {n}_2}+\frac{\partial \left({\rho}_l{W}_l^{\prime }{c}_l+{\rho}_r{W}_r^{\prime }{c}_r\right)}{\partial {n}_3}\right],\end{array}} $$
(52)
$$ {\displaystyle \begin{array}{l}{F}_3^{\ast }=\rho {U}^{\prime }{V}^{\prime }-\tau \rho \left[{A}_1\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^1\right\rangle \right.+{A}_2\left\langle {u^{\prime}}^2\right\rangle \left\langle {v^{\prime}}^1\right\rangle +{A}_3\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^2\right\rangle +{A}_4\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^1\right\rangle \\ {}\kern4em +\frac{1}{2}\left.{A}_5\left(\left\langle {u^{\prime}}^3\right\rangle \left\langle {v^{\prime}}^1\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^3\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^2\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^1\right\rangle \left\langle {\xi}^2\right\rangle \right)\right]\\ {}\kern4em -\tau \left[\frac{\partial \left({\rho}_l{V}_l^{\prime }{c}_l+{\rho}_r{V}_r^{\prime }{c}_r\right)}{\partial {n}_1}+\frac{\partial \left[\left({\rho}_l{V^{\prime}}_l^2+{p}_l\right){b}_l+\left({\rho}_r{V^{\prime}}_r^2+{p}_r\right){b}_r\right]}{\partial {n}_2}\right.\\ {}\kern6.5em \left.+\frac{\partial \left({\rho}_l{V_l}^{\prime }{W_l}^{\prime }{b}_l+{\rho}_r{V_r}^{\prime }{W_r}^{\prime }{b}_r\right)}{\partial {n}_3}\right],\end{array}} $$
(53)
$$ {\displaystyle \begin{array}{l}{F}_4^{\ast }=\rho {U}^{\prime }{W}^{\prime }-\tau \rho \left[{A}_1\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^1\right\rangle \right.+{A}_2\left\langle {u^{\prime}}^2\right\rangle \left\langle {w^{\prime}}^1\right\rangle +{A}_3\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^1\right\rangle +{A}_4\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^2\right\rangle \\ {}\kern4em +\frac{1}{2}\left.{A}_5\left(\left\langle {u^{\prime}}^3\right\rangle \left\langle {w^{\prime}}^1\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^2\right\rangle \left\langle {w^{\prime}}^1\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^3\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^1\right\rangle \left\langle {\xi}^2\right\rangle \right)\right]\\ {}\kern4em -\tau \left[\frac{\partial \left({\rho}_l{W}_l^{\prime }{c}_l+{\rho}_r{W}_r^{\prime }{c}_r\right)}{\partial {n}_1}\right.+\frac{\partial \left({\rho}_l{V_l}^{\prime }{W_l}^{\prime }{b}_l+{\rho}_r{V_r}^{\prime }{W_r}^{\prime }{b}_r\right)}{\partial {n}_2}\\ {}\kern6.5em \left.+\frac{\partial \left[\left({\rho}_l{W^{\prime}}_l^2+{p}_l\right){b}_l+\left({\rho}_r{W^{\prime}}_r^2+{p}_r\right){b}_r\right]}{\partial {n}_3}\right],\end{array}} $$
(54)
$$ {\displaystyle \begin{array}{l}{F}_5^{\ast }=\left(\rho E+p\right){U}^{\prime }-\frac{1}{2}\tau \rho \left\{{A}_1\left[\left\langle {u^{\prime}}^3\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^2\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^2\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {\xi}^2\right\rangle \right]\right.\\ {}\kern2.5em +{A}_2\left[\left\langle {u^{\prime}}^4\right\rangle +\left\langle {u^{\prime}}^2\right\rangle \left\langle {v^{\prime}}^2\right\rangle +\left\langle {u^{\prime}}^2\right\rangle \left\langle {w^{\prime}}^2\right\rangle +\left\langle {u^{\prime}}^2\right\rangle \left\langle {\xi}^2\right\rangle \right]\\ {}\kern2.5em +{A}_3\left[\left\langle {u^{\prime}}^3\right\rangle \left\langle {v^{\prime}}^1\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^3\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^2\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^1\right\rangle \left\langle {\xi}^2\right\rangle \right]\\ {}\kern2.5em +{A}_4\left[\left\langle {u^{\prime}}^3\right\rangle \left\langle {w^{\prime}}^1\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^2\right\rangle \left\langle {w^{\prime}}^1\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^3\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^1\right\rangle \left\langle {\xi}^2\right\rangle \right]\\ {}\kern2.5em +\frac{1}{2}{A}_5\left[\left\langle {u^{\prime}}^5\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^4\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^4\right\rangle +\left\langle {u^{\prime}}^1\right\rangle \left\langle {\xi}^4\right\rangle +2\left\langle {u^{\prime}}^3\right\rangle \left\langle {v^{\prime}}^2\right\rangle +2\left\langle {u^{\prime}}^3\right\rangle \left\langle {w^{\prime}}^2\right\rangle \right.\\ {}\kern6.5em +\left.\left.2\left\langle {u^{\prime}}^3\right\rangle \left\langle {\xi}^2\right\rangle +2\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^2\right\rangle \left\langle {w^{\prime}}^2\right\rangle +2\left\langle {u^{\prime}}^1\right\rangle \left\langle {v^{\prime}}^2\right\rangle \left\langle {\xi}^2\right\rangle +2\left\langle {u^{\prime}}^1\right\rangle \left\langle {w^{\prime}}^2\right\rangle \left\langle {\xi}^2\right\rangle \right]\right\}\\ {}\kern2em -\frac{1}{2}\tau \left\{\frac{\partial }{\partial {n}_1}\left\{{\rho}_l\left[{e}_l+\left({V^{\prime}}_l^2+{W^{\prime}}_l^2+\left(b-1\right){RT}_l\right){c}_l\right]\right.\right.\\ {}\kern8.5em +\left.{\rho}_r\left[{e}_r+\left({V^{\prime}}_r^2+{W^{\prime}}_r^2+\left(b-1\right){RT}_r\right){c}_r\right]\right\}\\ {}\kern5.5em +\frac{\partial }{\partial {n}_2}\left\{{\rho}_l{V}_l^{\prime}\left[{d}_l+\left({V^{\prime}}_l^2+{W^{\prime}}_l^2+\left(b+1\right){RT}_l\right){b}_l\right]\right.\\ {}\kern9em +\left.{\rho}_r{V}_r^{\prime}\left[{d}_r+\left({V^{\prime}}_r^2+{W^{\prime}}_r^2+\left(b+1\right){RT}_r\right){b}_r\right]\right\}\\ {}\kern5.5em +\frac{\partial }{\partial {n}_3}\left\{{\rho}_l{W}_l^{\prime}\left[{d}_l+\left({V^{\prime}}_l^2+{W^{\prime}}_l^2+\left(b+1\right){RT}_l\right){b}_l\right]\right.\\ {}\kern9em +\left.\left.{\rho}_r{W}_r^{\prime}\left[{d}_r+\left({V^{\prime}}_r^2+{W^{\prime}}_r^2+\left(b+1\right){RT}_r\right){b}_r\right]\right\}\right\}.\end{array}} $$
(55)

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 Ot), while that of the gas-kinetic BGK scheme is Ot2). 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.

3.3 Determination of collision time τ

Theoretically, the collision time τ in Eq. (1) is proportional to the physical viscosity

$$ \tau =\mu /p, $$
(56)

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:

$$ \tau =\frac{\mu }{p}+\frac{\mid {p}_l-{p}_r\mid }{p_l+{p}_r}\Delta t, $$
(57)

where Δt is the time step in the solution of Navier-Stokes equations, pl and pr 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.

3.4 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

$$ {g}^s=g\left[1+\left(1-\Pr \right)\mathbf{c}\cdot \mathbf{q}\left(\frac{c^2}{RT}-5\right)/\left(5 pRT\right)\right], $$
(58)

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

$$ \mathbf{q}=\frac{1}{2}\int \left(\mathbf{u}-\mathbf{U}\right)\left({\left(u-U\right)}^2+{\left(v-V\right)}^2+{\left(w-W\right)}^2+{\xi}^2\right) fd\Xi . $$
(59)

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].

$$ {F}_5^{new}={F}_5+\left(\frac{1}{\Pr }-1\right)\mathbf{q}\cdot {\mathbf{n}}_1, $$
(60)

where F5 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 F5, 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.

3.5 Computational sequence

In this section, the basic solution procedure of the current 3D GKFS is summarized as follows:

  1. (1)

    Firstly, we need to calculate the derivatives of conservative variables and reconstruct the initial conservative variables at two sides of cell interface.

  2. (2)

    Compute the unit vector in the normal direction n1 and in the tangential directions n2 and n3 of cell interface. Convert the velocities in the Cartesian coordinate system into the local coordinate system via Eq. (17).

  3. (3)

    Calculate the conservative variables at cell interface W by using Eqs. (33)–(37).

  4. (4)

    Calculate the vector (G1,  G2,  G3,  G4,  G5)T by using Eqs. (A.19)-(A.23) and further compute coefficients A1, A2, A3, A4, A5 by Eqs. (45)–(49).

  5. (5)

    Calculate the numerical flux \( {\mathbf{F}}_n^{\ast } \) by Eqs. (51)–(55).

  6. (6)

    Compute the heat flux q via Eq. (59), and make correction for energy flux by using Eq. (60).

  7. (7)

    Convert the numerical flux in the local coordinate system \( {\mathbf{F}}_n^{\ast } \) to the global Cartesian coordinate system Fn by using Eq. (22).

  8. (8)

    Once the fluxes at all cell interfaces are obtained, solve ordinary differential equation (Eq. (10)) by using time marching method. This step gives the conservative variables at cell centers at new time step.

  9. (9)

    Repeat steps (1) to (8) until convergence criterion is satisfied.

4 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 WL and WR in the reconstruction process. Specifically, WL and WR are computed by

$$ {\displaystyle \begin{array}{l}{\mathbf{W}}^L={\mathbf{W}}_c^L+{\Psi}_c^L\left({\mathbf{x}}_b-{\mathbf{x}}_c^L\right)\cdot \nabla {\mathbf{W}}_c^L,\\ {}{\mathbf{W}}^R={\mathbf{W}}_c^R+{\Psi}_c^R\left({\mathbf{x}}_b-{\mathbf{x}}_c^R\right)\cdot \nabla {\mathbf{W}}_c^R,\end{array}} $$
(61)

where \( {\mathbf{W}}_c^L \) and \( {\mathbf{W}}_c^R \) are the conservative flow variables at centers of the left and the right cells, respectively; \( \nabla {\mathbf{W}}_c^L \) and \( \nabla {\mathbf{W}}_c^R \) are their corresponding first-order derivatives. \( {\mathbf{x}}_c^L \), \( {\mathbf{x}}_c^R \) and xb are the coordinates of the left cell center, the right cell center and the midpoint of cell interface, respectively. \( {\Psi}_c^L \) and \( {\Psi}_c^R \) 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

$$ {\displaystyle \begin{array}{l}\rho \left(x,y,z\right)=1+0.2\sin \left(\pi \left(x+y+z\right)\right),\\ {}u\left(x,y,z\right)=1,v\left(x,y,z\right)=1,w\left(x,y,z\right)=1,p\left(x,y,z\right)=1.\end{array}} $$
(62)

The exact solutions under periodic boundary condition are

$$ {\displaystyle \begin{array}{l}\rho \left(x,y,z,t\right)=1+0.2\sin \left(\pi \left(x+y+z-3t\right)\right),\\ {}u\left(x,y,z,t\right)=1,v\left(x,y,z,t\right)=1,w\left(x,y,z,t\right)=1,p\left(x,y,z,t\right)=1.\end{array}} $$
(63)

Since this test case belongs to the inviscid flow, the collision time τ takes

$$ \tau =\varepsilon \Delta t+\frac{\mid {p}_l-{p}_r\mid }{p_l+{p}_r}\Delta t, $$
(64)

where ε = 0.01 is used. Numerical tests are conducted on the computational domain of [0, 2] × [0, 2] × [0, 2]. The uniform meshes with Δx = Δy = Δz = 2/N and N = 20, 40, 60, 80 are used. The L1 error of the density field at t = 2 is extracted and shown in Fig. 1. It can be seen that the GKFS is about the second order of accuracy in space.

Fig. 1
figure 1

L1 error of the density field for the advection of density perturbation

4.1 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

$$ {\displaystyle \begin{array}{l}{x}_i=0.5\left(1-\eta {\tan}^{-1}\left(\left(1-{\kappa}_i\right)\cdot \tan \left(1/\eta \right)\right)\right),\kern5em i\le \left(i\max +1\right)/2,\\ {}{x}_i=1.0-{x}_{i\max +1-i},\kern20.5em else.\end{array}} $$
(65)

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 stretching and is selected as 1.1 in this study. Similarly, the mesh point in the y- and z-directions 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].

Fig. 2
figure 2

Comparison of velocity profiles on the plane of z = 0.5 for 3D lid-driven cavity flow. Upper: Re = 100; Lower: Re = 400

Fig. 3
figure 3

Streamlines on three mid-planes for Re = 100 (left) and Re = 400 (right). a mid-plane of z = 0.5. b mid-plane of y = 0.5. c mid-plane of x = 0.5

4.2 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 = ρUD/μ, 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.

Fig. 4
figure 4

Partial view of computational mesh for flow past a sphere

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.

Table 1 Comparison of drag coefficient for flow past a stationary sphere

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 Ls increases with Reynolds number. Quantitative comparison between the present results of Ls 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 non-axisymmetric 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].

Fig. 5
figure 5

Streamlines at four different Reynolds numbers of 50, 100, 150 and 200 in the steady axisymmetric regime

Fig. 6
figure 6

Comparison of recirculation length Ls

Fig. 7
figure 7

Streamlines for flow past a stationary sphere at Re = 250 in the steady non-axisymmetric regime

4.3 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 free-stream Mach number is taken as M = 0.8395, the mean-chord based Reynolds number is chosen as Re = 11.72 × 106 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.

Fig. 8
figure 8

Pressure contours of flow around ONERA M6 wing

Fig. 9
figure 9

Comparison of pressure coefficient distributions at selected positions for ONERA M6 Wing

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 coefficient distribution at 65% span. It can be seen that the GKFS captures strong shock wave without any oscillation and the pressure coefficient distribution agrees well with the AUSM scheme [43].

Fig. 10
figure 10

Pressure contours and pressure coefficient distribution at 65% span for ONERA M6 Wing at M = 5

4.4 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 × 106 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 wing-body 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].

Fig. 11
figure 11

Pressure contours of DLR-F6 wing/body

Fig. 12
figure 12

Separation bubble on the intersection of wing and body (left: Vassberg et al. [46]; right: present)

Fig. 13
figure 13

Comparison of pressure coefficient distributions of DLR-F6 wing/body at different locations

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 × 106 and angle of attach α = 0. Table 2 shows the present results of force coefficients, including lift coefficient Cl, pressure drag coefficient Cd, p, friction drag coefficient Cd, f, total drag coefficient Cd and moment coefficient CM. 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].

Table 2 Comparison of force coefficients for DLR-F6 wing-body configuration

5 Conclusions

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.

Availability of data and materials

All data generated or analyzed during this study are included in this published article.

Abbreviations

2D:

Two-dimensional

3D:

Three-dimensional

KFVS:

Kinetic flux vector scheme

BGK:

Bhatnagar-Gross-Krook

FDS:

Flux difference splitting

CFD:

Computational fluid dynamics

GKFS:

Gas-kinetic flux solver

LU-SGS:

Lower-upper symmetric-Gauss-Seidel

References

  1. Sun Y, Shu C, Teo CJ et al (2015) Explicit formulations of gas-kinetic flux solver for simulation of incompressible and compressible viscous flows. J Comput Phys 300:492–519

    Article  MathSciNet  MATH  Google Scholar 

  2. Xu K (1998) Gas-kinetic schemes for unsteady compressible flow simulations. VKI for Fluid Dynamics Lecture Series

    Google Scholar 

  3. Su M, Xu K, Ghidaoui MS (1999) Low-speed flow simulation by the gas-kinetic scheme. J Comput Phys 150(1):17–39

    Article  MATH  Google Scholar 

  4. Xu K, Mao M, Tang L (2005) A multidimensional gas-kinetic BGK scheme for hypersonic viscous flow. J Comput Phys 203(2):405–421

    Article  MathSciNet  MATH  Google Scholar 

  5. Tian CT, Xu K, Chan KL et al (2007) A three-dimensional multidimensional gas-kinetic scheme for the Navier-Stokes equations under gravitational fields. J Comput Phys 226(2):2003–2027

    Article  MathSciNet  MATH  Google Scholar 

  6. Jiang J, Qian Y (2012) Implicit gas-kinetic BGK scheme with multigrid for 3D stationary transonic high-Reynolds number flows. Comput Fluids 66:21–28

    Article  MathSciNet  MATH  Google Scholar 

  7. Yang LM, Shu C, Wu J (2014) A simple distribution function-based gas-kinetic scheme for simulation of viscous incompressible and compressible flows. J Comput Phys 274:611–632

    Article  MathSciNet  MATH  Google Scholar 

  8. Li W, Kaneda M, Suga K (2014) An implicit gas kinetic BGK scheme for high temperature equilibrium gas flows on unstructured meshes. Comput Fluids 93:100–106

    Article  MathSciNet  MATH  Google Scholar 

  9. Li ZH, Zhang HX (2004) Study on gas kinetic unified algorithm for flows from rarefied transition to continuum. J Comput Phys 193(2):708–738

    Article  MATH  Google Scholar 

  10. Xu K, Huang JC (2010) A unified gas-kinetic scheme for continuum and rarefied flows. J Comput Phys 229(20):7747–7764

    Article  MathSciNet  MATH  Google Scholar 

  11. Guo Z, Xu K, Wang R (2013) Discrete unified gas kinetic scheme for all Knudsen number flows: low-speed isothermal case. Phys Rev E 88(3):033305

    Article  Google Scholar 

  12. Liu S, Yu P, Xu K et al (2014) Unified gas-kinetic scheme for diatomic molecular simulations in all flow regimes. J Comput Phys 259:96–113

    Article  MathSciNet  MATH  Google Scholar 

  13. Roe PL (1981) Approximate Riemann solvers, parameter vectors, and difference schemes. J Comput Phys 43(2):357–372

    Article  MathSciNet  MATH  Google Scholar 

  14. Engquist B, Osher S (1981) One-sided difference approximations for nonlinear conservation laws. Math Comput 36(154):321–351

    Article  MathSciNet  MATH  Google Scholar 

  15. Van Leer B (1997) Flux-vector splitting for the Euler equation. In: Hussaini MY, van Leer B, Van Rosendale J (eds) Upwind and high-resolution schemes. Springer, Heidelberg

    MATH  Google Scholar 

  16. Osher S, Chakravarthy S (1983) Upwind schemes and boundary conditions with applications to Euler equations in general geometries. J Comput Phys 50(3):447–481

    Article  MathSciNet  MATH  Google Scholar 

  17. Pullin DI (1980) Direct simulation methods for compressible inviscid ideal-gas flow. J Comput Phys 34(2):231–244

    Article  MATH  Google Scholar 

  18. Deshpande SM (1986) A second-order accurate kinetic-theory-based method for inviscid compressible flows. NASA STI/Recon Technical Report N 87

    Google Scholar 

  19. Perthame B (1992) Second-order Boltzmann schemes for compressible Euler equations in one and two space dimensions. SIAM J Numer Anal 29(1):1–19

    Article  MathSciNet  MATH  Google Scholar 

  20. Mandal JC, Deshpande SM (1994) Kinetic flux vector splitting for Euler equations. Comput Fluids 23(2):447–478

    Article  MathSciNet  MATH  Google Scholar 

  21. Chou SY, Baganoff D (1997) Kinetic flux-vector splitting for the Navier-Stokes equations. J Comput Phys 130(2):217–230

    Article  MATH  Google Scholar 

  22. Prendergast KH, Xu K (1993) Numerical hydrodynamics from gas-kinetic theory. J Comput Phys 109(1):53–66

    Article  MathSciNet  MATH  Google Scholar 

  23. Chae D, Kim C, Rho OH (2000) Development of an improved gas-kinetic BGK scheme for inviscid and viscous flows. J Comput Phys 158(1):1–27

    Article  MATH  Google Scholar 

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

    Article  MathSciNet  MATH  Google Scholar 

  25. Xu K, Sun Q, Yu P (2010) Valid physical processes from numerical discontinuities in computational fluid dynamics. Int J Hypersonics 1(3):157–172

    Article  Google Scholar 

  26. Bhatnagar PL, Gross EP, Krook M (1954) A model for collision processes in gases. I Small amplitude processes in charged and neutral one-component systems. Phys Rev 94(3):511–525

    Article  MATH  Google Scholar 

  27. Shakhov EM (1968) Generalization of the Krook kinetic relaxation equation. Fluid Dyn 3(5):95–96

    Article  MathSciNet  Google Scholar 

  28. Yoon S, Jameson A (1988) Lower-upper symmetric-Gauss-Seidel method for the Euler and Navier-Stokes equations. AIAA J 26(9):1025–1026

    Article  Google Scholar 

  29. Venkatakrishnan V (1995) Convergence to steady state solutions of the Euler equations on unstructured grids with limiters. J Comput Phys 118(1):120–130

    Article  MATH  Google Scholar 

  30. Pan L, Xu K (2020) High-order gas-kinetic scheme with three-dimensional WENO reconstruction for the Euler and Navier-Stokes solutions. Comput Fluids 198:104401

    Article  MathSciNet  MATH  Google Scholar 

  31. Shu C, Niu XD, Chew YT (2003) Taylor series expansion and least squares-based lattice Boltzmann method: three-dimensional formulation and its applications. Int J Mod Phys C 14(07):925–944

    Article  MATH  Google Scholar 

  32. Wu J, Shu C (2010) An improved immersed boundary-lattice Boltzmann method for simulating three-dimensional incompressible flows. J Comput Phys 229(13):5022–5042

    Article  MATH  Google Scholar 

  33. Wu J, Shu C (2009) Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications. J Comput Phys 228(6):1963–1979

    Article  MATH  Google Scholar 

  34. Wang Y, Shu C, Teo CJ et al (2015) An immersed boundary-lattice Boltzmann flux solver and its applications to fluid-structure interaction problems. J Fluids Struct 54:440–465

    Article  Google Scholar 

  35. Johnson TA, Patel VC (1999) Flow past a sphere up to a Reynolds number of 300. J Fluid Mech 378:19–70

    Article  Google Scholar 

  36. Kim J, Kim D, Choi H (2001) An immersed-boundary finite-volume method for simulations of flow in complex geometries. J Comput Phys 171(1):132–150

    Article  MathSciNet  MATH  Google Scholar 

  37. Wang XY, Yeo KS, Chew CS et al (2008) A SVD-GFD scheme for computing 3D incompressible viscous fluid flows. Comput Fluids 37(6):733–746

    Article  MathSciNet  MATH  Google Scholar 

  38. Gilmanov A, Sotiropoulos F, Balaras E (2003) A general reconstruction algorithm for simulating flows with complex 3D immersed boundaries on Cartesian grids. J Comput Phys 191(2):660–669

    Article  MATH  Google Scholar 

  39. Slater JW (2002) ONERA M6 Wing. http://www.grc.nasa.gov/WWW/wind/valid/m6wing/m6wing01/m6wing01.html. Accessed 30 Aug 2002

    Google Scholar 

  40. Spalart P, Allmaras S (1992) A one-equation turbulence model for aerodynamic flows. AIAA paper 92-0439

    Google Scholar 

  41. Yang LM, Shu C, Wu J (2015) A three-dimensional explicit sphere function-based gas-kinetic flux solver for simulation of inviscid compressible flows. J Comput Phys 295:322–339

    Article  MathSciNet  MATH  Google Scholar 

  42. Schmitt V, Charpin F (1979) Pressure distributions on the ONERA-M6-wing at transonic Mach numbers, experimental data base for computer program assessment. Report of the Fluid Dynamics Panel Working Group 04, AGARD AR 138, May 1979

    Google Scholar 

  43. Liou MS (2006) A sequel to AUSM, part II: AUSM+-up for all speeds. J Comput Phys 214:137–170

    Article  MathSciNet  MATH  Google Scholar 

  44. Frink NT (2006) 3rd AIAA CFD Drag Prediction Workshop. https://aiaa-dpw.larc.nasa.gov/Workshop3/workshop3.html. Accessed 29 Nov 2006

    Google Scholar 

  45. Morrison J (2006) https://dpw.larc.nasa.gov/DPW3/multiblock_Boeing_CGNS/F6wb_struc_boeing. Accessed 18 Jan 2006

  46. Vassberg J, Tinoco E, Mani M et al (2007) Summary of the third AIAA CFD drag prediction workshop. Paper presented at the 45th AIAA aerospace sciences meeting and exhibit, 1 Jan 2007

    Book  Google Scholar 

  47. Rumsey CL, Rivers SM, Morrison JH (2005) Study of CFD variation on transport configurations from the second drag-prediction workshop. Comput Fluids 34(7):785–816

    Article  MATH  Google Scholar 

  48. Yang LM, Shu C, Wu J (2014) A hybrid lattice Boltzmann flux solver for simulation of 3D compressible viscous flows. Paper presented at the eighth international conference on computational fluid dynamics, Chengdu, China, 24-29 July 2014

    Google Scholar 

Download references

Acknowledgments

This work is supported by National Natural Science Foundation of China (Grant Nos. 11772157 and 11832012).

Funding

National Natural Science Foundation of China (Grant Nos. 11772157 and 11832012).

Author information

Authors and Affiliations

Authors

Contributions

The contribution of the authors to this work is equivalent. All authors read and approved the final manuscript.

Corresponding author

Correspondence to L. M. Yang.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix

Appendix

1.1 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))

$$ g=\rho {\left(\frac{\lambda }{\pi}\right)}^{\frac{K+3}{2}}{e}^{-\lambda \left({\left(u-U\right)}^2+{\left(v-V\right)}^2+{\left(w-W\right)}^2+{\xi}^2\right)}. $$

Following the idea of [2], the notation for the moments of g is defined as

$$ \rho \left\langle \cdot \right\rangle =\int \left(\cdot \right) gdudvdwd\xi . $$
(A.1)

Then the general moment formualtion becomes

$$ \left\langle {u}^n{v}^m{w}^l{\xi}^p\right\rangle =\left\langle {u}^n\right\rangle \left\langle {v}^m\right\rangle \left\langle {w}^l\right\rangle \left\langle {\xi}^p\right\rangle . $$
(A.2)

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 un and ξn are

$$ \left\langle {u^{\prime}}^0\right\rangle =1, $$
(A.3)
$$ \left\langle {u^{\prime}}^1\right\rangle ={U}^{\prime }, $$
(A.4)
$$ \left\langle {u^{\prime}}^{(n)}\right\rangle ={U}^{\prime}\left\langle {u^{\prime}}^{\left(n-1\right)}\right\rangle +\frac{n-1}{2\lambda}\left\langle {u^{\prime}}^{\left(n\hbox{-} 2\right)}\right\rangle, $$
(A.5)

and

$$ \left\langle {\xi}^0\right\rangle =1, $$
(A.6)
$$ \left\langle {\xi}^2\right\rangle =\frac{K}{2\lambda }, $$
(A.7)
$$ \left\langle {\xi}^4\right\rangle =\frac{K^2+2K}{4{\lambda}^2}. $$
(A.8)

When the moments for un are calculated in the half space, the exponential function and the complementary error function appear in the formulation. Take notation of the integral from 0 to +∞ as 〈>0 and integral from −∞ to 0 as 〈<0, the moments become

$$ {a}_l={\left\langle {u^{\prime}}^0\right\rangle}_{>0}=\frac{1}{2}\mathit{\operatorname{erfc}}\left(-\sqrt{\lambda_l}{U}_l^{\prime}\right), $$
(A.9)
$$ {b}_l={\left\langle {u^{\prime}}^1\right\rangle}_{>0}={U}_l^{\prime }{a}_l+\frac{1}{2}\frac{e^{-{\lambda}_l{U_l^{\prime}}^2}}{\sqrt{{\pi \lambda}_l}}, $$
(A.10)
$$ {c}_l={\left\langle {u^{\prime}}^2\right\rangle}_{>0}={U}_l^{\prime }{b}_l+\frac{1}{2{\lambda}_l}{a}_l, $$
(A.11)
$$ {d}_l={\left\langle {u^{\prime}}^3\right\rangle}_{>0}={U}_l^{\prime }{c}_l+\frac{1}{\lambda_l}{b}_l, $$
(A.12)
$$ {e}_l={\left\langle {u^{\prime}}^4\right\rangle}_{>0}={U}_l^{\prime }{d}_l+\frac{3}{2{\lambda}_l}{c}_l, $$
(A.13)

and

$$ {a}_r={\left\langle {u^{\prime}}^0\right\rangle}_{<0}=\frac{1}{2}\mathit{\operatorname{erfc}}\left(\sqrt{\lambda_r}{U}_r^{\prime}\right), $$
(A.14)
$$ {b}_r={\left\langle {u^{\prime}}^1\right\rangle}_{<0}={U}_r^{\prime }{a}_r-\frac{1}{2}\frac{e^{-{\lambda}_r{U_r^{\prime}}^2}}{\sqrt{{\pi \lambda}_r}}, $$
(A.15)
$$ {c}_r={\left\langle {u^{\prime}}^2\right\rangle}_{<0}={U}_r^{\prime }{b}_r+\frac{1}{2{\lambda}_r}{a}_r, $$
(A.16)
$$ {d}_r={\left\langle {u^{\prime}}^3\right\rangle}_{<0}={U}_r^{\prime }{c}_r+\frac{1}{\lambda_r}{b}_r, $$
(A.17)
$$ {e}_r={\left\langle {u^{\prime}}^4\right\rangle}_{<0}={U}_r^{\prime }{d}_r+\frac{3}{2{\lambda}_r}{c}_r. $$
(A.18)

In order to compute the coefficients A1 to A5 via Eqs. (45)–(49), the values of G1 to G5 in Eq. (44) should be calculated in advance. The explicit formulations of G1 to G5 can be given as

$$ {G}_1=\frac{\partial \left({\rho}_l{b}_l+{\rho}_r{b}_r\right)}{\partial {n}_1}+\frac{\partial \left({\rho}_l{V}_l^{\prime }{a}_l+{\rho}_r{V}_r^{\prime }{a}_r\right)}{\partial {n}_2}+\frac{\partial \left({\rho}_l{W}_l^{\prime }{a}_l+{\rho}_r{W}_r^{\prime }{a}_r\right)}{\partial {n}_3}, $$
(A.19)
$$ {G}_2=\frac{\partial \left({\rho}_l{c}_l+{\rho}_r{c}_r\right)}{\partial {n}_1}+\frac{\partial \left({\rho}_l{V}_l^{\prime }{b}_l+{\rho}_r{V}_r^{\prime }{b}_r\right)}{\partial {n}_2}+\frac{\partial \left({\rho}_l{W}_l^{\prime }{b}_l+{\rho}_r{W}_r^{\prime }{b}_r\right)}{\partial {n}_3}, $$
(A.20)
$$ {\displaystyle \begin{array}{l}{G}_3=\frac{\partial \left({\rho}_l{V_l}^{\prime }{b}_l+{\rho}_r{V_r}^{\prime }{b}_r\right)}{\partial {n}_1}+\frac{\partial \left[\left({\rho}_l{V^{\prime}}_l^2+{p}_l\right){a}_l+\left({\rho}_r{V^{\prime}}_r^2+{p}_r\right){a}_r\right]}{\partial {n}_2}\\ {}\kern3.5em +\frac{\partial \left({\rho}_l{V_l}^{\prime }{W_l}^{\prime }{a}_l+{\rho}_r{V_r}^{\prime }{W_r}^{\prime }{a}_r\right)}{\partial {n}_3},\end{array}} $$
(A.21)
$$ {\displaystyle \begin{array}{l}{G}_4=\frac{\partial \left({\rho}_l{W_l}^{\prime }{b}_l+{\rho}_r{W_r}^{\prime }{b}_r\right)}{\partial {n}_1}+\frac{\partial \left({\rho}_l{V_l}^{\prime }{W_l}^{\prime }{a}_l+{\rho}_r{V_r}^{\prime }{W_r}^{\prime }{a}_r\right)}{\partial {n}_2}\\ {}\kern3.5em +\frac{\partial \left[\left({\rho}_l{W^{\prime}}_l^2+{p}_l\right){a}_l+\left({\rho}_r{W^{\prime}}_r^2+{p}_r\right){a}_r\right]}{\partial {n}_3},\end{array}} $$
(A.22)
$$ {\displaystyle \begin{array}{l}{G}_5=\frac{\partial }{\partial {n}_1}\Big\{{\rho}_l\left[{d}_l+\left({V^{\prime}}_l^2+{W^{\prime}}_l^2+\left(b-1\right){RT}_l\right){b}_l\right]\operatorname{}\\ {}\kern5.5em +\operatorname{}{\rho}_r\left[{d}_r+\left({V^{\prime}}_r^2+{W^{\prime}}_r^2+\left(b-1\right){RT}_r\right){b}_r\right]\Big\}\\ {}\kern2.5em +\frac{\partial }{\partial {n}_2}\Big\{{\rho}_l{V}_l^{\prime}\left[{c}_l+\left({V^{\prime}}_l^2+{W^{\prime}}_l^2+\left(b+1\right){RT}_l\right){a}_l\right]\operatorname{}\\ {}\kern5.5em +\operatorname{}{\rho}_r{V}_r^{\prime}\left[{c}_r+\left({V^{\prime}}_r^2+{W^{\prime}}_r^2+\left(b+1\right){RT}_r\right){a}_r\right]\Big\}\\ {}\kern2.5em +\frac{\partial }{\partial {n}_3}\Big\{{\rho}_l{W}_l^{\prime}\left[{c}_l+\left({V^{\prime}}_l^2+{W^{\prime}}_l^2+\left(b+1\right){RT}_l\right){a}_l\right]\operatorname{}\\ {}\kern5.5em +\operatorname{}{\rho}_r{W}_r^{\prime}\left[{c}_r+\left({V^{\prime}}_r^2+{W^{\prime}}_r^2+\left(b+1\right){RT}_r\right){a}_r\right]\Big\}.\end{array}} $$
(A.23)

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Sun, Y., Yang, L.M., Shu, C. et al. A three-dimensional gas-kinetic flux solver for simulation of viscous flows with explicit formulations of conservative variables and numerical flux. Adv. Aerodyn. 2, 13 (2020). https://doi.org/10.1186/s42774-020-00039-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s42774-020-00039-6

Keywords