Conserved discrete unified gas-kinetic scheme with unstructured discrete velocity space

TL;DR

This paper introduces an extended, unstructured discrete velocity space for the conserved discrete unified gas-kinetic scheme (C-DUGKS), significantly improving computational efficiency for multi-scale flow simulations across various regimes.

Contribution

It develops a non-isothermal C-DUGKS with unstructured velocity space, reducing computational cost and enhancing flexibility in simulating multi-scale flows.

Findings

Significant increase in computational efficiency using unstructured velocity space.

Accurate simulation results across a wide range of Knudsen numbers.

Validated effectiveness through diverse numerical tests.

Abstract

Discrete unified gas-kinetic scheme (DUGKS) is a multi-scale numerical method for flows from continuum limit to free molecular limit, and is especially suitable for the simulation of multi-scale flows, benefiting from its multi-scale property. To reduce integration error of the DUGKS and ensure the conservation property of the collision term in isothermal flow simulations, a Conserved-DUGKS (C-DUGKS) is proposed. On the other hand, both DUGKS and C-DUGKS adopt Cartesian-type discrete velocity space, in which Gaussian and Newton-Cotes numerical quadrature are used for calculating the macroscopic physical variables in low speed and high speed flows, respectively. While, the Cartesian-type discrete velocity space leads to huge computational cost and memory demand. In this paper, the isothermal C-DUGKS is extended to the non-isothermal case by adopting coupled mass and inertial energy distribution functions. Moreover, since the unstructured mesh, such as the triangular mesh in two dimensional case, is more flexible than the structured Cartesian mesh, it is introduced to the discrete velocity space of C-DUGKS, such that more discrete velocity points can be arranged in the velocity regions that enclose large number of molecules, and only a few discrete velocity points need to be arranged in the velocity regions with small amount of molecules in it. By using the unstructured discrete velocity space, the computational efficiency of C-DUGKS is significantly increased. A series of numerical tests in a wide range of Knudsen numbers, such as the Couette flow, lid-driven cavity flow, two-dimensional rarefied Riemann problem and the supersonic cylinder flows, are carried out to examine the validity and efficiency of the present method.

Tables

Table 1. Table 1: Convergence of C-DUGKS with both structured and unstructured velocity meshes (in seconds) for Courtte flow.

	Structured mesh (A)	Unstructured mesh (B)	ratio (A/B)
