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

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.


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-distributionfunction (TEDDF) model [11], and the Rykov (R) model [12][13][14]. 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 nonequilibrium 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 highand 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 [17][18][19][26][27][28].
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 A, we include the details on the Hermite polynomials and Hermite expansion. Appendix B contains the derivations of the requirements for the two reduced source terms. Appendix C documents some details on the Rykov model.

The reduced model Boltzmann system with source terms
The Boltzmann equation with an additional source term can be expressed as where f (x, ξ , t) is the particle distribution function, x = (x 1 , ..., x 3 ) 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 = f eq − f /τ . τ = μ/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. S f 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] 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 where = C v T is the internal energy per unit mass, C v 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 f eq . C v and the specific heat at constant pressure C p 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 C v = (K + 3)R/2 and C p = (K + 5)R/2, which implies that the specific heat ratio and thus the Prandtl number are γ = C p /C v = (K + 5)/(K + 3) and Pr = μC p /κ, 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 and, by comparing the energy moment of the Boltzmann equation with the macroscopic energy equation, the heat flux q can be determined as The physical conservative requirements can be expressed through the moments of the collision operator f , which reads Therefore, provided that the mass conservation and the momentum conservation laws are observed, we have the following basic requirements for the source term 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 = fdζ and h = ζ 2 fdζ . Therefore, the two model Boltzmann equations residing in lower dimensional space can be obtained In Eq. (8), the collision operators and source terms are where the equilibrium distribution functions g eq and h eq are , h eq = KRTg eq .
Based on Eq. (6), the conservation laws can be recasted in terms of the collision operators g and h , as follows, From Eq. (7), the two reduced source terms must satisfy the following requirements In addition, from Eq. (3), we find that the conservative variables can be evaluated as Moreover, from Eqs. (4) and (5), the viscous stress σ and the heat flux q become

