Skip to main content

Inverse design of mesoscopic models for compressible flow using the Chapman-Enskog analysis

Abstract

In this paper, based on simplified Boltzmann equation, we explore the inverse-design of mesoscopic models for compressible flow using the Chapman-Enskog analysis. Starting from the single-relaxation-time Boltzmann equation with an additional source term, two model Boltzmann equations for two reduced distribution functions are obtained, each then also having an additional undetermined source term. Under this general framework and using Navier-Stokes-Fourier (NSF) equations as constraints, the structures of the distribution functions are obtained by the leading-order Chapman-Enskog analysis. Next, five basic constraints for the design of the two source terms are obtained in order to recover the NSF system in the continuum limit. These constraints allow for adjustable bulk-to-shear viscosity ratio, Prandtl number as well as a thermal energy source. The specific forms of the two source terms can be determined through proper physical considerations and numerical implementation requirements. By employing the truncated Hermite expansion, one design for the two source terms is proposed. Moreover, three well-known mesoscopic models in the literature are shown to be compatible with these five constraints. In addition, the consistent implementation of boundary conditions is also explored by using the Chapman-Enskog expansion at the NSF order. Finally, based on the higher-order Chapman-Enskog expansion of the distribution functions, we derive the complete analytical expressions for the viscous stress tensor and the heat flux. Some underlying physics can be further explored using the DNS simulation data based on the proposed model.

1 Introduction

The Boltzmann equation is of vital importance in the kinetic theory of dilute gases [1]. It is noted that the construction of gas kinetic models has a long history. After the pioneering work of Wang-Chang and Uhlenbeck [2] on the eigenvalues and eigenfunctions of the linearized Boltzmann equation, Gross and Jackon [3] established the relations between the linearized Boltzmann equation and some models for monatomic gas by comparing the eigenvalue spectra of the collision operators. Later, Hanson and Morse [4] obtained the kinetic model equations for polyatomic gases based on Wang Chang-Uhlenbeck equation. The original collision operator in the Boltzmann equation is a complex integral term, which makes direct numerical simulation of the system very costly. The simplest choice is to replace the original collision operator with the Bhatnager-Gross-Krook (BGK) model [5]. It should be noted that the original collision operator is directly based on the physical description of molecule interactions while the BGK model describes the fact that the distribution of the molecules relaxes to the local equilibrium state through particle collisions without considering the detailed molecule interactions. It has long been recognized that such an approximation works well beyond its theoretical limit as long as the relaxation time can be made to capture the relevant physics [6, 7].

By applying the Chapman-Enskog expansion to the Boltzmann equation with the BGK collision operator, the NSF system can be recovered, but with a unit Prandtl number which does not obey the physical reality. Hence, some improved models have been developed to overcome this limitation from different physical considerations, such as the Shakhov (S) model [8], the ellipsoidal statistical (ES) model [9], the internal energy double-distribution-function (IEDDF) model [10], the total energy double-distribution-function (TEDDF) model [11], and the Rykov (R) model [1214]. Further, we notice that, for example, the ratio between the bulk to shear viscosity in S model is always less than 2/3. Therefore, for a fixed specific heat ratio, the S model cannot be used to investigate the physical effect of the bulk-to-shear viscosity ratio in compressible flows.

Some merits and drawbacks of these models are briefly discussed below. The S model may encounter a negative value of the particle distribution function because of the modified equilibrium distribution to accommodate arbitrary Prandtl number, while the ES model and the TEDDF model can avoid such unphysical deficiency. However, Chen et al. [15] showed that the S model may yield more accurate solutions than that from the ES model in the transition regime and they proposed a generalized model which combines the advantages of the S model and ES model. For both IEDDF and TEDDF models, two distribution functions are introduced with different relaxation times for the non-equilibrium part of the particle distribution function because the momentum and energy have different relaxation time scales during the collision process as suggested by Wood [16]. In the TEDDF model, spatial and time derivatives of the hydrodynamic velocity are not involved in the source terms while they are involved in the source terms of the IEDDF model which may introduce some numerical errors and lead to some unphysical phenomena in fluid systems containing large spatial gradients. In the R model, by considering the elastic and non-elastic particle collision processes, the hydrodynamic flow variables corresponding to the translational and rotational processes can be evaluated separately. Both the total internal energy and the total heat flux are the sum of the contributions from the two processes. The bulk-to-shear viscosity ratio can be modified in the R model through the ratio of the total number of translational and rotational collisions to that of rotational collisions.

During the past few decades, the BGK model has been widely used to simulate different flows such as homogeneous isotropic turbulence [17], turbulent channel flows [18] and multiphase flows [19], by different numerical approaches such as the lattice Boltzmann method (LBM) [20], the gas kinetic scheme (GKS) [21], the unified gas kinetic scheme (UGKS) [22], and the discrete unified gas kinetic scheme (DUGKS) [23, 24]. Recently, Liu et al. [25] claimed that the predictions based on the BGK model for highly nonequilibrium flows are only qualitatively correct in the transitional regime since the BGK model filters out the information of the detailed molecular-interaction processes. They compare the Boltzmann equation and its model equations through some test cases where the distribution functions are far from equilibrium. From these tests, they found that information contained in the nonequilibrium moments and the different relaxation rates of high- and low-speed molecules is essential in adjusting the behaviors of model collision terms. However, many existing works have shown that the BGK model is adequate in simulating many flows accurately for both the continuum and rarefied regimes [1719, 2628].

In this paper, we focus on the inverse design of the source term in the model Boltzmann equation for compressible flows. Following the work done by Guo et al. [24], an adjustable parameter representing the internal degree of freedom of molecules is introduced to the Maxwellian equilibrium distribution function. We will demonstrate that the two source terms in the two reduced model Boltzmann equations can be redesigned to attain the following objectives. First, the NSF system can be recovered in the continuum limit by applying the Chapman-Enskog analysis. Second, the model Boltzmann system can have flexible Prandtl number as well as adjustable bulk-to-shear viscosity ratio. Third, an arbitrary thermal source/sink term can be added to the internal energy equation.

The rest of the paper is organized as follows. In Section 2, the model Boltzmann equation with an additional source term is introduced. By introducing two reduced distribution functions, two reduced model Boltzmann equations are obtained. Some notations and conventions are given in Section 3. In Section 4, the structures of the Boltzmann equations are obtained by applying the first-order Chapman-Enskog expansion. Next, five requirements for the two reduced source terms are given in Section 5. In Section 6, we present one design for the two source terms by applying the Hermite expansion to the two source terms. In Section 7, we show that the S model, the TEDDF model as well as the R model are compatible with the five derived constraints. Then we discuss the derivation of the proper implementation of the hydrodynamic boundary conditions in Section 8. Next, we derive the complete analytical expression for the viscous stress and the heat flux based on the second-order Chapman-Enskog expansion in the following three sections. Major conclusions are drawn in Section 12. In Appendix 13, we include the details on the Hermite polynomials and Hermite expansion. Appendix 14 contains the derivations of the requirements for the two reduced source terms. Appendix 15 documents some details on the Rykov model.

2 The reduced model Boltzmann system with source terms

The Boltzmann equation with an additional source term can be expressed as

$$\begin{array}{@{}rcl@{}} \frac{\partial f}{\partial t}+\boldsymbol{\xi}\cdot \boldsymbol{\nabla} f+\boldsymbol{a}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}}f=\Omega_{f}+S_{f}, \end{array} $$
(1)

where f(x,ξ,t) is the particle distribution function, x=(x1,...,x3) is the spatial location, t is the time, ξ=(ξ1,...,ξ3) is the particle velocity in three dimensional space, and ζ=(ζ1,...,ζK) represents K-dimensional internal degree of freedoms. a represents the body force per unit mass. The single-relaxation-time Bhatnager-Gross-Krook (BGK) model [5] is used for the collision operator, i.e. Ωf=(feqf)/τ. τ=μ/p is the molecular relaxation time and μ is the shear viscosity which can be approximated by the hard-sphere model [24, 29] or Sutherland’s law [30, 31]. p=ρRT is the pressure for ideal gas. Sf is a source term to be designed, which will allow for modification of both the Prandtl number Pr as well as the bulk-to-shear viscosity ratio χ=μV/μ with μV being the bulk viscosity.

By assuming that the particle motion in ζ subspace is at local equilibrium, the local Maxwellian equilibrium distribution function can be written as [24]

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

where ρ is the density of the fluid, R is the specific gas constant, T is the temperature, and c=ξu is the peculiar velocity with u being the hydrodynamic velocity.

The conservative variables are defined as the moments of the particle distribution function

$$\begin{array}{@{}rcl@{}} \rho=\int fd\boldsymbol{\xi}d\boldsymbol{\zeta},~ \rho\boldsymbol{u}=\int\boldsymbol{\xi}fd\boldsymbol{\xi}d\boldsymbol{\zeta},~ \rho E=\frac{1}{2} \rho u^{2}+\rho \epsilon=\int \frac{\xi^{2}+\zeta^{2}}{2}fd\boldsymbol{\xi}d\boldsymbol{\zeta}, \end{array} $$
(3)

where ε=CvT is the internal energy per unit mass, Cv is the specific heat capacity at constant volume, and ρE is the total energy per unit volume which is the sum of the internal energy and the kinetic energy. All relations in Eq. (3) remain valid if f is replaced by feq. Cv and the specific heat at constant pressure Cp are determined by the number of the internal degrees of freedom, K, and the gas constant, R. By integrating the energy moment of the equilibrium distribution, we can obtain Cv=(K+3)R/2 and Cp=(K+5)R/2, which implies that the specific heat ratio and thus the Prandtl number are γ=Cp/Cv=(K+5)/(K+3) and Pr=μCp/κ, where κ is the thermal conductivity.

In addition, by comparing the first-order moment of the model Boltzmann equation with the Navier-Stokes equation, it can be shown that the viscous stress tensor σ is determined by the non-equilibrium part of the particle distribution function as

$$\begin{array}{@{}rcl@{}} \boldsymbol{\sigma}=-\int \boldsymbol{c} \boldsymbol{c}\left(f-f^{eq}\right)d\boldsymbol{\xi} d\boldsymbol{\zeta}, \end{array} $$
(4)

and, by comparing the energy moment of the Boltzmann equation with the macroscopic energy equation, the heat flux q can be determined as

$$\begin{array}{@{}rcl@{}} \boldsymbol{q}=\frac{1}{2}\int \boldsymbol{c}\left(c^{2}+\zeta^{2}\right)fd\boldsymbol{\xi} d\boldsymbol{\zeta}. \end{array} $$
(5)

The physical conservative requirements can be expressed through the moments of the collision operator Ωf, which reads

$$\begin{array}{@{}rcl@{}} \int\Omega_{f}d\boldsymbol{\xi} d\boldsymbol{\zeta}=0,~\int\boldsymbol{\xi}\Omega_{f} d\boldsymbol{\xi} d\boldsymbol{\zeta}=\boldsymbol{0},~ \int\frac{1}{2}\left(\xi^{2}+\zeta^{2}\right)\Omega_{f} d\boldsymbol{\xi} d\boldsymbol{\zeta}=0. \end{array} $$
(6)

Therefore, provided that the mass conservation and the momentum conservation laws are observed, we have the following basic requirements for the source term

$$\begin{array}{@{}rcl@{}} \int S_{f} d\boldsymbol{\xi} d\boldsymbol{\zeta}=0,~\int \boldsymbol{\xi} S_{f} d\boldsymbol{\xi} d\boldsymbol{\zeta}=\boldsymbol{0},~ \int\frac{1}{2}\left(c^{2}+\zeta^{2}\right)S_{f}d\boldsymbol{\xi} d\boldsymbol{\zeta}=-\rho\Lambda, \end{array} $$
(7)

where Λ represents a source term applied to the energy equation, an example of which is the thermal cooling function [30].

Physically, the evolution of the particle distribution function only depends on the particle velocity ξ. In order to remove the dependence of the passive variables and also reduce the computational cost in the practical implementation, two independent, reduced distribution functions g and h, residing in lower dimensional phase space, are introduced [24], namely, \(g=\int f d\boldsymbol {\zeta }\) and \(h=\int \zeta ^{2} fd\boldsymbol {\zeta }\). Therefore, the two model Boltzmann equations residing in lower dimensional space can be obtained

$$\begin{array}{@{}rcl@{}} \frac{\partial g}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla}{g}+\boldsymbol{a}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}}g=\Omega_{g}+S_{g},~ \frac{\partial h}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla}{h}+\boldsymbol{a}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}}h=\Omega_{h}+S_{h}. \end{array} $$
(8)

In Eq. (8), the collision operators and source terms are

$$\begin{array}{@{}rcl@{}} \Omega_{g}=\frac{g^{eq}-g}{\tau},~\Omega_{h}=\frac{h^{eq}-h}{\tau},~ S_{g}=\int S_{f} d\boldsymbol{\zeta},~S_{h}=\int \zeta^{2} S_{f} d\boldsymbol{\zeta}, \end{array} $$
(9)

where the equilibrium distribution functions geq and heq are

$$\begin{array}{@{}rcl@{}} g^{eq}=\frac{\rho}{(2\pi RT)^{3/2}}\exp\left[-\frac{c^{2}}{2RT}\right],~ h^{eq}=KRT g^{eq}. \end{array} $$
(10)

Based on Eq. (6), the conservation laws can be recasted in terms of the collision operators Ωg and Ωh, as follows,

$$\begin{array}{@{}rcl@{}} \int\Omega_{g}d\boldsymbol{\xi}=0,~\int\boldsymbol{\xi}\Omega_{g}d\boldsymbol{\xi}=\boldsymbol{0},~\int\left(\xi^{2}\Omega_{g}+\Omega_{h}\right)d\boldsymbol{\xi}=0. \end{array} $$
(11)

From Eq. (7), the two reduced source terms must satisfy the following requirements

$$\begin{array}{@{}rcl@{}} \int S_{g} d\boldsymbol{\xi}=0,~\int \boldsymbol{\xi} S_{g}d\boldsymbol{\xi}=\boldsymbol{0},~\int\frac{1}{2}\left(c^{2} S_{g}+S_{h}\right)d\boldsymbol{\xi}=-\rho\Lambda. \end{array} $$
(12)

In addition, from Eq. (3), we find that the conservative variables can be evaluated as

$$\begin{array}{@{}rcl@{}} \rho=\int gd\boldsymbol{\xi},~ \rho\boldsymbol{u}=\int\boldsymbol{\xi}gd\boldsymbol{\xi}, ~\rho E=\frac{1}{2}\int\left(\xi^{2} g+h\right)d\boldsymbol{\xi}. \end{array} $$
(13)

Moreover, from Eqs. (4) and (5), the viscous stress σ and the heat flux q become

$$\begin{array}{@{}rcl@{}} \boldsymbol\sigma=-\int\boldsymbol{c}\boldsymbol{c}\left(g-g^{eq}\right)d\boldsymbol{\xi},~\boldsymbol{q}=\frac{1}{2}\int \boldsymbol{c}\left(c^{2}g+h\right)d\boldsymbol{\xi}. \end{array} $$
(14)