Discrete velocity points	6400	1882	3.4
Convergence time (Kn=0.01)	137.0	40.3	3.4
Convergence time (Kn=$\text{0.2/}\sqrt{\pi }$)	36.3	9.3	3.9
Convergence time ((Kn=$\text{2/}\sqrt{\pi }$)	28.1	9.3	3.0
Convergence time ((Kn=$\text{20/}\sqrt{\pi }$)	271.7	96.1	2.8

Table 2. Table 2: Convergence time of C-DUGKS with both structured and unstructured velocity meshes (in seconds) for Kn=0.075 case of Cavity flow.

	Structured mesh (A)	Unstructured mesh (B)	ratio (A/B)
Discrete velocity points	2500	1192	2.1
Convergence time (Kn=0.075)	677.1	235.7	2.9

Table 3. Table 3: Convergence time of C-DUGKS with both structured and unstructured velocity meshes (in seconds) for Kn=1 and 10 cases of Cavity flow.

	Structured mesh (A)	Unstructured mesh (B)	ratio (A/B)
Discrete velocity points	10201	1192	8.6
Convergence time (Kn=1)	1537.3	189.3	8.1
Convergence time (Kn=10)	2319.1	189.1	12.3

Table 4. Table 4: Run time of C-DUGKS with both structured and unstructured velocity meshes (in seconds) for Riemann problem.

	Structured mesh (A)	Unstructured mesh (B)	ratio (A/B)
Discrete velocity points	10201	2788	3.7
Run time (Kn=27.7)	642.2	175.9	3.7

Table 5. Table 5: Convergence time of C-DUGKS with both structured and unstructured velocity meshes (in seconds) for hypersonic cylinder case.

	Structured mesh (A)	Unstructured mesh (B)	ratio (A/B)
Discrete velocity points	7921	3045	2.6
Convergence time (Kn=0.1)	45664.0	17596.4	2.6
Convergence time (Kn=1)	13352.4	5312.2	2.5

Equations

$$\frac{\partial f}{\partial t}+\mathbf{\xi }\cdot \nabla f=\Omega \equiv -\frac{f-f^s}{\tau },$$

$$f^s=f^{eq}\left[1+\left(1-\mathrm{Pr}\right)\frac{\mathbf{c}\cdot \mathbf{q}}{5pRT}\left(\frac{c^2+{\eta }^2}{RT}-5\right)\right]=f^{eq}+f_{\mathrm{Pr}},$$

$$f^{eq}=\frac{\rho }{{\left(2\pi RT\right)}^{\left(3+K\right)/2}}\exp \left(-\frac{c{}^2+{\eta }^2+{\zeta }^2}{2RT}\right),$$

$${\bm{W}}=\left({\begin{array}[]{*{20}{c}}\rho\\ {\rho{\bm{u}}}\\ {\rho E}\end{array}}\right)=\int{{\bm{\psi}}\left({{{\bm{\xi}}},{{\bm{\eta}}},{{\bm{\zeta}}}}\right)f}d{{\bm{\xi}}}d{{\bm{\eta}}}d{{\bm{\zeta}}},$$

$$\gamma =\frac{c_p}{c_V}=\frac{5+K}{3+K},$$

$$\mathbf{q}=\frac{1}{2}\int \mathbf{c}\left(c^2+{\eta }^2+{\zeta }^2\right)fd\mathbf{\xi }d\mathbf{\eta }d\mathbf{\zeta },$$

$$\mathbf{\tau }=\int \mathbf{cc}\left(f-f^{eq}\right)d\mathbf{\xi }d\mathbf{\eta }d\mathbf{\zeta }.$$

$$\mu ={\mu }_{ref}{\left(\frac{T}{T_{ref}}\right)}^{\omega },$$

$$\lambda =\frac{2\mu \left(5-2\omega \right)\left(7-2\omega \right)}{15\rho {\left(2\pi RT\right)}^{1/2}}.$$

$$Kn=\frac{\lambda }{L}=\frac{2\mu \left(5-2\omega \right)\left(7-2\omega \right)}{15\rho {\left(2\pi RT\right)}^{1/2}L}=\sqrt{\frac{\gamma }{\pi }}\frac{\sqrt{2}\left(5-2\omega \right)\left(7-2\omega \right)}{15}\frac{Ma}{\mathrm{Re}},$$

$$\left\{{\begin{array}[]{*{20}{c}}{g\left({{\bm{x}},{{\bm{\xi}}},t}\right)=\int{f\left({{\bm{x}},{{\bm{\xi}}},{{\bm{\eta}}},{{\bm{\zeta}}},t}\right)d{{\bm{\eta}}}d{{\bm{\zeta}}}}},\\ {h\left({\bm{x},\bm{\xi},t}\right)=\int{\left({{{\eta}}^{2}+{{\zeta}}^{2}}\right)f\left({{\bm{x}},{{\bm{\xi}}},{{\bm{\eta}}},{{\bm{\zeta}}},t}\right)d{{\bm{\eta}}}d{{\bm{\zeta}}}}.}\end{array}}\right.$$

$$\left\{{\begin{array}[]{*{20}{c}}{\frac{{\partial g}}{{\partial t}}+{\bm{\xi}}\cdot\nabla g=\Omega_{g}\equiv-\frac{{g-g^{s}}}{\tau}},\\ {\frac{{\partial h}}{{\partial t}}+{\bm{\xi}}\cdot\nabla h=\Omega_{h}\equiv-\frac{{h-h^{s}}}{\tau}},\end{array}}\right.$$

$$\left\{{\begin{array}[]{*{20}{c}}{g^{s}=\int{f^{s}\left({{\bm{x}},{{\bm{\xi}}},{{\bm{\eta}}},{{\bm{\zeta}}},t}\right)d{{\bm{\eta}}}d{{\bm{\zeta}}}}=g^{eq}+g_{\Pr}},\\ {h^{s}=\int{\left({\eta^{2}+\zeta^{2}}\right)f^{s}\left({{\bm{x}},{{\bm{\xi}}},{{\bm{\eta}}},{{\bm{\zeta}}},t}\right)d{{\bm{\eta}}}d{{\bm{\zeta}}}}=h^{eq}+h_{\Pr}},\end{array}}\right.$$

$$g^{eq}=\int f^{eq}d\mathbf{\eta }d\mathbf{\zeta }=\frac{\rho }{{\left(2\pi RT\right)}^{D/2}}\exp \left(-\frac{c^2}{2RT}\right),$$

$$g_{\mathrm{Pr}}=\int f_{\mathrm{Pr}}d\mathbf{\eta }d\mathbf{\zeta }=\left(1-\mathrm{Pr}\right)\frac{\mathbf{c}\cdot \mathbf{q}}{5pRT}\left(\frac{c^2}{RT}-D-2\right)g^{eq},$$

$$h^{eq}=\int \left({\eta }^2+{\zeta }^2\right)f^{eq}d\mathbf{\eta }d\mathbf{\zeta }=\left(K+3-D\right)RT{g}^{eq},$$

$$h_{\mathrm{Pr}}=\int \left({\eta }^2+{\zeta }^2\right)f_{\mathrm{Pr}}d\mathbf{\eta }d\mathbf{\zeta }=\left(1-\mathrm{Pr}\right)\frac{\mathbf{c}\cdot \mathbf{q}}{5pRT}\left[\left(\frac{c^2}{RT}-D\right)\left(K+3-D\right)-2K\right]RT{g}^{eq}.$$

$$\left\{{\begin{array}[]{*{20}{c}}{\rho=\int{gd{{\bm{\xi}}}}},\\ {\rho{\bm{u}}=\int{{{\bm{\xi}}}gd{{\bm{\xi}}}}},\\ {\rho E=\frac{1}{2}\int{\left({{{\xi}}^{2}g+h}\right)d{{\bm{\xi}}}}},\end{array}}\right.$$

$$\mathbf{q}=\frac{1}{2}\int \mathbf{c}\left(c^{2}g+h\right)d\mathbf{\xi },$$

$$\mathbf{\tau }=\int \mathbf{cc}\left(g-g^{eq}\right)d\mathbf{\xi }.$$

$$\left\{{\begin{array}[]{*{20}{c}}{\int{\Omega_{g}d{{\bm{\xi}}}=0}},\\ {\int{{{\bm{\xi}}}\Omega_{g}d{{\bm{\xi}}}=0}},\\ {\int{\left({{{\xi}}^{2}\Omega_{g}+\Omega_{h}}\right)d{{\bm{\xi}}}=0}}.\end{array}}\right.$$

$$\frac{\partial \varphi }{\partial t}+\mathbf{\xi }\cdot \nabla \varphi ={\Omega }_{\varphi }\equiv -\frac{\varphi -{\varphi }^s}{\tau }.$$

$${\varphi }_j^{n+1}\left(\mathbf{\xi }\right)-{\varphi }_j^n\left(\mathbf{\xi }\right)+\frac{\Delta t}{\left|V_j\right|}{F}_j^{n+1/2}\left(\mathbf{\xi }\right)=\frac{\Delta t}{2}\left[{\Omega }_j^{n+1}\left(\mathbf{\xi }\right)+{\Omega }_j^n\left(\mathbf{\xi }\right)\right],$$

$$F_j^{n+1/2}\left(\mathbf{\xi }\right)={\int }_{\partial V_j}\left(\mathbf{\xi }\cdot \mathbf{n}\right)\varphi \left(\mathbf{x},\mathbf{\xi },t_{n+1/2}\right)d\mathbf{S},$$

$${\varphi }_j^{n+1}\left(\mathbf{\xi }\right)={\left(1+\frac{\Delta t}{2{\tau }_j^{n+1}}\right)}^{-1}[{\varphi }_j^n\left(\mathbf{\xi }\right)-\frac{\Delta t}{\left|V_j\right|}{F}_j^{n+1/2}\left(\mathbf{\xi }\right)+$$

$$\frac{\Delta t}{2}\left(\frac{{\varphi }_j^{s,n+1}\left(\mathbf{\xi }\right)}{{\tau }_j^{n+1}}+\frac{{\varphi }_j^{s,n}\left(\mathbf{\xi }\right)-{\varphi }_j^n\left(\mathbf{\xi }\right)}{{\tau }_j^n}\right)].$$

$$\displaystyle\left({\begin{array}[]{*{20}{c}}{\rho_{j}^{n+1}}\\ {\left({\rho{\bm{u}}}\right)_{j}^{n+1}}\\ {\left({\rho E}\right)_{j}^{n+1}}\end{array}}\right)$$

$$\displaystyle=\left({\begin{array}[]{*{20}{c}}{\rho_{j}^{n}}\\ {\left({\rho{\bm{u}}}\right)_{j}^{n}}\\ {\left({\rho E}\right)_{j}^{n}}\end{array}}\right)-$$

$$\varphi \left({\mathbf{x}}_b,\mathbf{\xi },t_n+h\right)-\varphi \left({\mathbf{x}}_b-\mathbf{\xi }h,\mathbf{\xi },t_n\right)=\frac{h}{2}\left[\Omega \left({\mathbf{x}}_b,\mathbf{\xi },t_n+h\right)+\Omega \left({\mathbf{x}}_b-\mathbf{\xi }h,\mathbf{\xi },t_n\right)\right],$$

$$\begin{array}{c}\stackrel{-}{\varphi }=\varphi -\frac{h}{2}\Omega =\frac{2\tau +h}{2\tau }\varphi -\frac{h}{2\tau }{\varphi }^s\text{hfill},\\ \stackrel{-}{{\varphi }^{+}}=\varphi +\frac{h}{2}\Omega =\frac{2\tau -h}{2\tau }\varphi +\frac{h}{2\tau }{\varphi }^s\text{hfill},\end{array}$$

$$\stackrel{-}{\varphi }\left({\mathbf{x}}_b,\mathbf{\xi },t_n+h\right)=\stackrel{-}{{\varphi }^{+}}\left({\mathbf{x}}_b-\mathbf{\xi }h,\mathbf{\xi },t_n\right).$$

Full text

Keywords: unstructured discrete velocity space, Conserved-DUGKS, multi-scale flow, continuum flow, rarefied flow

1 Introduction

With the fast development of hypersonic vehicles and near-space aerocraft, numerical flow solvers are required to be applicable to flows in all flow regimes from continuum limit to free molecular limit. In multi-scale simulations, both numerical methods for continuum flow, such as computational fluid dynamics (CFD) methods based on Navier-Stokes (N-S) equations, and rarefied numerical methods, such as direct simulation Monte Carlo (DSMC) [1]method for Boltzmann equation, should face the problems of model inaccuracy and large computational cost, respectively. Although the CFD methods have achieved great success in continuum flow simulations, it can not be directly extended for rarefied flows with non-equilibrium effect [1, 2]. On the other hand, due to the limitation of cell size and time step, the computational cost of the DSMC method becomes huge in continuum and near-continuum flow simulations [2].

In recent years, the numerical schemes based on gas-kinetic theory have achieved a fast development [2, 3, 4, 5, 6, 7]. As the governing equation of the gas-kinetic theory, theoretically, the Boltzmann equation can be used for flow simulations in all flow regimes [8, 9]. The Discrete velocity method (DVM) [10, 11, 12] is one of the numerical methods for Boltzmann equation and its Bhatnagar-Gross-Krook-type (BGK-type) [13] model equations, in which a discrete velocity space is used. Since most DVM schemes are single-scale methods with split particle transport and collision processes as the DSMC method, their iteration time step should be less than the mean collision time, and their cell size should be less than the mean free path [14, 15]. To address this problem, some asymptotic preserving schemes [16, 17], such as the Unified Gas-Kinetic Scheme (UGKS) [6, 18, 19, 20, 21] and Discrete Unified Gas-Kinetic Scheme (DUGKS) [2, 3, 7], have been proposed. In the algorithm of UGKS, the particle transport and collision processes are coupled by the analytical solution of the BGK-type model equation for updating the discrete distribution functions, as a result, there is no restrictions of the cell size, and the time step can be chosen according to CFL condition. DUGKS has a similar physical process as the UGKS, except the way of calculating the distribution functions at the cell interface. As a results, the numerical flux of distribution functions at the cell interface is obtained from the time integral of the analytical solution of the BGK-type equation in the UGKS algorithm, while, in the algorithm of DUGKS, this numerical flux is obtained from the distribution functions at the half-time step, which can be determined by the numerical characteristic solution of the BGK-type equation. Being the same with the UGKS, the numerical flux of the DUGKS also couples the particle transport and collision processes, while its mathematical form is simple and easy to compute [2].

Recently, Liu et al. proposed a conserved-DUGKS (C-DUGKS) for isothermal micro-channel flows in all flow regimes [22]. Like the original DUGKS, the C-DUGKS also calculates the micro-flux at the cell interface by means of the characteristic-line theory using the auxiliary distribution functions, while the discrete distribution functions and macroscopic conservative variables (mass, momentum and energy) are updated simultaneously to ensure the conservation of numerical collision term in a time implicitly form which is similar to that of the UGKS method.

Up to now, both the DUGKS and the C-DUGKS use a Cartesian discrete velocity space with either Gauss-Hermite or Newton-Cotes numerical quadrature. For hypersonic rarefied gas flow simulations, a huge number of discrete velocity points should be adopted to capture the non-equilibrium distribution functions accurately, and the corresponding computational cost and memory demand is huge.

To improve the computational efficiency of the asymptotic preserving gas-kinetic methods such as the UGKS and the DUGKS, Zhu et al. and Pan et al. proposed implicit UGKS [23] and implitct DUGKS schemes [24] for non-equilibrium steady flow, respectively, and the multi-grid technology is used to accelerate the convergence [25]. Then, the corresponding schemes are extended to unsteady flow simulations further [26]. Up to now, the implicit scheme with multi-grid technology improves the computational efficiency significantly (up to 100 times faster than the explicit ones). Besides the implicit method with multi-grid technology, the velocity space adaptation [27, 28, 29, 30, 31] is another feasible method for acceleration propose. To improve the computational efficiency and reduce the memory consumption of the UGKS method, Yuan et al. introduced a unstructured discrete velocity space in the numerical scheme [32], and the corresponding technology is proposed to suppress the integration error on unstructured mesh. Since the unstructured mesh is flexible, it can be made dense in the regions easily where the value of distribution functions is large, and be made coarse when the value of distribution functions is small. As a result, a small amount of discrete velocity points are sufficient for capturing the non-equilibrium distribution functions.

In this paper, the C-DUGKS method is extended from isothermal and incompressible one to more general non-isothermal and compressible one by adopting the coupled mass and inertial energy distribution functions and the corresponding coupled BGK-type equation system. On the other hand, the unstructured velocity space is introduced to C-DUGKS to improve the computational efficiency, and the validity of the unstructured velocity space is further examined. The Shakhov model is chosen from BGK-type model equations, since its numerical behavior is better than the other models especially for high Mach number flows.

The rest of the paper is organized as follows: the present non-isothermal and compressible C-DUGKS based on the Shakhov model introduced in Section 2 along with the unstructured discrete velocity space. In Section 3, a series of numerical tests in different flow regimes, including the Couette flow, lid-driven cavity flow, two-dimensional rarefied Riemann problem and the supersonic cylinder flow, are used to examine the validity, accuracy and computational efficiency of the present method. The discussion and conclusion are in Section 4.

2 The conserved discrete unified gas-kinetic scheme

2.1 Shakhov model for Boltzmann equation

The starting point of the DUGKS is the BGK-type Boltzmann model equations that use simple relaxation operators instead of complicate Boltzmann collision term. The BGK equation corresponds to a fixed unity Prandtl (Pr) number [3]. To get the right Pr number, several modified models are proposed based on different physical considerations, such as the Shakhov model [33] and the ellipsoidal statistical model (ES model) [34] . The D-dimensional (D is from 1 to 3) Shakhov model equation used in this paper is in the form of

[TABLE]

where $f=f\left({{\bm{x}},{{\bm{\xi}}},{{\bm{\eta}}},{{\bm{\zeta}}},t}\right)$ is the velocity distribution functions for particles moving in D-dimensional physical space with velocity ${{\bm{\xi}}}=\left({\xi_{1},...,\xi_{D}}\right)$ at position ${\bm{x}}=\left({x_{1},...,x_{D}}\right)$ and time $t$. Here ${{\bm{\eta}}}=\left({\xi_{D+1},{\text{ }}.{\text{ }}.{\text{ }}.{\text{ }},\xi_{3}}\right)$ is the dummy velocity (with the degree of freedom $L=3-D$) consisting of the rest components of the translational velocity of particles in three-dimensional space; ${{\bm{\zeta}}}$ is a vector of $K$ elements representing the internal degree of freedom of molecules; $\Omega$ is the collision operator, $\tau$=$\mu$/$p$, $\mu$ is dynamic viscous coefficient, $p$ is static pressure, and $f^{s}$ is the Shakhov equilibrium distribution functions in the following form:

[TABLE]

where $f^{eq}$ is the Maxwellian distribution functions, ${\bm{c}}={{\bm{\xi}}}-{\bm{u}}$ is the peculiar velocity, ${\bm{u}}$ is the macroscopic flow velocity, ${{\bm{q}}}$ is the heat flux, $R$ is the specific gas constant, and $T$ is the temperature. The Maxwellian distribution functions is given by

[TABLE]

where $\rho$ is the density. In the gas-kinetic theory, the conservative flow variables are defined by the moments of the distribution functions as follows:

[TABLE]

where $\bm{\psi}=\left({{\text{1, }}{{\bm{\xi}}},\left({{\xi}^{2}+{\eta}^{2}+{\zeta}^{2}}\right)/2}\right)^{T}$ , $\rho E=\rho u^{2}/2+\rho\varepsilon$ is the total energy density, $\rho\varepsilon=\rho c_{V}T$ is the inertial energy density, $c_{V}=\left({3+K}\right)R/2$ is the specific heat capacity at constant volume. The specific heat ratio is

[TABLE]

where $c_{p}=\left({5+K}\right)R/2$ is the specific heat capacity at constant pressure. In the gas-kinetic theory, the heat flux is defined by

[TABLE]

and the stress tensor is defined by

[TABLE]

The expression of viscous coefficient $\mu$ is related to the law of intermolecular interactions [1, 3, 35]. For variable hard-sphere (VHS) molecules, the viscous coefficient is

[TABLE]

here $\omega$ is the heat index, which is 0.5 for HS model and 0.81 for Argon gas. The relation between the mean free path $\lambda$ and $\mu$ (VHS model) can be written as the follows:

[TABLE]

In the non-dimensional system, the Knudsen (Kn) number is defined as

[TABLE]

where $Ma=\left.{u}/\sqrt{\gamma RT}\right.$ and $\operatorname{Re}=\left.\rho{u}L/\mu\right.$ are the Mach number and Reynolds number, respectively, and L is the reference length.

2.2 Reduced distributions

The transport process of molecules depends on the D-dimensional particle velocity ${\bm{\xi}}$ only, and is irrelevant to ${\bm{\eta}}$ and ${\bm{\zeta}}$. To avoid the discretization of ${\bm{\eta}}$ and ${\bm{\zeta}}$, two reduced distribution functions [36] are adopted in the present numerical scheme as follows:

[TABLE]

Given Eq.(1) and Eq.(11), the following BGK-type equation system can be obtained,

[TABLE]

where the reduced Shakhov equilibrium distribution functions $g^{s}$ and $h^{s}$ are given by

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

By substituting Eq.(11) into Eq.(4), (6), (7), it can be obtained that

[TABLE]

[TABLE]

[TABLE]

Using the reduced distribution functions, the conservation property of Shakhov collision operator can also be written as:

[TABLE]

2.3 The Conserved DUGKS

As can be seen from Eq.(12), $g$ and $h$ have the same updating rule, therefor they can be replaced by a new symbol $\phi$ in the algorithm for simplicity, and the governing equation can be written as:

[TABLE]

Integrating Eq.(22) in control volume $V_{j}$ from time $t_{n}$ to $t_{n+1}{\text{ = }}t_{n}{\text{ + }}\Delta t$, it can obtain that

[TABLE]

where $\left|V_{j}\right|$ is the volume of $V_{j}$ and $F^{n+1/2}\left({{\bm{\xi}}}\right)$ is the microflux across the cell interface given by

[TABLE]

where $\partial V_{j}$ is the surface of the control volume $V_{j}$ and ${\bm{n}}$ is the external normal unit vector of $d\bm{S}$ (differential of $V_{j}$). It should be noted that a trapezoidal rule is used for the time discretization of the collision term.

In order to avoid calculating the collision term at $t_{n+1}$ time on the right side of Eq.(23), DUGKS introduces two new auxiliary distribution functions. While, the conservative variables of C-DUGKS are updated first [22], so that the equilibrium distribution functions at $t_{n+1}$ can be obtained, then the implicit evolution equation can be transformed into the following explicit one,

[TABLE]

Similar to the UGKS, conservative variables ${\bm{W}_{j}}^{n+1}$ need to be updated before calculating the equilibrium distribution functions $\phi_{j}^{s,n+1}\left({\bm{\xi}}\right)$ and the relaxation time $\tau_{j}^{n+1}$ in Eq.(25). Given Eq.(18), Eq.(21) and Eq.(23), ${\bm{W}_{j}}^{n+1}$ can be updated using the following equations:

[TABLE]

Eq.(25) and Eq.(26) are the evolution equations for the microscopic distribution functions and the macroscopic conservative variables, respectively. Once the distribution functions $\phi\left(\bm{x},{\bm{\xi}},t_{n+1/2}\right)$ at the half-time step and at the cell interface is obtained, the whole scheme can be established.

In order to obtain $\phi\left(\bm{x},\bm{\xi},t_{n+1/2}\right)$ , the Shakhov model equation Eq.(22) is integrated along the characteristic-line (in the direction of particle velocity) from $t_{n}$ to $t_{n+1/2}$. As shown in Fig. 1, the characteristic ends at the midpoint of the cell interface. This process can be expressed by the following equation:

[TABLE]

where ${\bm{x}}_{b}$ is the midpoint of the cell interface, and $h=\Delta t/2$ is the half-time step. By introducing the following auxiliary distribution functions:

[TABLE]

Eq.(27) can be written as

[TABLE]

That means as long as $\mathop{\phi^{+}}\limits^{-}\left({{\bm{x}}_{b}-{\bm{\xi}}h,{\bm{\xi}},t_{n}}\right)$ is known, $\mathop{\phi}\limits^{-}\left({{\bm{x}}_{b},{\bm{\xi}},t_{n}+h}\right)$ can be directly abtained. Since the collision operator fulfills the conservation property (Eq.(21)), the macroscopic variables at the cell interface at half-time step can be obtained from the auxiliary distribution functions as follows:

[TABLE]

[TABLE]

At the half-time step, the Shakhov equilibrium distribution functions $\phi^{s}\left({{\bm{x}}_{b},{\bm{\xi}},t_{n}+h}\right)$ at the cell interface can be obtained by calculating the macroscopic conservative variables ${\bm{W}}\left({{\bm{x}}_{b},{\bm{\xi}},t_{n}+h}\right)$ and the heat flux ${\bm{q}}\left({{\bm{x}}_{b},{\bm{\xi}},t_{n}+h}\right)$ at the cell interface. Thus, the distribution functions at the cell interface can be recovered through the following equation:

[TABLE]

Since both distribution functions and macroscopic variables at the cell interface at half-step are constructed from $\mathop{\phi}\limits^{-}\left({{\bm{x}}_{b},{\bm{\xi}},t_{n}+h}\right)$, and $\mathop{\phi}\limits^{-}\left({{\bm{x}}_{b},{\bm{\xi}},t_{n}+h}\right)$ is equal to $\mathop{\phi^{+}}\limits^{-}\left({{\bm{x}}_{b}-{\bm{\xi}}h,{\bm{\xi}},t_{n}}\right)$ (Eq.(29)). Therefore, the calculation of $\mathop{\phi^{+}}\limits^{-}\left({{\bm{x}}_{b}-{\bm{\xi}}h,{\bm{\xi}},t_{n}}\right)$ is very important in the present scheme. By using the Taylor expansion at the center of the control volume, $\mathop{\phi^{+}}\limits^{-}\left({{\bm{x}}_{b}-{\bm{\xi}}h,{\bm{\xi}},t_{n}}\right)$ can be obtained from the following reconstruction:

[TABLE]

where $V_{c}$ stands for control volume which is centered at point $C$ (Fig. 1). If ${\bm{\xi}}\cdot{\bm{n}}_{b}\geqslant 0$, point $C$ is $P$ (the center the left cell) in Fig. 1, otherwise point $C$ is $Q$ (the center of the right cell). $\nabla\mathop{\phi^{+}}\left({{\bm{x}}_{c},{\bm{\xi}},t_{n}}\right)$ is the gradient of the auxiliary distribution functions at point $C$, which is calculated by the least square method in this paper, and $\Psi\left({{\bm{x}}_{c},{\bm{\xi}},t_{n}}\right)$ is the gradient limiter used to suppress the numerical oscillations. The Venkatakrishnan limiter [37] is chosen in this paper.

The time step of the present explicit scheme is determined by the following Courant-Friedrichs-Lewy (CFL) condition,

[TABLE]

where $\alpha$ is the CFL number, and $\Delta x$ is the minimum grid spacing.

Fig. 2 shows the flow chart of the present compressible and non-isothermal C-DUGKS. The key point of the present scheme is calculating the distribution functions at the cell interface at the half-time step.

2.4 Unstructured discrete velocity space

By integrating the distribution functions in the continuous velocity space, the conservative variables can be obtained from Eq.(18). In numerical simulations, when discrete velocity space is adopted, the continuous integration in Eq.(18) is replaced by the following numerical quadrature:

[TABLE]

Up to now, the DUKGS generally uses Gauss-Hermite numerical quadrature for low-speed continuum flows and Newton-Cotes numerical quadrature for high-speed rarefied flows. These two numerical quadrature are conducted on structured Cartesian mesh in velocity space. Therefore, when the flow Mach number is large, the amount of discrete velocity points need to be extremely large. Actually, according to Eq.(35), the velocity mesh should be dense only in the regions where the value of distribution function is large, and can be coarse in the regions with small distribution function.

To reduce the amount of discrete velocity points and improve the computational efficiency, Yuan et al. introduced the unstructured discrete velocity into the UGKS method, since the unstructured mesh (the triangular mesh for example) is more flexible than the structured one [32]. Fig. 3 illustrates a simple comparison between the structured Cartesian mesh and unstructured triangular mesh for discrete velocity space about a Maxwellian distribution with zero macroscopic speed. The triangular mesh can be dense at the zero velocity point where the distribution functions is the maximum, and gradually become coarse when the molecular speed is large (far from the zero velocity point in velocity space). Obviously, the unstructured discrete velocity space is flexible. It can refine and coarsen the mesh when needed. This makes the unstructured mesh a efficient way in organizing the discrete velocity space.

In the unstructured discrete velocity space, the discrete points and their weights are the central point and the area of mesh cells, respectively. Theoretically, the accuracy of this type of integration is slightly lower than that of Gauss-Hermite and Newton-Cotes ones on structured Cartesian mesh. While, since the unstructured mesh is flexible, its integration accuracy is often higher than that of the structured one with the same amount of discrete points. In this paper, the unstructured discrete velocity space is introduced to the C-DUGKS for the first time to improve the computational efficiency and reduce memory cost. Since DUGKS has a similar physical process as UGKS, it is easy to believe that unstructured discrete velocity space will work well in C-DUGKS. Several examples will be carried out to explore the performance of present method in the next section.

3 Numerical experiments

A series of test cases are carried out using the present C-DUGKS with both structured Cartesian discrete velocity space and unstructured triangular discrete velocity space. The accuracy and numerical efficiency of the two methods are compared. The results of UGKS and DSMC are chosen as the benchmark solutions.

The first test case is the Couette flow with Kn number $0.01$, ${\text{0}}{\text{.2/}}\sqrt{\pi}$, ${\text{2/}}\sqrt{\pi}$ and ${\text{20/}}\sqrt{\pi}$. The second test case is the lid-driven cavity flow with Kn number $0.075$, $1$ and $10$. The third test case is the two-dimensional rarefied Riemann problem. The last test case is the supersonic cylinder flow with Kn number $0.1$ and $1$. For simulation of steady flow in this paper, the relative errors of macroscopic variables are calculated every 1000 steps. When the relative error is less than $\varepsilon$ as the following,

[TABLE]

it is considered that the flow has reached the stable state. In this paper, $\varepsilon$ is chosen as ${10^{-6}}$ for all steady test cases. All numerical simulations are conducted on a compute node with 24 cores, and the CPU information is Intel(R) Xeon(R) CPU E5-2685 v3 @ 2.6GHz.

3.1 Couette flow

The Couette flow is driven by two parallel plates whose distance is $H$ as shown in Fig. 4. The upper and lower plates move in the opposite directions with a magnitude $U_{w}$. In this case, the physical space between the two plates is discretized into 100 cells. The working gas is Argon ($R=208J/kg/K$), and HS model is adopted ($\omega=0.5$). Taking $H$ as the reference length, four cases with Kn number $0.01$, ${\text{0}}{\text{.2/}}\sqrt{\pi}$, ${\text{2/}}\sqrt{\pi}$ and ${\text{20/}}\sqrt{\pi}$ are carried out. The flow field is in stationary initially. Both the initial temperature and wall temperatures are $273K$. The velocity of plates are $U_{w}=0.5\sqrt{RT}=119.15m/s$. For simplicity, dimensionless physical variables are used. The reference temperature and velocity are $T_{ref}=273K$ and $U_{ref}=\sqrt{2RT_{ref}}=337.00m/s$, respectively. The dimensionless density and temperature of the initial flow are $\rho_{0}=1.0$ and $T_{0}=1.0$ respectively, and the magnitude of the dimensionless velocity of both plates are $U=0.3535$.

Similar to Ref. [23], $80\times 80$ uniform discrete velocity points in a range of $\left[{-4,\;4}\right]\times\left[{-4,\;4}\right]$ are chosen as the structured discrete velocity space. And the unstructured discrete velocity space shown in Fig. 5 is adopted. The range of this unstructured velocity space is a circle surface that is centered at $\left({0,0}\right)$ with a radius of 4. The number of discrete velocity points is $1,681$. Although different types of discrete velocity spaces are chosen, they can use the same CFL number, since their minimum spatial mesh size and maximum discrete velocity are the same.

Fig. 6 shows the horizontal velocity profile in the vertical direction predicted by C-DUGKS with both structured and unstructured velocity spaces, along with the benchmark UGKS data. Since the vertical velocity is almost zero for Couette flow, the discrete velocity spaces in the vertical direction is set coarse, and it is refined in the horizontal direction. It is found that the results predicted by the UGKS and the present C-DUGKS are quite consistent with each other. Since the unstructured mesh is flexible, 1882 discrete velocity points are sufficient in this case, while the number of discrete velocity points are 6400 (80*80) for the structured velocity space. Table 1 is the convergence time of C-DUGKS with both structured and unstructured velocity space. Since the computation time is proportional to the number of discrete velocity points, the utilization of the unstructured velocity space leads to a 3 to 4 times speedup.

3.2 Cavity flow

Two-dimensional lid-driven cavity flow is one of the classic cases to test numerical scheme, especially their performance on viscous effect. Here, cavity flow cases with wide range of Kn numbers (0.075, 1 and 10) are carried out. The length of the four walls of the cavity (shown in Fig. 7) are L. The top wall is moving from left to right with a constant velocity $U_{w}$, and the other walls are static. The working gas is Argon, and the HS molecular model is chosen. The temperatures of all walls are 273K. Initially, the flow field is in stationary with a temperature of 273K. The sonic speed and Mach number of the moving wall are $a=\sqrt{\gamma RT}=307.64m/s$ and $M\text{a}=0.16$, respectively. In numerical simulation, $T_{ref}=227.50K$ and $U_{ref}=\sqrt{2RT_{ref}}=307.64m/s$ are set as the reference temperature and velocity, respectively. Therefore, the dimensionless density and temperature of the initial flow field are $\rho_{0}=1.0$ and $T_{0}=1.2$ respectively. The dimensionless velocity of the top wall is $U_{w}=0.16$. The present C-DUGKS with both structured and unstructured velocity space are used for simulation, and the UGKS data is used as the benchmark solution. The spatial mesh for cavity (2460 cells) is shown in Fig. 8(a). In order to capture the boundary layer, a body-fitted rectangular mesh is used near the wall, while unstructured triangular mesh is used inside the cavity. For unstructured discrete velocity spaces, the velocity mesh (1192 points) in Fig. 8(b) is chosen for all Kn numbers, which is centered at $\left({0,0}\right)$ with a radius of 4. For structured discrete velocity spaces, a structured and uniform mesh of 50x50 points is chosen for Kn=0.075, while the mesh with 101x101 points is chosen for Kn=1 and 10. The ranges of both structured velocity space are $\left[{-4,\;4}\right]\times\left[{-4,\;4}\right]$.

The results of cavity flow at Kn = 0.075 and 10 are presented in Figs. 9 and 10, where the temperature contours, the streamline of heat flux, and the velocities along the central horizontal and central vertical lines are illustrated. It can be seen that the results from C-DUGKS with both structured and unstructured velocity space are in good agreement with the UGKS. In addition, a close observation on Fig. 10(a) shows that the result of C-DUGKS with the structured mesh (10201 discrete velocity points) has a nonphysical ray effect at the bottom-left corner of the cavity, i.e. the temperature contour is jittery and discontinuous. On the contrary, the result of C-DUGKS with the unstructured velocity mesh (only 1192 discrete velocity points) is consistent with the benchmark solution with no ray effect. Since the speed of the top moving wall is only 0.16 (near the origin of the velocity space), and the temperature is around 1.0 in the whole flow field, then particles are concentrated near the origin of the velocity space. Since the unstructured mesh is flexible and can arrange sufficient discrete velocity points near the velocity origin, then the distribution functions is well captured, and the ray effect is avoided. Moreover, the unstructured velocity mesh can be set coarse in the regions far from the velocity origin, then the total number of discrete velocity points is still much less than that of the structured and uniform velocity mesh. Table 2 and table 3 show the convergence time taken by C-DUGKS with both structured and unstructured mesh. It is found that the convergence time of C-DUGKS with the unstructured mesh is much less than that of C-DUGKS with the structured velocity space (with a speedup of 8 to 12), benefiting from use of the unstructured velocity mesh.

3.3 Two-dimensional rarefied Riemann problem

In this section, the two-dimensional rarefied Riemann problem is used to test the validity of the present C-DUGKS method. And the computational domain is divided into quadrants with different macroscopic variables as follows:

[TABLE]

The working gas is composed by diatomic molecule (HS model) with $\gamma=1.4$. The viscous coefficient is $\mu_{0}=10$. Taking the physical variables in the first quadrant as the reference values, then Kn is 27.7. The computational domain $\left[-{{\text{ 0.5 }},\;0.5}\right]\times\left[-{{\text{ 0.5 }},\;0.5}\right]$ is discretized into a Cartesian mesh with 60x60 cells. The unstructured discrete velocity spaces used in C-DUGKS (2788 discrete velocity points) is shown in Fig. 11, which is centered at $\left({0,0}\right)$ with a radius of 7. While the structured discrete velocity spaces used in C-DUGKS has 101x101 discrete velocity points in the range of $\left[{-{\text{7}},\;{\text{7}}}\right]\times\left[{-{\text{7}},\;{\text{7}}}\right]$.

The density, temperature, velocity magnitude $\sqrt{u_{x}^{2}+u_{y}^{2}}$, and streamlines at t = 0.15 are shown in Fig. 12. The results of C-DUGKS with both structured and unstructured mesh match well with the results from the solution of the collisionless Boltzmann equation. Table 4 is number of discrete velocity points and the computation time of the C-DUGKS with both structured and unstructured mesh, respectively. It is found that the speed-up ratio is the same with the ratio of the mesh numbers (the number of structured velocity mesh over the number of unstructured velocity mesh). Although the velocity space is unstructured, the algorithm of C-DUGKS only use the information of the coordinates of the discrete velocity points (the centers of the mesh cells) and their weights (the areas (2-D) of the mesh cells), other information, such as the face-node connection, face-cell connection, etc. are not used. As a results, the efficiencies of structured and unstructured mesh are almost the same if used as the discrete velocity space. Moreover, since the unstructured mesh is more flexible, the number of discrete velocity points can be reduced significantly, leading to a large speed-up ratio.

3.4 Supersonic flow over a cylinder

To examine the accuracy and efficiency of the present C-DUGKS method for predicting high-speed non-equilibrium flows, the supersonic flow over a circular cylinder are conducted. The working gas is Argon for which the VHS molecular model is chosen. The velocity and temperature of the free stream are 1538.18m/s and 273K, respectively, yielding a free stream Mach number equal to 5. The temperature on the wall of cylinder is a constant at 273K. The above setting is identical to that in Ref. [23]. The radius $R$ of the cylinder is taken as the reference length. Two Kn numbers 0.1 and 1.0 in transition regime are chosen for simulation. The computational domain of physical space, which is a circle with a radius of 15R, is discretized using a mesh of 64x61 cells. In order to capture the boundary layer correctly, the mesh near the cylinder wall is refined in the normal direction. The height of the first layer near the solid boundary is set 0.005 and 0.01 for Kn=0.1 and 1, respectively, since the thickness of boundary layer decreases as Kn number decreases. During the simulation, the reference temperature and velocity are $T_{ref}=273K$ and $U_{ref}=\sqrt{2RT_{ref}}=337.00m/s$, respectively. Therefore, the dimensionless free-stream temperature and the wall temperature are both 1.0, the free-stream density is 1.0 and free-stream velocity is 4.56. The dimensionless density and temperature of the initial flow field are $\rho_{0}=1.0$ and $T_{0}=1.0$ respectively. The unstructured velocity mesh (3045 points) used in C-DUGKS is shown in Fig. 13, which is centered at $\left({0,\;0}\right)$ with a radius of 10. Note that, since the free-stream velocity is 4.56, the unstructured velocity mesh is refined near both the origin $\left({0,\;0}\right)$ and the point $\left({4.56,\;0}\right)$. For structured velocity mesh, 89x89 (7921) points is used, and the range of discrete velocity is $\left[{-{\text{10}},\;{\text{10}}}\right]\times\left[{-{\text{10}},\;{\text{10}}}\right]$.

The flow variables along the stagnation line in front of the cylinder in Kn = 0.1 case are illustrated in Fig. 14. The results of C-DUGKS with both structured and unstructured mesh match well the UGKS and DSMC results in Ref. [21]. In Fig. 14(c), the temperature profile predicted by both the UGKS and the C-DUGKS deviate slightly from the DSMC result in the front of the bow shock. This deviation is due to the chosen of Shakhov model equation [2], and the same phenomenon can be found in the profiles of the normal shock structures [38]. Fig. 15 are the profiles of the heat flux, pressure and shear stress on the surface of a cylinder. It is found that the results of C-DUGKS match well with those of UGKS and DSMC. For Kn=1 case, the flow variables along the stagnation line in front of the cylinder are illustrated in Fig. 16; the heat flux, pressure and shear stress on the surface of a cylinder are illustrated in Fig. 17. Being the same with the Kn=0.1 case, the C-DUGKS results for Kn=1 case match well with those from the UGKS and DSMC methods. As the increasing of Kn number, the non-equilibrium effect becomes stronger, and the deviations of C-DUGKS and UGKS temperature profiles from the DSMC one becomes more obvious. Table 5 is the convergence time of C-DUGKS with both structured and unstructured mesh. It is found that the speed-up ratios benefiting from the unstructured velocity space are 2.6 and 2.5 for Kn 0.1 and 1, respectively, which are almost the same as the ratio of discrete velocity points.

4 Conclusion

In this paper, the coupled mass and energy distributions are introduced into the isothermal C-DUGKS, expanding its scope of simulation to non-isothermal and compressible flows in all flow regimes. Moreover, to relieve the dimensional crisis adjoint to discrete velocity space of the C-DUGKS method, a unstructured velocity space with triangular elements, which is developed for the UGKS method earlier, are introduced in this paper. The validity of the present method is proved by a series of test cases with a wide range of Kn numbers, including the Couette flow, lid-driven cavity flow, two-dimensional rarefied Rieman problem, and hypersonic flow over a circular cylinder. Especially, the validity of the present method for compressible and non-isothermal flows are verified by its precise prediction of two-dimensional rarefied Riemann problem and hypersonic cylinder flows. One important advantage of the structured uniform mesh is that high order integration methods can be used, leading to less discrete points. While, since the present unstructured velocity mesh can adjust the mesh density wherever the physical problems needs, the number discrete velocity used can be further reduced, and is significantly less than that of the structured mesh with high order integration methods. This improvement leads to a large speed-up in the two-dimensional cases of this paper. For high Mach number, this speed-up ratio will further increase, since the unstructured mesh can be made dense easily for high speed molecules concentrated in a small region of velocity space (such as hypersonic inflow).

Acknowledgements

The authors thank Prof. Kun Xu in Hong Kong University of Science and Technology and Porf. Zhaoli Guo in Huazhong University of Science and Technology for discussions of the UGKS, the DUGKS and multi-scale flow simulations. Jianfeng Chen thanks Dr. Yong Wang in Northwestern Polytechnical University for his help in developing C-DUGKS solver on Code-Satrune. Jianfeng Chen thanks Dr. Ruifeng Yuan in Northwestern Polytechnical University for useful discussions on the unstructured velocity space.

This work has been financially supported by National Natural Science Foundation of China (No. 11472219), 111 project of China (No. B17037) and National Key Laboratory Foundation (No. G2018KC0118D).