1 Introduction

Accurate physical models and efficient numerical methods for the solution to the Boltzmann equation are required for predicting the non-equilibrium transport phenomena of the rarefied gas flows encountered in high-altitude aerodynamics, vacuum technology, and micro-electromechanical systems (MEMS). Deterministic solvers of the Boltzmann equation^[1] have been widely developed due to their efficiency in resolving near continuum flows, micro-scale flows, and unsteady flows^[2], as well as their directness in coupling the continuum computational fluid dynamics (CFD) and deterministic structural, thermal, and electrostatic modelling^[3]. The deterministic approaches are attributed to the discrete ordinate method (DOM)^[4], which also corresponds to the discrete velocity method, in which the control equations are firstly discretized in the velocity space and cast into a system of partial differential equations depending on physical coordinates and time only. Then, the discrete equations are solved using the traditional CFD methods. Recently, this method is attractive for the simulation of a flow surrounding a re-entry flying object at the hypersonic speed in transitional regimes^[5].

The finite difference method (FDM) and finite volume method (FVM) have been widely developed for the deterministic solutions to the Boltzmann equation or kinetic models^[6-10]. The Runge-Kutta discontinuous Galerkin (RKDG) method is another high-order discrete formalism based on the finite element method (FEM)^[11]. This method is well suited for the solutions to the time-dependent hyperbolic equations. One of the most attractive features of this method is its ability to naturally obtain fluxes both in the interior and at the boundaries of a domain with the same high-order accuracy^[12-13]. This is especially important for applications of nonequilibrium flow analysis where stress and heat flux distributions on the solid boundaries are the main quantities of interest. Other advantages include easy formulations on irregular meshes, straightforward implementation of boundary conditions, as well as efficient parallel computation due to compactness of the scheme^[14]. This method has been successfully used to model a wide range of physical phenomena such as fluid dynamics, electromagnetics, aeroacoustics, and advection of contaminants. Recently, the RKDG method is used to solve the Boltzmann equations for the simulation of rarefied gaseous flow^[15-20].

The RKDG method resolves solutions in the piecewise finite element space combining with an explicit-multistep Runge-Kutta time marching. Borrowing from the FVM and FDM, the numerical fluxes are evaluated at edges of the physical elements. This approach is a shockcapture scheme, which may automatically capture discontinuities in solutions even if the initial conditions are smooth. It has been proven that this method is energy stable so that it could directly be used to solve control equations with smooth solutions or solutions with weak discontinuities^[21]. However, when strong shock waves appear in the flow field, the approximate solutions may exhibit spurious oscillations, which generate non-physical data and lead to instability. Usually, some forms of limiters are used to deal with the non-physical oscillations in the presence of strong discontinuities. Treated as a post processor of the approximate solutions, a limiter uses a new polynomial with the same order of degree to replace the old one which is deemed to contain oscillations. Several limiters have been designed, such as the total variation diminishing (TVD) limiters^[22], the total variation bounded (TVB) limiters^[23], the moment-based limiters^[24-25], the weighted essentially non-oscillatory (WENO) limiters^[26-27], and recently the bound-preserving limiters^[28-29].

To achieve the non-oscillatory property for strong shocks is one of the bottlenecks to the development of the discontinuous Galerkin (DG) based method. To our best knowledge, none of the previous research works about the RKDG solvers of the Boltzmann equations emphasized on the stability of the scheme and were specially designed for hypersonic rarefied gaseous flows. In this work, we develop a stable deterministic solver for the Boltzmann kinetic model equations. Three essences determine the stability as follows: (ⅰ) the accuracy of the velocity discretization, (ⅱ) the conservation of the collision terms, (ⅲ) a limiter. Based on the maximum-principlesatisfying scheme for scalar conservation laws^[29], a positivity-preserving limiter is designed for the kinetic models with the source term. The remainder of the paper is organized as follows. In Section 2, the basic numerical method is described with emphasis on the velocity discretization and conservative calculation of the collision relaxation term. In Section 3, the positivity-preserving limiter is illustrated in detail. The normal shock with the Mach number 2.05, the hypersonic flow about a two-dimensional (2D) cylinder with Mach numbers 6.0 and 12.0, and an unsteady shock tube flow are used to verify the proposed solver in Section 4. Finally, the paper is closed with some conclusions in Section 5.