3 Notations and conventions

For convenience, two time derivatives are introduced

$$\begin{array}{@{}rcl@{}} &&\frac{D}{Dt} \equiv \frac{\partial}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla} +\boldsymbol{a}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}},~~~ \frac{d}{dt} \equiv \frac{\partial}{\partial t}+\boldsymbol{u}\cdot\boldsymbol{\nabla}, \end{array} $$
(15)

where D/Dt is the time derivative along the phase-space trajectory of a particle subjected to a body force a per unit mass and d/dt is the rate of change of a physical quantity along the path of a fluid element in the physical space. Three variables including the time t, the spacial location x, and the particle velocity ξ are assumed to be independent when these time derivatives act on the distribution functions g(x,ξ,t) and h(x,ξ,t).

In addition, S=(uT+u)/2 is the strain rate tensor and Ω=(uTu)/2 is the rotation tensor. 𝜗=·u is the velocity divergence (also known as dilatation).

The Newtonian constitutive relation for the viscous stress σ(NS) and the Fourier’s law for the heat flux q(NS) are, respectively,

$$\begin{array}{@{}rcl@{}} &&\boldsymbol\sigma^{(NS)}=2\mu\left(\boldsymbol{S}-\frac{1}{3}\vartheta\boldsymbol{I}\right)+\mu_{V}\vartheta\boldsymbol{I},~\boldsymbol{q}^{(NS)}=-\kappa\boldsymbol{\nabla} T. \end{array} $$
(16)

4 The structure of the particle distribution functions

In the continuum limit, the relaxation time τ, when normalized by the acoustic time scale l0/c0, is proportional to the Knudsen number, where l0 is a system length scale and c0 is the speed of the sound at a reference temperature T0. Therefore, τ may be taken as a small parameter in the Boltzmann equation. At the level of NSF equations, terms higher than O(τ) in the distribution functions can be neglected.

The derivatives of the equilibrium distribution function geq will be multiplied by τ to form the O(τ) terms in the distribution functions. Therefore, we only need to evaluate them to O(1). Direct evaluation yields the derivatives of the equilibrium distribution function geq as

$$\begin{array}{@{}rcl@{}} &&\frac{\partial g^{eq}}{\partial t} =\left[\frac{1}{\rho}\frac{\partial \rho}{\partial t} +\left(\frac{c^{2}}{2RT}-\frac{3}{2}\right)\frac{1}{T}\frac{\partial T}{\partial t}+\frac{\partial \boldsymbol{u}}{\partial t}\cdot\frac{\boldsymbol{c}}{RT}\right]g^{eq},\\ &&\boldsymbol{\nabla} g^{eq} =\left[\frac{1}{\rho}\boldsymbol{\nabla}\rho +\left(\frac{c^{2}}{2RT}-\frac{3}{2}\right)\frac{1}{T}\boldsymbol{\nabla} T+\boldsymbol{\nabla}\boldsymbol{u}\cdot\frac{\boldsymbol{c}}{RT}\right]g^{eq},\\ &&\boldsymbol{\nabla}_{\boldsymbol{\xi}}g^{eq}=-\frac{\boldsymbol{c}}{RT}g^{eq}. \end{array} $$
(17)

The three coefficients in three derivatives are found to be polynomials of the peculiar velocity c and are related to the time t and spatial location x through the relation c=ξu.

By employing the Euler equations in Appendix 14 to replace the time derivatives of the hydrodynamic variables with spatial derivatives in geq/t, we obtain the expression for Dgeq/Dt and Dheq/Dt to the leading order, as

$$\begin{array}{@{}rcl@{}} &&\frac{Dg^{eq}}{Dt}=Gg^{eq}+O(\tau),~\frac{Dh^{eq}}{Dt}=\left(G+\Phi_{1}\right)h^{eq}+O(\tau), \end{array} $$
(18)

where G=G1+G2+G3. The coefficients are given explicitly as

$$\begin{array}{@{}rcl@{}} &&G_{1}=\left(\frac{c^{2}}{2RT}-\frac{5}{2} \right)\boldsymbol{c}\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right),\\ &&G_{2}=\frac{\boldsymbol{c} \cdot\boldsymbol{S}\cdot \boldsymbol{c}}{RT}-\frac{1}{K+3} \left(\frac{c^{2}}{RT}+K\right)\vartheta,\\ &&G_{3}=-\left(\frac{c^{2}}{2RT}-\frac{3}{2}\right) \frac{\Lambda}{C_{v}T},\\ &&\Phi_{1}=\boldsymbol{c}\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right)-\frac{2}{K+3}\vartheta-\frac{\Lambda}{C_{v}T}. \end{array} $$
(19)

Therefore, up to the order O(τ) in the Chapman-Enskog expansion [1], the structure of the distribution functions g and h can be obtained and they are

$$\begin{array}{@{}rcl@{}} g=(1-\tau G)g^{eq}+\tau S_{g}+O\left(\tau^{2}\right),~ h=\left(1-\tau G-\tau \Phi_{1}\right)h^{eq}+\tau S_{h}+O\left(\tau^{2}\right). \end{array} $$
(20)

5 Five basic requirements for the two source terms

Based on the structure of the distribution function, we shall now propose five basic requirements for the two source terms, Sg and Sh. The requirements are given as follows and the details for their derivations are included in Appendix 14. If these five requirements are satisfied up to the order of O(τ), then the NSF system can be recovered by applying the Chapman-Enskog expansion to the two model Boltzmann equations.

The first and second requirements come from the continuity and the momentum equations

$$ \int S_{g} d\boldsymbol{\xi}=0,~ \int \boldsymbol{c} S_{g}d\boldsymbol{\xi}=\int \boldsymbol{\xi} S_{g} d\boldsymbol{\xi}=\boldsymbol{0}. $$
(21a)

The third requirement is used to modify the bulk viscosity in the viscous stress term, and it is

$$ \int\boldsymbol{\xi} \boldsymbol{\xi} S_{g}d\boldsymbol{\xi}=\int \boldsymbol{c}\boldsymbol{c} S_{g}d\boldsymbol{\xi}=-\left(\chi-\frac{2K}{3(K+3)}\right)p\vartheta\boldsymbol{I}-\frac{p\Lambda}{C_{v}T}\boldsymbol{I}. $$
(21b)

The fourth requirement follows from the energy equation and it is

$$ \int S_{h}d\boldsymbol{\xi}=-\int \xi^{2} S_{g}d\boldsymbol{\xi}-2\rho\Lambda. $$
(21c)

The fifth requirement is expressed as

$$ \int \boldsymbol{c} S_{h}d\boldsymbol{\xi}=\frac{2\left(1-Pr\right)\boldsymbol{q}^{(NS)}}{\tau}-\int \boldsymbol{c} c^{2} S_{g}d\boldsymbol{\xi}, $$
(21d)

which is needed to modify the heat flux and thus the resulting Prandtl number.

As a result of the design constraints, Eqs. (21a)–(21d), the model Boltzmann equation will yield the following NSF system

$$\begin{array}{@{}rcl@{}} &&\frac{\partial\rho}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u})=0,\\ &&\frac{\partial (\rho \boldsymbol{u})}{\partial t}+\boldsymbol{\nabla}\cdot(\rho \boldsymbol{u}\boldsymbol{u})=-\boldsymbol{\nabla} p+ \boldsymbol{\nabla}\cdot\boldsymbol\sigma^{(NS)}+\rho \boldsymbol{a},\\ &&\frac{\partial\left(\rho E\right)}{\partial t}+\nabla\cdot(\rho E \boldsymbol{u})= -\boldsymbol{\nabla}\cdot\boldsymbol{q}^{(NS)}-\boldsymbol{\nabla}\cdot(p\boldsymbol{u}) +\boldsymbol{\nabla}\cdot\left(\boldsymbol\sigma^{(NS)}\cdot\boldsymbol{u}\right) +\rho \boldsymbol{a}\cdot \boldsymbol{u}-\rho \Lambda. \end{array} $$
(22)

6 A possible design of the two reduced source terms

There are many possible ways to design the specific form of the two source terms. By applying the Hermite expansion to the two source terms, a new mesoscopic model with both adjustable Prandtl number and bulk viscosity is proposed next. Any reasonable design of the two source terms should satisfy the five basic requirements presented in Eqs. (21a) to (21d) in the continuum limit.

Due to the desire to keep the order of Gauss-Hermite quadrature as low as feasible in the numerical implementation, we further require \(\int \boldsymbol {\xi \xi \xi } S_{g} d\boldsymbol {\xi }=\boldsymbol {0}\). Then, using Eq. (21a), we have \(\int S_{g}\mathscr {H}^{(3)}\left (\boldsymbol {\xi },T_{0}\right)d\boldsymbol {\xi }=\boldsymbol {0}\). The Eqs. (21a)–(21b) can also be written as

$$\begin{array}{@{}rcl@{}} &&\int S_{g}\mathscr{H}^{(0)}\left(\boldsymbol{\xi},T_{0}\right)d\boldsymbol{\xi}=0,~ \int S_{g}\mathscr{H}^{(1)}\left(\boldsymbol{\xi},T_{0}\right) d\boldsymbol{\xi}=\boldsymbol{0},\\ && \int S_{g}\mathscr{H}^{(2)}\left(\boldsymbol{\xi},T_{0}\right)d\boldsymbol{\xi}=-\left(\chi-\frac{2K}{3(K+3)}\right)\frac{p\vartheta}{RT_{0}}\boldsymbol{I}-\frac{p}{RT_{0}}\frac{\Lambda}{C_{v}T}\boldsymbol{I}, \end{array} $$
(23)

where T0 is a reference temperature and the velocity ξ are scaled with \(\sqrt {RT_{0}}\).

Therefore, by using Eq. (23) and keeping the Hermite expansion (see Appendix 13) of the source term Sg up to the third-order, we obtain

$$\begin{array}{@{}rcl@{}} &&S_{g}\left(\boldsymbol{x},\boldsymbol{\xi},t\right) = \frac{1}{2!}\omega\left(\boldsymbol{\xi},T_{0}\right)\boldsymbol{a}^{(2)}(\boldsymbol{x},t):\mathscr{H}^{(2)}(\boldsymbol{\xi},T_{0})\\ && =-\omega\left(\boldsymbol{\xi},T_{0}\right)\left[\left(\chi-\frac{2K}{3(K+3)}\right)\frac{p\vartheta}{2RT_{0}}+\frac{p}{2RT_{0}}\frac{\Lambda}{C_{v}T}\right]\left(\frac{\xi^{2}}{RT_{0}}-3\right). \end{array} $$
(24)

From Eqs. (21a)–(21b), we can obtain

$$\begin{array}{@{}rcl@{}} &&\int \boldsymbol{c}c^{2}S_{g}d\boldsymbol{\xi}=5\left(\chi-\frac{2K}{3(K+3)}\right)p\vartheta\boldsymbol{u}+\frac{10}{K+3}\rho\Lambda\boldsymbol{u}. \end{array} $$
(25)

Combination of Eqs. (21c), (21d) and (25) yields

$$\begin{array}{@{}rcl@{}} &&\int\boldsymbol{\xi} S_{h}d\boldsymbol{\xi}=\frac{2(1-Pr)\boldsymbol{q}^{(NS)}}{\tau}+\boldsymbol{u}\int S_{h}d\boldsymbol{\xi}-\int \boldsymbol{c} c^{2}S_{g}d\boldsymbol{\xi}\\ &&=\frac{2(1-Pr)\boldsymbol{q}^{(NS)}}{\tau}-2\left(\chi-\frac{2K}{3(K+3)}\right)p\vartheta\boldsymbol{u} -\frac{2(K+5)}{K+3}{\rho\Lambda \boldsymbol{u}}. \end{array} $$
(26)

Therefore, we obtain

$$\begin{array}{@{}rcl@{}} &&\int S_{h}\mathscr{H}^{(1)}\left(\boldsymbol{\xi},T_{0}\right)d\boldsymbol{\xi}\\ &&=\frac{2(1-Pr)\boldsymbol{q}^{(NS)}}{\tau\sqrt{RT_{0}}} -2\left(\chi-\frac{2K}{3(K+3)}\right)\frac{p\vartheta\boldsymbol{u}}{\sqrt{RT_{0}}}-\frac{2(K+5)}{K+3}\frac{\rho\Lambda\boldsymbol{u}}{\sqrt{RT_{0}}}. \end{array} $$
(27)

Eqs. (21b) and (21c) together imply that

$$\begin{array}{@{}rcl@{}} &&\int S_{h}\mathscr{H}^{(0)}\left(\boldsymbol{\xi},T_{0}\right)d\boldsymbol{\xi} =3\left(\chi-\frac{2K}{3(K+3)}\right)p\vartheta-\frac{2K}{K+3}\rho\Lambda. \end{array} $$
(28)

Combination of Eqs. (27) and (28) yields one design for the source term Sh

$$\begin{array}{@{}rcl@{}} S_{h}=\omega\left(\boldsymbol{\xi},T_{0}\right)\left[ \begin{aligned} &\frac{2(1-Pr)\boldsymbol{q}^{(NS)}\cdot\boldsymbol{\xi}}{\tau RT_{0}}+\left(3-2\frac{\boldsymbol{u}\cdot\boldsymbol{\xi}}{RT_{0}}\right)\left(\chi-\frac{2K}{3(K+3)}\right)p\vartheta\\ &-2\left(\frac{K}{K+3}+\frac{K+5}{K+3}\frac{\boldsymbol{u}\cdot\boldsymbol{\xi}}{RT_{0}}\right)\rho\Lambda \end{aligned} \right]. \end{array} $$
(29)

In low Mach number thermal flows, the second-order central finite difference scheme is adequate for calculating ·u in Eqs. (24) and (29). In addition, two different methods can be used to evaluate the heat flux q in Eq. (29). One way is to first obtain T by the second-order central finite difference scheme, then it follows that q=−κT. Another way is to evaluate the heat flux through the velocity moment of distribution functions. For example, under the framework of DUGKS proposed by Guo et al. [24], the heat flux term at the cell center and the cell interface can be calculated by \(\boldsymbol {q}=\frac {2\tau }{2\tau +\Delta t Pr}\frac {1}{2}\int \boldsymbol {c}\left (c^{2} \tilde g+\tilde h\right) d\boldsymbol {\xi }\) and \(\boldsymbol {q}=\frac {2\tau }{2\tau +(\Delta t/2) Pr}\frac {1}{2}\int \boldsymbol {c}\left (c^{2} \bar {g}+\bar {h}\right) d\boldsymbol {\xi }\) after the distribution functions \(\tilde g, \tilde h\) and \(\bar g, \bar h\) are updated. No obvious difference can be observed for global turbulence statistics when the velocity divergence in the source terms is evaluated by either a second-order or fourth-order finite difference scheme in compressible isotropic turbulence [32, 33]. Some improved shock capturing schemes with low dissipation and high-order accuracy could also be tried to evaluate ·u when the turbulent Mach number is further increased.

