An efficient high-order gas-kinetic scheme with hybrid WENO-AO method for the Euler and Navier-Stokes solutions

The two-dimensional Boltzmann-BGK equation is written as

where $\mathbf{u}=\left(u,v\right)$ is the particle velocity, $\tau$ is the collision time. $f$ is the gas distribution function and $g$ is the equilibrium state [28].

The equilibrium state $g$ is a Maxwellian distribution,

where $\rho$ is the density and $\lambda=m/2kT$, where $m$ is the molecular mass, $k$ is the Boltzmann constant, and $T$ is the temperature. $U$ and $V$ are the macroscopic velocities in the x-direction and y-direction. For a 2D flow, the internal variable $\xi$ include the random particle motion in the z-direction, and the total number of degrees of freedom $K=\left(5-3\gamma\right)/\left(\gamma-1\right)+1$. $\gamma$ is the specific heat ratio. In the equilibrium state, the internal variable ${{\xi}^{2}}=\xi_{1}^{2}+\xi_{2}^{2}+\cdots+\xi_{K}^{2}$. The relation between the conservative variables $\left(\rho,\rho U,\rho V,\rho E\right)$ and the distribution function $f$ is

where ${{\psi}_{\alpha}}={{\left({{\psi}_{1}},{{\psi}_{2}},{{\psi}_{3}},{{\psi}_{4}}\right)}^{T}}={{\left(1,u,v,\frac{1}{2}\left({{u}^{2}}+{{v}^{2}}+{{\xi}^{2}}\right)\right)}^{T}}$, $d\Xi=dudvd{{\xi}_{1}}\ldots d{{\xi}_{K}}$. The collision term satisfies the following compatibility condition

Based on the Chapman-Enskog expansion for Boltzmann-BGK equation [29, 30], the gas distribution function in the continuum regime can be expanded as

where ${{D}_{u}}=\partial/\partial t+u\cdot\nabla$. By truncating on different orders of $\tau$, the corresponding macroscopic equations can be derived [31]. For the Euler equations, the zeroth order truncation is taken, i.e. $f=g$. For the Navier-Stokes equations, when the first order truncation is used, the distribution function is

The solution $f$ of the Boltzmann-BGK equation at a cell interface ${{x}_{j+1/2}}$ is

where ${x}^{\prime}={{x}_{i+1/2}}-u\left(t-{t}^{\prime}\right)$ and ${y}^{\prime}=y-v\left(t-{t}^{\prime}\right)$ are the particle trajectory, ${{x}_{i+1/2}}=0$ is the location of the cell interface. ${{f}_{0}}$ is the initial gas distribution function $f$ at the beginning of each time step $\left(t=0\right)$. The integral solution simulates the physical process from the free transport of particles in the kinetic scale physics ${{f}_{0}}$ to the evolution of hydrodynamic flow in the integral of $g$ term. The flow behavior at the cell interface depends on the ratio of the time step to the local particle collision time $\Delta t/\tau$.

The initial gas distribution function ${{f}_{0}}$ with considering the first order non-equilibrium part is expressed as follows:

The relation between the terms of the ${{f}_{0}}$ and macroscopic variables ${Q}$ is

where $k=l,r$ and $\left\langle\ldots\right\rangle$ are the moments of the equilibrium $g$ and defined by

The coefficients can be determined by the spatial derivatives of macroscopic flow variables and the compatibility condition as follows

where $a\equiv\left(\partial g/\partial x\right)/g={{g}_{x}}/g,A\equiv\left(\partial g/\partial t\right)/g={{g}_{t}}/g$, After obtaining all the required terms, substitute these terms into Eq. (6) and solve it can obtain the gas distribution function on the interface [2].

where $H(u)$ is the Heaviside function. Derivations related to gas kinetic scheme can be found in  [2].

Finally, the gas-kinetic numerical fluxes in the x-direction on the cell interface can be computed as

By integrating the above equation in the time interval $\left[{{t}_{n}},\ {{t}_{n}}+\Delta t\right]$, the total mass, momentum, and energy transport can be written as

More details can be found in  [2].

There are many ways to solve ODE Eq. (1). Here we define

The two-stage fourth-order time marching scheme is used to solve the initial value problem Eq. (1) and is written as

where $\partial\mathcal{L}\left(\mathbf{W}\right)/\partial t$ is the time derivative of the spatial operator. If the time derivatives of ${{\mathcal{L}}^{\left(n\right)}}$ up to $\left(n-1\right)$ can be given, a $n$th-order time marching scheme can be constructed directly. Usually only a few low-order derivative is used to solve nonlinear system problems, such as $\mathcal{L}$ for approximating the Riemann solver. For hyperbolic conservation law, the two-stage fourth-order temporal discretization is derived in  [10] and ${{\mathcal{L}}^{\left(1\right)}}$ for the generalized Riemann problem (GRP) solver. In the establishment of high-order gas-kinetic schemes, ${{\mathcal{L}}^{\left(1\right)}}$ and ${{\mathcal{L}}^{\left(2\right)}}$ are often used by HGKS for 2nd-order GKS flux functions and 3rd-order GKS flux functions respectively [32, 33, 34]. However, higher-order derivatives are prohibited due to the complexity of the HGKS flux function. Another multi-stage multi-derivative method (MSMD) was proposed with HGKS and continues to be used. Similar to the Runge-Kutta (RK) time marching scheme, the MSMD introduce middle stages. The update at ${{t}^{n+1}}$ becomes a linear combination of $\mathcal{L}$ and their derivatives in multiple stages. Because of the use of $\mathcal{L}$ and their derivatives, MSMD can achieve the required time accuracy with fewer middle stages than RK method. This paper mainly focuses on the two-stage fourth-order temporal discretization scheme (S2O4), and more HGKS with MSMD method is detailed and analyzed in  [12].

