A Hybrid Discontinuous Galerkin Scheme for Multi-scale Kinetic Equations

TL;DR

This paper introduces a hybrid discontinuous Galerkin method for multi-scale kinetic equations, enabling model selection and efficient numerical flux computation, demonstrated through various numerical simulations.

Contribution

It presents a novel multi-dimensional hybrid DG scheme utilizing moment realizability matrices for model selection and flux approximation in multi-scale kinetic problems.

Findings

Effective model selection indicator developed

Numerical flux constructed for asymptotic fluid limit

Numerical simulations validate method performance

Abstract

We develop a multi-dimensional hybrid discontinuous Galerkin method for multi-scale kinetic equations. This method is based on moment realizability matrices, a concept introduced by D. Levermore, W. Morokoff and B. Nadiga for one dimensional problem. The main issue addressed in this paper is to provide a simple indicator to select the most appropriate model and to apply a compact numerical scheme to reduce the interface region between different models. We also construct a numerical flux for the fluid model obtained as the asymptotic limit of the flux of the kinetic equation. Finally we perform several numerical simulations for time evolution and stationary problems.

Equations

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial f}{\partial t}+{\bf v}\cdot\nabla_{\bf x}f\,=\,\frac{1}{\varepsilon}\,\mathcal{Q}(f),\\ \\ f(0,{\bf x},{\bf v})\,=\,f_{0}({\bf x},{\bf v}),\end{array}\right.$$

$$\mathcal{Q}(f)(\mathbf{v})={\int }_{{\mathbb{R}}^3}{\int }_{{\mathbb{S}}^2}B(|\mathbf{v}-{\mathbf{v}}_{\star }|,\cos \theta )(f_{\star }^{\prime }{f}^{\prime }-f_{\star }f)d\sigma d{\mathbf{v}}_{\star },$$

$${\mathbf{v}}^{\prime }=\frac{\mathbf{v}+{\mathbf{v}}_{\ast }}{2}+\frac{|\mathbf{v}-{\mathbf{v}}_{\ast }|}{2}\sigma ,{\mathbf{v}}_{\ast }^{\prime }=\frac{\mathbf{v}+{\mathbf{v}}^{\ast }}{2}-\frac{|\mathbf{v}-{\mathbf{v}}_{\ast }|}{2}\sigma .$$

$${\mathcal{Q}}(f)({\bf v})\,=\,\nu(\rho,T)\big{(}\mathcal{M}[f]-f),$$

$${\mathcal{Q}}(f)({\bf v})\,=\,\nu(\rho,T)\big{(}\mathcal{G}[f]-f),$$

$$\mathcal{G}[f]=\frac{\rho }{\sqrt{\det (2\pi \mathcal{T})}}\exp \left(-\frac{(\mathbf{v}-\mathbf{u}){\mathcal{T}}^{-1}(\mathbf{v}-\mathbf{u})}{2}\right),$$

$$\mathcal{T}:=(1-\beta )T\mathbf{I}+\beta \Theta ,\rho \Theta ={\int }_{{\mathbb{R}}^3}(\mathbf{v}-\mathbf{u})\otimes (\mathbf{v}-\mathbf{u})fd\mathbf{v},$$

$$0<\mathrm{Pr}=\frac{1}{1-\beta }\le +\infty .$$

$${\int }_{{\mathbb{R}}^3}\mathcal{Q}(f)(\mathbf{v})m(\mathbf{v})d\mathbf{v}=0_{{\mathbb{R}}^5}^T,$$

$$\mathcal{M}[f]\,:=\,\frac{\rho}{(2\pi T)^{3/2}}\exp\Big{(}-\frac{|{\bf v}-{\bf u}|^{2}}{2T}\Big{)}.$$

$$\rho ={\int }_{{\mathbb{R}}^3}f(\mathbf{v})d\mathbf{v},\mathbf{u}=\frac{1}{\rho }{\int }_{{\mathbb{R}}^3}\mathbf{v}f(\mathbf{v})d\mathbf{v},T=\frac{1}{3\rho }{\int }_{{\mathbb{R}}^3}|\mathbf{v}-\mathbf{u}{|}^{2}f(\mathbf{v})d\mathbf{v}.$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}\rho+\nabla_{\bf x}\cdot(\rho{\bf u})\,=\,0,\\ \\ \displaystyle\partial_{t}(\rho{\bf u})+\nabla_{\bf x}\cdot\left(\int_{\mathbb{R}^{3}}{\bf v}\otimes{\bf v}\,f({\bf v})\,d{\bf v}\right)\,=\,{\bf 0}_{\mathbb{R}^{3}},\\ \\ \displaystyle\partial_{t}E+\nabla_{\bf x}\cdot\left(\int_{\mathbb{R}^{3}}\frac{1}{2}|{\bf v}|^{2}\,{\bf v}\,f({\bf v})\,d{\bf v}\right)\,=\,0.\end{array}\right.$$

$$\mathbf{U}:=(\rho ,\rho \mathbf{u},E{)}^{\text{T}}={\int }_{{\mathbb{R}}^3}m(\mathbf{v})f(\mathbf{v})d\mathbf{v},$$

$$f^{\varepsilon}({\bf v})\,:=\,\mathcal{M}\big{[}1+\varepsilon\,g^{(1)}+\varepsilon^{2}\,g^{(2)}+\cdots\big{]},$$

$${\int }_{{\mathbb{R}}^3}{g}^{(i)}(\mathbf{v})m(\mathbf{v})d\mathbf{v}=0_{{\mathbb{R}}^5}^T,i=1,2,\cdots .$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}\rho+\nabla_{\bf x}\cdot(\rho{\bf u})\,=\,0,\\ \\ \displaystyle\partial_{t}(\rho{\bf u})+\nabla_{\bf x}\cdot\left(\rho\,{\bf u}\otimes{\bf u}+\rho\,T({\bf I}+\bar{{\bf A}}^{\varepsilon})\right)\,=\,{\bf 0}_{\mathbb{R}^{3}},\\ \\ \displaystyle\partial_{t}E+\nabla_{\bf x}\cdot\Bigg{(}\frac{1}{2}\rho\,|{\bf u}|^{2}\,{\bf u}+\rho\,T\,\Big{(}\frac{3+2}{2}{\bf I}+\bar{\bf A}^{\varepsilon}\Big{)}{\bf u}+\rho\,T^{3/2}\bar{\bf B}^{\varepsilon}\Bigg{)}\,=\,0,\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\bar{\bf A}^{\varepsilon}\,:=\,\frac{1}{\rho}\int_{\mathbb{R}^{3}}{\bf A}({\bf V})f^{\varepsilon}({\bf v})d{\bf v},\\ \\ \displaystyle\bar{\bf B}^{\varepsilon}\,:=\,\frac{1}{\rho}\int_{\mathbb{R}^{3}}{\bf B}({\bf V})f^{\varepsilon}({\bf v})d{\bf v},\end{array}\right.$$

$${\bf A}({\bf V})={\bf V}\otimes{\bf V}-\frac{|{\bf V}|^{2}}{3}{\bf I},\quad{\bf B}({\bf V})=\frac{1}{2}\big{[}|{\bf V}|^{2}-(3+2)\big{]}{\bf V}$$

$$\mathbf{V}:=\frac{\mathbf{v}-\mathbf{u}}{\sqrt{T}}.$$