7 An examination of three existing mesoscopic models in our design framework

7.1 The Shakhov model

In this section, we will prove that the well-known Shakhov model [8, 24] is compatible with our inverse-design here. Through the Chapman-Enskog analysis, it can be verified that the ratio of bulk viscosity to shear viscosity in the Shakhov model is μV/μ=2K/[3(K+3)]. In addition, no cooling function is considered, namely, Λ=0. Therefore, the five general requirements (Eqs. (21a) to (21d)) for the two source terms can be reduced as

$$\begin{array}{@{}rcl@{}} &&\int S_{g}d\boldsymbol{\xi}=0,~\int \boldsymbol{c} S_{g}d\boldsymbol{\xi}=\boldsymbol{0},~\int \boldsymbol{cc} S_{g}d\boldsymbol{\xi}=\boldsymbol{0},~\int S_{h} d\boldsymbol{\xi}=0,\\ &&\int \boldsymbol{c} c^{2} S_{g}d\boldsymbol{\xi}+\int \boldsymbol{c} S_{h}d\boldsymbol{\xi}=\frac{2(1-Pr)\boldsymbol{q}^{(NS)}}{\tau}. \end{array} $$
(30)

By using Eq. (30), the source term Sg can be assumed as

$$\begin{array}{@{}rcl@{}} &&S_{g}=\frac{1}{3!}\omega(\boldsymbol{c}, T) \boldsymbol{a}^{(3)}(\boldsymbol{x},t):\mathscr{H}^{(3)}(\boldsymbol{c}, T), \end{array} $$
(31)

where ω(c,T) is the peculiar-velocity-based weighting function, a(3)(x,t) is the coefficient and \(\mathscr {H}^{(3)}(\boldsymbol {c}, T)\) is the third-order Hermite polynomial.

Note that \(\boldsymbol {a}^{(3)}(\boldsymbol {x},t)\equiv \int S_{g}\mathscr {H}^{(3)}(\boldsymbol {c}, T)d\boldsymbol {\xi }\) is also symmetrical with respect to the components of c because \(\mathscr {H}^{(3)}(\boldsymbol {c}, T)\) remains unchanged under the permutation operation, the simplest way is to assume that the coefficient a(3) takes the form \(a_{ijk}^{(3)}=A_{i} \delta _{jk}+A_{j} \delta _{ki}+A_{k}\delta _{ij}\), where A(x,t)=(A1,A2,A3) is a vector coefficient to be determined. Then, Eq. (31) gives

$$\begin{array}{@{}rcl@{}} &&S_{g}=\frac{1}{2}\omega(\boldsymbol{c}, T)A_{i} \mathscr{H}_{ijj}^{(3)}\left(\boldsymbol{c}, T\right)=\frac{1}{2}\omega(\boldsymbol{c}, T)\boldsymbol{A}\cdot\frac{\boldsymbol{c}}{\sqrt{RT}} \left(\frac{c^{2}}{RT}-5\right), \end{array} $$
(32)

Moreover, it can be shown that \(\int \boldsymbol {c} c^{2}S_{g}d\boldsymbol {c}=5\boldsymbol {A} (RT)^{3/2}\).

Similarly, the source term Sh can be designed as

$$\begin{array}{@{}rcl@{}} &&S_{h}=\omega(\boldsymbol{c},T)\boldsymbol{B}\cdot \mathscr{H}^{(1)}(\boldsymbol{c},T)+\frac{1}{2}\omega(\boldsymbol{c}, T)\boldsymbol{C}\cdot\frac{\boldsymbol{c}}{\sqrt{RT}} \left(\frac{c^{2}}{RT}-5\right). \end{array} $$
(33)

Therefore, the coefficient A and B should satisfy the following relation,

$$\begin{array}{@{}rcl@{}} &&5(RT)^{3/2}\boldsymbol{A}+\sqrt{RT}\boldsymbol{B}=\frac{2(1-Pr)\boldsymbol{q}^{(NS)}}{\tau}, \end{array} $$
(34)

and C is a vector coefficient to be determined.

If the coefficients A,B and C are chosen specifically as

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{A}=\frac{2}{5}\frac{(1-Pr)\boldsymbol{q}^{(NS)}}{\tau(RT)^{3/2}},~ \boldsymbol{B}=\boldsymbol{0},~ \boldsymbol{C}=\frac{2K}{5}\frac{(1-Pr)\boldsymbol{q}^{(NS)}}{\tau(RT)^{1/2}}, \end{array} $$
(35)

then Eq. (34) is satisfied and the two source terms Sg and Sh are given as

$$\begin{array}{@{}rcl@{}} &&S_{g}=\frac{1-Pr}{\tau}\frac{\boldsymbol{c}\cdot\boldsymbol{q}^{(NS)}}{5pRT}\left(\frac{c^{2}}{RT}-5\right)g^{eq}, ~S_{h}=KRTS_{g}. \end{array} $$
(36)

In the Shakhov model, the source term Sf corresponding to the original distribution function f is given by

$$\begin{array}{@{}rcl@{}} &&S_{f}=\frac{f^{Pr}}{\tau}=\frac{1-Pr}{\tau}\frac{\boldsymbol{c}\cdot \boldsymbol{q}^{(NS)}}{5pRT}\left(\frac{c^{2}}{RT}-5\right)f^{eq}. \end{array} $$
(37)

Substitution of Eq. (37) into Eq. (9) yields the same results given in Eq. (36). Therefore, the Shakhov model is consistent with five constraints here.

7.2 The total energy double-distribution-function model

The total energy double-distribution-function model (TEDDF) is originally proposed by Guo et al. [11] and then generalized by Liu et al. [34] to simulate thermal compressible flows. The TEDDF model is briefly introduced as follows. From the original distribution function f, two new distribution functions g and b are introduced, namely, \(g=\int f d\boldsymbol {\zeta }\) and \(b=(1/2)\int \left (\xi ^{2}+\zeta ^{2}\right)f d\boldsymbol {\zeta }\). Therefore, this kinetic system can be expressed as

$$\begin{array}{@{}rcl@{}} &&\frac{\partial g}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla} g+\boldsymbol{a}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}}g =\frac{g^{eq}-g}{\tau_{g}},\\ &&\frac{\partial b}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla} b+\boldsymbol{a}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}}b=\frac{b^{eq}-b}{\tau_{b}}+\left(\boldsymbol{\xi}\cdot\boldsymbol{u}-\frac{1}{2}u^{2}\right)\frac{g-g^{eq}}{\tau_{bg}}+g\boldsymbol{\xi}\cdot\boldsymbol{a}, \end{array} $$
(38)

where the relaxation times are τ=τg=μ/p, τb=τg/Pr and τbg=τbτg/(τgτb). The bulk viscosity of this model is shown to be 2Kμ/[3(K+3)]. The equilibrium beq can be written as beq=(ξ2geq+heq)/2. From Eq. (38), we find that the relation 2b=ξ2g+h holds. The expressions for the hydrodynamic variables are the same as those given in Eqs. (13) and (14). Using g and h, the kinetic system Eqs. (38) can be rewritten as

$$\begin{array}{@{}rcl@{}} &&\frac{\partial g}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla} g+\boldsymbol{a}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}}g=\Omega_{g},\\ &&\frac{\partial h}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla} h+\boldsymbol{a}\cdot\boldsymbol{\nabla}_{\boldsymbol{\xi}}h=\Omega_{h}-(1-Pr)\left(c^{2}\Omega_{g}+\Omega_{h}\right). \end{array} $$
(39)

Therefore, the two source terms are

$$\begin{array}{@{}rcl@{}} &&S_{g}=0,~S_{h}=-(1-Pr)\left(c^{2}\Omega_{g}+\Omega_{h}\right). \end{array} $$
(40)

Eq. (40) indicates that the source terms in the TEDDF model can be expressed in terms of the collision operators.

By noticing that the conservation law for the internal energy, \(\int \left (c^{2}\Omega _{g}+\Omega _{h}\right)d\boldsymbol {\xi }=0\), and the heat flux can be expressed as the moments of the collision operators, \(\boldsymbol {q}=-(1/2)\tau \int \boldsymbol {c}\left (c^{2}\Omega _{g}+\Omega _{h}\right)d\boldsymbol {\xi }\), we find that the five general conditions given by Eqs. (21a) – (21d) are satisfied.

Therefore, we conclude that the TEDDF model is also a special design of the two source terms. Although two relaxation times are used to modify the Prandtl number, the TEDDF model is equivalent to a mesoscopic model with a single relaxation time.

7.3 The Rykov model

The well-known Rykov model for diatomic gases with rotational degrees of freedom is originally obtained by Rykov [12, 13]. Recently, Wu et al. [14] has generalized this model to polyatomic gases. The elastic and non-elastic collision processes are considered respectively in this model. By integrating the particle distribution function f with respect to the rotational energy e, the following two-equation kinetic system can be established.

$$\begin{array}{@{}rcl@{}} &&\frac{\partial f_{0}}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla} f_{0} =\frac{1}{\tau Z}\left(f_{0}^{r}-f_{0}\right) +\frac{1}{\tau}\left(1-\frac{1}{Z}\right)\left(f_{0}^{t}-f_{0}\right),\\ &&\frac{\partial f_{1}}{\partial t}+\boldsymbol{\xi}\cdot\boldsymbol{\nabla} f_{1} =\frac{1}{\tau Z}\left(f_{1}^{r}-f_{1}\right) +\frac{1}{\tau}\left(1-\frac{1}{Z}\right)\left(f_{1}^{t}-f_{1}\right), \end{array} $$
(41)

where the equilibrium distribution functions corresponding to the elastic and nonelastic processes are

$$\begin{array}{@{}rcl@{}} &&f_{0}^{r}=f_{M}(T)\left[1-\omega_{0} \boldsymbol{q}^{t} \cdot \boldsymbol{a}(T)\right], ~f_{0}^{t}=f_{M}(T_{t})\left[1-\boldsymbol{q}^{t}\cdot \boldsymbol{a}(T_{t})\right],\\ &&f_{1}^{r}=RTf_{0}^{r}+\omega_{1}RTf_{M}(T)(1-\delta)\frac{\boldsymbol{q}^{r} \cdot \boldsymbol{c}}{pRT},\\ &&f_{1}^{t}=RT_{r}f_{0}^{t}+RT_{r} f_{M}(T_{t}) (1-\delta) \frac{\boldsymbol{q}^{r}\cdot \boldsymbol{c}}{p_{t} R T_{r}},\\ &&f_{M}(T)=\frac{\rho}{(2\pi RT)^{3/2}}\exp\left(-\frac{c^{2}}{2RT}\right),~ \boldsymbol{a}(T)=\frac{2}{15}\frac{\boldsymbol{c}}{pRT}\left(\frac{5}{2}-\frac{c^{2}}{2RT}\right). \end{array} $$
(42)

Here f0 is the velocity distribution function and f1 is the distribution for rotational energy. \(f_{0}^{t}\) and \(f_{0}^{r}\) denote the equilibrium distributions of the elastic and nonelastic collision processes for f0, respectively. Similarly, \(f_{1}^{t}\) and \(f_{1}^{r}\) denote the equilibrium distributions of the elastic and nonelastic collision processes for f1, respectively. fM is the Maxwellian equilibrium distribution function. Tt is the translational temperature corresponding to the translational degrees of freedom of particles, Tr is the rotational temperature corresponding to the rotational degrees of freedom, and T is the total temperature in the local equilibrium state. qt is the translational heat flux and qr is the rotational heat flux. The total heat flux is decomposed as q=qt+qr.

The relaxation time τ is related to the shear viscosity μ and pressure p through the relation τ=μ/p with p=ρRT. Physically, the relaxation time τ is related to the translational temperature Tt instead of the rotational temperature Tr. Therefore, in analogy to the case of a monatomic gas, the following assumption is used, τ=μt/pt, μt=μ(Tt) and pt=ρRTt, pr=ρRTr. Z indicates the ratio of the total number of translational and rotational collisions to that of rotational collisions. When Z goes to infinity, the Rykov model can reduce to the Shakhov model for monatomic gas without energy exchange between translational and rotational motions. δ=(μt/ρ)D, where D is the gas self-diffusion coefficient. ω0 and ω1 are two parameters which can be selected to achieve proper relaxation of the heat fluxes.

The hydrodynamic variables are defined by the following relationships.

$$\begin{array}{@{}rcl@{}} &&\rho=\int f_{0}d\boldsymbol{\xi},~\rho\boldsymbol{u}=\int \boldsymbol{\xi} f_{0}d\boldsymbol{\xi},\\ &&\frac{3}{2}\rho R T_{t}=\frac{1}{2}\int c^{2} f_{0}d\boldsymbol{\xi},~\rho R T_{r}=\int f_{1}d\boldsymbol{\xi},\\ &&\frac{5}{2}\rho RT=\frac{3}{2} \rho RT_{t} +\rho R T_{r}= \frac{1}{2}\int (c^{2} f_{0}+2f_{1}) d\boldsymbol{\xi},\\ &&\boldsymbol{q}^{t}=\int \frac{1}{2}\boldsymbol{c} c^{2} f_{0} d\boldsymbol{\xi},~ \boldsymbol{q}^{r}=\int \boldsymbol{c} f_{1}d\boldsymbol{\xi},~\boldsymbol{q}=\boldsymbol{q}^{t}+\boldsymbol{q}^{r}. \end{array} $$
(43)

In order to use our general results, we first notice K=2, Λ=0 in this case. Then we introduce two new distribution functions, g=f0 and h=2f1.

Two new collision operators are defined as

$$\begin{array}{@{}rcl@{}} &&\Omega_{g}=\frac{1}{\tau}\left[f_{M}(T)-g\right],~\Omega_{h}=\frac{1}{\tau}\left[2RTf_{M}(T)-h\right]. \end{array} $$
(44)

Two new source terms are given by

$$\begin{array}{@{}rcl@{}} &&S_{g}=\frac{1}{\tau}\left[\frac{1}{Z}f_{0}^{r}+\left(1-\frac{1}{Z}\right)f_{0}^{t}-f_{M}(T)\right],\\ &&S_{h}=\frac{1}{\tau}\left[\frac{2}{Z}f_{1}^{r}+2\left(1-\frac{1}{Z}\right)f_{1}^{t}-2RT f_{M}(T)\right]. \end{array} $$
(45)

The expressions for the hydrodynamic variables in Eq. (43) can be rewritten in terms of g and h, which are found to be the same as Eqs. (13) and (14). From Eqs. (43) and (44), we confirm that the newly defined collision operators, Ωg and Ωh, still satisfy the conservative requirements in Eq. (11). Furthermore, it can be shown that Sg and Sh indeed satisfy five basic requirements in Eqs. (21a) – (21d). The details of proof are provided in Appendix 15 in which we can confirm that the bulk-to-shear viscosity ratio χ is determined by the collision ratio Z through χ=4Z/15 in the continuum limit assuming that both τ and τZ are small in Chapman-Enskog analysis. By contrast, both the S model and TEDDF model give a constant bulk-to-shear viscosity ratio 4/15 for diatomic gas.