2 Kinetic models and basic RKDG method 2.1 Boltzmann kinetic model equations

Due to the presence of the collision integral term, the Boltzmann equation is difficult to solve with analytical or deterministic methods, and therefore the kinetic model equation is usually used as an alternative control system. In the kinetic model, the original Boltzmann equation is reduced to a relatively simpler partial differential equation by introducing the collision relaxation term, which is much easier but can reproduce the important kinetic properties of the gaseous flows. The two frequently used forms, i.e., the Bhatnagar-Gross-Krook (BGK)^[30] and the ellipsoidal statistical BGK (ES-BGK)^[31] models, are used. The kinetic model equation is written as

(1)

where f is the distribution function of molecular velocities, which is defined so that f(t, r, c) drdc is the number of molecules at the time instant t, with velocity components lying within the limits c and c + dc, and spatial coordinates lying within the limits r and r + dr. f_E denotes a suitable distribution under an equilibrium. In the standard BGK model which corresponds to the gas Prandtl number Pr=1 for monatomic gases, f_E is the local Maxwell distribution, while in the ES-BGK model, f_E is the local anisotropic Gaussian. υ is the collision frequency given as

(2)

where p is the local gas pressure, and μ is the viscosity coefficient. An arbitrary viscosity law could be specified, and a power-law coefficient μ=μ_ref(T/T_ref)^ω is used with the reference properties μ_ref, T_ref, and ω given by Bird^[32]. The Prandtl number is a free parameter allowing the ES-BGK collision model to reproduce both the viscosity and thermal conductivity corresponding to an arbitrary Pr^[33]. Therefore, the BGK model is a special case of the ES-BGK for Pr=1.

In order to reduce the round-off error during the numerical procedure, all the variables and functions need to be non-dimensionalized. Let L, n₀, and T₀ be reference dimensional values of the length, number density, and temperature. Then, the reference velocity is , the reference time is t₀=L/u₀, the reference pressure is p₀=n₀kT₀, and the reference viscosity is n₀kT₀L/u₀, where R and k are the gas constant and the Boltzmann constant, respectively. Finally, the velocity distribution functions including the local equilibrium are non-dimensionalized by n₀/(2RT₀)^3/2. The non-dimensional macroscopic gas flow parameters are obtained by integrating the distribution functions over the velocity phase,

(3)

In the remainder of the paper, the tilde symbols denoting the non-dimensional variables will be omitted.

2.2 Discrete velocity model

The velocity distribution function is continuous in the velocity and physics spaces. The DOM method chooses a set of values of the velocity and interpolates the distribution function in terms of its values corresponding to the discrete velocities, i.e., . Therefore, the kinetic model equation is replaced by a system of differential equations for the function f^j(t, r), which is dependent on physical ordinates and time only. The differential equations can be solved numerically by a certain CFD method. Then, the macroscopic parameters calculated from the exact integration over the velocity space are approximated from some quadrature. Both the Cartesian and spherical meshes have been developed for the velocity discretization^[20]. The Cartesian mesh adopts the composite Newton-Cotes formulas, thus it has finite limits which must be chosen carefully to ensure that the effects of the tails of the distribution function are negligible, while the spherical mesh involves the Gaussian-Laguerre quadrature and has no bound limitation. The accuracy of the velocity discretization highly depends on the localization of the distributions. Figure 1 shows the non-dimensionalized x-velocity distribution functions of upstream and downstream for normal shock waves with the Mach numbers 1.55 and 9.0. The values are obtained according to an argon gas with the upstream density 1.067 ×10^-4 kg/m³ and the temperature 300 K. Therefore, for hypersonic flows, as the Mach number increases, the appropriate velocity range could become quite wide, and a large number of discrete velocities are needed to meet the accuracy requirement. Besides, the distribution function could change sharply in a small limit of the molecular velocity, which infers that the discrete velocities might become very dense. Although the GaussianLaguerre quadrature gives an optimal approximation, the quadrature points and weights are very difficult to calculate as the number of quadrature points increases. Thus, the composite Newton-Cotes formulae are more convenient^[6]. In this work, the second-order mid-point rule and the fourth-order Simpson's rule, i.e., the three-point quadrature^[34] is used.