Notations and conventions
For convenience, two time derivatives are introduced 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 = ∇u T + ∇u /2 is the strain rate tensor and = ∇u T − ∇u /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, 4 The structure of the particle distribution functions In the continuum limit, the relaxation time τ , when normalized by the acoustic time scale l 0 /c 0 , is proportional to the Knudsen number, where l 0 is a system length scale and c 0 is the speed of the sound at a reference temperature T 0 . 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 g eq 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 g eq as 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 B to replace the time derivatives of the hydrodynamic variables with spatial derivatives in ∂g eq /∂t, we obtain the expression for Dg eq /Dt and Dh eq /Dt to the leading order, as where G = G 1 + G 2 + G 3 . The coefficients are given explicitly as 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 (20)

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, S g and S h . The requirements are given as follows and the details for their derivations are included in Appendix B. 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 The third requirement is used to modify the bulk viscosity in the viscous stress term, and it is The fourth requirement follows from the energy equation and it is The fifth requirement is expressed as 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

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 ξξξ S g dξ = 0. Then, using Eq. (21a), can also be written as where T 0 is a reference temperature and the velocity ξ are scaled with √ RT 0 . Therefore, by using Eq. (23) and keeping the Hermite expansion (see Appendix A) of the source term S g up to the third-order, we obtain From Eqs. (21a)-(21b), we can obtain Combination of Eqs. (21c), (21d) and (25) yields Therefore, we obtain Eqs. (21b) and (21c) together imply that Combination of Eqs. (27) and (28) 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 secondorder 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 q = 2τ 2τ + tPr 1 2 c c 2g +h dξ and q = 2τ 2τ +( t/2)Pr 1 2 c c 2ḡ +h dξ after the distribution functionsg,h andḡ,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.

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 ]. 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 S g dξ = 0, cS g dξ = 0, ccS g dξ = 0, S h dξ = 0, By using Eq. (30), the source term S g can be assumed as where ω(c, T) is the peculiar-velocity-based weighting function, a (3) (x, t) is the coefficient and H (3) (c, T) is the third-order Hermite polynomial. Note that a (3) T)dξ is also symmetrical with respect to the components of c because H (3) (c, T) remains unchanged under the permutation operation, the simplest way is to assume that the coefficient a (3) takes the form a is a vector coefficient to be determined. Then, Eq. (31) gives Moreover, it can be shown that cc 2 S g dc = 5A(RT) 3/2 . Similarly, the source term S h can be designed as Therefore, the coefficient A and B should satisfy the following relation, and C is a vector coefficient to be determined. If the coefficients A, B and C are chosen specifically as then Eq. (34) is satisfied and the two source terms S g and S h are given as In the Shakhov model, the source term S f corresponding to the original distribution function f is given by 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.

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 = fdζ and b = (1/2) ξ 2 + ζ 2 fdζ . Therefore, this kinetic system can be expressed as 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 b eq can be written as b eq = ξ 2 g eq + h eq /2. From Eq. (38), we find that the relation 2b = ξ 2 g + 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 ∂g ∂t Therefore, the two source terms are 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, c 2 g + h dξ = 0, and the heat flux can be expressed as the moments of the collision operators, q = −(1/2)τ c c 2 g + h dξ , 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.

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.
where the equilibrium distribution functions corresponding to the elastic and nonelastic processes are , a(T) = 2 15 c pRT Here f 0 is the velocity distribution function and f 1 is the distribution for rotational energy. f t 0 and f r 0 denote the equilibrium distributions of the elastic and nonelastic collision processes for f 0 , respectively. Similarly, f t 1 and f r 1 denote the equilibrium distributions of the elastic and nonelastic collision processes for f 1 , respectively. f M is the Maxwellian equilibrium distribution function. T t is the translational temperature corresponding to the translational degrees of freedom of particles, T r is the rotational temperature corresponding to the rotational degrees of freedom, and T is the total temperature in the local equilibrium state. q t is the translational heat flux and q r is the rotational heat flux. The total heat flux is decomposed as q = q t + q r .
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 T t instead of the rotational temperature T r . Therefore, in analogy to the case of a monatomic gas, the following assumption is used, τ = μ t /p t , μ t = μ(T t ) and p t = ρRT t , p r = ρRT r . 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.
In order to use our general results, we first notice K = 2, = 0 in this case. Then we introduce two new distribution functions, g = f 0 and h = 2f 1 . Two new collision operators are defined as Two new source terms are given by 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 S g and S h indeed satisfy five basic requirements in Eqs. (21a) -(21d). The details of proof are provided in Appendix C 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 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.

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 G 1 (ξ ), G 2 (ξ ), G 3 (ξ ) and 1 (ξ ) in Eq. (19) can be rewritten in terms of the particle velocity ξ .
where φ = g or h. Obviously, the coefficients satisfy the relations A φ (ξ ) = A φ (−ξ ) and B φ (ξ ) = −B φ (−ξ ). They can be expressed explicitly as follows, 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 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 S g given in Eq. (24) and S h given in Eq. (29), we have Consider the three-dimensional isothermal flow in the incompressible limit with constant temperature T 0 . 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 S g = 0. Then the Eq. (52) with φ = g and β = 1 can be simplified as In the lattice Boltzmann method [20], we first introduce a transformation as 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, where α denotes the directions of the discrete velocities ξ α and W α denotes the corresponding weight. After some reorganization, the final result is If we only keep the first term in Eq. (58), we can obtaiñ 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 u wall = 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.

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 singlerelaxation-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 Dg eq /Dt and Dh eq /Dt to the order of O(τ ), where G and 1 have been given above, L and 2 are given as By applying the Chapman-Enskog expansion, we can obtain the structure of the particle distribution function as

Viscous stress tensor up to O τ 2
When the local Knudsen number becomes finite, additional contributions from the nonequilibrium 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 singlerelaxation-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.
After some computation and simplification, we obtain the complete expression for the viscous stress tensor σ as where the material derivative of the strain rate tensor dS/dt can be derived from the Euler equations, which reads 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 S g = 0 and the following approximate result obtained by Chen et al. [7] can be reproduced from the Eq. (64), namely, from the Grad's 13 moment equations. These expressions could potentially give guidelines to the functional form for Reynolds stress modeling in compressible turbulence.

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 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 where the time and spatial derivatives of the relaxation time are given by 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.

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 bulkto-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.