In addition, the Prandtl number can be identified as Pr=7Rμ/[2(κt+κr)], where the tranlational and rotational thermal conductivity coefficients κt and κr as well as the total thermal conductivity coefficient are shown to be κt=15Rμt/(4A),κr=Rμt/B and κ=κt+κr, where A=1+0.5(1−ω0)/Z and B=δ+(1−δ)(1−ω1)/Z.

Therefore, we have proved that the Rykov model is compatible with our design in the continuum limit. The Rykov model is constructed from physical viewpoint with a broad range of bulk-to-shear viscosity ratios. In contrast, inverse design approach presented in this paper is directly based on the Chapman-Enskog analysis in order to obtain feasible model for compressible flows. The underlying assumptions for our model are (a) τ is small and (b) τχ(τZ) is small, which originates from the premise of the Chapman-Enskog analysis that the non-equilibrium part and the source terms only serve as small corrections to the equilibrium distribution function. As a result, the bulk-to-shear viscosity ratio should be limited such that the premise for the Chapman-Enskog expansion can be preserved.

8 Implementation of macroscopic hydrodynamic and thermodynamic boundary conditions

When numerically implementing mesoscopic methods based on a model Boltzmann equation, a challenge is to properly determine the unknown distribution functions near a solid boundary, such that the resulting scheme is fully consistent with the NSF system near the boundary. Since the NSF system is derived from the Chapman-Enskog expansion up to O(τ), it follows that the proper implementation of the boundary condition should be based on a consistent application of the Chapman-Enskog expansion up to O(τ). In the literature, this requirement is often not checked and thus not met rigorously, leading to degradation of the accuracy of a mesoscopic method. Furthermore, for thermal or compressible flows, as will be shown below, the implementations of velocity and temperature boundary conditions, at the level of the distribution functions, could be coupled. Source terms could also affect the implementation details. Such fine points are not fully realized in the literature. Below we shall explore the relations between the components of the distribution functions (typically the distribution functions between two opposite particle velocity directions after the particle velocity space is discretized).

By using the relation c=ξu, the expression for G1(ξ),G2(ξ),G3(ξ) and Φ1(ξ) in Eq. (19) can be rewritten in terms of the particle velocity ξ.

$$\begin{array}{@{}rcl@{}} &&G_{1}(\boldsymbol{\xi})=G_{11}(\boldsymbol{\xi})+G_{12}(\boldsymbol{\xi}),\\ &&G_{11}(\boldsymbol{\xi})=\left(\frac{\xi^{2}\boldsymbol{\xi}}{2RT}-\frac{5}{2}\boldsymbol{\xi}+\frac{\boldsymbol{\xi}\cdot\boldsymbol{u} \boldsymbol{u}}{RT}+\frac{u^{2}\boldsymbol{\xi}}{2RT}\right)\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right),\\ &&G_{12}(\boldsymbol{\xi})=\left(-\frac{u^{2} \boldsymbol{u}}{2RT}+\frac{5}{2}\boldsymbol{u}-\frac{\boldsymbol{u}\cdot\boldsymbol{\xi}\boldsymbol{\xi}}{RT}-\frac{\xi^{2}\boldsymbol{u}}{2RT}\right)\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right). \end{array} $$
(46)
$$\begin{array}{@{}rcl@{}} &&G_{2}(\boldsymbol{\xi})=G_{21}(\boldsymbol{\xi})+G_{22}(\boldsymbol{\xi}),\\ &&G_{21}(\boldsymbol{\xi})=2\left(-\frac{\boldsymbol{\xi}\cdot\boldsymbol{S}\cdot\boldsymbol{u}}{RT}+\frac{1}{K+3}\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT}\vartheta\right),\\ &&G_{22}(\boldsymbol{\xi})=\frac{\boldsymbol{\xi}\cdot\boldsymbol{S}\cdot\boldsymbol{\xi}}{RT}+\frac{\boldsymbol{u}\cdot\boldsymbol{S}\cdot\boldsymbol{u}}{RT}-\frac{1}{K+3}\left(\frac{\xi^{2}+u^{2}}{RT}+K\right)\vartheta. \end{array} $$
(47)
$$\begin{array}{@{}rcl@{}} &&G_{3}(\boldsymbol{\xi})=G_{31}(\boldsymbol{\xi})+G_{32}(\boldsymbol{\xi}),\\ &&G_{31}(\boldsymbol{\xi})=\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT}\frac{\Lambda}{C_{v}T},~ G_{32}(\boldsymbol{\xi})=\left(-\frac{\xi^{2}+u^{2}}{2RT}+\frac{3}{2}\right)\frac{\Lambda}{C_{v}T}. \end{array} $$
(48)
$$\begin{array}{@{}rcl@{}} &&\Phi_{1}(\boldsymbol{\xi})=\Phi_{11}(\boldsymbol{\xi})+\Phi_{12}(\boldsymbol{\xi}),\\ &&\Phi_{11}(\boldsymbol{\xi})=\boldsymbol{\xi}\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right),~ \Phi_{12}(\boldsymbol{\xi})=-\boldsymbol{u}\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right)-\frac{2}{K+3}\vartheta-\frac{\Lambda}{C_{v} T}. \end{array} $$
(49)

Obviously, we have Gi1(ξ)=−Gi1(−ξ),Gi2(ξ)=Gi2(−ξ) and Φ11(ξ)=−Φ11(−ξ),Φ12(ξ)=Φ12(−ξ),(i=1,2,3).

From Eq. (20), we have

$$\begin{array}{@{}rcl@{}} &&\phi(\boldsymbol{\xi})=\left(A_{\phi}(\boldsymbol{\xi})-B_{\phi}(\boldsymbol{\xi})\right)\phi^{eq}(\boldsymbol{\xi})+\tau S_{\phi}(\boldsymbol{\xi})+O\left(\tau^{2}\right),\\ &&\phi(-\boldsymbol{\xi})=\left(A_{\phi}(\boldsymbol{\xi})+B_{\phi}(\boldsymbol{\xi})\right)\phi^{eq}(-\boldsymbol{\xi})+\tau S_{\phi}(-\boldsymbol{\xi})+O\left(\tau^{2}\right), \end{array} $$
(50)

where ϕ=g or h. Obviously, the coefficients satisfy the relations Aϕ(ξ)=Aϕ(−ξ) and Bϕ(ξ)=−Bϕ(−ξ). They can be expressed explicitly as follows,

$$\begin{array}{@{}rcl@{}} &&A_{g}(\boldsymbol{\xi})=1-\tau G_{12}(\boldsymbol{\xi})-\tau G_{22}(\boldsymbol{\xi})-\tau G_{32}(\boldsymbol{\xi}),\\ &&B_{g}(\boldsymbol{\xi})=\tau G_{11}(\boldsymbol{\xi})+\tau G_{21}(\boldsymbol{\xi})+\tau G_{31}(\boldsymbol{\xi}),\\ &&A_{h}(\boldsymbol{\xi})=1-\tau G_{12}(\boldsymbol{\xi})-\tau G_{22}(\boldsymbol{\xi})-\tau G_{32}(\boldsymbol{\xi})-\tau\Phi_{12}(\boldsymbol{\xi}),\\ &&B_{h}(\boldsymbol{\xi})=\tau G_{11}(\boldsymbol{\xi})+\tau G_{21}(\boldsymbol{\xi})+\tau G_{31}(\boldsymbol{\xi})+\tau \Phi_{11}(\boldsymbol{\xi}). \end{array} $$
(51)

If the particle distribution function ϕ(−ξ) is already known, then the particle distribution function ϕ(ξ) in the opposite direction can be obtained in the following way. From Eq. (50), we obtain a generalized bounce back scheme

$$\begin{array}{@{}rcl@{}} &&\phi(\boldsymbol{\xi})-\beta\phi(-\boldsymbol{\xi})\\ &&=A_{\phi}(\boldsymbol{\xi})\left(\phi^{eq}(\boldsymbol{\xi})-\beta\phi^{eq}(-\boldsymbol{\xi})\right) -B_{\phi}(\boldsymbol{\xi})\left(\phi^{eq}(\boldsymbol{\xi})+\beta\phi^{eq}(-\boldsymbol{\xi})\right)\\ &&+\tau\left(S_{\phi}(\boldsymbol{\xi})-\beta S_{\phi}(-\boldsymbol{\xi})\right)+O\left(\tau^{2}\right), \end{array} $$
(52)

where β is a coefficient to be determined. Specially, we can choose β=1 or β=−1 in real implementation. For this purpose, we have to evaluate the sum or difference of the equilibriums and source terms.

The sum and difference of the source terms depend on the specific form used. For Sg given in Eq. (24) and Sh given in Eq. (29), we have

$$\begin{array}{@{}rcl@{}} &&S_{g}(\boldsymbol{\xi})+S_{g}(-\boldsymbol{\xi})\\ &&=-2\omega(\boldsymbol{\xi},T_{0})\left[\left(\chi-\frac{2K}{3(K+3)}\right)\frac{p\vartheta}{2RT_{0}} +\frac{p}{2RT_{0}}\frac{\Lambda}{C_{v} T}\right]\left(\frac{\xi^{2}}{RT_{0}}-3\right),\\ &&S_{g}(\boldsymbol{\xi})-S_{g}(-\boldsymbol{\xi})=0, \end{array} $$
(53)
$$\begin{array}{@{}rcl@{}} &&S_{h}(\boldsymbol{\xi})+S_{h}(-\boldsymbol{\xi}) =\omega(\boldsymbol{\xi},T_{0})\left[6\left(\chi-\frac{2K}{3(K+3)}\right)p\vartheta -\frac{4K}{K+3}\rho\Lambda\right],\\ &&S_{h}(\boldsymbol{\xi})-S_{h}(-\boldsymbol{\xi})\\ &&=\omega(\boldsymbol{\xi},T_{0})\left[ \begin{aligned} &\frac{4(1-Pr)\boldsymbol{q}^{(NS)}\cdot\boldsymbol{\xi}}{\tau RT_{0}}-4\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\left(\chi-\frac{2K}{3(K+3)}\right)p\vartheta\\ &-4\frac{K+5}{K+3}\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\rho\Lambda \end{aligned} \right]. \end{array} $$
(54)

Consider the three-dimensional isothermal flow in the incompressible limit with constant temperature T0. The internal degree of freedom is K=0 and the bulk viscosity is μV=0. No thermal cooling function is applied, i.e. Λ=0. The source term Sg=0. Then the Eq. (52) with ϕ=g and β=1 can be simplified as

$$\begin{array}{@{}rcl@{}} &&g(\boldsymbol{\xi})-g(-\boldsymbol{\xi})\\ &&=2\omega(\boldsymbol{\xi},T_{0})\rho\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)+\frac{1}{3}\omega(\boldsymbol{\xi},T_{0})\rho\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)\left[\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)^{2}-3\frac{u^{2}}{RT_{0}}\right] \\ &&+4\omega(\boldsymbol{\xi},T_{0})\rho \frac{\boldsymbol{u}\cdot (\tau\boldsymbol{S})\cdot \boldsymbol{\xi}}{RT_{0}} -2\omega(\boldsymbol{\xi},T_{0})\rho\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)\frac{\boldsymbol{\xi}\cdot(\tau\boldsymbol{S})\cdot\boldsymbol{\xi}}{RT_{0}}\\ &&+2\omega(\boldsymbol{\xi},T_{0})\rho\tau\vartheta\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right) +O\left(\tau Ma^{4}, Ma^{4}, \tau^{2}, \tau^{2} Ma^{2}\right). \end{array} $$
(55)

In the lattice Boltzmann method [20], we first introduce a transformation as

$$\begin{array}{@{}rcl@{}} \tilde{g}(\boldsymbol{x},\boldsymbol{\xi}, t)=g(\boldsymbol{x},\boldsymbol{\xi}, t)+\frac{\Delta t}{2\tau}\left(g(\boldsymbol{x},\boldsymbol{\xi}, t)-g^{eq}(\boldsymbol{x},\boldsymbol{\xi}, t)\right)-\frac{\Delta t}{2}\frac{\boldsymbol{a}\cdot \boldsymbol{c}}{RT_{0}}g^{eq}(\boldsymbol{x},\boldsymbol{\xi}, t), \end{array} $$
(56)

where Δt is the time step.

Next, in order to use the Gauss-Hermite quadrature for the evaluation of the integrals, we introduce another transformation as,

$$\begin{array}{@{}rcl@{}} \tilde{g}\left(\boldsymbol{x},\boldsymbol{\xi}_{\alpha},t\right)=\frac{W(\boldsymbol{\xi}_{\alpha})}{\omega(\boldsymbol{\xi}_{\alpha},T_{0})}\tilde{g}\left(\boldsymbol{x},\boldsymbol{\xi}_{\alpha},t\right), \end{array} $$
(57)

where α denotes the directions of the discrete velocities ξα and Wα denotes the corresponding weight.

After some reorganization, the final result is

$$\begin{array}{@{}rcl@{}} &&\tilde{g}\left(\boldsymbol{x},\boldsymbol{\xi}_{\alpha},t\right)-\tilde{g}\left(\boldsymbol{x},-\boldsymbol{\xi}_{\alpha},t\right)\\ &&=2\rho W_{\alpha}\frac{\left(\boldsymbol{u}-\frac{\Delta t}{2}\boldsymbol{a}\right)\cdot\boldsymbol{\xi}_{\alpha}}{RT_{0}}+\frac{1}{3}\rho W_{\alpha} \left(\frac{\boldsymbol{\xi}_{\alpha}\cdot\boldsymbol{u}}{RT_{0}}\right)\left[\left(\frac{\boldsymbol{\xi}_{\alpha}\cdot\boldsymbol{u}}{RT_{0}}\right)^{2}-3\frac{u^{2}}{RT_{0}}\right]\\ &&+\frac{2\tau+\Delta t}{2\tau}\rho W_{\alpha}\left[ \begin{aligned} &4 \frac{\boldsymbol{u}\cdot (\tau\boldsymbol{S})\cdot \boldsymbol{\xi}_{\alpha}}{RT_{0}} -2\left(\frac{\boldsymbol{\xi}_{\alpha}\cdot\boldsymbol{u}}{RT_{0}}\right)\frac{\boldsymbol{\xi}_{\alpha}\cdot(\tau\boldsymbol{S})\cdot\boldsymbol{\xi}_{\alpha}}{RT_{0}} +2\tau\vartheta\left(\frac{\boldsymbol{\xi}_{\alpha}\cdot\boldsymbol{u}}{RT_{0}}\right) \end{aligned} \right]\\ &&-\frac{\Delta t}{2}\rho W_{\alpha}\frac{\boldsymbol{a}\cdot \boldsymbol{\xi}_{\alpha}}{RT_{0}}\left[\left(\frac{\boldsymbol{\xi}_{\alpha}\cdot\boldsymbol{u}}{RT_{0}}\right)^{2}-\frac{u^{2}}{RT_{0}}\right]+\Delta t\rho W_{\alpha} \frac{\boldsymbol{a}\cdot \boldsymbol{u}}{RT_{0}}\frac{\boldsymbol{\xi}_{\alpha}\cdot \boldsymbol{u}}{RT_{0}}\\ &&+O\left(\tau Ma^{4}, Ma^{4}, \tau^{2}, \tau^{2} Ma^{2}\right). \end{array} $$
(58)