2.3 DG formulation and conservative discretization of collision term

The control system which is discrete in the velocity space will be solved using the RKDG method on the 2D unstructured mesh. The 2D computational domain is partitioned with triangulations {Δ_i}_i=1^M, and the approximate solutions f^j(t, r) are sought in the finite element space of piecewise polynomials P^k(Δ_i) of degree at most k defined on Δ_i. Then, the method is uniformly (k + 1) th-order accurate^[11]. In the DG method, the solution as well as the test function space is given by {φ(x, y)|_Δi∈ P^k(Δ_i)}. A local orthogonal basis^[27] is used here, and the approximate solutions are expressed as

(4)

where F_l^j(t) is the degree of freedom, and the first one F₀^j is the average value on the local triangular cell. In order to determine the approximate solution, the standard technique of the finite element method is used to obtain the weak formulations,

(5)

where w^l is the diagonal component of the mass matrix, and h_e^j is the numerical flux. Since the values of f^j at the edges of the triangles calculated from the adjacent cells are not required equally, the fluxes borrowed from the FVM and FDM are calculated based on Riemann solvers. In this work, the simplest upwind scheme is used, and we write it as

(6)

Here, n_e is the outward unit normal to the edge, f_{e, int}^j is the approximate solution obtained from the interior of the triangle Δ_i, and f_{e, ext}^j is the one obtained from the exterior of Δ_i. It is easy to show that the upwind flux satisfies the conservativity and consistency,

(7)

Moreover, five different types of boundary conditions, including the symmetric boundary, the specular-diffuse moving wall with the given accommodation coefficient, the periodic boundaries, the far-pressure inlet/outlet boundaries, and supersonic inlet/outlet boundaries, are incorporated to specify the boundary values of f_ext^[20].

The resulting system of ordinary differential equations, i.e., Eq. (5), is discretized in time by a special class of explicit TVD Runge-Kutta schemes^[35]. For the DG discretization in the P^k space, the Runge-Kutta time marching is required with at least (k+1) th-order accuracy^[11]. At each intermediate time step, the discrete local equilibrium distribution f_E^j is specified such that the mass, momentum and energy conservations are enforced in the collision relaxation term. The equilibrium distribution which satisfies the minimization of entropy has the form^[36],

(8)

In order to be consistent with the weak formulation of the DG method and to retain the high-order accuracy, the unknown coefficient as is also sought in the finite element space,

(9)

The 4(K + 1) or 7(K + 1) unknown coefficients A_s^l are found from the weak formulations of the conservative law in the collision relaxation term. As an example, the weak formulation of the mass conservation is written as

(10)

The yielded non-linear equations are solved iteratively using the Newton's method. More details can be found in Ref. [20]. In this way, the discrete collision term does not give rise to any source or sink of mass, momentum or energy. The conservative nature with respect to the collision relaxation term is very essential to iteration stability, especially for flow at the low Kn number^[37-38].

3 Courant-Friedrichs-Lewy (CFL) condition and positivity-preserving limiter for hypersonic rarefied flow

The RKDG scheme is a shock-capture scheme, which may automatically capture the discontinuities in solutions even if the initial conditions are smooth. This method is energy stable so that it could directly be used to solve control equations with smooth solutions or solutions with weak discontinuities^[21]. However, for solutions with strong shocks, the RKDG scheme will generate spurious oscillations near discontinuities. Because of oscillations, the numerical solutions to the velocity distribution functions may become negative, which is either non-physical or could lead to instability and non-convergence of iterations. Usually, some forms of limiters are used to keep the scheme stable. As a post processor of the approximate solutions, a limiter uses a new polynomial to replace the old one which is deemed to contain oscillations or non-physical values. In order to obtain reasonable results, the new polynomials are required to keep the accuracy and conservativity of the solutions. The most widely used limiter is the so-called TVB limiter^[11]. However, the TVB limiters involve a parameter which should be chosen by the user for different test cases. Therefore, these limiters may take effect in certain smooth cells and destroy accuracy in these cells. The WENO limiters^[15-26] can maintain the high-order accuracy even if they take effect in smooth cells. However, these limiters need to use the information from not only the immediate neighboring cells but also neighbors' neighbors, making its implementation complicated in multi-dimensions^[26]. In Refs.^[28]-^[29], special limiters have been developed to preserve the maximum principle for DG and FVM schemes solving the 2D scalar conservation law on triangular meshes and to preserve positivity for DG schemes solving one-dimensional compressible Euler equations involving source terms. This kind of limiters can preserve strict bounds as well as maintain high-order accuracy. Besides, it is very easy to be implemented. With the same methodology, a positivity-preserving limiter to solve the Boltzmann kinetic model equations on arbitrary triangular meshes is designed, which is independent of the discrete collision operator.