Eq. (10) provides a time-dependent gas distribution function. And the flux for the macroscopic flow variables can be calculated by Eq. (12). The flux of the two-stage fourth-order GKS in the time interval $\left[{{t}_{n}},\ {{t}_{n}}+\Delta t\right]$ needs the ${{F}_{i\pm 1/2}}\left(\mathbf{W}\right)$ and ${{\partial}_{t}}{{F}_{i\pm 1/2}}\left(\mathbf{W}\right)$ at both ${{t}_{n}}$ and ${{t}_{*}}={{t}_{n}}+\Delta t/2$. The followings mainly introduce the two-stage algorithm.
(1) With the reconstruction at ${{t}_{n}}$, the x-direction flux ${{\mathbb{F}}_{i+1/2,j}}\left({{\mathbf{W}}^{n}},\Delta t\right)$, ${{\mathbb{F}}_{i+1/2,j}}\left({{\mathbf{W}}^{n}},\Delta t/2\right)$ and y-direction flux ${{\mathbb{G}}_{i,j+1/2}}\left({{\mathbf{W}}^{n}},\Delta t\right)$ and ${{\mathbb{G}}_{i,j+1/2}}\left({{\mathbf{W}}^{n}},\Delta t/2\right)$ in the time interval $\left[{{t}_{n}},\ {{t}_{n}}+\frac{\Delta t}{2}\right]$ can be evaluated by Eq. (12).
(2) In the time interval $\left[{{t}_{n}},\ {{t}_{n}}+\Delta t\right]$, the flux is expanded as the linear form as

The terms $F_{{}_{i+1/2,j}}\left({{\mathbf{W}}^{n}},{{t}_{n}}\right)$ and $\partial tF_{{}_{i+1/2,j}}\left({{\mathbf{W}}^{n}},{{t}_{n}}\right)$ can be determined as

The terms $F_{{}_{i+1/2,j}}\left({{\mathbf{W}}^{n}},{{t}_{n}}\right)$ and $\partial tF_{{}_{i+1/2,j}}\left({{\mathbf{W}}^{n}},{{t}_{n}}\right)$ can be computed by solving this linear system as

Similar to $F_{{}_{i+1/2,j}}\left({{\mathbf{W}}^{n}},{{t}_{n}}\right)$ and $\partial tF_{{}_{i+1/2,j}}\left({{\mathbf{W}}^{n}},{{t}_{n}}\right)$, the $G_{{}_{i,j+1/2}}\left({{\mathbf{W}}^{n}},{{t}_{n}}\right)$ and $\partial tG_{{}_{i,j+1/2}}\left({{\mathbf{W}}^{n}},{{t}_{n}}\right)$ can be computed from ${{\mathbb{G}}_{i,j+1/2}}\left({{\mathbf{W}}^{n}},\Delta t\right)$ and ${{\mathbb{G}}_{i,j+1/2}}\left({{\mathbf{W}}^{n}},\Delta t/2\right)$.
(3) Update $\mathbf{W}_{ij}^{*}$ at ${{t}_{*}}={{t}_{n}}+\Delta t/2$ by

Again through the step (1) and step (2) to compute the middle stage value, we can get $\partial tF_{{}_{i+1/2,j}}\left({{\mathbf{W}}^{*}},{{t}_{*}}\right)$ and $\partial tG_{{}_{i,j+1/2}}\left({{\mathbf{W}}^{*}},{{t}_{*}}\right)$ in the time interval $\left[{{t}_{*}},\ {{t}_{*}}+\Delta t\right]$.
(4) The numerical fluxes $\mathcal{F}_{i+1/2,j}^{n}$ and $\mathcal{G}_{i+1/2,j}^{n}$ can be computed by

Update $\mathbf{W}_{ij}^{n+1}$ by

In summary, for most of the tests in the third section, the flux function from the distribution function Eq. (10) is typically used. As analyzed in [11], in smooth region the Euler equation with leading error $\mathcal{O}\left({{\left(\Delta x\right)}^{5}},{{\left(\Delta t\right)}^{4}}\right)$ is solved in the S2O4 scheme, and the NS equations with the error $\mathcal{O}\left({{\left(\Delta x\right)}^{5}},\tau\Delta t\right)$. Normally, the condition $\tau\ll\Delta t$ is satisfied for many flow computations.

In this following, we present the main idea about the hybrid WENO scheme. We describe how to determine the point-wise values and derivatives required on the interface in the present hybrid WENO scheme. Next, we extend the HKGS with hybrid WENO method to the two-dimensional case (the process for extending to the three-dimensional case is straightforward). Finally, we provide a detailed analysis of why the improved HGKS can reduce computational costs.