If we only keep the first term in Eq. (58), we can obtain

$$\begin{array}{@{}rcl@{}} \tilde{g}\left(\boldsymbol{x},\boldsymbol{\xi}_{\alpha},t\right)-\tilde{g}\left(\boldsymbol{x},-\boldsymbol{\xi}_{\alpha},t\right)=2\rho W_{\alpha}\frac{\left(\boldsymbol{u}-\frac{\Delta t}{2}\boldsymbol{a}\right)\cdot\boldsymbol{\xi}_{\alpha}}{RT_{0}}+O\left(\tau^{2}, \tau Ma, Ma^{2}\right). \end{array} $$
(59)

We note that the body force enters the implementation of the bounce-back scheme, which is not well documented in the literature. Furthermore, it must be cautioned that Eq. (59) is not fully consistent with the Chapman-Enskog expansion of the NSF system as the O(τ) terms in Eq. (58) are not included. Luckily, in the special case of no-slip boundary uwall=0, the O(τ) terms in Eq. (58) will disappear. Note that the source term and velocity could all enter the implementation of the thermal boundary conditions.

9 High-order structure of the distribution functions

The NSF equations, which are based on the continuum hypothesis, have been widely used in understanding flow behaviors in many natural and engineering problems. However, in some cases such as microchannel flows [35], compressible turbulence [30] and space vehicles in low earth orbits [36], the local Knudsen number may be finite such that the flow may lie in the continuum-transition regime locally. Therefore, the NSF equations are not adequate to capture the finite Knudsen number effect while the Boltzmann equation can describe the flows in all Kn number regimes.

In order to quantitatively estimate the departure from the local thermodynamic equilibrium and study the extended hydrodynamics, the second-order Chapman-Enskog expansion of the particle distribution function is desired, which results in the so-called “Burnett equations” [37]. The Burnett equations have been derived from the original Boltzmann equation by applying the Chapman-Enskog expansion [1] or the Grad’s 13 moment equations [38] by the iteration approach [39]. However, these theoretical results are seldom compared with those using the single-relaxation-time BGK model. In addition, detailed derivations are less reported or the final results are not presented in a general form. In this section, we will derive the structure of the distribution functions up to the order O(τ2). Then, the complete analytical expressions for the viscous stress tensor and the heat flux are obtained in the subsequent sections. Furthermore, by comparing our results with those from Grad’s 13 moment equations, it is found that the mathematical form of the viscous stress tensor and the heat flux can be fully determined in the single-relaxation-time BGK model. The difference from the literature in the coefficients could be attributed to different relaxation rates to the local equilibrium for different moments used in the literature.

By using the NSF equations, we obtain the expression for Dgeq/Dt and Dheq/Dt to the order of O(τ),

$$\begin{array}{@{}rcl@{}} \frac{Dg^{eq}}{Dt}=\left(G+L\right)g^{eq}+O\left(\tau^{2} \right),~ \frac{D h^{eq}}{Dt}=\left(G+\Phi_{1}+L+\Phi_{2}\right)h^{eq}+O\left(\tau^{2}\right), \end{array} $$
(60)

where G and Φ1 have been given above, L and Φ2 are given as

$$\begin{array}{@{}rcl@{}} &&L=\frac{\boldsymbol{c}}{\rho RT}\cdot \left(\boldsymbol{\nabla}\cdot\boldsymbol\sigma^{(NS)}\right)+\left(\frac{c^{2}}{2RT}-\frac{3}{2}\right)\frac{1}{\rho C_{v} T}\left(\boldsymbol\sigma^{(NS)}:\boldsymbol{S}-\boldsymbol{\nabla}\cdot\boldsymbol{q}^{(NS)}\right),\\ &&\Phi_{2}=\frac{1}{\rho C_{v} T}\left(\boldsymbol\sigma^{(NS)}:\boldsymbol{S}-\boldsymbol{\nabla}\cdot\boldsymbol{q}^{(NS)}\right). \end{array} $$
(61)

By applying the Chapman-Enskog expansion, we can obtain the structure of the particle distribution function as

$$\begin{array}{@{}rcl@{}} &&g=g^{eq}-\tau \frac{Dg^{eq}}{Dt}+\tau\frac{D}{Dt}\left(\tau\frac{Dg^{eq}}{Dt}\right)+\tau S_{g}-\tau\frac{D(\tau S_{g})}{Dt}+O\left(\tau^{3}\right)\\ &&=(1-\tau G)g^{eq}+\tau S_{g}\\ &&\hspace{0.4cm}+\left(-\tau L+\tau\frac{D(\tau G)}{Dt}+\tau^{2} G^{2}\right)g^{eq}-\tau\frac{D(\tau S_{g})}{Dt}+O\left(\tau^{3}\right),\\ &&h=h^{eq}-\tau \frac{Dh^{eq}}{Dt}+\tau\frac{D}{Dt}\left(\tau\frac{Dh^{eq}}{Dt}\right)+\tau S_{h}-\tau\frac{D(\tau S_{h})}{Dt}+O\left(\tau^{3}\right)\\ &&=\left(1-\tau G-\tau \Phi_{1}\right)h^{eq}+\tau S_{h} -\tau (L+\Phi_{2})h^{eq}+\tau\frac{D\left[\tau \left(G+\Phi_{1}\right)\right]}{Dt}h^{eq}\\ &&\hspace{0.4cm}+\tau^{2}\left(G+\Phi_{1}\right)^{2}h^{eq}-\tau\frac{D(\tau S_{h})}{Dt}+O\left(\tau^{3}\right). \end{array} $$
(62)

10 Viscous stress tensor up to O(τ 2)

When the local Knudsen number becomes finite, additional contributions from the non-equilibrium part of the distribution function result in the high-order components of the viscous stress tensor. Agarwal et al. [39] and Struchtrup [40] derived the viscous stress tensor up to the second-order for Maxwell molecules from the Grad’s 13 moment equations by the series expansion in terms of the shear viscosity. Chen et al. [7] obtained the expression of viscous stress tensor up to the second-order based on the single-relaxation-time BGK model. However, they mainly focus on the incompressible limit and the terms proportional to the density gradient, the temperature gradient and the velocity divergence have been neglected in their derivation. By making an analogy between the turbulent fluctuations and microscale thermal fluctuations, they show that the Reynolds stress obtained by the BGK-Boltzmann equation has model coefficients similar to some existing turbulence models. They also claimed that the turbulence phenomenon such as the secondary flow structures and rapid distortion processes [41] can be better understood according to the kinetic theory. As an extension of Chen’s work, the complete form of the viscous stress tensor will be derived up to O(τ2) using the single-relaxation-time BGK model considering the internal degree of freedom of molecules. Moreover, this new result will be compared to that obtained by Agarwal et al. [39] for Maxwell molecules.

The general expression of the viscous stress tensor is given as follows.

$$\begin{array}{@{}rcl@{}} &&\boldsymbol\sigma=-\int \boldsymbol{c}\boldsymbol{c}\left(g-g^{eq}\right)d\boldsymbol{\xi}\\ &&=\tau\int \boldsymbol{c}\boldsymbol{c}G g^{eq}d\boldsymbol{\xi}-\tau\int \boldsymbol{c}\boldsymbol{c} S_{g} d\boldsymbol{\xi}+\tau\int \boldsymbol{c}\boldsymbol{c} L g^{eq}\boldsymbol{\xi}-\tau\int \boldsymbol{c}\boldsymbol{c} \frac{D(\tau G)}{Dt}g^{eq}d\boldsymbol{\xi}\\ &&-\int \boldsymbol{c}\boldsymbol{c} \tau^{2} G^{2} g^{eq}d\boldsymbol{\xi}+\tau\int \boldsymbol{c}\boldsymbol{c} \frac{D(\tau S_{g})}{Dt}d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{3}\right). \end{array} $$
(63)

After some computation and simplification, we obtain the complete expression for the viscous stress tensor σ as

$$\begin{array}{@{}rcl@{}} &&\boldsymbol\sigma=\boldsymbol\sigma^{(NS)}-\frac{2\tau}{K+3}\left(\boldsymbol{\nabla}\cdot\boldsymbol{q}^{(NS)}\right)\boldsymbol{I}\\ &&\hspace{0.2cm}-\tau^{2}pRT\left\{ \begin{aligned} & \left(\frac{1}{\tau}\boldsymbol{\nabla}\tau\right)\left(\frac{1}{T}\boldsymbol{\nabla} T\right) +\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\left(\frac{1}{\tau}\boldsymbol{\nabla}\tau\right)\\ &+\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\cdot\left(\frac{1}{\tau}\boldsymbol{\nabla}\tau\right)\boldsymbol{I} +\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\\ &+\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)+\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\boldsymbol{I}\\ &+\left\vert\frac{1}{T}\boldsymbol{\nabla} T\right\vert^{2}\boldsymbol{I} +2\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\left(\frac{1}{T}\boldsymbol{\nabla} T\right) +\left(\frac{1}{T}\boldsymbol{\nabla}^{2} T\right)\boldsymbol{I}+2\left(\frac{1}{T}\boldsymbol{\nabla\nabla}T\right)\\ \end{aligned} \right\}\\ &&\hspace{0.2cm}+\tau^{2} p\left\{ \begin{aligned} &-2\left(\frac{1}{\tau}\frac{d\tau}{dt}\right)\boldsymbol{S}+\frac{2}{K+3} \left(\frac{1}{\tau}\frac{d\tau}{dt}\right)\vartheta\boldsymbol{I} +\frac{2}{K+3}\left[\chi-\frac{2(K+6)}{3(K+3)}\right]\vartheta^{2} \boldsymbol{I}\\ &+\frac{8}{K+3} \vartheta \boldsymbol{S}-2\frac{d\boldsymbol{S}}{dt}+\frac{2}{K+3}\frac{d\vartheta}{dt}\boldsymbol{I}\\ &+\frac{4}{K+3}(\boldsymbol{S}:\boldsymbol{S})\boldsymbol{I}-4 \boldsymbol{S}\cdot\boldsymbol{S}+2 \left(\boldsymbol{S}\cdot\boldsymbol{\Omega}-\boldsymbol{\Omega}\cdot\boldsymbol{S}\right)\\ \end{aligned} \right\}\\ &&\hspace{0.2cm}+\tau^{2} p\left\{ \begin{aligned} &4\frac{\Lambda}{C_{v} T}\boldsymbol{S}+\frac{d}{dt}\left(\frac{\Lambda}{C_{v} T}\right)\boldsymbol{I}+\left(\frac{1}{\tau}\frac{d\tau}{dt}\right)\frac{\Lambda}{C_{v} T}\boldsymbol{I}\\ &-\frac{4}{K+3}\frac{\Lambda}{C_{v} T}\vartheta\boldsymbol{I}-\left(\frac{\Lambda}{C_{v} T}\right)^{2} \boldsymbol{I}\\ \end{aligned} \right\}\\ &&\hspace{0.2cm}+\tau\int \boldsymbol{c}\boldsymbol{c}\frac{D(\tau S_{g})}{Dt}d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{3}\right), \end{array} $$
(64)

where the material derivative of the strain rate tensor dS/dt can be derived from the Euler equations, which reads

$$\begin{array}{@{}rcl@{}} &&\frac{d\boldsymbol{S}}{dt}=-\left(\boldsymbol{S}\cdot\boldsymbol{S}+\boldsymbol{\Omega}\cdot\boldsymbol{\Omega}\right)\\ &&\hspace{0.2cm}+RT\left[\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)+\frac{1}{2}\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\left(\frac{1}{T}\boldsymbol{\nabla} T\right)+\frac{1}{2}\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\right]\\ &&\hspace{0.2cm}-\frac{1}{\rho}\boldsymbol{\nabla}\boldsymbol{\nabla} p +\frac{1}{2}\left(\boldsymbol{\nabla}\boldsymbol{a}+\boldsymbol{\nabla} \boldsymbol{a}^{T}\right)+\boldsymbol{O}(\tau). \end{array} $$
(65)

Furthermore, by setting K=0, χ=0, Λ=0 and neglecting all the terms proportional to the velocity divergence 𝜗, the density gradient ρ and the temperature gradient T, it is found that Sg=0 and the following approximate result obtained by Chen et al. [7] can be reproduced from the Eq. (64), namely,

$$\begin{array}{@{}rcl@{}} \boldsymbol\sigma \approx 2\mu \boldsymbol{S} +\frac{4\mu\tau}{3}\boldsymbol{S}:\boldsymbol{S}\boldsymbol{I}-4\tau^{2}p \boldsymbol{S}\cdot\boldsymbol{S}+2\tau^{2} p \left(\boldsymbol{S}\cdot\boldsymbol{\Omega}-\boldsymbol{\Omega}\cdot\boldsymbol{S}\right)-2\tau p\frac{d\left(\tau\boldsymbol{S}\right)}{dt}. \end{array} $$
(66)

For Maxwell molecules of which the shear viscosity is linearly proportional to the temperature, we have dμ/dT=μ/T. Therefore, the relations (1/τ)τ=−(1/ρ)ρ and (1/τ)dτ/dt=−(1/ρ)dρ/dt hold. Using the continuity equation, it follows that (1/τ)dτ/dt=𝜗. Then, using our notations, the results obtained by Agarwal et al. [39] can be rewritten as

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{\sigma}=2\mu\left(\boldsymbol{S}-\frac{1}{3}\vartheta\boldsymbol{I}\right)-\frac{10}{9}\frac{\mu^{2}}{p}\vartheta^{2}\boldsymbol{I}+\frac{2}{3}\frac{\mu^{2}}{p}\boldsymbol{S}\vartheta\\ &&+4\frac{\mu^{2}}{p}(\boldsymbol{S}:\boldsymbol{S})\boldsymbol{I}-4\frac{\mu^{2}}{p}\boldsymbol{S}\cdot\boldsymbol{S} +2\frac{\mu^{2}}{p}(\boldsymbol{S}\cdot\boldsymbol{\Omega}-\boldsymbol{\Omega}\cdot\boldsymbol{S})-2\frac{\mu^{2}}{p}\frac{d}{dt}\left(\boldsymbol{S}-\frac{1}{3}\vartheta\boldsymbol{I}\right)\\ &&+\frac{\mu^{2}}{\rho T^{2}}\vert\boldsymbol{\nabla}T\vert^{2}\boldsymbol{I}-3\frac{\mu^{2}}{\rho T^{2}}\left(\boldsymbol{\nabla}T\right)\left(\boldsymbol{\nabla}T\right)+\frac{\mu^{2}}{\rho T}\left(\nabla^{2}{T}\right)\boldsymbol{I} -3\frac{\mu^{2}}{\rho T}\boldsymbol{\nabla}\boldsymbol{\nabla}T+\boldsymbol{O}\left(\tau^{3}\right). \end{array} $$
(67)