The Boltzmann kinetic model equations are scalar conservative laws with source terms. First of all, applying the first-order forward Euler time discretization to the first weak formulations (see Eq.(5)), i.e., the control equations for cell averages f_Δi, we have

(11)

Here, the superscript n indicates the time step, and the discrete velocity index j has been omitted. For the P^k methods, the quadrature rules for the edges must be exact for polynomials of degree 2k + 1^[11], so that the (k + 1)-point Gaussian rule is used to approximate the edge integral ∫_edΓ. Then, Eq.(11) becomes

(12)

where ω_β is the quadrature weight of the Gauss rule on the interval with , and s_e is the length of the edge e. In order to rewrite the right hand side of Eq.(12) as a monotonically increasing function of some point values of fⁿ under a certain CFL condition, a special rule is introduced to discretize the triangular integral^[29]. The rule combines with the 3-point Gauss-Lobatto rule and the (k + 1)-point Gaussian rule, which is exact for a polynomial f(x, y) if its degree is not larger than k. Given the points and weights of the Gauss-Lobatto rule on the interval as and the ones of Gaussian rule as , the triangle quadrature points indicated by area coordinates and their relative weights W_δ are listed in Table 1. The quadrature points contain 3(k+1) points on the triangle edges which are coincident with the ones of Gaussian quadrature rule, and another 3(k+1) points on the interior of the triangle. Then, the cell average value Q_Δin is expressed as

(13)

where Q_{int, γ}ⁿ is the point value of the interior, and is the weight.

With the conservativity of the flux, the flux integral is decomposed as

(14)

Then, substituting Eqs. (13) and (14) into Eq. (12), we obtain

(15)

where υ_{e, β}, υ_{int, γ}, f_{e, β}^E, and f_{int, γ}^E are the point values of the collision frequency and the equilibrium distribution, and they are all positive values due to the ways through which they are calculated. The functions H_{e, β} are

(16)

Taken H_{2, β} as an example to discuss its property, the function is rewritten as

(17)

Therefore, under the following condition:

(18)

the function H_{2, β} is a monotonically increasing function of , and . Similarly, under the condition, i.e.,

(19)

the functions H_{1, β} and H_{3, β} are also monotonically increasing functions of the point values.

Finally, writing the right-hand side of Eq. (15) as a function of all the point values of fⁿ, this function is monotonically increasing with respect to each argument under Eqs. (18) and (19) plus (1 -Δtυ_{int, γ}) > 0. Therefore, if all the point values fⁿ(xδ) at the time step tⁿ are non-negative, the average value will still be positive at the next time step tⁿ⁺¹. The similar result holds for other monotone fluxes^[29].

Now, we obtain the sufficient conditions to preserve the positivity of the distribution functions. Firstly, the time interval should satisfy the CFL condition,

(20)

where , and υ_max is the maximum collision frequency, and therefore the time interval must be smaller than the minimum mean collision time. Secondly, the values of the distribution functions on the quadrature points x_δ should be positive, which are achieved by a linear scaling limiter^[29]. At each time step, the new polynomial f_Δi^new(x, y) is constructed to replace the solution f_Δiⁿ(x, y), which is defined as

(21)

with

This limiter retains the accuracy and conservativity of the solutions^[39]. Since the high-order Runge-Kutta time discretization is a convex combination of the forward Euler scheme, the full scheme with the high-order time discretization still satisfies the positivity preserving property on the cell average values^[29].

4 Results and discussion