$$\left\{\begin{array}[]{l}\displaystyle\bar{\bf A}^{\varepsilon}_{Euler}\,:=\,\frac{1}{\rho}\int_{\mathbb{R}^{3}}{\bf A}({\bf V})\mathcal{M}({\bf v})d{\bf v}\,=\,{\bf 0}^{T}_{\mathbb{M}^{3}},\\ \\ \displaystyle\bar{\bf B}^{\varepsilon}_{Euler}\,:=\,\frac{1}{\rho}\int_{\mathbb{R}^{3}}{\bf B}({\bf V})\mathcal{M}({\bf v})d{\bf v}\,=\,{\bf 0}^{T}_{\mathbb{R}^{3}},\end{array}\right.$$

$$\left\{\begin{array}[]{l}\partial_{t}\rho+\nabla_{\bf x}\cdot(\rho{\bf u})\,=\,0,\\ \\ \partial_{t}(\rho{\bf u})+\nabla_{\bf x}\cdot\left(\rho\,{\bf u}\otimes{\bf u}+\rho\,T\,{\bf I}\right)\,=\,{\bf 0}_{\mathbb{R}^{3}},\\ \\ \partial_{t}E+\nabla_{\bf x}\cdot\Big{(}{\bf u}\big{(}E+\rho\,T\big{)}\Big{)}\,=\,0.\end{array}\right.$$

$${\partial }_t\mathcal{M}+\mathbf{v}\cdot {\nabla }_{\mathbf{x}}\mathcal{M}={\mathcal{L}}_{\mathcal{M}}{g}^{(1)}+\mathcal{O}(\epsilon ),$$

$$\partial_{t}\mathcal{M}+{\bf v}\cdot\nabla_{\bf x}\mathcal{M}\,=\,\mathcal{M}\Big{[}\partial_{t}\rho+{\bf v}\cdot\nabla_{\bf x}\rho+\frac{1}{\sqrt{T}}\left({\bf V}\cdot\partial_{t}{\bf u}+{\bf V}\otimes{\bf v}:\nabla_{\bf x}{\bf u}\right)+\frac{1}{2T}(|{\bf V}|^{2}-3)(\partial_{t}T+{\bf v}\cdot\nabla_{\bf x}T)\Big{]}.$$

$${\mathcal{L}}_{\mathcal{M}}{g}^{(1)}=\mathcal{M}\left[\mathbf{A}(\mathbf{V}):\mathbf{D}(\mathbf{u})+2\mathbf{B}(\mathbf{V})\cdot \frac{{\nabla }_{\mathbf{x}}T}{\sqrt{T}}\right],$$

$$\mathbf{D}(\mathbf{u}):={\nabla }_{\mathbf{x}}\mathbf{u}+({\nabla }_{\mathbf{x}}\mathbf{u}{)}^T-\frac{2}{3}({\nabla }_{\mathbf{x}}\cdot \mathbf{u})\mathbf{I}.$$

$${\rm ker}\,\mathcal{L}_{\mathcal{M}}\,=\,{\rm Span}\left\{\frac{1}{\rho},\frac{{\bf V}}{\rho},\frac{1}{2\rho}\big{(}|{\bf V}|^{2}-3\big{)}\right\},$$

$$g^{(1)}=\mathcal{L}^{-1}_{\mathcal{M}}\big{(}\mathcal{M}{\bf A}({\bf V})\big{)}:{\bf D}({\bf u})+2\mathcal{L}^{-1}_{\mathcal{M}}\big{(}\mathcal{M}{\bf B}({\bf V})\big{)}\cdot\frac{\nabla_{\bf x}T}{\sqrt{T}}.$$

$$\left\{\begin{array}[]{l}\displaystyle\bar{\bf A}^{\varepsilon}_{NS}\,:=\,\frac{1}{\rho}\int_{\mathbb{R}^{3}}{\bf A}({\bf V})\mathcal{M}({\bf v})\big{[}1+\varepsilon g^{(1)}({\bf v})\big{]}d{\bf v}\,=\,-\varepsilon\frac{\mu}{\rho\,T}{\bf D}({\bf u}),\\ \\ \displaystyle\bar{\bf B}^{\varepsilon}_{NS}\,:=\,\frac{1}{\rho}\int_{\mathbb{R}^{3}}{\bf B}({\bf V})\mathcal{M}({\bf v})\big{[}1+\varepsilon g^{(1)}({\bf v})\big{]}d{\bf v}\,=\,-\varepsilon\frac{\kappa}{\rho\,T^{3/2}}\nabla_{\bf x}\,T.\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\mu\,:=\,-T\int_{\mathbb{R}^{3}}\mathcal{M}({\bf v}){\bf A}({\bf V}):\mathcal{L}^{-1}_{\mathcal{M}}\big{(}\mathcal{M}{\bf A})({\bf v})d{\bf v},\\ \\ \displaystyle\kappa\,:=\,-T\int_{\mathbb{R}^{3}}\mathcal{M}({\bf v}){\bf B}({\bf v})\cdot\mathcal{L}^{-1}_{\mathcal{M}}\big{(}\mathcal{M}{\bf B})({\bf v})d{\bf v}.\end{array}\right.$$

$$\left\{\begin{array}[]{l}\displaystyle\partial_{t}\rho+\nabla_{\bf x}\cdot\,(\rho\,{\bf u})\,=\,0,\\ \\ \,\displaystyle\partial_{t}(\rho\,{\bf u})+\nabla_{\bf x}\cdot\,(\rho\,{\bf u}\otimes{\bf u}+\rho\,T\,{\bf I})\,=\,\varepsilon\,\nabla_{\bf x}\cdot({\bf\sigma}({\bf u})),\\ \\ \,\displaystyle\partial_{t}E+\nabla_{\bf x}\cdot\,({\bf u}\,(E+\rho\,T))\,=\,\varepsilon\,\nabla_{\bf x}\cdot({\bf\sigma}({\bf u})\cdot{\bf u}+{\bf q}),\end{array}\right.$$

Full text

Francis Filbet

Institut de Mathématiques de Toulouse, UMR5219 Université de Toulouse; CNRS & IUF

UPS, F-31400, Toulouse, France

Tao Xiong

School of Mathematical Sciences, Fujian Provincial Key Laboratory of Mathematical Modeling

and High-Performance Scientific Computing, Xiamen University

Xiamen, Fujian, 361005, P.R. China

1. 1 Introduction

2. 2 Hydrodynamic limit and domain decomposition

3. 2.1 The hydrodynamic limit

4. 2.2 Domain decomposition

5. 3 Hybrid discontinuous Galerkin scheme

6. 3.1 Preliminary

7. 3.2 Discontinuous Galerkin for the kinetic equation

8. 3.3 Discontinuous Galerkin for the compressible Navier-Stokes equations

9. 3.4 Interface coupling condition

10. 4 Numerical examples

11. 4.1 Flow caused by evaporation and condensation.

12. 4.2 2D Riemann problem.

13. 4.3 2D ghost effect.

14. 5 Conclusion

15. A BGK equation and Chu reduction

1. Introduction

Many physical problems, such as micro-electro-mechanical systems, involve boundary layers or transitional regimes which cannot be solved by using standard fluid models. Hence a kinetic model is needed to accurately describe complex phenomena occurring around the boundary. However, a kinetic description is computationally very expensive to simulate and we desire to use it only locally in space. An interesting approach is to design hybrid kinetic/fluid schemes, with domain decomposition which can automatically and accurately identify fluid and kinetic zones. In this case we only solve the kinetic model in a minimized area around the kinetic boundary layers or shocks, while take advantage of the low computational cost of numerical methods for the fluid system elsewhere.

Taking into account of collisions, the multi-scale kinetic equation for the description of particle dynamics in a dilute gas is given by

[TABLE]

with ${\bf x}\in\Omega\subset\mathbb{R}^{d_{x}},{\bf v}\in\mathbb{R}^{3}$. The open set $\Omega$ is a bounded Lipschitz-continuous domain in $\mathbb{R}^{d_{x}}$ and supplied with some boundary conditions, whereas the particle distribution function $f:=f(t,{\bf x},{\bf v})$ depends on time and phase space variables $({\bf x},{\bf v})\in\Omega\times\mathbb{R}^{3}$ and the initial datum $f_{0}$ is a non-negative function. The parameter $\varepsilon>0$ is the dimensionless Knudsen number, which is defined as the ratio of the mean free path of particles over a typical length scale such as the size of the spatial domain. It measures the rarefaction of the gas, namely, the gas is in a rarefied or kinetic regime if $\varepsilon\sim 1$ and in a dense or fluid regime if $\varepsilon\ll 1$. The collision operator ${\mathcal{Q}}(f)$ may be given by the full Boltzmann operator [13, 14]

[TABLE]

where we used the shorthand $f=f(v)$, $f_{*}=f(v_{*})$, $f^{{}^{\prime}}=f(v^{\prime})$, $f_{*}^{{}^{\prime}}=f(v_{*}^{{}^{\prime}})$. The velocities of the colliding pairs $(v,v_{*})$ and $(v^{\prime},v^{\prime}_{*})$ are related by

[TABLE]

The collision kernel $B$ is a non-negative function which by physical arguments of invariance only depends on $|{\bf v}-{\bf v}_{*}|$ and $\cos\theta={\widehat{\bf u}}\cdot\sigma$, where ${\widehat{\bf u}}=({\bf v}-{\bf v}_{*})/|{\bf v}-{\bf v}_{*}|$. A simpler choice is the BGK operator [9] given by

[TABLE]

or a modified ES-BGK operator [3]

[TABLE]

where $\nu$ is a collision frequency. In the ES-BGK operator, $\mathcal{G}(f)$ is a Gaussian defined as

[TABLE]

with a corrected stress tensor

[TABLE]

where ${\bf I}$ is the identity matrix. The parameter $\beta\in(-\infty,1)$ is used to modify the value of the Prandtl number through the formula

[TABLE]

The correct Prandtl number for a monoatomic gas of hard spheres is equal to $2/3$ for $\beta=-1/2$, whereas the classical BGK operator has a Prandtl number equal to $1$ for $\beta=0$. We refer to [25] for more discussions about these collision operators. These Boltzmann-like collision operators share the fundamental properties of conserving mass, momentum and energy, that is

[TABLE]

where $m({\bf v})=(1,{\bf v},\frac{1}{2}|{\bf v}|^{2})^{T}$ and super index $T$ denotes the transpose of the corresponding vector. The equilibrium of the collision operator, when ${\mathcal{Q}}(f)=0$, is given by the local Maxwellian distribution function $\mathcal{M}[f]$, which is defined as :

[TABLE]

The density $\rho$, mean velocity ${\bf u}$ and temperature $T$ are macroscopic moments of the distribution function $f$, which can be computed as

[TABLE]

There are already several works about hybrid methods for multi-scale kinetic equations in the literature, which mostly rely on the same domain decomposition technique. I. Boyd, G. Chen and G. Candler [11] first proposed a macroscopic criterion based on the local Knudsen number, they pass from a hydrodynamic description to a kinetic one when the quantity is below a (problem-dependent) threshold. It was practically used by V. Kolobov et. al. [30] with a discrete velocity model of direct numerical solution (DNS) for the Boltzmann equation and a gas-kinetic scheme for the hydrodynamic part, and recently by P. Degond and G. Dimarco [21] with a Monte Carlo solver for the kinetic equation and a finite volume method for the macroscopic ones. Another hydrodynamic breakdown indicator based on the viscous and heat flux of the Navier-Stokes equations is introduced by S. Tiwari [40] and has been used with various deterministic kinetic and hydrodynamic solvers [22, 41, 42, 1]. Recently, F. Filbet and T. Rey [25] proposed a domain decomposition indicator based on moment realizability matrices, a concept first introduced by D. Levermore, W. Morokoff and B. Nadiga [31]. They have shown that criteria from/to hydrodynamic to/from kinetic are both needed, and finite volume schemes for both the Boltzmann equation and the hydrodynamic equations are explored. This criterion has lately been used by T. Xiong and J.-M. Qiu [44] to form a hierarchy high order discontinuous Galerkin schemes for the BGK equation under a micro-macro decomposition framework.

However, this criterion requires some computational effort since it is based on the comparison of local eigenvalues on each cell, hence the cost is not negligible in high dimension. Furthermore, high order finite volume schemes often need wide stencils. This may easily be done in 1D in space [25], however, it is not very convenient or robust to be generalized to multi-dimensional problems, as it would require each region as wide as the stencil along each space dimension. It will cause some trouble when switching from one region to the other dynamically with time evolution.

The aim of this paper is to develop a hybrid discontinuous Galerkin scheme, which extends the work in [25] in high dimension. The well recognized discontinuous Galerkin method [29], due to its advantages of compactness, $h$-$p$ adaptivity, high efficiency on parallelization, and flexibility on complicated geometries, has been widely applied to physical and engineering problems. For compressible Euler or Navier-Stokes systems, F. Bassi and S. Rebay first introduced a discontinuous Galerkin method with macroscopic numerical flux such as Godunov flux for the advection and a primal formulation with centered numerical flux for the diffusion [7, 6, 8]. It has provided a good framework for solving the compressible Navier-Stokes system and many works are followed, e.g. a positivity preserving high order discontinuous Galerkin scheme recently proposed by X. Zhang [45]. Some other type of discontinuous Galerkin schemes can be referred to [33, 34, 18] and many references therein. For the kinetic Boltzmann, ellipsoidal statistical BGK equations, A. Alexeenko et al. [2] proposed a high order explicit Runge-Kutta discontinuous Galerkin method, with Newton’s iteration solving the nonlinear collisional source term. It has recently been generalized to 2D in space by W. Su et al. in [39]. On the other hand, to take care of the stiff collisional term and avoid a stringent time step size, T. Xiong et al. [43] proposed a high order asymptotic preserving nodal discontinuous Galerkin method coupled with an implicit-explicit scheme for the BGK equation. There are also some other works which would like to build a connection from the kinetic equation to the hydrodynamic system, and discontinuous Galerkin methods for compressible Navier-Stokes equations with kinetic flux-vector splitting are proposed, such as the kinetic flux-vector splitting flux following Chapman-Enskog velocity distribution function [19, 16], or gas-kinetic schemes mimicking a Hilbert expansion [32, 35], etc.

In this paper, we propose a simplified indicator, based on [25], to determine the kinetic and fluid regions and develop a discontinuous Galerkin method, which is well suited to couple the kinetic and fluid model at the interface. Due to the compactness of the discontinuous Galerkin method, the interface coupling condition for two different regions in the hybrid discontinuous Galerkin scheme only requires consistent numerical fluxes defined at the cell interfaces. For the domain decomposition, it allows isolated hydrodynamic or kinetic cells, which will be very convenient especially when extended to high-dimensional space problems and/or on unstructured meshes. We will solve the kinetic equation in the kinetic region with a discontinuous Galerkin scheme under the asymptotic preserving framework [24], while solving the compressible Navier-Stokes equations, which can be (formally) obtained by Chapman-Enskog expansion [13, 15, 17, 5], by following the discontinuous Galerkin scheme of F. Bassi and S. Rebay [6]. We apply an upwind numerical flux for the kinetic solver, while to be consistent, especially in order to match the numerical fluxes at the cell interface between two different regions, we define a flux for the hydrodynamic solver obtained as the asymptotic limit of the kinetic flux. Let us emphasize that this particular hydrodynamic flux has already been introduced in [19] as a kinetic flux splitting scheme. Furthermore, we will simplify the computation of the moment realizability criterion for domain decomposition, by extracting some main derivatives in the moment realization matrix. For illustration in our numerical example, we take the BGK equation and reduce the 3D in velocity by the technique of Chu reduction [20]. Our new proposed hybrid discontinuous Galerkin scheme will be more robust on $h$-$p$ adaptivity, parallelization, and on high dimensional problems with unstructured meshes. We will perform some numerical tests on some physical relevant problems with shocks or boundary layers, such as flow caused by evaporation and condensation, 2D Riemann problem as well as 2D ghost effect. The results will demonstrate the efficiency and effectiveness of our proposed approach.

The rest of the paper is organized as follows. In Section 2, we recall the hydrodynamic limit of the multi-scale kinetic equation based on Chapman-Enskog expansion, and describe the domain decomposition from moment realizability criterion. In Section 3, a discontinuous Galerkin scheme for the kinetic equation and a discontinuous Galerkin scheme for the compressible Navier-Stokes equations will be formed, and consistent numerical fluxes for the interface coupling condition between two different regions are stated. Numerical tests are followed in Section 4. Conclusion and future work are made in Section 5.

2. Hydrodynamic limit and domain decomposition

In this section, we will recall the hydrodynamic limit of the kinetic equation (1.1). By formally doing the Chapman-Enskog expansion, the equation (1.1), to the zeroth order limit as $\varepsilon\rightarrow 0$ will converge to the compressible Euler equations, while to the first order limit when $0<\varepsilon\ll 1$, it is going to be the compressible Navier-Stokes equations. In the following, we will briefly review the derivation of compressible Euler and Navier-Stokes limit [25, 5], and describe the domain decomposition criteria to divide the domain into kinetic and hydrodynamic regions, which will be used to define our hybrid scheme.

2.1. The hydrodynamic limit

According to the conservation property (1.5), by integrating in the velocity space, without any closure, we have

[TABLE]

Let ${\bf U}$ denote the conservative macroscopic components of density, momentum and energy, which correspond to the first three moments of the distribution function $f$:

[TABLE]

where the energy is $E=\frac{1}{2}\rho|{\bf u}|^{2}+\frac{3}{2}\rho\,T$. $f$ can be approximated by the Chapman-Enskog expansion, which is

[TABLE]

where $\mathcal{M}$ is in short of $\mathcal{M}[f]$ and the fluctuations satisfy

[TABLE]

Substituting $f^{\varepsilon}({\bf v})$ into (2.1), we get

[TABLE]

where the traceless matrix $\bar{\bf A}^{\varepsilon}\in\mathbb{M}^{3}$ and vector $\bar{\bf B}^{\varepsilon}\in\mathbb{R}^{3}$ are

[TABLE]

with

[TABLE]

and

[TABLE]

In the zeroth order limit when $f^{\varepsilon}({\bf v})=\mathcal{M}$, since the matrix $\bar{\bf A}^{\varepsilon}$ is traceless, and due to $\bar{\bf B}^{\varepsilon}$ only involving odd, centered moments of $f$, we obtain

[TABLE]

hence (2.3) gives us the compressible Euler system

[TABLE]

In the first order limit, by plugging (2.2) into (2.3), since the Maxwellian distribution $\mathcal{M}$ is an equilibrium of the collision operator ${\mathcal{Q}}(\mathcal{M})=0$, the fluctuation $g^{(1)}$ is given by

[TABLE]

where $\mathcal{L}_{\mathcal{M}}$ is the linearized collision operator around the Maxwellian distribution $\mathcal{M}$. Besides, we also have

[TABLE]

Replacing the time derivatives by spatial ones from (2.3), and dropping all terms of order $\varepsilon$ in (2.7), after some computations, we get

[TABLE]

where ${\bf A}$ and ${\bf B}$ are defined in (2.4) and the deformation tensor ${\bf D}({\bf u})$ is given by

[TABLE]

The linear operator $\mathcal{L}_{\mathcal{M}}$ is invertible on the orthogonal of its kernel [25], where its kernel is

[TABLE]

and ${\bf A}({\bf V}),\,{\bf B}({\bf V})\,\in\,{\rm ker}^{\perp}\left(\mathcal{L}_{\mathcal{M}}\right)$. From (2.8), it yields

[TABLE]

Plugging this expression into (2.4), it gives

[TABLE]

The viscosity and thermal conductivity coefficients are given by

[TABLE]

Substituting (2.11) into (2.3), it gives us the compressible Navier-Stokes system

[TABLE]

where ${\bf\sigma}:=-\mu\,{\bf D}({\bf u})$ is the viscosity tensor and ${\bf q}:=-\kappa\nabla_{\bf x}T$ is the heat flux.

For the Boltzmann equation in the hard sphere case, $\mu$ and $\kappa$ can be expressed as [26]

[TABLE]

for some positive constants $\mu_{0}$ and $\kappa_{0}$. In the BGK case, $\mu$ and $\kappa$ are related to the macroscopic quantities of $\rho$ and $T$ by [38]

[TABLE]

and the collision frequency can be taken as $\nu=\frac{2}{\sqrt{\pi}}\rho$ [4].

The compressible Navier-Stokes equations (2.13) can be written in a vector form as

[TABLE]

with the advection flux ${\bf F}^{a}({\bf U})$ and the viscous diffusion flux ${\bf F}^{d}({\bf U},\nabla_{\bf x}{\bf U})$ to be

[TABLE]

Moreover, if we define

[TABLE]

it becomes

[TABLE]

The flux function ${\bf F}({\bf U},\nabla_{\bf x}{\bf U})$ corresponds to the kinetic description (2.1) by taking the first order distribution function $f^{\varepsilon}({\bf v})$ (2.2), that is

[TABLE]

Formally by letting $\varepsilon\rightarrow 0$, (2.18) corresponds to the compressible Euler equations

[TABLE]

2.2. Domain decomposition

We now follow the domain decomposition indicator as defined in [25], to divide the computation domain into kinetic and hydrodynamic regions, where in the kinetic region (mainly around shocks or boundary layers) the kinetic equation (1.1) is used, otherwise in the hydrodynamic region, the low computational cost compressible Navier-Stokes equations (2.13) or (2.18) as the hydrodynamic equation will be solved. We review the main principle on how to define the criteria, however we would like to simplify the tedious computation of the eigenvalues for the moment realizability matrix.

From fluid to kinetic, the hydrodynamic breakdown criterion is defined by the moment realizability matrix. The deviation of the eigenvalue for the moment realizability matrix measures the appropriateness of the corresponding hydrodynamic limit. The compressible Euler system has all eigenvalues to be $1$ with an identity moment realizability matrix. Here since we use the compressible Navier-Stokes equations as the fluid model, we only concern about the deviation of the eigenvalues between the first order compressible Navier-Stokes equations and the second order Burnett equations. We recall the moment realization matrix from [25], which is given by

[TABLE]

For compressible Navier-Stokes equations, $\bar{{\bf A}}^{\varepsilon}_{NS}$ and $\bar{{\bf B}}^{\varepsilon}_{NS}$ are defined in (2.11), while for the second order Burnett equations, $\bar{{\bf A}}^{\varepsilon}_{Burnett}$ and $\bar{{\bf B}}^{\varepsilon}_{Burnett}$ are [25]

[TABLE]

Noticing that for the main components $\bar{{\bf A}}^{\varepsilon}$ and $\bar{{\bf B}}^{\varepsilon}$ in (2.21), the major differences between $\bar{{\bf A}}^{\varepsilon}_{NS}$ and $\bar{{\bf A}}^{\varepsilon}_{Burnett}$ as well as $\bar{{\bf B}}^{\varepsilon}_{NS}$ and $\bar{{\bf B}}^{\varepsilon}_{Burnett}$ are the $\mathcal{O}(\varepsilon^{2})$ terms appeared in (2.22). As an indicator, we may only need to extract the crucial spatial derivatives appeared in those terms and roughly measure their magnitudes, to simplify the tedious computations of all terms, thus avoiding too much computational cost on the domain decomposition indicator along this side. From (2.14) and (2.15), either $\mu/\sqrt{T}$ and $\kappa/\sqrt{T}$ for the Boltzmann operator or $\nu\mu/(\rho\,T)$ and $\nu\kappa/(\rho\,T)$ for the BGK operator are $\mathcal{O}(1)$ if we assume macroscopic quantities $\rho$ and $T$ are $\mathcal{O}(1)$, a new indicator from fluid to kinetic we are now using is defined as

[TABLE]

We take the fluid model of compressible Navier-Stokes equations to be not appropriate at some point $(t,{\bf x})$ if

[TABLE]

$\lambda_{\varepsilon^{2}}(t,{\bf x})$ is measured at the cell center in the numerical section and we take $\eta_{0}=10^{-3}$.

From kinetic to fluid, the criterion is that a kinetic description at some point $(t,{\bf x})$ corresponds to a hydrodynamic equations if

[TABLE]

where $\delta_{0}$ is a small parameter and we take $\delta_{0}=10^{-3}$ and the function $f_{2}=\mathcal{M}[1+\varepsilon g^{(1)}]$ corresponds to the Chapman-Enskog expansion (2.2) up to first order.

3. Hybrid discontinuous Galerkin scheme

In this section, we will propose our hybrid discontinuous Galerkin scheme for the multi-scale kinetic equations, namely we define either an approximation of the kinetic equation (1.1) according to the value of the indicator (2.23) or another one of the fluid equations (2.18) in the region where (2.24) is satisfied. Due to the compactness of discontinuous Galerkin schemes, only consistent numerical fluxes need to be defined at the cell interface between two regions, which is also important to ensure mass conservation.

We adopt the discontinuous Galerkin scheme for spatial discretization, as it is more convenient than a finite volume scheme used in [25], especially when extending to high dimensions in space and in the domain decomposition. Indeed, the finite volume method requires either the kinetic or fluid region to be at least as wide as the stencil of the scheme. Besides, the discontinuous Galerkin scheme is $h$-$p$ adaptive and very flexible to nonuniform meshes, making it more suitable to physical problems with boundary layers (see Section 4).

Our hybrid discontinuous Galerkin scheme can be defined for the multi-scale kinetic equation with either the full Boltzmann collision operator or the (ES-) BGK operator. The Boltzmann operator (1.2) in a bounded domain can be discretized by the method introduced in [23] and references therein, and in an asymptotic preserving framework can be penalized by the BGK operator [24]. We will only take the BGK operator (1.3) in the kinetic equation and describe the discontinuous Galerkin scheme for the BGK equation (1.1). The ES-BGK operator (1.4) can be done similarly and we can use the technique in [24] for the Boltzmann operator. The discontinuous Galerkin scheme for the hydrodynamic equations (2.18) is followed. Only the discontinuous Galerkin scheme in 2D in space is presented, the 1D case can be easily deduced from the 2D ones. We take the kinetic flux-vector splitting [19] to define the hydrodynamic numerical flux at the cell interface between two hydrodynamic cells or two cells between two regions. A nodal discontinuous Galerkin scheme [28] is implemented to produce the numerical results.

For the velocity space, we will take a large enough cut-off domain $\Omega_{v}=[-V_{c},V_{c}]$ and discretize it uniformly with $N_{v}$ points, $\{v^{j}\}_{j=1}^{N_{v}}$, along each $v$ direction. We integrate the velocity space by mid-point rule, which is spectrally accurate for smooth functions with periodic boundary conditions or compact supports [12], as there is no differential operator acting on ${\bf v}$. However, we note that the domain cut-off and the discretization in ${\bf v}$ will lead to conservation of mass, momentum and energy not exactly but approximately up to the integral error along ${\bf v}$ direction. In the following, for easy presentation, we still keep ${\bf v}$ to be continuous and only discuss the discretizations in space and in time.

3.1. Preliminary

We consider the space domain in 2D to be a rectangular $\Omega_{\bf x}=[a,\,b]\times[c,\,d]$, and divide it by $a=x_{1/2}<x_{{3}/{2}}<\cdots<x_{N_{x}+1/2}=b$ and $c=y_{1/2}<y_{{3}/{2}}<\cdots<y_{N_{y}+1/2}=d$. Let $I_{i,j}=I_{i}\times I_{j}=[{x_{i-\frac{1}{2}}},{x_{i+\frac{1}{2}}}]\times[{y_{j-\frac{1}{2}}},{y_{j+\frac{1}{2}}}]$ denote an element with its length $\Delta x_{i}={x_{i+\frac{1}{2}}}-{x_{i-\frac{1}{2}}}$ and $\Delta y_{j}={y_{j+\frac{1}{2}}}-{y_{j-\frac{1}{2}}}$. Let $h_{x}=\max_{i=1}^{N_{x}}\Delta x_{i}$, $h_{y}=\max_{j=1}^{N_{y}}\Delta y_{j}$ and $h=\max(h_{x},\,h_{y})$. Given any non-negative integer vector ${\bf K}=(K_{1},K_{2})$, we define a finite dimensional discrete space,

[TABLE]

Its vector version is denoted as

[TABLE]

For simplicity, we take the 2D local space $Q^{\bf K}(I_{i,j})$ as a tensor product space of $P^{K_{1}}(I_{i})\otimes P^{K_{2}}(I_{j})$, where $P^{K}(I)$ consists of polynomials of degree at most $K$ on $I$. Functions in $\mathcal{Z}_{h}^{\bf K}$ are piecewise defined and may be discontinuous across cell interfaces. The left and right limits of a function $u\in\mathcal{Z}_{h}^{\bf K}$ at the interface $(x_{i+\frac{1}{2}},y)$ along $y$ direction are denoted as $u(x_{i+\frac{1}{2}}^{\pm},y)=\lim_{\varepsilon\rightarrow\pm 0}u(x_{i+\frac{1}{2}}+\varepsilon,y)$, similarly for $u(x,y^{\pm}_{j+\frac{1}{2}})$.

We will use a nodal basis to represent functions in the discrete space $P^{K}(I)$, and approximate integrals by numerical quadratures. Specifically, we choose the local nodal basis (also called Lagrangian basis) $\{\phi_{i}^{k}(x)\}_{k=0}^{K}$ associated with the $K+1$ Gaussian quadrature points $\{x^{k}_{i}\}^{K}_{k=0}$ on $I_{i}$, defined as below

[TABLE]

where $\delta_{kk^{\prime}}$ is the Kronecker delta function. We denote by $\{\omega_{k}\}^{K}_{k=0}$ the corresponding quadrature weights on the reference element $(-{1}/{2},{1}/{2})$. For the two dimensional local nodal basis, we choose a tensor product of the one dimensional local nodal basis (3.3) but associated with $K_{1}+1$ and $K_{2}+1$ Gaussian quadrature points on $I_{i}$ and $I_{j}$ along $x$ and $y$ directions respectively.

3.2. Discontinuous Galerkin for the kinetic equation

For the kinetic equation (1.1) with the BGK operator, we take an implicit-explicit (IMEX) time discretization as in [24, 25]. Considering the continuous problem in phase space $({\bf x},{\bf v})$, a first order implicit/explicit scheme is defined as follows to pass from time level $t^{n}$ to $t^{n+1}$

[TABLE]

The implicit Maxwellian $\mathcal{M}({\bf v},{\bf U}^{n+1})$ is first computed from the macroscopic quantity ${\bf U}^{n+1}$, where

[TABLE]

For each ${\bf v}$, if we define

[TABLE]

in vector form we have

[TABLE]

Following this strategy and for simplicity keeping ${\bf v}\in\mathbb{R}^{3}$, a two dimensional in space discontinuous Galerkin scheme for (3.7) is defined as follows: for given $f^{n}_{h}({\bf v})$, we find $f^{n+1}_{h}({\bf v})$, by computing a discrete approximation ${\bf R}_{h}^{n+1}$, which solves the following problem : for any $\zeta\in\mathcal{Z}_{h}^{\bf K}$ and for $1\leq i\leq N_{x}$, $1\leq j\leq N_{y}$,

[TABLE]

then we compute ${\bf U}^{n+1}_{h}$ given for any $\beta\in\mathcal{Z}_{h}^{\bf K}$ and for $1\leq i\leq N_{x}$, $1\leq j\leq N_{y}$, by

[TABLE]

Finally the discrete distribution function $f_{h}^{n+1}$ is defined such that for any $\alpha\in\mathcal{Z}_{h}^{\bf K}$ and for $1\leq i\leq N_{x}$, $1\leq j\leq N_{y}$,

[TABLE]

where $\widetilde{vf}$ is an upwind numerical flux along its direction,

[TABLE]

whereas $f^{\pm}$ are the left and right limits of $f^{n}_{h}$ at the cell interface of two adjacent cells respectively.

For a nodal discontinuous Galerkin scheme with Lagrangian bases, the scheme (3.2) becomes

[TABLE]

whereas ${\bf U}$ is given by

[TABLE]

and the discrete approximation $f_{h}^{n+1}$ is defined by

[TABLE]

where $f^{n}_{i,k_{1},j,k_{2}}({\bf v}),{\bf R}^{n+1}_{i,k_{1},j,k_{2}}({\bf v}),{\bf U}^{n+1}_{i,k_{1},j,k_{2}}$ with subindex $(i,k_{1},j,k_{2})$ are the corresponding numerical values at the $(k_{1},k_{2})$-th Gaussian quadrature point in cell $I_{i,j}$, for $0\leq k_{1}\leq K_{1}$ and $0\leq k_{2}\leq K_{2}$.

3.3. Discontinuous Galerkin for the compressible Navier-Stokes equations

We will follow F. Bassi and S. Rebay [6] to define a discontinuous Galerkin scheme for the compressible Navier-Stokes equations (2.18). Letting $\nabla_{\bf x}{\bf U}:=({\bf S}_{1},{\bf S}_{2})$, with an Euler forward time discretization, the discontinuous Galerkin scheme for (2.18) is defined as follows: we seek ${\bf U}^{n}_{h}\in{{\bf Z}}_{h}^{\bf K}$, such that for any $\eta({\bf x})\in\mathcal{Z}_{h}^{{\bf K}}$ and $1\leq i\leq N_{x}$, $1\leq j\leq N_{y}$, we have

[TABLE]

whereas $\nabla_{\bf x}{\bf U}^{n}_{h}=({\bf S}^{n}_{1,h},{\bf S}^{n}_{2,h})$ is such that ${\bf S}^{n}_{1,h},{\bf S}^{n}_{2,h}\in{{\bf Z}}_{h}^{\bf K}$ for any $\eta({\bf x})\in\mathcal{Z}_{h}^{{\bf K}}$ and $1\leq i\leq N_{x}$, $1\leq j\leq N_{y}$, we have

[TABLE]

and

[TABLE]

where $\hat{\bf U}$ is taken to be a central flux

[TABLE]

along $x$ or $y$ direction. ${\bf U}^{\pm}$ are the left and right limits of ${\bf U}^{n}_{h}$ at the cell interface. Finally to construct a nodal discontinuous Galerkin scheme, we apply a quadrature formula as it has been done in the previous section.

The flux $\hat{\bf F}$ is defined in the next subsection. An important issue is to construct a numerical flux which is consistent with the kinetic flux in the asymptotic limit $\varepsilon\rightarrow 0$ in order to avoid spurious oscillations of order $\varepsilon$ or $\varepsilon^{2}$.

3.4. Interface coupling condition

For a discontinuous Galerkin scheme, due to its compactness, we only need to pass an interface coupling condition between two cells of different regions, that is, to define consistent numerical fluxes at such a cell interface for both the kinetic equation and the hydrodynamic system. Let us describe how to define the numerical fluxes $\widetilde{v_{1}f}$ and $\hat{{\bf F}}_{1}$ in the $x$ direction, while $\widetilde{v_{2}f}$ and $\hat{{\bf F}}_{2}$ in the $y$ direction can be defined similarly.

We assume $I_{i,j}$ is a kinetic cell and $I_{i+1,j}$ is a hydrodynamic cell, and the cell interface is at $\{x_{i+\frac{1}{2}}\}\times[{y_{j-\frac{1}{2}}},{y_{j+\frac{1}{2}}}]$:

• •

For the discontinuous Galerkin scheme of the kinetic equation (3.2)-(3.2), an upwind flux (3.11) in (3.2) is needed,

[TABLE]

at each Gaussian quadrature point within $[{y_{j-\frac{1}{2}}},{y_{j+\frac{1}{2}}}]$. For $v_{1}\geq 0$, $f^{-}$ in cell $I_{i,j}$ can be defined from $f^{n}_{h}$, while $f^{+}$ for $v_{1}<0$ in cell $I_{i+1,j}$ cannot. In this case, we take $f^{n}_{h}$ in cell $I_{i+1,j}$ to be the first order truncated distribution functions as defined in (2.2) with $g^{(1)}$ defined in (2.10), so that $f^{+}$ can be recovered from $v_{1}$, ${\bf U}(x^{+}_{i+\frac{1}{2}},\cdot)$ and ${\bf S}(x^{+}_{{i+\frac{1}{2}}},\cdot)$, where ${\bf S}=\nabla_{\bf x}{\bf U}$.

• •

For the discontinuous Galerkin scheme of the fluid equations (3.3)-(3.3), for conservation, the flux function $\hat{{\bf F}}_{1}$ in (3.3) is defined as the integral of the kinetic upwind flux (3.15) on ${\bf v}$,

[TABLE]

Based on $f^{\pm}$ in (3.15), it can be obtained.

The other way when $I_{i,j}$ is a hydrodynamic cell and $I_{i+1,j}$ is a kinetic cell can be defined symmetrically. The hydrodynamic flux $\hat{\bf F}_{1}(x_{i+\frac{1}{2}},\cdot)$ at the interface of two cells inside the hydrodynamic region is defined in the same form as (3.16), with $f^{\pm}$ both taken as the first order truncated distribution function (2.2) in cell $I_{i,j}$ and $I_{i+1,j}$ respectively, which is known as the kinetic flux-vector splitting flux [19].

This choice is particularly important to avoid that some spurious oscillations appear when the coupling between kinetic and fluid regions occur. Furthermore, this numerical flux is consistent with the continuous flux given in (2.19). Indeed, it is given by,

[TABLE]

When we consider the BGK equation, we have

[TABLE]

where ${\bf A}$ and ${\bf B}$ are given in (2.4). The resulting integrals on $v_{1}$, either on $[0,\infty]$ or on $[-\infty,0]$, can be expressed in terms of the Gauss error function

[TABLE]

and some of its related functions. We can get explicit expressions in terms of $\textrm{erfc}(x)$ to avoid integration on $v_{1}$. We mainly need to do some integrals in the following form

[TABLE]

for $0\leq n\leq 3$ and $0\leq m\leq 3$. First by letting

[TABLE]

It can be transformed to

[TABLE]

Except the coefficient $\frac{\rho}{2}$, it is in the form of

[TABLE]

which is the Gauss error function $\textrm{erfc}(s)$ in case of $n=0$ and $m=0$. Expanding on $(a\,z+b)^{n}$ for $n=0,1,2,3$, what we need are

[TABLE]

for $1\leq n\leq 6$. From integration by parts, they are

[TABLE]

For the integral

[TABLE]

on the lower half plane, by letting $v=z\,\sqrt{2T}+u$, it can be transformed to

[TABLE]

replacing $z$ by $-z$, it becomes to

[TABLE]

except the coefficient ${\rho}/{2}$, it is in the form of (3.17).

4. Numerical examples

In this section, we will test the hybrid discontinuous Galerkin scheme with first order time discretizations for some physical relevant problems. For simplicity and reduce the high dimension in velocity, we take the reduced BGK system described in Appendix A as our kinetic equation. No limiters on the discontinuous Galerkin method are applied. We take a second order discontinuous Galerkin method for 1D problem, while for 2D problems specify it in the examples. We take a cut-off domain with $V_{c}=8$ and discretize it with $N_{v}=32$ uniform points along each direction if not specified.

4.1. Flow caused by evaporation and condensation.

In this example, we consider a rarefied gas flow caused by evaporation and condensation between two parallel planes with condensed phases [4]. We consider two cases: one is a weak evaporation and condensation with phase conditions to be

[TABLE]

the other is a strong evaporation and condensation

[TABLE]

The boundary conditions are fixed wall temperature and pressure, and we take the initial conditions to be a linear function connecting the wall temperature or pressure. We take a nonuniform mesh with $N_{x}/4$ cells on the width of $0.05$ at the boundary, while with $N_{x}/2$ cells in the center on a width of $0.9$. For the hybrid scheme, we force the cells within a width of $0.1$ at the boundary always to be in the kinetic region.

In Figure 4.1, we show the pressure and temperature with $N_{x}=40$ for the weak evaporation and condensation (4.1), at $t=100$ for $\varepsilon=10^{-1},10^{-2},10^{-3}$ respectively. The results are comparable to those in [4, Figure 2] obtained with a much refined mesh $N_{x}=600$ for $10^{-3}\leq\varepsilon<1$. The mean velocity is nearly a constant up to the spatial discretization error and is omitted. For this example, our hybrid method can well match the full kinetic scheme and identify the flat part of the steady state solution in the middle to be in the hydrodynamic region, while the boundary layers are in the kinetic region. We notice that even for a relatively large $\varepsilon=10^{-1}$, the hybrid scheme has a result very comparable to the full kinetic scheme. Furthermore we compare the CPU cost between the hybrid scheme and the full kinetic scheme. We take a uniform mesh with $N_{x}=40$ and time step $1/5000$. We run the code up to $t=40$ for three times. For the full kinetic scheme, the averaged CPU cost (real time) is $2$ minutes and $45$ seconds, while the hybrid scheme costs about $1$ minute $50$ seconds for $\varepsilon=10^{-2}$ and $1$ minute $33$ seconds for $\varepsilon=10^{-3}$, so the hybrid scheme saves about $1/3$ and $2/5$ of the CPU cost for $\varepsilon=10^{-2}$ and $\varepsilon=10^{-3}$ respectively. We note that for the nonuniform mesh, since we put half of mesh points around the kinetic boundary, the CPU costs for the hybrid scheme and the full kinetic scheme are almost the same. However, we would emphasize that here we use the reduced BGK system as the kinetic model, the savings of a hybrid scheme would be more significant if we consider the kinetic equation with the full Boltzmann collision operator and full velocity in 3D. Nevertheless the hybrid scheme seems to be only very efficient when the kinetic region (and kinetic cells) is very locally, which may usually happen when $\varepsilon$ is relatively small, e.g., $\varepsilon<10^{-1}$. Besides, for relatively large $\varepsilon$, e.g., $\varepsilon\geq 10^{-1}$, very less hydrodynamic cells can be identified for a hybrid scheme. It will be more efficient to directly use a full kinetic scheme, to avoid the computation of the domain indicator. In the following, we will use $\varepsilon=10^{-1}$ as a reference and only consider our hybrid scheme for $\varepsilon<10^{-1}$.

Next, in Figure 4.2, we show the pressure, mean velocity and temperature with $N_{x}=40$ for the strong evaporation and condensation (4.2), at $t=100$ for $\varepsilon=10^{-1},10^{-2},10^{-3}$ respectively. The mean velocity is scaled by a factor of $-{1}/{\sqrt{2}}$. The results are comparable to those in [4, Figure 11] for the results on a much refined nonuniform mesh. This example is more demanding, as it requires a good resolution for small values like $\varepsilon=10^{-3}$ in order to well capture the boundary layer on both sides. We can observe that for large $\varepsilon=10^{-1}$, the hybrid scheme automatically becomes a full kinetic scheme, whereas when $\varepsilon$ goes to zero, the numerical scheme well identifies the hydrodynamic region away from the boundary layers. In conclusion, the above two cases show the capability of our hybrid scheme, especially the good capturing of our defined domain decomposition indicator.

4.2. 2D Riemann problem.

In this example, we now consider a 2D Riemann problem for polytropic gas with initial datum in four quadrants [36],

[TABLE]

here we have $\gamma=5/3$.

To avoid numerical oscillations from high order discontinuous Galerkin discretization for the shock problem (usually limiters are needed in order to control numerical oscillations), we take first order discontinuous Galerkin discretizations along both $x$ and $y$ directions. The mesh size is uniform and $N_{x}=N_{y}=80$. For the velocity, we take a cut-off domain $\Omega_{v}=[-8,8]$ with $N_{v}=64$ along each direction, due to discontinuous macroscopic density and temperature in the Maxwellian distribution function. We run the code with MPI parallelization along $y$ direction with $4$ processors.

In Figure 4.3, we first show the density profiles with $\varepsilon=10^{-2}$, at time $t=0.01,0.2,0.35$ for the hybrid scheme and the full kinetic scheme, as well as the domain indicators for the hybrid scheme. By comparing the contour lines, we can see that the hybrid scheme and the full kinetic scheme have almost the same results. However, due to a little large viscosity of $\varepsilon=10^{-2}$, the four shock lines are smeared as time evolves. Our domain decomposition indicator can well capture the kinetic region around the smeared shock lines. The corresponding temperature profiles are displayed in Figure 4.4, similarly almost the same results for the hybrid scheme and the full kinetic scheme can be observed.

In Figures 4.5 and 4.6, we now show the density and temperature profiles with $\varepsilon=10^{-3}$, at time $t=0.01,0.2,0.35$ for the hybrid scheme and the full kinetic scheme, as well as the domain indicators for the hybrid scheme. Similarly the hybrid scheme has almost the same results as the full kinetic one. However, due to small viscosity and sharp shock lines, our domain indicator tightly follows the moving of the shock lines. This problem well demonstrates the good performance of our hybrid scheme and the domain indicator in the 2D case.

Concerning the CPU cost for the 2D problem, we run the code with MPI parallelization on 4 processors up to $t=0.02$ with a time step $\Delta t\,=\,2.5\,\times\,10^{-5}$ and $800$ time steps for $\varepsilon=10^{-2}$, and $t=0.2$ with a time step $\Delta t=10^{-4}$ with $2000$ time steps for $\varepsilon=10^{-3}$ respectively. We run the code three times, and the averaged CPU cost (real time) for the hybrid scheme is about $7$ minutes $33$ seconds for $\varepsilon=10^{-2}$ and $15$ minutes and $4$ seconds for $\varepsilon=10^{-3}$, while a full kinetic scheme costs about $14$ minutes and $54$ seconds for $\varepsilon=10^{-2}$ and $37$ minutes for $\varepsilon=10^{-3}$, so that the hybrid scheme saves about $1/2$ for $\varepsilon=10^{-2}$ and $3/5$ for $\varepsilon=10^{-3}$ of the CPU cost. For the case of $\varepsilon=10^{-2}$, due to large smearing of shock profile, more and more kinetic cells are identified as time goes on, see Figure 4.3. Our hybrid scheme although still is effective on capturing the kinetic cells, but is getting less efficient. However, we would also mention that there are some communication costs from MPI parallelization, generally the hybrid scheme may save even more as compared to the full kinetic scheme.

4.3. 2D ghost effect.

On the basis of kinetic theory, the heat-conduction equation from the stationary state solution of the classical Navier-Stokes equations is not suitable for describing the temperature field of a gas in the continuum limit in an infinite domain without flow at infinity, where the flow vanishes in this limit. By the asymptotic theory, as the Knudsen number of the system approaches zero, the temperature field should be obtained by the kinetic equation, this phenomenon is called the ghost effect [23, 37, 10]. In this example, we will numerically study this effect based on our hybrid scheme.

We consider a rarefied gas between two parallel plane walls at $x=0$ and $x=1$. The walls both have a common periodic temperature distribution $T_{w}$

[TABLE]

and move in a common small mean velocity $u_{w}$ of order $\varepsilon$

[TABLE]

We take a uniform mesh with $N_{x}=N_{y}=40$ and time step $\Delta t=1/5000$. For velocity, since the solution is smooth and does not vary too much, we take a cut-off domain $\Omega_{v}=[-8,8]$ and discretize it with $N_{v}=16$ points along each direction. We apply a second order nodal discontinuous Galerkin scheme and force the grids inside the width of $0.1$ along $x$ direction around the walls always to be kinetic. We show the results (rotated by $90^{\circ}$) for $\varepsilon=0.02$ in Figure 4.7 for the isothermal lines, the mean velocity field as well as the domain indicators for the hybrid scheme, at time $t=10,80$ respectively. The solution at $t=80$ is assumed to be close to a steady state. We can observe that the solution at $t=80$ is similar to the results in [37] which is obtained from solving the BGK equation by a finite difference scheme, and also very similar to the results obtained from a $2d_{x}\times 2d_{v}$ Boltzmann system [23].

From the domain indicators at different time in Figure 4.7, we can find that the indicator well defines the kinetic region where the solution has large variants (left and right sides) due to the effect from the moving walls, while other parts are almost in the hydrodynamic region (note that top and bottom are forced to be in the kinetic region). This example with long time simulation well demonstrates the good property of our hybrid scheme, in which we have almost only solved the expensive kinetic equation in the middle region where it is necessary and use the cheap hydrodynamic equations elsewhere. Due to the compactness of our discontinuous Galerkin scheme, we may observe some isolated hydrodynamic or kinetic cells, which show the robustness of the hybrid discontinuous Galerkin method in the multi-dimensional case, while this cannot be easily achieved by a finite volume scheme due to the requirement of wide stencils.

5. Conclusion

In this paper, we developed a hierarchical of hybrid discontinuous Galerkin scheme for some physically relevant problems in both 1D and 2D space dimensions. The compact discontinuous Galerkin scheme has shown its robustness on the domain decomposition and $h$-$p$ adaptivity on capturing the boundary layers. Although only a first order time discretization is used in the paper for some long time simulations, it can be easily generalized to high order in time based on some high order implicit-explicit time discretizations. Extensions to unstructured meshes on more complicated geometries and to ES-BGK or Boltzmann collision operators will be investigated in our future work.

Acknowledgement

T. Xiong acknowledges support from the Marie Skłodowska-Curie Individual Fellowships H2020-MSCA-IF-2014 of the European Commission, under the project HNSKMAP 654175. This work has also been supported by NSFC grant 11601455, NSAF grant U1630247, NSF grant of Fujian Province 2016J05022 and the Fundamental Research Funds for the Central Universities No. 20720160009.

Appendix A BGK equation and Chu reduction

We consider three dimensional problems in velocity, but to save the computational cost, we take the BGK equation and use the technique of Chu reduction [20]. It has also recently been used by M. Groppi, G. Russo and G. Stracquadanio in [27]. The main idea is that under suitable geometry symmetry assumptions, we can integrate the distribution function $f$ in the velocity phase space and introduce several auxiliary variables, to transform the BGK equation (1.1) in 3D in velocity into a system of equations in 2D or 1D in velocity, such that the system has the same dimension in space and in velocity. We refer to [20] for more details.

Chu reduction

We start with the BGK equation in 3D in both space and velocity, which reads

[TABLE]

${\bf x}=(x,y,z)$, ${\bf v}=(v_{1},v_{2},v_{3})$ and $\mathcal{M}[f]$ is the Maxwellian defined as

[TABLE]

where ${\bf u}=(u_{1},u_{2},u_{3})$ is the macroscopic mean velocity. As explained in [20, 27], we consider the BGK equation to problems with axial symmetry, where all transverse spatial gradients vanish and the gas is drifting only in the axial directions. In such cases, $f(t,{\bf x},{\bf v})$ depends on the full velocity, that is, molecular trajectories are 3D, however all transverse components of the macroscopic velocity ${\bf u}$ vanish.

We first consider problems are axial symmetry with respect to two axises $(x_{1},x_{2})$, so that the problem does not depend on $x_{3}$ and the macroscopic velocity ${\bf u}$ vanishes at $u_{3}=0$. For simplicity, let us now denote $\tilde{\bf x}=(x_{1},x_{2})$, $\tilde{\bf v}=(v_{1},v_{2})$, $\tilde{\bf u}=(u_{1},u_{2})$ and remove the dependence of $f$ on $x_{3}$. We introduce two new unknowns

[TABLE]

which are integrations on $v_{3}$ and from assumption we notice that

[TABLE]

Multiplying (A.1) by $1$ and $v_{3}^{2}/2$, integrating on $v_{3}$, it yields a system for the new unknown vector ${\bf g}=(g_{1},g_{2})$, coupled with suitable initial conditions

[TABLE]

The new BGK system (A.3) describes a relaxation process towards the vector function $(\mathcal{M}[g]_{1},\mathcal{M}[g]_{2})$, which is a Chu reduction of the Maxwellian (A.2) from 3D in velocity to 2D in velocity and has the form

[TABLE]

where

[TABLE]

We only list the argument $\tilde{\bf v}$ in $\mathcal{M}[g]_{1}(\tilde{\bf v})$ to emphasize it is $\tilde{\bf v}$-dependent, while $(t,\tilde{\bf x})$ are omitted. The macroscopic moments $\rho$, $\tilde{\bf u}$ and $T$ of $f$ can be given in terms of ${\bf g}$ as

[TABLE]

If we further assume the model is axial symmetry only with respect to axis $x_{1}$, and $f$ only depends on $x_{1}$, we can transform (A.1) into a system of 1D problem in both $x$ and $v$. The procedure is similar, that is, we introduce

[TABLE]

multiplying (A.1) by $1$ and $(v_{2}^{2}+v_{3}^{2})/2$, integrating with respect to $v_{2}$ and $v_{3}$, it yields a system for the unknown vector ${\bf h}=(h_{1},h_{2})$, coupled with suitable initial conditions

[TABLE]

The vector $(\mathcal{M}[h]_{1},\mathcal{M}[h]_{2})$ has the form

[TABLE]

where

[TABLE]

The macroscopic moments $\rho$, $u_{1}$ and $T$ of $f$ can be given in terms of ${\bf h}$

[TABLE]

Flux relation

Now we establish the relation between the hydrodynamic equations (2.18) and the reduced BGK systems (A.4) in 1D and (A.3) in 2D in space respectively. For the hybrid discontinuous Galerkin scheme proposed in the following section, due to its compactness, we only need to define consistent numerical fluxes at the interface of two cells in the hydrodynamic region and the kinetic region respectively. Here we would like to borrow the form (2.19) and provide explicit formulae for the flux functions in (2.18) to the unknowns of the reduced BGK systems (A.4) and (A.3).

First for the 1D BGK system, the deformation tensor is simply

[TABLE]

The first order truncated distribution function $f^{\varepsilon}({\bf v})$ (2.2) becomes

[TABLE]

where ${\bf V}=(V_{1},V_{2},V_{3}):=\frac{{\bf v}-{\bf u}}{\sqrt{T}}$. Corresponding to the reduced BGK system (A.4) in 1D, we have

[TABLE]

The subindex “T” denotes truncations for the corresponding unknown functions. The nonzero first component of the flux function (2.19) becomes

[TABLE]

which will be used to define the numerical hydrodynamic flux at the interface of two cells between two regions. Here $m(v_{1})=(1,v_{1},v_{1}^{2}/2)^{T}$ and ${\bf e}_{3}=(0,0,1)^{T}$.

For the 2D reduced BGK system (A.3), the truncated distribution function $f_{T}({\bf v})$ is

[TABLE]

so that

[TABLE]

The flux function (2.19) in 2D becomes

[TABLE]

in this equation we denote $\tilde{\bf v}=(v_{1},v_{2})$ and $m(\tilde{\bf v})=(1,v_{1},v_{2},(v_{1}^{2}+v_{2}^{2})/2)^{T}$.

Discontinuous Galerkin scheme

The discontinuous Galerkin scheme (3.2)-(3.2) can be easily adapted to the reduced BGK system (A.3), which is for any $\zeta\in\mathcal{Z}_{h}^{\bf K}$

[TABLE]

where $\widetilde{v\,{\bf g}}$ is an upwind numerical flux. Then, ${\bf U}^{n+1}_{h}$ is given for any $\beta\in\mathcal{Z}_{h}^{\bf K}$

[TABLE]

where the matrix ${\bf\Phi}$ is given by

[TABLE]

Finally, for any $\alpha\in\mathcal{Z}_{h}^{\bf K}$ we have

[TABLE]

with ${\bf g}({\bf v}):=(g_{1},g_{2})^{\text{T}}$ and $\mathcal{\bf M}[{\bf g}]:=(\mathcal{M}[g]_{1},\mathcal{M}[g]_{2})^{\text{T}}$, whereas ${\bf R}({\bf v})={\bf g}({\bf v})-\Delta t\,{\bf v}\cdot\nabla_{\bf x}\,{\bf g}({\bf v})$.