Moreover, Eq. (64) can be simplified as

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{\sigma}=2\mu\left(\boldsymbol{S}-\frac{1}{3}\vartheta\boldsymbol{I}\right)-\frac{2}{9}\frac{\mu^{2}}{p}\vartheta^{2}\boldsymbol{I}+\frac{2}{3}\frac{\mu^{2}}{p}\boldsymbol{S}\vartheta\\ &&+\frac{4}{3}\frac{\mu^{2}}{p}(\boldsymbol{S}:\boldsymbol{S})\boldsymbol{I}-4\frac{\mu^{2}}{p}\boldsymbol{S}\cdot\boldsymbol{S} +2\frac{\mu^{2}}{p}(\boldsymbol{S}\cdot\boldsymbol{\Omega}-\boldsymbol{\Omega}\cdot\boldsymbol{S})-2\frac{\mu^{2}}{p}\frac{d}{dt}\left(\boldsymbol{S}-\frac{1}{3}\vartheta\boldsymbol{I}\right)\\ &&+\frac{3}{2}\frac{\mu^{2}}{\rho T^{2}}\vert\boldsymbol{\nabla}T\vert^{2}\boldsymbol{I}-2\frac{\mu^{2}}{\rho T^{2}}\left(\boldsymbol{\nabla}T\right)\left(\boldsymbol{\nabla}T\right)+\frac{3}{2}\frac{\mu^{2}}{\rho T}\left(\nabla^{2}{T}\right)\boldsymbol{I} -2\frac{\mu^{2}}{\rho T}\boldsymbol{\nabla}\boldsymbol{\nabla}T+\boldsymbol{O}\left(\tau^{3}\right). \end{array} $$
(68)

From Eqs. (67) and (68), we observe that these two expressions share identical mathematical form up to the order O(τ2) except for values of some coefficients. It is observed that the nonlinear terms S·S,S·ΩΩ·S and S𝜗 are exactly identical. Besides, the material derivative term d(S−(1/3)𝜗I)/dt is also the same and the terms related to the temperature gradient and temperature diffusion are also very close to each other. The sign of corresponding coefficients is also the same, which implies that the negative or positive contribution to the viscous stress tensor can be qualitatively determined based on the BGK collision model. Moreover, it is found that the viscous stress tensor can be changed by the body force effect included in the material derivative term d(S−(1/3)𝜗I)/dt at the second-order expansion but not at the first-order. Therefore, we conclude that although the BGK model only uses single relaxation time to characterize the relaxation process to the local equilibrium without considering rigorous collision interaction details, all the dominant terms in the viscous stress tensor can be recovered compared to those obtained from the Grad’s 13 moment equations. These expressions could potentially give guidelines to the functional form for Reynolds stress modeling in compressible turbulence.

11 Heat flux up to O(τ 2)

Based on the second-order Chapman-Enskog expansion of the distribution functions, the analytical expression for the heat flux q is given by

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{q}=\frac{1}{2}\int \boldsymbol{c} \left(c^{2} g+h\right)d\boldsymbol{\xi}\\ &&=\frac{1}{2}\int \boldsymbol{c} c^{2}\left[g^{eq}-\tau Gg^{eq}+\tau S_{g}\right]d\boldsymbol{\xi}+\frac{1}{2}\int \boldsymbol{c} \left[h^{eq}-\tau(G+\Phi_{1})h^{eq}+\tau S_{h}\right]d\boldsymbol{\xi}\\ &&+\frac{1}{2}\int \boldsymbol{c} c^{2} \left[-\tau Lg^{eq}+\tau\frac{D(\tau G)}{Dt}g^{eq}+\tau^{2}G^{2}g^{eq}-\tau\frac{D\left(\tau S_{g}\right)}{Dt}\right]d\boldsymbol{\xi}\\ &&+\frac{1}{2}\int \boldsymbol{c}\left\{ \begin{aligned} &-\tau(L+\Phi_{2})h^{eq}+\tau\frac{D(\tau G)}{Dt}h^{eq}+\tau^{2} (G+\Phi_{1})^{2} h^{eq}\\ &+\tau\frac{D(\tau \Phi_{1})}{Dt}h^{eq}-\tau\frac{D\left(\tau S_{h}\right)}{Dt} \end{aligned} \right\}d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{3}\right). \end{array} $$
(69)

Noticing Eqs. (21a)–(21d) and Eqs. (19), (61), all the integrals in Eq. (69) can be evaluated term by term. After some reorganization, we obtain

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{q}=\boldsymbol{q}^{(NS)}-\frac{1}{2}(K+5)\tau RT\boldsymbol{\nabla}\cdot\boldsymbol\sigma^{(NS)}\\ &&+\tau^{2}\rho(RT)^{2}\left\{ \begin{aligned} &\frac{1}{2}(K+5)\left(\frac{1}{\tau}\frac{d\tau}{dt}\right)\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\\ &+(K+7) \boldsymbol{S}\cdot\left(\frac{1}{\tau}\boldsymbol{\nabla}\tau\right)-\frac{K+7}{K+3}\vartheta\left(\frac{1}{\tau}\boldsymbol{\nabla}\tau\right) \end{aligned} \right\}\\ &&+\tau^{2}\rho (RT)^{2}\left\{ \begin{aligned} &2\left(K+7\right) \boldsymbol{S}\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right)-\frac{3K+19}{K+3}\vartheta\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\\ &+(K+5)\boldsymbol{S}\cdot\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)-\frac{K+5}{K+3}\vartheta\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\\ & +(K+7) \boldsymbol{\nabla}\cdot\boldsymbol{S}-2\frac{K+6}{K+3}\boldsymbol{\nabla}\vartheta+\left(K+5\right) \boldsymbol{\Omega}\cdot\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\\ \end{aligned} \right\}\\ &&+\tau^{2}\rho (RT)^{2}\left\{ \begin{aligned} &-\frac{3}{2}(K+5)\boldsymbol{\nabla}\left(\frac{\Lambda}{C_{v} T}\right)-(K+5) \left(\frac{1}{\tau}\boldsymbol{\nabla}\tau\right)\frac{\Lambda}{C_{v} T}\\ &-\frac{1}{2}(K+5) \left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\frac{\Lambda}{C_{v} T} -\frac{5}{2}(K+5)\left(\frac{1}{T}\boldsymbol{\nabla} T\right)\frac{\Lambda}{C_{v} T} \end{aligned} \right\}\\ &&-\frac{1}{2}\int \boldsymbol{c} c^{2}\tau\frac{D}{Dt}(\tau S_{g})d\boldsymbol{\xi}-\frac{1}{2}\int \boldsymbol{c}\tau\frac{D}{Dt}\left(\tau S_{h}\right)d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{3}\right), \end{array} $$
(70)

where the time and spatial derivatives of the relaxation time are given by

$$\begin{array}{@{}rcl@{}} &&\frac{1}{\tau}\frac{d\tau}{dt}=\left(\gamma-\frac{2}{K+3}\frac{T}{\mu}\frac{d\mu}{dT}\right)\vartheta +\left(1-\frac{T}{\mu}\frac{d\mu}{dT}\right)\frac{\Lambda}{C_{v} T}+O(\tau),\\ &&\frac{1}{\tau}\boldsymbol{\nabla}\tau=-\left(1-\frac{T}{\mu}\frac{d\mu}{dT}\right)\frac{1}{T}\boldsymbol{\nabla} T-\frac{1}{\rho}\boldsymbol{\nabla}\rho. \end{array} $$
(71)

The results in Eq. (70) are briefly discussed here. The first term is the Fourier’s law. The second term is determined by the divergence of the viscous stress tensor. The third term is caused by the variation of the particle relaxation time in both space and time. The fourth term is composed of the coupling terms between the strain rate, rotation rate, temperature gradient and density gradient as well as the divergence of the strain rate. The fifth term represents the contributions from the terms relevant to the thermal energy source. The last two integrals depend on the specific form of the source terms Sg and Sh used in different models.

Similar to what we have done for the viscous stress tensor, by setting K=0, χ=0, Λ=0 and Sg=0, a comparison would also be performed for heat flux for Maxwell molecules. The result obtained by Agarwal et al. [39] can be reformulated as

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{q}=-\frac{15}{4}\mu R\boldsymbol{\nabla} T+15\frac{\mu^{2}}{\rho}\boldsymbol{S}\cdot\left(\frac{1}{T}\boldsymbol{\nabla}T\right) -3\frac{\mu^{2}}{\rho}\boldsymbol{S}\cdot\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\\ &&-\frac{25}{8}\frac{\mu^{2}}{\rho}\vartheta\left(\frac{1}{T}\boldsymbol{\nabla}T\right) +\frac{\mu^{2}}{\rho}\vartheta\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)+3\frac{\mu^{2}}{\rho}\boldsymbol{\nabla}\cdot\boldsymbol{S}\\ &&-\frac{19}{4}\frac{\mu^{2}}{\rho}\boldsymbol{\nabla}\vartheta +\frac{45}{4}\frac{\mu^{2}}{\rho}\boldsymbol{\Omega}\cdot\left(\frac{1}{T}\boldsymbol{\nabla}T\right)+\boldsymbol{O}\left(\tau^{3}\right). \end{array} $$
(72)

Correspondingly, Eq. (70) can be simplified as

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{q}=-\frac{15}{4}\mu R\boldsymbol{\nabla} T+9\frac{\mu^{2}}{\rho}\boldsymbol{S}\cdot\left(\frac{1}{T}\boldsymbol{\nabla}T\right) -2\frac{\mu^{2}}{\rho}\boldsymbol{S}\cdot\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)\\ &&-\frac{13}{6}\frac{\mu^{2}}{\rho}\vartheta\left(\frac{1}{T}\boldsymbol{\nabla}T\right) +\frac{2}{3}\frac{\mu^{2}}{\rho}\vartheta\left(\frac{1}{\rho}\boldsymbol{\nabla}\rho\right)+2\frac{\mu^{2}}{\rho}\boldsymbol{\nabla}\cdot\boldsymbol{S}\\ &&-\frac{7}{3}\frac{\mu^{2}}{\rho}\boldsymbol{\nabla}\vartheta +5\frac{\mu^{2}}{\rho}\boldsymbol{\Omega}\cdot\left(\frac{1}{T}\boldsymbol{\nabla}T\right)-\frac{1}{2}\int\boldsymbol{c}\tau\frac{D(\tau S_{h})}{Dt}d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{3}\right). \end{array} $$
(73)

Again, Eqs. (72) and (73) share the same mathematical form and the same sign for each contribution up to the order O(τ2). In our model, Sh is mainly designed to modify the Prandtl number and thermal energy source. Note that we keep the term relevant to Sh in Eq. (73) but it can be evaluated once the specific form of Sh is given.

12 Conclusions

In this paper, an inverse design approach of mesoscopic models for compressible flows in continuum or near-continuum regime has been explored based on the Chapman-Enskog analysis. The design began with a model Boltzmann equation in a high dimensional phase space and with an undetermined source term. Then two reduced model Boltzmann equations in seven dimensional phase space are introduced, each containing a source term. First, it is found that there are many possible ways to design the source terms in order to recover the NSF system in the continuum limit, as long as five newly-derived requirements for the two source terms are met. These source terms allow for flexible Prandtl number, bulk-to-shear viscosity ratio, and a thermal energy source/sink term.

Second, based on the Hermite expansion, we have provided one design for the two source terms. This newly introduced model has been utilized to simulate decaying compressible isotropic turbulence [32] and forced compressible isotropic turbulence [33], achieving global turbulence statistics in excellent agreement with those based on solving the NSF system [29, 30]. Our model can be viewed as a Boltzmann-equation based mesoscopic solver for compressible flows that solves the NSF system.

Third, three well accepted kinetic models, namely, the Shakhov model, the total energy double-distribution-function model, and the Rykov model, have been shown to be compatible with five basic constraints derived from the Chapman-Enskog analysis. For Rykov model, translational and rotational temperatures are introduced, which allows for larger bulk-to-shear viscosity ratio compared to other models. Similarly, in our model, the bulk-to-shear viscosity ratio should be limited such that the premise of Chapman-Enskog expansion is satisfied.

Furthermore, by applying the first-order Chapman-Enskog expansion to the distribution functions, we discuss the structures of the distribution functions and the implementation of bounce back boundary conditions. These results can be used to improve the implementation of hydrodynamic boundary conditions in terms of the distribution functions, namely, constructing the missing distributions from the known distribution near a solid boundary, in both laminar and turbulent flows.

Finally, although present model is mainly designed for continuum or near-continuum compressible flows, with BGK approximation, we derive the complete analytical expressions for the viscous stress tensor and the heat flux based on the second-order Chapman-Enskog expansion of the distribution functions, generalizing the previous results in the incompressible limit. These new results have been compared with those obtained from Grad’s 13 moment equations, which demonstrates that the final structure of the viscous stress tensor and heat flux can be fully determined by the single-relaxation-time BGK model except for differences in some coefficients. We believe that high-order effects in compressible turbulence could be partially captured by the BGK-Boltzmann equation, which certainly deserves further investigation. It would be desirable to explore underlying physics associated with the second-order terms especially in compressible turbulence, in the future using DNS data. The second-order terms may also provide a way to assess the difference between NSF flows and the flows governed by the model Boltzmann equation.

13 Hermite polynomials and hermite expansion

The n-th order Hermite polynomial is defined by [42, 43],

$$\begin{array}{@{}rcl@{}} \mathscr{H}^{(n)}({\boldsymbol{\xi},T_{0}})\equiv\left(\sqrt{RT_{0}}\right)^{n}\frac{(-1)^{n}}{\omega({\boldsymbol{\xi},T_{0}})}\boldsymbol{\nabla}_{\boldsymbol{\xi}}^{n}\omega({\boldsymbol{\xi},T_{0}}),~n=0,1,2,\cdots, \end{array} $$
(74)

where \(\boldsymbol {\nabla }_{\boldsymbol {\xi }}^{n}=\boldsymbol {\nabla }_{\boldsymbol {\xi }}\boldsymbol {\nabla }_{\boldsymbol {\xi }}\cdots \boldsymbol {\nabla }_{\boldsymbol {\xi }}\) implies that \(\mathscr {H}^{(n)}({\boldsymbol {\xi },T_{0}})\) is a symmetrical tensor of rank-n. The weighting function is ω(ξ,T0)=(2πRT0)D/2 exp(−ξ2/2RT0).

From the 4th-order Hermite expansion of geq and the 2nd-order Hermite expansion of heq [32], we obtain