In this section, we provide numerical results to demonstrate the performance of the proposed scheme for the solutions of hypersonic rarefied gaseous flow. The density and temperature profiles of the flow fields are obtained using the second-and third-order RKDG schemes. The solutions are compared with those obtained from experimental data, direct simulation Monte Carlo (DSMC) method as well as analytical calculations. The numerical schemes are implemented in C++, and MPI is used for the parallel version of the code. The three cases are all involved in the near continuum region, where the strong shock wave may occur. The strength of the discontinuity would be gradually weakened as the Knudsen number increases^[40]. As a result, the troubled cell seldom appears in the large Kn number flow, and the limiter would not be used then.

4.1 Normal shock wave: accuracy of positivity-preserving limiter

We first test the accuracy of the positivity-preserving limiter. The results presented below are computed from a normal shock wave in the argon gas at the Mach number Ma₁=2.05. Under this Mach number condition, the RKDG schemes are stable even without limiters. The same parameters as those in the experiment^[41] are used here with the upstream density ρ₁=1.067 ×10^-4 kg/m³ and the temperature T₁=300 K, corresponding to the mean free path of the hard sphere model as λ₁=1.098 ×10^-3 m. The non-dimensionalized boundary conditions obtained from the Rankine-Hugoniot relation are shown in Table 2. The reference speed is c₀=. The second-order and third-order RKDG methods are used to solved the 2D/threevelocity (3V) ES-BGK model on the 2D spatial mesh with the uniform triangulation. The spatial mesh and boundary conditions are shown in Fig. 2. The length of the domain is 40λ₁. The top and bottom boundaries are taken symmetrical. The left boundary is the hypersonic inlet with the upstream condition, and the right boundary is a specular wall moving with the down stream. At the beginning, the x < 0 region is initialized with the upstream condition, while the x > 0 region is initialized with the downstream one. The velocity grid of 13 ×13 ×13 discrete velocities with bounds in the x-direction and in the other two directions, and the second-order quadrature rule are used.

When the L₂ norm of the change in velocity distributions at each time step is less than 10^-5, the steady state solutions are assumed to have been reached.

The obtained results are compared with the DSMC solutions, experimental data, and analytical results. The DSMC simulation uses 300 cells. The average number of molecules per cell is about 50. About 30 000 iterations with a time step of 7.5 ×10^-8 s are needed to reach the steady state. The macro-parameters are sampled over another 50 000 steps. One of the important properties of a shock wave with in a monatomic gas is the overshoot of temperature associated with the longitudinal component of thermal velocities T_x, which could be larger than the gas temperature behind the front of the shock due to the non-equilibrium in translational energies of longitudinal and transversal directions. T_x is related to the number density n as^[42]

(22)

Here, we use this parameter as an analytical criterion to demonstrate the errors of the RKDG solutions.

The relative L₁ and L_∞ errors of the overshoot temperature for the RKDG method with the positivity-preserving limiter compared with the original RKDG method without a limiter are shown in Table 3. The downstream parameters are used to determine the time step Δt. For the methods without limiter, the L₁ and L_∞ errors of the overshoot temperature of the second-order discontinuous Galerkin discretization (DG-2) results on 128 triangles reduce to 5.46×10^-4 and 5.37 ×10^-3, respectively, while the third-order discontinuous Galerkin discretization (DG-3) method obtains the results of the same order accuracy on just 32 triangles. For each case under the same computing conditions, the RKDG scheme combining with the positivity-preserving limiter uses a bit fewer iterations to get the steady flow field, and the errors of the solutions are at the same level as the ones from the method without limiter. Therefore, the limiter keeps the designed order of accuracy. Figure 3 shows the normalized density and temperature . The solid lines are the profiles of DG-3 solution with limiter on 32 triangles, the dash-dot lines are the density distributions from the experiment, and the dots illustrate the DSMC results. The proposed scheme captures the normal shock structure very well.

4.2 Normal shock wave: compared with TVB and WENO limiters

Then, we use the same normal shock wave case for comparison of the performance of the positivity-preserving limiter with the widely used TVB and WENO limiters on the simulation of rarefied gaseous flow. In any cell which is deemed to contain spurious oscillation, i.e., the troubled cell, the DG polynomial is replaced by a new one of the same degree which is less oscillatory than the old one. Different limiters compute this new polynomial in different ways. The TVB limiters are local slope-limiting operators. Here, we implement the one described in Ref. [11]. The TVB limiters involve a problem-dependent parameter M, related to the value of the second derivative of the exact solution near smooth extrema^[43]. Numerical tests obtained with various values of M=10, M=1, M=0.1, and M=0.01 are presented. Generally, the smaller M is, the more troubled cells are. The WENO limiter is based on the WENO methodology for finite volume schemes, which reconstructs the polynomials in the troubled cell using information of neighboring cells. The usual WENO method, which involves the non-linear combination of a set of polynomials based on the cell averages of neighboring cells to directly reconstruct the moments^[26], is used. In order to minimize the influence of accuracy in smooth regions, a trouble indicator is also needed. The TVB indicator and the so-called KrivodonovaXin-Remacle-Chevaugeon-Flaherty (KXRCF) shock detection^[44] are the best choices^[45]. Here, we use the KXRCF indicator in the WENO limiter. In our calculations, the troubled cell indicators are performed for every discrete velocity distribution function during running. Once the threshold is satisfied, the approximation of the distribution function would be reconstructed. Both the TVB and WENO limiters need information from neighboring cells. Ghost cells are built for cells on boundaries, on which the molecular velocity distribution functions are defined by agreeing with the boundary conditions.

The shock wave is still in the argon gas at Ma₁=2.05. The second-order RKDG method combining with the positivity-preserving limiter, TVB limiter, or WENO limiter is used to solve the 2D/3V ES-BGK kinetic model on the 2D spatial mesh with 64 uniform triangles. Other conditions have been described in the previous section. All the cases are run with the same time step of Δt=2.03 ×10^-8 s according to Eq. (20). However, it should be mentioned that it is not necessary for TVB and WENO limiters to follow this CFL condition. Each case is run for 20 000 iterations. Figure 4 shows the residuals, i.e.,

(23)

where f_{Δi, new}^j is the distribution obtained at the current time step, and f_{Δi, old}^j is the one obtained at the previous time step. At the first 5 000 iterations, the behaviors of the solvers are similar to residuals reduced from about 0.2 to 10^-3. Thereafter, the residuals separate a lot. At steps 20 000, the residual of positivity-preserving limiter decreases below 10^-5, while the ones of TVB limiter with M=10, M=1, M=0.1, and M=0.01, as well as the one of WENO limiter vibrate at about 1.2 ×10^-5, 7.0 ×10^-5, 1.0 ×10^-5, 2.0 ×10^-4, and 3.0 ×10^-4, respectively. Figure 5 shows the selected normalized density and temperature profiles at the 15 035th iteration when the residual of positivity-preserving limiter decreases to 10^-5. Table 4 gives the results of residual at the 20 000th iteration, the maximum number of troubled cells, and errors of T_x at the 15 035th iteration for different limiters. It shows that the reconstruction using information from neighborhood brings additional round-off errors to the residual. There exists a best choice of the value M for the TVB limiter, with which the solver could reproduce an accurate solution. However, the best M is problem-dependent. Up to now, it cannot be predicted before calculation. The results of TVB limiter coincide well with the ones obtained from positivity-preserving limiter. The WENO limiter produces results with notable discrepancy.

The comparison indicates that the accuracy of TVB and WENO limiters is highly dependent on problems. The WENO limiter has the largest error mainly due to that the new polynomial of WENO is a non-linear combination ofⁿine linear polynomials constructed using the averages on the troubled cell and on nine neighbouring cells. The combination coefficients are determined based on the smoothness of the linear polynomials. As a result, the new polynomial takes more information from the smoother polynomials. In this case, the gradient only exists along the axial direction. Therefore, the reconstruction of polynomial is mainly based on the information from the transversal direction. As a result, ladder-shaped profiles appear in the density and temperature. Compared with the TVB and WENO limiters, the positivity-preserving limiter has the advantage of easy implementation. It also possesses the advantages as follows: (ⅰ) The troubled-cell detection automatically embeds into the reconstruction procedures, and no additional operator and problem-dependent parameters are required; (ⅱ) The reconstructed polynomials are only dependent on the target piecewise functions, which could minimize the influence of accuracy of the original RKDG method.

4.3 2D supersonic flow past cylinder