$$\begin{array}{@{}rcl@{}} &&g^{eq}(\boldsymbol{\xi})+g^{eq}(-\boldsymbol{\xi})=2\rho\omega(\boldsymbol{\xi},T_{0})\times\\ &&\left\{ \begin{aligned} &1+\frac{1}{2}\left[\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)^{2}-\frac{u^{2}}{RT_{0}}+\left(\theta-1\right)\left(\frac{\xi^{2}}{RT_{0}}-3\right)\right]\\ &+\frac{1}{24}\left[ \begin{aligned} &\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)^{4}-6\frac{u^{2}}{RT_{0}}\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)^{2}+3\left(\frac{u^{2}}{RT_{0}}\right)^{2}\\ &+6\left(\theta-1 \right)\left[\left(\frac{\xi^{2}}{RT_{0}}-7\right)\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)^{2}+\frac{u^{2}}{RT_{0}}\left(5-\frac{\xi^{2}}{RT_{0}}\right)\right]\\ &+3\left(\theta-1\right)^{2}\left[\left(\frac{\xi^{2}}{RT_{0}}\right)^{2}-10\frac{\xi^{2}}{RT_{0}}+15\right] \end{aligned} \right] \end{aligned} \right\}\\ &&+\text{high\ order\ terms}\\ &&g^{eq}(\boldsymbol{\xi})-g^{eq}(-\boldsymbol{\xi})=2\rho\omega(\boldsymbol{\xi},T_{0})\times\\ &&\left\{{\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}}+\frac{1}{6}\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}} \left[\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)^{2}-3\frac{u^{2}}{RT_{0}}+3\left(\theta-1\right)\left(\frac{\xi^{2}}{RT_{0}}-5\right)\right]\right\}\\ &&+\text{high~order~terms}. \end{array} $$
(75)
$$\begin{array}{@{}rcl@{}} &&h^{eq}(\boldsymbol{\xi})+h^{eq}(-\boldsymbol{\xi})=2\omega(\boldsymbol{\xi},T_{0})K\rho RT\times\\ &&\left\{1+\frac{1}{2}\left[\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right)^{2}-\frac{u^{2}}{RT_{0}}+\left(\theta-1\right)\left(\frac{\xi^{2}}{RT_{0}}-3\right)\right] \right\}+\text{high order terms},\\ &&h^{eq}(\boldsymbol{\xi})-h^{eq}(-\boldsymbol{\xi}) =2\omega(\boldsymbol{\xi},T_{0})K\rho RT\left(\frac{\boldsymbol{\xi}\cdot\boldsymbol{u}}{RT_{0}}\right) +\text{high order terms}, \end{array} $$
(76)

where θ=T/T0 is the normalized temperature.

14 Derivations of the requirements for the source terms

The Euler equations can be obtained by assuming that g=geq+O(τ) when evaluating the viscous stress and the heat flux. This leads to σO(τ) and qO(τ). Therefore, the Euler equations are

$$\begin{array}{@{}rcl@{}} &&\frac{\partial \rho}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u})=0,~ \rho\left(\frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\cdot\boldsymbol{\nabla} \boldsymbol{u}\right)=-\boldsymbol{\nabla} p+\rho \boldsymbol{a}+\boldsymbol{O}(\tau),\\ &&\rho C_{v}\frac{dT}{dt}=-p\boldsymbol{\nabla}\cdot\boldsymbol{u}-\rho \Lambda+O(\tau). \end{array} $$
(77)

Taking the first-order moment of the Boltzmann equation for g, we have

$$\begin{array}{@{}rcl@{}} &&\frac{\partial (\rho \boldsymbol{u})}{\partial t}+\boldsymbol{\nabla}\cdot(\rho\boldsymbol{u}\boldsymbol{u})=-\boldsymbol{\nabla} p+ \boldsymbol{\nabla}\cdot\boldsymbol\sigma+\rho \boldsymbol{a}+\int\boldsymbol{\xi} S_{g}d\boldsymbol{\xi}. \end{array} $$
(78)

By using the Eq. (20), we can make a closure of the viscous stress,

$$\begin{array}{@{}rcl@{}} &&\boldsymbol\sigma=-\int\boldsymbol{cc}\left(g-g^{eq}\right)d\boldsymbol{\xi}\\ &&=\tau\int \boldsymbol{cc} G_{2} g^{eq}d\boldsymbol{\xi}+\tau\int \boldsymbol{cc} G_{3} g^{eq}d\boldsymbol{\xi} -\tau\int \boldsymbol{cc}S_{g}d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{2}\right). \end{array} $$
(79)

All the integrals in the RHS of Eq. (79) can be evaluated term by term,

$$\begin{array}{@{}rcl@{}} &&\tau\int \boldsymbol{cc}G_{2} g^{eq}d\boldsymbol{\xi}=\tau p\left[2\boldsymbol{S}-\frac{2}{K+3}\vartheta\boldsymbol{I}\right],~ \tau\int \boldsymbol{cc} G_{3}g^{eq}d\boldsymbol{\xi}=-\frac{p\tau\Lambda}{C_{v} T}\boldsymbol{I}. \end{array} $$
(80)

Substitution of Eq. (80) into Eq. (79) yields

$$\begin{array}{@{}rcl@{}} &&\boldsymbol\sigma=2\mu\left(\boldsymbol{S}-\frac{1}{3}\vartheta\boldsymbol{I}\right)+\mu_{V}\vartheta\boldsymbol{I}-\left(\chi-\frac{2K}{3(K+3)}\right)\mu\vartheta\boldsymbol{I}\\ &&\hspace{0.5cm}-\frac{\tau p\Lambda}{C_{v}T}\boldsymbol{I}-\tau\int \boldsymbol{c}\boldsymbol{c}S_{g}d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{2}\right). \end{array} $$
(81)

Hence, in order to recover the Newtonian constitutive law (see Eq. (16)), the Eq. (21b) must be satisfied.

Similarly, combining the second-order moment of the Boltzmann equation for g and the zeroth-order moment of the Boltzmann equation for h yields the following equation,

$$\begin{array}{@{}rcl@{}} \frac{\partial(\rho E)}{\partial t}+\boldsymbol{\nabla}\cdot\left(\boldsymbol{q}+\rho E \boldsymbol{u}+p \boldsymbol{u}-\boldsymbol\sigma\cdot\boldsymbol{u}\right) =\rho \boldsymbol{a}\cdot\boldsymbol{u}+\frac{1}{2}\int\left(\xi^{2}S_{g}+S_{h}\right)d\boldsymbol{\xi}. \end{array} $$
(82)

By using the Eqs. (20), we can make a closure of the heat flux term,

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{q}=\frac{1}{2}\int \boldsymbol{c}\left(c^{2} g+h\right)d\boldsymbol{\xi}\\ &&=-\frac{1}{2}\tau\int \boldsymbol{c} c^{2} G g^{eq}d\boldsymbol{\xi} -\frac{1}{2}\tau\int \boldsymbol{c} G h^{eq}d\boldsymbol{\xi}-\frac{1}{2}\tau\int \boldsymbol{c}\Phi_{1} h^{eq}d\boldsymbol{\xi}\\ &&\hspace{0.4cm}+\frac{1}{2}\tau\int \boldsymbol{c} c^{2} S_{g} d\boldsymbol{\xi} +\frac{1}{2}\tau\int \boldsymbol{c} S_{h} d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{2}\right). \end{array} $$
(83)

Because of \(\int \boldsymbol {c} c^{2} G_{1} g^{eq}d\boldsymbol {\xi }=5\rho (RT)^{2}\left ((1/T)\boldsymbol {\nabla } T\right)\) and \(\int \boldsymbol {c} c^{2} G_{2} g^{eq}d\boldsymbol {\xi }=\int \boldsymbol {c} c^{2} G_{3} g^{eq}d\boldsymbol {\xi }=\boldsymbol {0}\), we obtain

$$\begin{array}{@{}rcl@{}} -\frac{1}{2}\tau\int \boldsymbol{c} c^{2} G g^{eq}d\boldsymbol{\xi}=-\frac{5}{2}\tau\rho(RT)^{2}\left(\frac{1}{T}\boldsymbol{\nabla} T\right). \end{array} $$
(84)

Because of \(\int \boldsymbol {c} G_{1} g^{eq}d\boldsymbol {\xi }=\int \boldsymbol {c} G_{2} g^{eq}d\boldsymbol {\xi }=\int \boldsymbol {c} G_{3} g^{eq}d\boldsymbol {\xi }=\boldsymbol {0}\), we have \(\tau \int \boldsymbol {c} G h^{eq}d\boldsymbol {\xi }=\boldsymbol {0}\). Further, we have

$$\begin{array}{@{}rcl@{}} -\frac{1}{2}\tau\int \boldsymbol{c} \Phi_{1} h^{eq}d\boldsymbol{\xi}=-\frac{k}{2}\tau\rho(RT)^{2}\left(\frac{1}{T}\boldsymbol{\nabla} T\right). \end{array} $$
(85)

By substituting of Eqs. (84) and (85) into Eq. (83), we obtain

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{q}=\boldsymbol{q}^{(NS)}-\left(1-Pr\right)\boldsymbol{q}^{(NS)}+\frac{1}{2}\tau\int \boldsymbol{c} c^{2} S_{g} d\boldsymbol{\xi}+\frac{1}{2}\tau\int \boldsymbol{c} S_{h} d\boldsymbol{\xi}+\boldsymbol{O}\left(\tau^{2}\right). \end{array} $$
(86)

Therefore, we have derived the fifth requirement for the source term in Eq. (21d).

15 Details in the derivations of the Rykov model

Here we prove that the two source terms in Eq. (45) should satisfy the five general requirements. From the Euler equations for the Rykov model, the time derivative for the translational temperature Tt and rotational temperature Tr are

$$\begin{array}{@{}rcl@{}} &&\frac{\partial T_{t}}{\partial t}=-\boldsymbol{u}\cdot\boldsymbol{\nabla} T_{t}-\frac{2}{3}T_{t}\vartheta +\frac{1}{\tau Z}\left(T-T_{t}\right)+O(\tau),\\ &&\frac{\partial T_{r}}{\partial t}=-\boldsymbol{u}\cdot\boldsymbol{\nabla} T_{r}+\frac{1}{\tau Z}\left(T-T_{r}\right)+O(\tau). \end{array} $$
(87)

The first requirement for the source term is satisfied because of

$$\begin{array}{@{}rcl@{}} \int S_{g}d\boldsymbol{\xi}=\frac{1}{\tau} \left[\frac{1}{Z}\int f_{0}^{r}d\boldsymbol{\xi}+\left(1-\frac{1}{Z}\right)\int f_{0}^{t}d\boldsymbol{\xi}-\int f_{M}(T)d\boldsymbol{\xi}\right]=0, \end{array} $$
(88)

The second requirement for the source term is satisfied because of

$$\begin{array}{@{}rcl@{}} \int\boldsymbol{\xi} S_{g}d\boldsymbol{\xi}=\frac{1}{\tau} \left[\frac{1}{Z}\int \boldsymbol{c} f_{0}^{r}d\boldsymbol{\xi}+\left(1-\frac{1}{Z}\right)\int \boldsymbol{c} f_{0}^{t}d\boldsymbol{\xi}-\int \boldsymbol{c} f_{M}(T)d\boldsymbol{\xi}\right]=\boldsymbol{0}. \end{array} $$
(89)

According to Eqs. (88) and (89), we obtain \(\int \boldsymbol {\xi \xi } S_{g} d\boldsymbol {\xi }=\int \boldsymbol {cc}S_{g} d\boldsymbol {\xi }\).

Since \(\int \boldsymbol {cc} f_{0}^{r} d\boldsymbol {\xi }=p\boldsymbol {I}\), \(\int \boldsymbol {cc} f_{0}^{t} d\boldsymbol {\xi }=p_{t} \boldsymbol {I}\) and \(\int \boldsymbol {c} \boldsymbol {c} f_{M}(T) d\boldsymbol {\xi }=p\boldsymbol {I}\), therefore,

$$\begin{array}{@{}rcl@{}} \int \boldsymbol{cc} S_{g}d\boldsymbol{\xi}=\frac{1}{\tau} \left[ \frac{1}{Z}p\boldsymbol{I}+\left(1-\frac{1}{Z}\right)p_{t} \boldsymbol{I}-p\boldsymbol{I}\right]=-\frac{1}{\tau}\left(1-\frac{1}{Z}\right)\left(p-p_{t}\right)\boldsymbol{I}. \end{array} $$
(90)

From Eq. (87), we have

$$\begin{array}{@{}rcl@{}} &&T-T_{t} =\frac{2}{5}\left(T_{r}-T_{t}\right) =\frac{4}{15}\tau ZT_{t}\vartheta+O(\tau^{2}), \end{array} $$
(91)

Therefore, we obtain

$$\begin{array}{@{}rcl@{}} &&p-p_{t}=\rho R(T-T_{t})=\frac{4}{15}Z\mu_{t}\vartheta+O(\tau^{2}). \end{array} $$
(92)

Substituting Eq. (92) into Eq. (90), we have

$$\begin{array}{@{}rcl@{}} &&\int \boldsymbol{\xi\xi} S_{g}d\boldsymbol{\xi}=\int \boldsymbol{cc} S_{g}d\boldsymbol{\xi}=-\left(\frac{4}{15}Z-\frac{4}{15}\right)p\vartheta\boldsymbol{I}+O(\tau). \end{array} $$
(93)

Therefore, the third requirement for the source term is proved. The ratio of bulk-to-shear viscosity χ=4Z/15 in Chapman-Enskog analysis.

From Eq. (90), we have

$$\begin{array}{@{}rcl@{}} &&\int \xi^{2} S_{g} d\boldsymbol{\xi}=\int c^{2} S_{g}d\boldsymbol{\xi}=-\frac{3}{\tau}\left(1-\frac{1}{Z}\right)\left(p-p_{t}\right). \end{array} $$
(94)

Using the relation \(p-p_{r}=-\frac {3}{2}\left (p-p_{t}\right)\), we have

$$\begin{array}{@{}rcl@{}} &&\int S_{h} d\boldsymbol{\xi}=\frac{2}{\tau}\left[\frac{1}{Z}\int f_{1}^{r}d\boldsymbol{\xi}+\left(1-\frac{1}{Z}\right)\int f_{1}^{t}d\boldsymbol{\xi}-RT\int f_{M}(T)d\boldsymbol{\xi}\right]\\ &&=\frac{2}{\tau}\left[\frac{1}{Z}p+\left(1-\frac{1}{Z}\right)p_{r}-p\right] =\frac{3}{\tau}\left(1-\frac{1}{Z}\right)\left(p-p_{t}\right). \end{array} $$
(95)

From Eqs. (94) and (95), the fourth requirement is satisfied.

We note that the following integrals can be carried out directly.

$$\begin{array}{@{}rcl@{}} &&\int \boldsymbol{c} f_{1}^{t}d\boldsymbol{\xi}=\left(1-\delta\right)\boldsymbol{q}^{r},~\int \boldsymbol{c}f_{1}^{r}d\boldsymbol{\xi}=\omega_{1}(1-\delta)\boldsymbol{q}^{r},~\int \boldsymbol{c} f_{M}(T)d\boldsymbol{\xi}=\boldsymbol{0},\\ &&\int \boldsymbol{c} c^{2} f_{0}^{t}d\boldsymbol{\xi}=\frac{2}{3}\boldsymbol{q}^{t},~\int \boldsymbol{c} c^{2} f_{0}^{r}d\boldsymbol{\xi}=\frac{2}{3}\omega_{0} \boldsymbol{q}^{t},~\int \boldsymbol{c} c^{2} f_{M}(T)d\boldsymbol{\xi}=\boldsymbol{0}. \end{array} $$
(96)