In this test, we consider flow past a cylinder of the radius 0.04 m. The free stream is argon with ρ_∞=1.95 ×10^-4 kg/m³ and T_∞=200 K for the density and the temperature, respectively. The variable hard sphere (VHS) model is used to calculate the viscosity, which gives a Knudsen number 0.005 at infinity. We choose such a small Kn to test the stability of the schemes, since flow can generate a strong bow shock passing the cylinder, which is generally regarded as discontinuity in continuum flow.

Two different Mach numbers Ma_∞=6.0 and Ma_∞=12.0 are considered to illustrate the capability of the proposed scheme for flow of the large Mach number, where the velocity of free stream is u_∞=1 580.8 m/s and u_∞=3 161.6 m/s, and the temperature of cylinder wall is set as T_w=500 K and T_w=1 000 K, respectively. The thermal accommodation coefficient is 1.0. The second-order RKDG method with the positivity-preserving limiter is used to solve the BGK model equation. In the physical space, the unstructured triangular mesh with the local refinement is used to partition the computation domain, which is restricted to the upstream flow. 2 313 triangles and 2 203 triangles are generated for cases of Ma_∞=6.0 and Ma_∞=12.0. The left boundary is hypersonic inlet, the bottom one is symmetrical, while the top and right ones are hypersonic outlet. In the velocity space, the case of Ma_∞=6.0 uses discrete velocity points of 35 ranging from -11.07 to 11.07 (non-dimensionalized values) in each velocity direction, the case of Ma_∞=12.0 uses discrete velocity points of 65 in each direction, while the ranges of c_x and c_y are from -18.86 to 18.86 and the range of cz is from -21.23 to 21.23. For this problem, there exists a high temperature region near the stagnation point, where the distribution function has a flat shape with a wide spread, and a low temperature region, where the distribution function has a high peak with a narrow spread. Therefore, the discrete velocity space has to be accurate enough to recover all the transport phenomena. According to the CFL condition, the time step is chosen as Δt=9.02 ×10^-9 s and Δt=6.13 ×10^-9 s for Mach numbers 6.0 and 12.0, respectively. The computational domain is initialized by the free stream condition, and the time convergence criterion is that the L₂ norm of residual reduces by a factor of 104. The obtained results are compared with the DSMC solutions, which are computed on the same domain with the uniform cell size of Δx=2 ×10^-4 m and the time step of Δt=8 ×10^-8 s. About 6 ×10⁶ particles and 60 000 iterations are used, and the macroscopic parameters are sampled over the last 30 000 steps.

Figures 6 and 7 show the density and temperature contours of DG-2 solutions compared with the DSMC results. A strong bow shock is yielded. The flow field structures from both the methods display rather striking similarity. The stagnation line profiles of density, temperature, and velocity for the case Ma_∞=12.0 are illustrated in Fig. 8. In general, good agreement between DG-2 computations and DSMC solutions in the density and velocity can be observed with the maximum error about 5.9%. The larger deviation appears at the bow shock front in temperature. It is due to the use of the BGK kinetic model, which cannot give correct viscosity and conductivity coefficients at the same time. For steady hypersonic flow in the continuum regime, the DG-2 method consumes much more CPU time than the DSMC method. There are probably three reasons as follows: (ⅰ) The large mean collision frequency results in a small time step. (ⅱ) The low Kn number leads to a long time from the initial condition to converged flow. (ⅲ) The standard Cartesian grid brings plenty of velocity nodes. However, we emphasize that this test is used to illustrate stability of the proposed scheme for flow with strong discontinuity, which may occur in the high density and velocity regime. Such regime would not be widespread in rarefied gas flow, but is very important in some cases, e.g., the vacuum plume with the near continuum hypersonic flow around nozzle exit. Therefore, the positivity-preserving limiter is essential to keep stability in calculation of such cases. In order to reduce the computational cost, a local refined discrete velocity grid^[5] or automatic adaption in the physical mesh^[40] and the time step need be considered, especially for the large density difference occurring in flow. That is what we will do in our future works.

4.4 Unsteady shock tube

The third test is a well-known one-dimensional unsteady flow problem, Riemann problem^[46], which treats the development of flow due to initial discontinuity. Removing the diaphragm separating the gas in the two reservoirs of the zero velocity with potentially different pressures, densities and temperatures, results in a characteristic wave system consisting of a shock wave, an expansion fan, and a contact discontinuity. The inviscid shock tube problem can be solved exactly using the gas dynamic theory. The test case chosen here was used by Laney^[47] to evaluate and compare a number of different CFD schemes, with the initial density and pressure ratios of ρ_L/ρ_R=8 and p_L/p_R=10. We use argon gas on both the left and right sides of the tube, and assume an initial pressure of p_L=1.013 25×10³ Pa and temperature of T_L=273 K on the left side.

At an elapsed time of t=6.5 ×10^-6 s, the results from the seconder-order RKDG method solving BGK model equations, are compared with the results from the DSMC method. Both the simulations are carried out on the 2D domain of the size 0.01 m×0.000 05 m with 200 ×1 uniform cells, and the initial discontinuity is located at x=0.005 m. Particularly, reflecting boundaries are set on all sides. Each simulation uses a time step interval of Δt=1.0 ×10^-8 s, which is about 80% of the mean collision time on the left side. For the RKDG calculation, the discrete velocities of 10 ranging from -3.54 to 3.54 in each velocity direction are used. For the DSMC simulation, in order to reduce the sample noise, the initial population is about 32 000 particles per cell on the left side and 4 000 particles per cell on the right side. Calculations are performed on 8 CPUs, the RKDG method takes about 0.067 h, and the DSMC method takes about 0.1h. Exact Riemann results are also used for comparison. Gas density and temperature profiles are shown in Fig. 9, with values normalized by the properties on the left side. The solid lines are the analytical results, dots illustrate the solutions from the proposed scheme, and the dashed lines with the statistical scatter are the ones from the DSMC method. The RKDG results agree very well with the results from the particle method, and both agree reasonably well with the profiles from the exact Riemann solver. The deviations at the bends of the profiles are mainly due to the gas rarefaction, which has weaken the strength of fluctuations in gas properties. For the analytical solution, it is assumed that the gas flow follows a Maxwellian distribution, and its evolution is ruled by Euler's equations. However, in rarefied gas, there is not enough molecular collisions to reach equilibrium. Therefore, the thickness of the shock wave becomes wider, and the contact discontinuity gradually fades out. The higher the rarefaction is, the larger the deviation is^[48].

5 Conclusions

A stable RKDG scheme which strictly preserves positivity of the solutions, has been designed to solve the Boltzmann kinetic model equations for hypersonic rarefied gaseous flow. Stability is kept by accuracy of the velocity discretization, the conservation of the collision terms and a limiter. For hypersonic flow, the appropriate velocity ranges could become quite wide, and a large number of discrete velocities are needed to meet the accuracy requirement. Therefore, the equally-spaced composite Newton-Cotes formulae are more convenient to approximate the velocity integral with appropriate velocity bounds. For the calculation of collision term, the equilibrium velocity distribution function is estimated using a discontinuous conservative discretization method, which enforces a weak conservation of mass, momentum, and energy. Based on the first-order forward Euler time discretization and a special triangle quadrature rule, the sufficient conditions which keep positivity of the cell average solutions of the velocity distribution functions are obtained. These sufficient conditions require that the time step is smaller than the minimum local collision time, and that the values of velocity distribution functions on the quadrature points are always non-negative. The later requirement is forced by a linear limiter. Since the high-order Runge-Kutta time discretization is a convex combination of Euler forward, the full scheme with the high-order time discretization still satisfies the positivity-preserving property on the cell average values. Verification of the scheme has been performed by comparison with the experimental data, DSMC and analytical solutions for a normal shock wave at the Mach number 2.05, hypersonic flow at the Mach numbers 6.0 and 12.0 passing a 2D cylinder, and an unsteady shock tube flow. Results show that, the scheme is stable and accurate to capture the shock structures in steady and unsteady hypersonic rarefied gaseous flow. Compared with the widely used TVB and WENO limiters, the current limiter has the advantage of easy implementation and the ability of minimizing the influence of accuracy of the original RKDG method. It would not destroy the accuracy of the evolution of the Boltzmann model equations, which are valid for the entire range of rarefaction.

Acknowledgements This work is supported by the National Supercomputer Center in Tianjin. All the numerical tests are run on the TH-I cluster.