Hence, we have

$$\begin{array}{@{}rcl@{}} &&\int \boldsymbol{c} S_{h}d\boldsymbol{\xi}=\frac{2}{\tau}\left[\boldsymbol{q}^{r}-\left(\delta+\frac{1}{Z}\left(1-\omega_{1}\right)\left(1-\delta\right) \right)\boldsymbol{q}^{r}\right],\\ &&\int \boldsymbol{c} c^{2} S_{g} d\boldsymbol{\xi}=\frac{2}{\tau}\left[\boldsymbol{q}^{t}-\frac{2}{3}\left(1+0.5\frac{1-\omega_{0}}{Z}\right)\boldsymbol{q}^{t}\right]. \end{array} $$
(97)

By applying the Chapman-Enskog expansion, we can prove that the heat fluxes are given by

$$\begin{array}{@{}rcl@{}} &&\boldsymbol{q}^{t}=-\kappa^{t}\boldsymbol{\nabla} T_{t}+\boldsymbol{O}(\tau^{2})=-\kappa^{t}\boldsymbol{\nabla} T+\boldsymbol{O}(\tau^{2}),\\ &&\boldsymbol{q}^{r}=-\kappa^{r}\boldsymbol{\nabla} T_{r}+\boldsymbol{O}(\tau^{2})=-\kappa^{r}\boldsymbol{\nabla} T+\boldsymbol{O}(\tau^{2}),\\ &&\boldsymbol{q}=-\kappa\boldsymbol{\nabla} T+\boldsymbol{O}(\tau^{2}), \end{array} $$
(98)

Combining Eqs. (97) and (98) gives

$$\begin{array}{@{}rcl@{}} &&\int \boldsymbol{c} S_{h}d\boldsymbol{\xi}+\int \boldsymbol{c} c^{2} S_{g} d\boldsymbol{\xi}=\frac{2}{\tau}\left(1-Pr\right)\boldsymbol{q}+\boldsymbol{O}(\tau), \end{array} $$
(99)

where the Prandtl number is Pr=μCp/κ=7Rμ/2κ. Therefore, the fifth requirement is satisfied.

Availability of data and materials

Not applicable.

Abbreviations

NSF:

Navier-Stokes-Fourier

CFD:

Computational fluid dynamics

BGK:

Bhatnagar-Gross-Krook

S:

Shakhov

ES:

Ellipsoidal statistical

IEDDF:

Internal energy double-distribution-function

TEDDF:

Total energy double-distribution-function

R:

Rykov

LBM:

Lattice Boltzmann method

GKS:

Gas kinetic scheme

UGKS:

Unified gas kinetic scheme

DUGKS:

Discrete unified gas kinetic scheme

EOS:

Equation of state

DNS:

Direct numerical simulation

References

  1. Chapman S, Cowling TG (1970) The mathematical theory of non-uniform gases. Cambridge University Press, Cambridge, UK.

    MATH  Google Scholar 

  2. Wang-Chang CS, Uhlenbeck GE (1951) Transport phenomena in polyatomic gases. University of Michigan Engineering Research Rept. No. CM-681.

  3. Gross EP, Jackson EA (1959) Kinetic models and the linearized Boltzmann equation. Phys Fluids 2:432–441. https://doi.org/10.1063/1.1724415.

    Article  MathSciNet  Google Scholar 

  4. Hanson FB, Morse TF (1967) Kinetic models for a gas with internal structure. Phys Fluids 10:345–353. https://doi.org/10.1063/1.1762114.

    Article  Google Scholar 

  5. Bhatnagar PL, Gross EP, Krook MK (1954) A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems. Phys Rev E 94:511–525. https://doi.org/10.1103/PhysRev.94.511.

    Article  Google Scholar 

  6. Chen H, Kandasamy S, Orsazg S, Shock R, Succi S, Yakhot V (2003) Extended Boltzmann kinetic equation for turbulent flows. Science 301:633–636. https://doi.org/10.1126/science.1085048.

    Article  Google Scholar 

  7. Chen H, Orsazg SA, Staroselsky IY, Succi S (2004) Expanded analogy between Boltzmann kinetic theory of fluids and turbulence. J Fluid Mech 519:301–314. https://doi.org/10.1017/S0022112004001211.

    Article  MathSciNet  Google Scholar 

  8. Shakhov EM (1968) Generalization of the Krook kinetic relaxation equation. Fluid Dyn 3:95–96. https://doi.org/10.1007/bf01029546.

    Article  Google Scholar 

  9. Holway LH (1966) New statistical models for kinetic theory: methods of construction. Phys Fluids 9:1658–1673. https://doi.org/10.1063/1.1761920.

    Article  Google Scholar 

  10. He X, Chen S, Doolen GD (1998) A novel thermal model for the lattice Boltzmann method in incompressible limit. J Comput Phys 1:282–300. https://doi.org/10.1006/jcph.1998.6057.

    Article  MathSciNet  Google Scholar 

  11. Guo ZL, Zheng C, Shi B, Zhao TS (2007) Thermal lattice Boltzmann equation for low Mach number flows: Decoupling model. Phys Rev E 75:036704. https://doi.org/10.1103/PhysRevE.75.036704.

    Article  Google Scholar 

  12. Rykov VA (1975) A model kinetic equation for a gas with rotational degrees of freedom. Fluid Dyn 10:959–966. https://doi.org/10.1007/BF01023275.

    Article  Google Scholar 

  13. Rykov VA (1978) Macroscopic description of the motions of a gas with rotational degrees of freedom. Fluid Dyn 13:144–147. https://doi.org/10.1007/BF01094479.

    Article  Google Scholar 

  14. Wu L, White C, Scanlon TJ, Zhang JMRY (2015) A kinetic model of the Boltzmann equation for non-vibrating polyatomic gases. J Fluid Mech 763:24–50. https://doi.org/10.1017/jfm.2014.632.

    Article  MathSciNet  Google Scholar 

  15. Chen S, Xu K, Cai Q (2015) A comparison and unification of ellipsoidal statistical and Shakhov BGK models. Adv Appl Math Mech 7:245–266. https://doi.org/10.4208/aamm.2014.m559.

    Article  MathSciNet  Google Scholar 

  16. Woods LC (1993) An introduction to the kinetic theory of gases and magnetoplasmas. Oxford University Press, Oxford, UK.

    MATH  Google Scholar 

  17. Wang P, Wang L-P, Guo ZL (2016) Comparison of the lattice Boltzmann equation and discrete unified gas-kinetic scheme methods for direct numerical simulation of decaying turbulent flows. Phys Rev E 94:043304. https://doi.org/10.1103/PhysRevE.94.043304.

    Article  MathSciNet  Google Scholar 

  18. Bo YT, Wang P, Guo ZL, Wang L-P (2017) DUGKS simulations of three-dimensional Taylor-Green vortex flow and turbulent channel flow. Comput Fluids 155:9–21. https://doi.org/10.1016/j.compfluid.2017.03.007.

    Article  MathSciNet  Google Scholar 

  19. Yang Z, Zhong C, Zhuo C (2019) Phase-field method based on discrete unified gas-kinetic scheme for large-density-ratio two-phase flows. Phys Rev E 99:043302. https://doi.org/10.1103/physreve.99.043302.

    Article  MathSciNet  Google Scholar 

  20. Chen S, Doolen GD (1998) Lattice Boltzmann method for fluid flows. Annu Rev Fluid Mech 30:329–364. https://doi.org/10.1146/annurev.fluid.30.1.329.

    Article  MathSciNet  Google Scholar 

  21. Xu K (2001) A gas-kinetic BGK scheme for the Navier-Stokes equations and its connection with artificial dissipation and Godunov method. J Comput Phys 171:289–335. https://doi.org/10.1006/jcph.2001.6790.

    Article  MathSciNet  Google Scholar 

  22. Xu K, Huang J-C (2010) A unified gas-kinetic scheme for continuum and rarefied flows. J Comput Phys 229:7747–7764. https://doi.org/10.1016/j.jcp.2010.06.032.

    Article  MathSciNet  Google Scholar 

  23. Guo ZL, Xu K, Wang R (2013) Discrete unified gas kinetic scheme for all Knudsen number flows: low-speed isothermal case. Phys Rev E 88:033305. https://doi.org/10.1103/physreve.88.033305.

    Article  Google Scholar 

  24. Guo ZL, Xu K, Wang RJ (2015) Discrete unified gas kinetic scheme for all Knudsen number flows. II. Thermal compressible case. Phys Rev E 91:033313. https://doi.org/10.1103/physreve.91.033313.

    Article  Google Scholar 

  25. Liu S, Zhong C (2014) Investigation of the kinetic model equations. Phys Rev E 89:033306. https://doi.org/10.1103/PhysRevE.89.033306.

    Article  Google Scholar 

  26. Peng C, Geneva N, Guo ZL, Wang L-P (2017) Direct numerical simulation of turbulent pipe flow using the lattice Boltzmann method. J Comput Phys 357:16–42. https://doi.org/10.1016/j.jcp.2017.11.040.

    Article  MathSciNet  Google Scholar 

  27. Frapolli N, Chikatamarla SS, Karlin IV (2015) Entropic lattice Boltzmann model for compressible flows. Phys Rev E 92:061301. https://doi.org/10.1103/physreve.92.061301.

    Article  Google Scholar 

  28. Zhu L, Guo Z, Xu K (2016) Discrete unified gas kinetic scheme on unstructed meshes. Comput Fluids 127:211–225. https://doi.org/10.1016/j.compfluid.2016.01.006.

    Article  MathSciNet  Google Scholar 

  29. Samtaney R, Pullin DI, Kosović B (2001) Direct numerical simulation of decaying compressible turbulence and shocklet statistics. Phys Fluids 13:1415–1430. https://doi.org/10.1063/1.1355682.

    Article  Google Scholar 

  30. Wang JC, Wang L-P, Xiao ZL, Shi YP, Chen SY (2010) A hybrid numerical simulation of isotropic compressible turbulence. J Comput Phys 229:5257–5279. https://doi.org/10.1016/j.jcp.2010.03.042.

    Article  Google Scholar 

  31. Sutherland W (2009) The viscosity of gases and molecular force. Philos Mag 36:507–531. https://doi.org/10.1080/14786449308620508.

    Article  Google Scholar 

  32. Chen T, Wen X, Wang L-P, Guo Z, Wang J, Chen S (2020) Simulation of three-dimensional compressible decaying isotropic turbulence using a redesigned discrete unified gas kinetic scheme. Phys Fluids 32:125104. http://doi.org/10.1063/5.0029424.

    Article  Google Scholar 

  33. Chen T, Wen X, Wang L-P, Guo Z, Wang J, Chen S (2020) Simulation of three-dimensional forced compressible isotropic turbulence using a redesigned discrete unified gas kinetic scheme. J Comput Phys.

  34. Liu H, Kong M, Chen Q, Zheng L, Cao Y (2018) Coupled discrete unified gas kinetic scheme for the thermal compressible flows in all Knudsen number regimes. Phys Rev E 98:053310. https://doi.org/10.1103/PhysRevE.98.053310.

    Article  MathSciNet  Google Scholar 

  35. Xu K, Li Z (2004) Microchannel flow in slip regime: gas-kinetic BGK-Burnett solutions. J Fluid Mech 513:87–110. https://doi.org/10.1017/S0022112004009826.

    Article  MathSciNet  Google Scholar 

  36. Ivanov MS, Gimelshein SF (1998) Computational hypersonic rarefied flows. Annu Rev Fluid Mech 30:469–505. https://doi.org/10.1146/annurev.fluid.30.1.469.

    Article  MathSciNet  Google Scholar 

  37. Burnett D (1935) The distribution of velocities and mean motion in a slight nonuniform gas. Proc Lond Math Soc 39:385. https://doi.org/10.1112/plms/s2-39.1.385.

    Article  Google Scholar 

  38. Grad H (1949) On the kinetic theory of rarefied gases. Commun Pur Appl Math 2:331–407. https://doi.org/10.1002/cpa.3160020403.

    Article  MathSciNet  Google Scholar 

  39. Agarwal RK, Yun K-Y, Balakrishnan R (2001) Beyond Navier-Stokes: Burnett equations for flows in the continuum-transition regime. Phys Fluids 13:3061–3085. http://doi.org/10.1063/1.1397256.

    Article  Google Scholar 

  40. Struchrup H (2004) Some remarks on the equations of Burnett and Grad In: Transport in Transition Regimes, 265–276. https://doi.org/10.1007/978-1-4613-0017-5_17.

  41. Yakhot V, Orszag SA, Thangam S, Gatski TB, Speziale CG (1992) Development of turbulence models for shear flows by a double expansion technique. Phys Fluids 4:1510–1520. https://doi.org/10.1063/1.858424.

    Article  MathSciNet  Google Scholar 

  42. Grad H (1949) Note on N-dimensional hermite polynomials. Commun Pur Appl Math 2:325–330. https://doi.org/10.1002/cpa.3160020402.

    Article  MathSciNet  Google Scholar 

  43. Shan X, Yuan X-F, Chen H (2006) Kinetic theory representation of hydrodynamics: A way beyond the Navier-Stokes equation. J Fluid Mech 550:413–441. https://doi.org/10.1017/S0022112005008153.

    Article  MathSciNet  Google Scholar 

Download references

Acknowledgements

Computing resources are provided by the Center for Computational Science and Engineering of Southern University of Science and Technology and by National Center for Atmospheric Research (CISL-UDEL0001).

Funding

This work has been supported by the U.S. National Science Foundation (CNS-1513031, CBET-1706130), the National Natural Science Foundation of China (91852205, 91741101 & 11961131006), the National Numerical Wind Tunnel program, Guangdong Provincial Key Laboratory of Turbulence Research and Applications (2019B21203001), and Shenzhen Science & Technology Program (Grant No. KQTD20180411143441009).

Author information

Authors and Affiliations

Authors

Contributions

Tao Chen developed the derivations, drafted and edited the manuscript. Lian-Ping Wang conceptualized the methodology and reviewed the derivations, edited the manuscript, provided suggestions for improving the manuscript, acquired funding for this research, and served as the corresponding author. Jun Lai checked the derivations. Shiyi Chen participated in funding acquisition and supervised this research. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Lian-Ping Wang.

Ethics declarations

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

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

Rights and permissions

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

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Chen, T., Wang, LP., Lai, J. et al. Inverse design of mesoscopic models for compressible flow using the Chapman-Enskog analysis. Adv. Aerodyn. 3, 5 (2021). https://doi.org/10.1186/s42774-020-00059-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s42774-020-00059-2

Keywords