Simulation of three-dimensional compressible decaying isotropic turbulence using a redesigned discrete unified gas kinetic scheme

In this paper, a new mesoscopic approach with both the adjustable Prandtl number and the ratio of bulk to shear viscosity has been developed to simulate three-dimensional compressible decaying homogeneous isotropic turbulence under the framework of discrete unified gas kinetic scheme (DUGKS). In the new approach, two reduced model Boltzmann equations with newly designed source terms are solved. In the continuum limit, the Navier–Stokes–Fourier system can be recovered by applying the Chapman–Enskog analysis. A three-dimensional DUGKS code has been developed, incorporating the fifth-order weighted essentially non-oscillatory scheme to better reconstruct the particle distribution functions at the cell interfaces. In addition, a new lattice velocity model with 77 discrete particle velocities is applied to ensure that the accuracy of the Gauss–Hermite quadrature is up to the ninth-order, and as such, the heat flux can be accurately evaluated. To validate our code, we simulate two cases with different initial turbulent Mach numbers and Taylor microscale Reynolds numbers. The simulation results converge with the increase in resolution and agree well with the results from the literature. As a direct application of our DUGKS, we briefly study the influence of bulk viscosity on turbulence statistics and flow structures. Our results show that the DUGKS is a reliable tool for simulating compressible decaying isotropic turbulence at low and moderate turbulent Mach numbers. More parametric studies are needed in the future to further explore the full capabilities of this specific mesoscopic method. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0029424., s

only treat the continuum flows correctly and cannot handle local flow regimes with high non-equilibrium effects. In addition, it is not feasible to derive a simple system of closed partial differential equations such as NSF equations for flows with locally varying Kn numbers. On the other hand, conceptually, the Boltzmann equation can not only describe the continuum regime but can also deal with non-continuum flows. This is one of the reasons that motivated the developments of mesoscopic computational fluid dynamics (CFD) methods based on the Boltzmann equations during the last three decades.
There are mainly two numerical approaches to solve the Boltzmann equation. The first one is called the direct simulation Monte Carlo (DSMC) method, 4 which is mainly used to simulate highspeed rarefied gas flows. However, the particle transport and collision processes are decoupled in the DSMC method, which implies that the grid spacing and the time step should be smaller than the mean free path and the relaxation time (equal to the ratio of the shear viscosity to the hydrodynamic pressure), respectively. In the local regions of nearly continuum flow in a multiscale physical problem, the computational cost for DSMC will increase enormously. In addition, for flows with low speed and small temperature variation, the statistical noise in the DSMC method will be on the order of the flow variables and would thus contaminate the simulation results. 5 The second approach that is used to solve the Boltzmann equation is the Discrete Velocity Methods (DVMs). [6][7][8] Most DVMs are essentially single-scale methods with particle collision and transport processes decoupled, which leads to the requirement that the grid spacing and time step should also be restricted by the mean free path and the relaxation time. The well-known lattice Boltzmann method (LBM) can be viewed as a special type of the discrete velocity formulation. However, in the LBM, the particle transport and collision are properly coupled. The LBM is mainly used to simulate incompressible or weakly compressible continuum flows under the low Mach number approximation. [9][10][11] Recently, a unified gas kinetic scheme (UGKS) [12][13][14][15] has been developed to simulate all Knudsen number flows. The UGKS solves the hydrodynamic equations and the Boltzmann equation simultaneously. The evolution of kinetic variables and the evolution of hydrodynamic variables are coupled through the collision term and the cell interface fluxes in the UGKS. The UGKS can be viewed as an extension of the gas kinetic scheme (GKS), 16 which is mainly used to simulate continuum flows. A kinetic scheme is an asymptotic preserving (AP) 17 scheme if (a) the time step is not restricted by the particle collision time in the continuum regime and (b) it can recover the limiting hydrodynamic equations in the continuum limit. In the UGKS, the particle transport and collision processes are coupled due to the implicit treatment of the collision term such that the time step is no longer limited by the relaxation time. The AP property makes the UGKS a true multiscale scheme for simulating different flow regimes.
Guo et al. 17,18 proposed another unified gas kinetic scheme, namely, the discrete unified gas kinetic scheme (DUGKS) by combining the advantages of LBM and UGKS. Compared to the UGKS, the DUGKS is a simpler multiscale kinetic method. In the UGKS, the reconstruction of the cell-interface particle distribution function is based on the analytical solution of the Boltzmann equation along the characteristic line. However, in the DUGKS, the flux across the cell interface is evaluated from the particle distribution functions at the half time step, and the particle distribution functions are obtained by using the characteristic method instead of using the analytical solution used in the UGKS. In the DUGKS, the Boltzmann equation is solved using an accurate finite-volume formulation coupled tightly with the particle transport and collision processes. Thus, the time step in the DUGKS is not limited by the particle collision time for all Kn numbers. Compared to the LBMs, the DUGKS can more easily incorporate irregular meshes and incorporate different discrete velocity models as the space grid is no longer coupled with the lattice velocity set. In addition, the DUGKS can ensure a second-order accuracy in both space and time. Due to the need of the reconstruction of the particle distribution function, the DUGKS is slightly more dissipative but with a better numerical stability than the LBM. The DUGKS has been applied to simulate single-phase decaying homogeneous isotropic turbulence, 19 wall-bounded turbulent channel flow, 20 3D lid-driven cavity flow, 21 microchannel gas flows, 22 and immiscible two-phase flows. 23, 24 Liu et al. 25 combined the total energy double-distributionfunction (TEDDF) model 26 with the DUGKS to simulate thermal compressible flows in all flow regimes. For compressible flows, the Cartesian discrete velocity space may require a much larger degree of freedoms, and a higher order of the quadrature may also be needed; these will influence the accuracy of hydrodynamic variables and numerical instability. Recently, Chen et al. 27 combined the unstructured discrete velocity space 28 with the DUGKS in 2D. The unstructured velocity mesh adjusts the density of discrete velocities according to the local relative magnitude of the distribution function in the velocity space. Although it may be hard to optimize, it is usually better than a uniform velocity mesh. However, the DUGKS has not yet been applied to simulate three-dimensional compressible turbulence.
Compressible turbulence is of fundamental importance to many natural and engineering applications, such as solar winds, high-temperature reactive flows, and design of aircrafts. With the development of computational techniques, direct numerical simulation (DNS) has been applied as an independent research tool to probe complex flow structures in both incompressible turbulence and compressible turbulence. For compressible flows, shocklets may develop randomly due to local compressibility, causing local discontinuities as the turbulent Mach number increases, which can be a great challenge to the DNS.
The previous relevant DNS studies of compressible turbulence are briefly summarized here. Samtaney et al. 29 studied the decaying compressible turbulence and shocklet statistics at the initial turbulent Mach numbers of Mat 0 = 0.1-0.5 and the Taylor-Reynolds number of Re λ0 = O(50-100) using a tenth-order compact finite difference scheme. A shocklet extraction algorithm is also developed to extract and quantify the shocklet statistics from the DNS data. Honein and Moin 30 presented a numerical scheme with high entropy conservation and numerical stability to simulate compressible isotropic turbulence at the turbulent Mach number of Mat 0 = 0.3. Frapolli et al. 31 overcame the low Mach-number limit of the traditional LBM and proposed an entropic lattice Boltzmann model (ELBM) to simulate compressible transonic and supersonic flows using a D3Q7 3 lattice velocity model with temperature-dependent weights. In their work, they first simulated the flow around a NACA0012 airfoil immersed in a supersonic flow at Ma = 1.4 and Re = 3 × 10 6 . Then, they considered two  38 proposed a high-order gas kinetic scheme (HGKS) to simulate supersonic CDHIT. Their simulation covers the moderate Reynold numbers ranging from 10 to 72 and the Mach numbers ranging from 0.5 to 1.2. They also investigated the turbulent Mach number effect and the Reynolds number effect on flow statistical quantities. Two relevant comments can be made here in terms of DNS of compressible flows. From a physical standpoint, the local Knudsen number Kn may not be small when the local Mach number is high enough to induce local shocklets. In other words, the local flow state is not uniform in the continuum flow regime. Therefore, the NSF system may fail to describe all local non-equilibrium effects, while the DUGKS can capture the physics in all flow regimes. In the hypersonic argon flow past over a sphere, the local Mach number in the wake region can still reach around 1.5. 39 From a numerical perspective, the NSF system needs a low-order treatment to capture shocklets and, at the same time, high-order treatment to capture turbulence fluctuations. Furthermore, the NSF system handles waves related to the hydrodynamics, while the DUGKS, instead, handles movements of particles, not the waves directly. The wave dynamics is obtained as the moments in the particle velocity space. We believe that the discontinuity at the particle level is much weaker than the discontinuity in hydrodynamics. For this reason, as a multiscale method, it may be stated that the DUGKS is not required to incorporate drastically different treatments for different flow regimes.

Physics of Fluids
In this paper, we apply the DUGKS combined with a fifth-order WENO scheme [40][41][42][43] to simulate CDHIT. A newly designed BGK model is utilized, allowing the variation of both the Prandtl number and the ratio of bulk to shear viscosity. The NSF system can be recovered in the continuum limit by applying the Chapman-Enskog expansion. 44 The rest of this paper is organized as follows. In Sec. II, the original BGK-Shakhov model is briefly summarized, and a new model with an adjustable Prandtl number and a viscosity ratio is introduced. Through the Chapman-Enskog analysis, the NSF system for compressible flows can be recovered, and the relevant dimensionless parameters are briefly summarized. In Sec. III, we present the implementation details of the DUGKS under our altered design. Important statistical quantities about the CDHIT are reviewed in Sec. IV. A three-dimensional DUGKS code is developed to simulate CDHIT, and the results are discussed in Sec. V. Two simulation cases are considered, and the results are carefully compared with the literature data to validate our approach. Different WENO schemes and interfacial treatments are also compared to examine the sensitivity of flow statistics on implementation details. Finally, we investigate the effect of bulk to shear viscosity ratio on the statistics of turbulence by using our DUGKS approach. Major conclusions are summarized in Sec. VII, along with some additional discussions. In Appendixes A-E, we include the details on Hermite polynomials and our newly designed novel E 9 3,77 discrete particle velocity model (Appendix A), Hermite expansion of the equilibrium, the implementation of the fifth-order WENO scheme (Appendix B), the Chapman-Enskog analysis of our redesigned model (Appendix C), explicit expressions of two reduced distribution functions (Appendix D). The evolution equation of mean square velocity divergence is also derived in detail (Appendix E).
By assuming that the particle motion in the (η, ζ) subspace is at local equilibrium, the Shakhov equilibrium distribution function is expressed as where f eq is the Maxwellian distribution function, ρ is the density of fluid, R is the specific gas constant, T is the temperature, c = ξ − u is the particle thermal-fluctuation velocity with respect to the hydrodynamic velocity u, Pr is the Prandtl number, and q is the heat flux. The extra term f Pr is designed to adjust the heat flux and thus the Prandtl number. Without f Pr , the Prandtl number would be 1.
The equilibrium distribution implies that the ideal-gas equation of state (EOS), i.e., p = ρRT, is satisfied. The conservative variables are defined as the moments of the particle distribution function, Here, ρE is referred to as the total energy and ρe = ρCvT is the internal energy per unit volume. All relations in Eq. (3) remain valid if f is replaced by f S . With the above-mentioned equilibrium distribution, it can be shown that the specific heat capacity at constant volume Cv and that at constant pressure Cp are related to K and R as Cv = (3 + K)R/2 and Cp = (5 + K)R/2. Therefore, the specific heat ratio γ = Cp/Cv = (5 + K)/(3 + K).
The viscous stress tensor σ and the heat flux q are determined by The shear viscosity is determined by the intermolecular interactions and molecular thermal motions. For ideal gases, the shear viscosity increases with temperature, which introduces an additional effect of thermal field on the hydrodynamic velocity field. Two wellknown models are widely used in the existing literature. The first one is called the hard-sphere (HS) or variable hard-sphere (VHS) model, which takes the form of power law. 2,29 The second one is known as the Sutherland's law, 46 which for air can be written as It has been verified from the experimental data that the HS model has a maximum relative error of about 5% for extreme cases of T/T 0 → 0.55 and T/T 0 → 3, while the Sutherland's law has a maximum relative error of 2.0% at T/T 0 = 3 and less than 0.52% at T/T 0 = 0.55. 35 Therefore, the latter is used in our simulations. The evolution of the particle distribution function depends only on the particle velocity in D-dimensional space. To reduce the computational cost, two reduced distribution functions residing in the lower dimensional phase-space [of dimension (2D)] are therefore introduced, 18 Thus, the evolution equation for g and h can be obtained from Eq. (1) to give where the reduced Shakhov equilibrium functions are g S = g eq + g Pr and h S = h eq + h Pr . g eq , h eq and g Pr , h Pr can be obtained by integrating f eq and f Pr over the (η, ζ) space, respectively. 18 Two important observations are in order. The first concerns the computational cost of the model. Through the Chapman-Enskog analysis, we can observe that when using the BGK-Shakhov model, a Gauss-Hermite quadrature of up to the twelth-order accuracy is needed to compute heat flux accurately as compared to the eigthorder accuracy without the Shakhov correction, which implies that, rigorously speaking, more discrete particle velocities should be used and that the computational cost will be increased. Therefore, an alternative model for a flexible Prandtl number will be considered instead. The second concerns the physical bulk viscosity in the model. It can be shown from the Chapman-Enskog expansion that the ratio between the bulk viscosity μV and the shear viscosity μ is χ ≡ μV /μ = 2/D − 2/(3 + K), which is always less than 2/D. Therefore, for a fixed specific heat ratio, the original DUGKS cannot be used to investigate the physical effect of the bulk-to-shear viscosity ratio. Moreover, in the previous simulations done by using traditional CFD methods, the Stokes assumption with zero bulk viscosity is applied, and the specific heat ratio γ is equal to 1.4. 35 However, in the original DUGKS, setting the specific heat ratio γ = 1.4 means that the internal degree of freedom is K = 2. This will result in a non-zero bulk viscosity μV = 4μ/15. This difference in bulk viscosities in the BGK-Shakhov model and in the simulations reported in the literature makes precise inter-comparison impossible. To resolve these two issues together, we shall propose next a new model, which can have a flexible Prandtl number and an arbitrary ratio between bulk viscosity and shear viscosity; then, only a lattice velocity model with at least eigth-order quadrature accuracy is needed to correctly compute the heat flux, the highest order moment encountered in the Chapman-Enskog analysis for the NSF system.

B. A new model with adjustable Prandtl number and the ratio of bulk to shear viscosity
We shall begin with the same two reduced Boltzmann equations but instead with to-be-designed source terms Sg and S h ,

ARTICLE
scitation.org/journal/phf where the equilibrium distribution functions g eq and h eq are simply the Maxwellian The two source terms Sg and S h are designed using an inverse design approach, as presented in Ref. 47. The design process is used to first find five integral constraints for the two source terms necessary to reproduce the NSF hydrodynamic system; then, a specific choice for Sg and S h is derived based on the properties of the Hermite polynomials. The design is not unique. The results of our specific design, as explained in detail in Ref. 47, are where the weighting function ω(ξ, , and T 0 is the reference temperature (see Appendix A for details). The velocity divergence is evaluated by either a secondorder or fourth-order finite difference scheme for the present simulations, which can be further optimized.
The expressions for evaluating conservative variables, heat flux, and the viscous stress tensor are given as By using the Chapman-Enskog expansion, one can recover the following constitutive relation for the viscous stress tensor σ: where S = (∇u + ∇u T )/2 is the strain rate tensor and I is the unit tensor. In addition, the closure model for heat flux vector q is consistent with Fourier's law, where the thermal conductivity coefficient is κ = μCp/Pr.
Applying the Chapman-Enskog analysis (see Appendix C), we can prove that the above-mentioned model can yield the following compressible Navier-Stokes-Fourier (NSF) system up to the order of O(τ): The following characteristic scales 35 are chosen to nondimensionalize the above-mentioned NSF system. We choose a reference length scale L 0 , reference velocity scale U 0 , reference time scale L 0 /U 0 , reference temperature T 0 , reference density ρ 0 , reference speed of the sound c 0 = √ γRT 0 , reference energy per unit volume ρ 0 U 2 , reference pressure p 0 = ρ 0 c 2 0 /γ = ρ 0 RT 0 , reference viscosity μ 0 , and reference thermal conductivity κ 0 . Then, the dimension analysis shows that the system is governed by five dimensionless parameters: the Reynolds number Re = ρ 0 UL/μ 0 , Prandtl number Pr = μ 0 cp/κ 0 , Mach number Ma = U 0 /c 0 , specific heat ratio γ (or internal degree of freedom K), and ratio between bulk viscosity and shear viscosity χ. If for a given fluid, γ, χ, and Pr are assumed to be known constants, the evolution of the system is then governed by two dimensionless parameters, i.e., the Reynolds number and the Mach number. For CDHIT, the initial flow Taylor microscale Reynolds number Re λ0 and the turbulent Mach number Mat 0 , along with the initial energy spectrum, are often used to describe the system.

III. THE DUGKS APPROACH
The implementation details about the DUGKS approach are described below. For convenience, we can rewrite the two Boltzmann equations in a unified form, where ϕ = g or h and the new collision operators are defined as As a finite-volume scheme, the computational domain is divided into many subcells Vj, with the cell centers denoted by xj. Integrating Eq. (15) over the control volume Vj from tn to tn +1 (=tn + Δt) and using the midpoint rule for the integration of linear convective term and the trapezoidal rule for the combined collision operators, we can obtaiñ where the flux across the cell interface at the half time step (tn + Δt/2) is Physics of Fluids ARTICLE scitation.org/journal/phf |Vj| and ∂Vj denote the cell volume and the surface area of Vj, and n is the outward normal vector of the cell interface ∂Vj. tn +1/2 = tn + s represents the half time step with s = Δt/2. Two new distribution functionsφ andφ are introduced in order to convert the implicit time integration scheme to an explicit scheme, namely, In the implementation, we track the distribution functionφ instead of the original one.
Once the distribution functionφ is obtained, the conservative variables at the cell center can be computed as In terms of the transformed distributions, the heat flux q and the viscous stress tensor σ at the cell center can be readily obtained as The surface flux integral at each interface is evaluated via the midpoint rule. In order to evaluate the flux term in Eq. (17) at the half time step, the original distribution function at the half time step is needed at the cell interfaces. This can be done by integrating Eq. (15) along the characteristic line for a half time step s = Δt/2 with the ending point located at the center of the cell interface x b . By using the trapezoidal rule for the collision term, we have Once again, two new transformed distribution functionsφ = ϕ −(s/2)Ω ϕ andφ + = ϕ+(s/2)Ω ϕ are introduced to remove the implicity. Then, we haveφ(x b , ξ, tn + s) =φ + (x b − ξs, ξ, tn). By applying the Taylor expansion, we can obtain where . For compressible flows, the proper treatments of the local valuē ϕ + (x b , ξ, tn) and local gradient σ b are essential. First, we implement the fifth-order WENO-JS scheme 40 to reconstruct the distribution functionφ + (x b , ξ, tn) at the cell interface by using the cell-averaged valuesφ n j (ξ) and enforcing the upwind rule according to the direction of a given discrete particle velocity. The WENO-Z 41,42 and WENO-M 43 schemes are also tested, and the results are compared with those from the WENO-JS scheme. See Appendix B for more details of these reconstruction schemes. Second, we use van Leer limiter 48 to obtain the slope σj = ∇φ + (xj, ξ, tn) at the cell center; then, it follows that the gradient at the cell interface can be approximated by σ b = (σj + σj +1 )/2. An alternative way to evaluate the gradient σ b is to use the finite difference schemes of different orders. To demonstrate the effect of different finite difference schemes used in the simulation, we will compare flow statistics from the second-order and the fourth-order central finite difference schemes in the second (higher Mach number) simulation case. If a shocklet is located near the cell interface, then σ b could have a jump across the shocklet surface. In order to deal with this kind of discontinuity, σ b is approximated by the gradient in the nearest upstream cell center or evaluated by the cell-averaged values in several upstream cells according to the direction of a given particle velocity. Finally,φ(x b , ξ, tn + s) can be obtained according to Eq. (21). In order to transform back to the original distributions, we first need to evaluate the conservative variables at the cell interface. They can be computed as The heat flux can be evaluated by Therefore, the equilibrium distribution function ϕ eq (x b , ξ, tn + s) can be computed based on the hydrodynamic variables obtained at the cell interface at the half time step. Then, the original distribution function can be obtained as Finally, the flux across the cell interface J n+1/2 ϕ (ξ) can be obtained by using Eq. (17).
Due to the nature of the linear transformation, the following two useful relations can be derived easily: At last, the distribution functionφ can be updated through Eq. (16). In the DUGKS, the time step is governed by the condition imposed on the Courant-Friedrichs-Lewy (CFL) number, Δt = CFL ⋅ Δxmin/(|u|max + |ξ|max), where Δxmin is the minimal grid spacing, |u|max is the local maximum hydrodynamic fluid velocity, and |ξ|max is the local maximum discrete velocity. A new discrete particle velocity model (E 9 3,77 ) with 77 discrete particle velocities and ninth-order Gauss-Hermite quadrature has been used in the implementation. See Appendix A for details.

IV. PROBLEM DESCRIPTION OF COMPRESSIBLE DECAYING HOMOGENEOUS ISOTROPIC TURBULENCE (CDHIT)
The dynamic process is more complicated in a compressible turbulence than that in an incompressible turbulence, due to stronger coupling effects between the hydrodynamic flow quantities and thermodynamic quantities. With the increase in the turbulent Mach number, such coupling process will become increasingly stronger and finally leads to emergence of the shocklets. 35 We consider CDHIT without an external driving force. The flow domain is a cubic box with the side length of L = 2π. The periodic boundary conditions and the uniform Cartesian grids are applied in all three directions. An initial divergence-free velocity field u(t = 0) is generated by a Gaussian random field with a prescribed kinetic energy spectrum E(k, t = 0) in the spectral space. Then, the initial velocity field in the physical space is obtained by carrying out the inverse Fourier transformation. The initial kinetic energy spectrum is specified as 29 where A 0 is the constant, k is the wavenumber, and k 0 is the peak wavenumber.
The relevant flow statistical quantities are summarized below. The root-mean-square velocity is defined as urms = ⟨u ⋅ u⟩ 1/2 / √ 3, where the symbol ⟨⋯⟩ denotes the spatial average over the whole domain at a given time. The Taylor microscale and Taylor microscale Reynolds numbers are defined as The turbulent kinetic energy per unit mass and per unit volume is defined as where E(k, t) and Eρ(k, t) are the three-dimensional kinetic energy spectra corresponding to K(t) and Kρ(t). The longitudinal integral length scale and the large eddy turnover time are given by The average viscous dissipation rate and the Kolmogorov length scale are defined as The turbulent Mach number Mat and the local Mach number Ma loc are, respectively, defined as The root mean square dilatation, an important property for compressible flow, is defined as In addition, the skewness and flatness of the longitudinal velocity derivatives are where the velocity derivatives are evaluated by first transforming the velocity field into the spectral space, performing the gradient operations, and then transforming back to the physical space using the inverse Fourier transform. In our simulation, the initial density ρ 0 , pressure p 0 , and temperature T 0 are set as constants. The initial turbulent kinetic energy K 0 , the root mean square velocity urms 0 , dissipation rate ε 0 , Kolmogorov length scale η 0 , integral length scale L f 0 , and large eddy turnover time Te 0 can be derived from the initial energy spectrum to give Therefore, the initial Taylor microscale Reynolds number Re λ0 and the initial turbulent Mach number Mat 0 can be computed as The time will be normalized by the initial large eddy turnover time Te 0 , and the subscript "0" of Te 0 will be omitted in all the subsequent figures for simplicity.

V. SIMULATION RESULTS
A three-dimensional DUGKS code is developed to simulate CDHIT by using the newly developed mesoscopic model. To

A. The low turbulent Mach number case
A case with a relatively low initial turbulent Mach number and a low initial Reynolds number (Mat 0 = 0.3, Re λ0 = 30) is first simulated to validate our three-dimensional DUGKS code. The bulk-to-shear viscosity ratio is set as 0 or 4/15, and the specific heat ratio γ is 1.4. The constant A 0 of the prescribed spectrum is 3.74 × 10 −4 , and the corresponding peak wave number is k 0 = 4. The turbulent kinetic energy K(t) is normalized by its initial value K 0 , and three quantities representing normalized root mean square thermodynamic fluctuations are defined as follows:

ARTICLE scitation.org/journal/phf
where the specific volume of the fluid is defined as V = 1/ρ and V 0 = 1/ρ 0 . The physical quantity with a subscript rms denotes the root mean square value of the physical quantity. The CFL number is set to 0.5. In this case, the large eddy turnover time Te = 2 √ 3/(k 0 M t0 ) in order to be consistent with the definition in the literature for the postprocessing. 30 The fifth-order WENO schemes are used to reconstruct the particle distribution functions at the cell interface.
First, several low resolution simulations are conducted. With 64 3 uniform grids as that used in the literature, Fig. 1 shows the results obtained by the DUGKS combined with different WENO schemes. The results of WENO-Z (p = 1) and WENO-Z (p = 2) almost overlap with those of WENO-M. The simulation results are compared to the data given by Honein and Moin, 30 Liao et al., 32 and Subbareddy and Candler. 33 Overall, all the DUGKS results match well with the literature data even at such a low resolution. The results obtained by the WENO-JS scheme are a little bit lower than other DUGKS results because WENO-JS is more dissipative than other WENO schemes. In contrast, the results obtained by Subbareddy and Candler. 33 are lower than other results for normalized thermodynamic fluctuation quantities after two large eddy turnover times.
Next, Fig. 2 shows the evolution of the velocity derivative skewness and flatness with 512 3 uniform grids. Their values approach almost a constant after one large eddy turnover time, indicating that the flow becomes physically developed. Then, the time evolution of normalized turbulent kinetic energy K(t)/K 0 and those normalized thermodynamic fluctuation quantities are carefully compared with the existing data from the literature in Fig. 3. Overall, the DUGKS results agree well with the reference data from the literature. We note that the reference results in the literature are obtained with only 64 3 uniform grid points, while the DUGKS results are obtained with higher resolution. With the grid resolution 512 3 , the results for thermodynamic fluctuation quantities almost overlap with those obtained by Honein and Moin 30 who used an entropy conservation sixth-order compact finite difference scheme. However, they do not provide enough data points in the peak region for comparison. Therefore, we do not know the exact position of the peak. It is observed that the peak position could be different when using different numerical treatments. For example, near the peak value of the curve, our results are almost identical with those obtained by Subbareddy and Candler 33 who used a low dissipation, finite-volume scheme. In contrast, the GKS results obtained by Liao et al. 32 are lower in the peak region. Moreover, the results obtained by different bulk-to-shear viscosity ratios only exhibit minor difference because the compressible effect is not so strong, and the difference in the bulk-to-shear viscosity ratio is not large enough in the present case. Still, we can see that the curve with χ = 4/15 is a little bit lower than that with χ = 0 for the thermodynamic flow quantity because a small increase in the bulk-to-shear viscosity ratio can decrease the compressibility effect slightly.
In the simulation of incompressible isotropic turbulence, Wang et al. 19 have shown that the DUGKS can resolve the flow when kmaxη > 3 in their simulated cases, where kmax is the maximum wave number and η is the Kolmogorov length scale. However, for compressible isotropic turbulence, the approximate value for kmaxη is still unknown. The grid convergence study is done with the different uniform grids N 3 = 64 3 , 128 3 , 256 3 , 512 3 , and the simulation results are shown in Fig. 4. The corresponding initial values of kmaxη are 1.40, 2.88, 5.84, and 11.69, respectively. The DUGKS result with the lowest resolution 64 3 is reasonable but does not coincide very well with the reference data, especially for the kinetic energy. With the increase in the resolution, the results gradually converge to the data points given by Honein and Moin. 30 It is clear that the result with  18 Therefore, in all the subsequent sections, the CFL number is fixed at 0.5.

B. The higher turbulent Mach number case
To further validate the capability of the present DUGKS approach to simulate locally subsonic and supersonic turbulent flows, we consider a case with a higher initial turbulent Mach number of Mat 0 = 0.488 and a Taylor microscale Reynolds number of Re λ0 = 175, previously studied by Samtaney et al. 29 and Frapolli et al. 31 The initial divergence-free velocity field is generated by an energy spectrum given in Eq. (26) with the constant of A 0 = 0.011 and the peak wavenumber of k 0 = 4. Two different grid resolutions, 256 3 and 512 3 , are adopted for the newly developed mesoscopic model. The corresponding initial kmaxη 0 are 2.43 and 4.89, respectively. The CFL number is set to 0.5. The ratio of bulk to shear viscosity is 0 as in the simulations from the literature. First, in Fig. 6, we display the time evolution of velocity derivative skewness and the root mean square value of velocity divergence. Results are compared with those from a NSF solver using a hybrid scheme 35 and the results from a tenth-order compact finite difference scheme. 29 Overall, our results are in good agreement with the literature data. Based on the time evolution of the velocity derivative skewness, we conclude that the flow has not evolved into a physical flow before t/Te = 1.5, during which the specific random flow initialization may affect the flow statistics. Different numerical treatments and grid resolutions will also introduce numerical dissipations at different levels, which can lead to the slight deviations near the peak and valley regions, but these regions occur during the transition from the initial random flow to the physical CDHIT flow. For these reasons, slight deviations at the regions near the peak and valley points for the skewness and the peak regions for the root mean square velocity divergence should not be taken too seriously. The evolution of the rms velocity divergence for t/Te > 1.5 is in excellent agreement with the literature results. As illustrated in Fig. 6(b), even for the results between the work of Samtaney et al. 29 and Wang et al., 35 some difference can be observed in the neighborhood of the peak value of the curves. In the following discussions, we will further confirm our statement here by simulating the same case with different WENO schemes and different treatments for interfacial derivatives of distribution functions at the cell interface. Next, we compare the probability density function (PDF) of the local Mach number obtained by using the DUGKS approach with those obtained by Samtaney et al. 29 and Frapolli et al. 31 In Fig. 7, the DUGKS results and the benchmark data are compared with each other at the dimensionless time t/Te = 1.56. We provide both the PDFs with and without the logarithmic y coordinate for comparison. In Fig. 7(a), we found that the DUGKS results agree very well with the benchmark data except for a minor difference in the region at the very end of the tail, perhaps due to the insufficient data samples. In Fig. 7(b), it is also observed that the DUGKS results overlap with the benchmark data for almost all Mach number regions.

Physics of Fluids
A comparison of the PDF of normalized root mean square velocity divergence is shown in Fig. 8. We observe that the DUGKS results with 256 3 and 512 3 uniform girds are in excellent agreement with the data points taken from the work of Samtaney et al. 29 In contrast, the DUGKS result with 512 3 uniform grids exhibits a longer tail with large negative velocity divergence than that with 256 3 grids, although the probability of the large negative value is relatively small. The longer tail indicates that a high resolution simulation can better capture extreme turbulent fluctuations of different scales in local compression regions. Figure 9 displays the PDF of the local Mach number and the PDF of normalized root mean square velocity divergence obtained using the third-order WENO-JS scheme and the fifth-order WENO-JS scheme to reconstruct the particle distribution functions at the cell interface. Again, the DUGKS results are in good agreement with the benchmark data. In Fig. 9(a), only a slight deviation can be observed near the maximum Mach number. The results with the fifth-order WENO-JS scheme appear much closer to the benchmark data near the maximum Mach number when comparing with the third-order WENO-JS scheme. In Fig. 9(b), different treatments only show slight difference in the tails of the PDF. We choose to use the fifth-order WENO-JS scheme in the following simulations. To further validate our code, Fig. 10 shows the time history of the normalized turbulent Mach number and the normalized turbulent kinetic energy. The DUGKS results are in excellent agreement with the benchmark data. We also observe that the local maximum Mach number can reach values as high as 1.5, which is consistent with the observation reported by Frapolli et al. 31 It would be interesting to examine the distribution of the local Knudsen number, Kn ≡ λ/l, where λ is the molecular mean free path and l is the local characteristic length scale. The value of l could be defined based on a hydrodynamic field variable, so the choice is not unique. Although the local Knudsen number can have different values due to different definitions of l, the following discussion can provide us a better understanding of the local flow state. Based on the density field and the temperature field, we can define three local Knudsen numbers, namely, Knρ ≡ λ|∇ρ|/ρ, Knpress ≡ λ|∇p|/p, and Kn T ≡ λ|∇T|/T, where the mean free path is given by λ = (μ/p) √ πRT/2 for air. 17,49 The evolution of PDFs of Knρ and Knp is shown in Fig. 11. The peak value of the PDF increases, and the corresponding Knudsen number decreases over time. As a result, the flow state appears to be within the continuum regime completely as the flow evolves with time. We observe that the local Knudsen number can be larger than 0.001 at earlier times. In the region with a high local Knudsen number, the flow may not belong to the continuum flow regime so that the NSF system cannot provide a satisfying description to the local flow state, especially near or inside the shocklets. However, for most flow regions, the NSF system can be accurate such that the NSF-based solver and the Boltzmann-based solver can yield identical results. Figure 12 shows two snapshots of contours of normalized velocity divergence ϑ/ϑrms and local Mach number Ma loc at z/L = 0.502 and t/Te = 1.56. Figure 13 displays slices of contours of

and
Ma loc are all positively correlated with the local Knudsen numbers although the correlation coefficients are relatively small. These can be explained from two aspects. First, the negatively skewed PDF of the velocity divergence demonstrates that the very strong compression motions associated with the shocklets exist more frequently than the expansion motions in the whole flow field. Generally speaking, these strong compression motions could more effectively decrease the local length scale and thus increase the local Knudsen number. However, local expansion motions can also decrease the local density, increase the local density gradient, and thus increase the local Knudsen number, which is more pronounced when the expansion region is surrounded by some local compression regions (high density regions with negative velocity divergence

C. Comparison of different treatments of interfacial derivatives
In the simulations mentioned above, the interfacial distribution functionφ + (x b , ξ, tn) is reconstructed from the cell-averaged valueφ n j (ξ) by applying the fifth-order WENO-JS scheme. Based on the known distribution functionsφ + (x b , ξ, tn) andφ n j (ξ), different numerical methods can be applied to evaluate the interfacial gradient σ b . In addition to the van leer limiter (method A) mentioned above, we also test the following the second-order (method B) and the fourth-order (method C) central finite difference schemes in our simulation. Note that for smooth flows without shocklets, σ b should be continuous at the cell interface, while a jump of σ b may occur across a strong shocklet surface in high turbulent Mach number flows. For the latter cases, the interfacial gradient σ b can be approximated by the velocity gradient in the nearest upstream cell center obtained by using the van Leer limiter according to the direction of a given particle velocity (method D). The WENO-JS scheme and WENO-Z scheme are combined with methods A-D, respectively. The parameter p in the nonlinear weights of the WENO-Z scheme is equal to one. The grid resolution is 256 3 .
For instance, we consider the implementation for the cell interfaces facing the x direction. In the normal direction of the cell interface, the second-order and fourth-order central FD schemes can be written as In the tangential direction of the cell interface, the second-order and fourth-order central FD schemes can be expressed as where the subscripts i, j, k denote the position of the function f and Δx, Δy, and Δz represent the grid spacings in x, y, z, respectively. All comparisons of methods A-D are performed at t/Te = 1.56. Figures 15 and 16 show the PDFs of the local Mach number and normalized velocity divergence. The results obtained from different methods are in good agreement with each other, and only small difference can be seen at the tails, which represent local highly compressed or expanded regions. Therefore, the PDFs of the local Mach number and normalized velocity divergence are not sensitive to the methods used to evaluate the interfacial gradient σ b . Figures 17 and 18 display the time history of velocity derivative skewness and root mean square velocity divergence, respectively. As demonstrated in Fig. 17, the results obtained from methods A-D agree well with each other for both the WENO-JS and WENO-Z schemes, which indicates that the skewness is also not sensitive to the methods used. In Fig. 18(a), with the WENO-JS scheme, the DUGKS result from method A coincides well with the result from the hybrid scheme for most regions but has a lower peak value. In addition, the DUGKS results obtained from methods B-D are slightly lower than that obtained from method A. However, in Fig. 18(b), with the WENO-Z scheme, the DUGKS results obtained from methods B and C agree well with the hybrid scheme except for the peak region of the curve. In contrast, the DUGKS results obtained from methods A and D have a slightly higher value than that from the hybrid scheme for most regions but can reach the same level in the peak region. For this 3D compressible case, no analytical solution can be derived, and the truth benchmark can never be known thoroughly. If the traditional and mesoscopic CFD methods yield comparable results, then we can have more confidence in our simulation results. Moreover, we compare the instantaneous density and temperature profiles (not shown here). The results obtained from different methods also match very well with each other. From the above-mentioned results, we can conclude that different treatments at the cell interface can yield very similar results.

D. Comparison of different WENO schemes
In this subsection, we reconstruct the particle distribution functions at the cell interface by using WENO-JS, WENO-Z, and WENO-M schemes, respectively. The parameter p in the WENO-Z scheme is taken as 1 and 2 in the simulation. The details and some discussion on these three WENO schemes are included in Appendix C. As discussed above, method A is still used to compute σ b at the cell interface. A 256 3 uniform grid is applied.
The PDFs of the local Mach number and normalized velocity divergence are plotted in Figs. 19(a) and 19(b). The DUGKS results obtained from different WENO schemes agree well with each other and the reference data. Figure 20(a) shows the evolution of velocity derivative skewness. Only small differences can be observed between these four curves. The evolution of root mean square velocity divergence is also shown in Fig. 20(b). The results obtained from two WENO-Z schemes and WENO-M scheme are almost identical, but their magnitudes are slightly higher than those from the hybrid scheme and WENO-JS scheme for most regions although they share almost the same peak value. This is because the WENO-Z scheme and WENO-M scheme are more accurate and less dissipative than the WENO-JS scheme, particularly near the critical points. We also confirm that the results do not show visible difference with the small parameter of ε 0 = 10 −5 and 10 −6 (not shown here). The comparison leads to the conclusion that the PDFs of the local Mach number, normalized velocity divergence, and velocity derivative skewness are not very sensitive to the details of the WENO schemes applied. As to the root mean square velocity divergence, the results are still reasonable. Different dissipation levels of different WENO schemes contribute to the small difference. We also compare the instantaneous density and temperature profiles (not shown here). No obvious difference could be observed, which implies that the instantaneous flow fields are also not sensitive to the WENO schemes applied. The Stokes hypothesis (i.e., zero bulk viscosity) is extensively applied to the NSF-based simulation of viscous compressible flows, in particular to gas dynamics problems. 29,35 The compressible turbulence is not expected to be affected by the bulk viscosity significantly when the bulk viscosity is of the order of the shear viscosity. 50 Therefore, the simulations with and without the Stokes assumption will yield almost the same results. As a direct application of our newly developed model, the effect of the bulk to shear viscosity ratio on PDFs of some physical quantities is explored by using our DNS data. Three viscosity ratios χ = 0, 30, and 100 are considered, and some preliminary results are presented next.
First, the PDFs of normalized velocity divergence and local Mach number are shown in Figs. 21(a) and 21(b) at t/Te = 1.56. When computing the PDFs, we first find the maximum (Amax) and minimum (A min ) values of the field; then, the bin width is defined as (Amax − A min )/N bin . The number of bins, N bin , is chosen to be 400. Oscillations of the PDF tails reflect the inadequate number of samples used there. As illustrated in Fig. 21(a), we observe that the left tail with large negative values becomes shorter with the increase in the bulk to shear viscosity ratio, which is qualitatively consistent with those reported by Chen et al. 37 and Pan and Johnsen. 50 This is because the increase in the bulk to shear viscosity ratio can damp the compressibility. In Fig. 21(b), the maximum Mach number also decreases with the increase in the bulk to shear viscosity ratio, which further confirms the previous observation. An explanation for the bulk viscosity effect is given as follows. For simplicity, neglecting the spatial variation of the shear and bulk viscosities in the viscous term, the spatially averaged evolution equation for ⟨ϑ 2 ⟩/2 in a periodic domain can be derived from the momentum equation, yielding where the minus sign in the last term implies that the bulk viscosity always reduces the velocity divergence magnitude in both the expansion and compression regions. This equation can be used to further analyze different contributions to the change in the mean square dilatation, which is beyond the scope of this paper. The derivation details are included in Appendix E. Next, Figs. 22(a) and 22(b) present the PDFs of normalized vorticity magnitude ω/ωrms and the normalized vortex stretching term Φ = (ω ⋅ S ⋅ ω)/(⟨ω 2 ⟩(S : S) 1/2 ) at t/Te = 1.56. Different from the previous ones, it is found that the PDFs of normalized vorticity magnitude and the normalized vortex stretching term exhibit little sensitivity to the ratio of bulk to shear viscosity. The PDFs are almost identical for different bulk-to-shear viscosity. This is also consistent with the findings by Pan and Johnsen. 50 Another evidence can be found in Fig. 23 that displays the vortices extracted by Q-criterion 51 at t/Te = 1.56 for χ = 0, 30, respectively. We observe that the vortex structures are almost unaffected by the bulk-to-shear viscosity ratio. This can be understood from two aspects. On the one hand, the vorticity dynamics is dominated by the vortex stretching and sheardensity gradient effects and could be partly influenced by the bulk viscosity-density gradient effect. 50 On the other hand, increasing the bulk-to-shear viscosity ratio can attenuate the coupling between the solenoidal mode and the compressive mode such that the solenoidal-velocity-dominated vortical motion is not evidently altered by the bulk viscosity.
To have an intuitive understanding of effect of the bulk to shear viscosity ratio, the isosurfaces of the normalized velocity divergence ϑ/ϑrms with χ = 0, 30 are also visualized in Figs. 24(a) and 24(b) at t/Te = 1.56, respectively. The red isosurfaces (ϑ/ϑrms = −2) denote the compression regions, while the blue ones (ϑ/ϑrms = 2) represent the expansion regions. We observe that more blob-like structures emerge, and the sheet-like structures is less visible with the increase in the bulk to shear viscosity ratio. The increase in the bulk to shear viscosity ratio damps the compressible effect in the local strong compression regions and reduces the probability of the existence of the shocklets, which finally drives the flow toward the incompressible state.

VI. COMPUTATIONAL COST
We ran both the DUGKS code and the hybrid compact FD-WENO code for 10 000 time steps at the 256 3 resolution on the Taiyi cluster provided by the Center for Computational Science and Engineering at the Southern University of Science and Technology. Taiyi is equipped with 815 nodes (2.4 GHz Intel Xeon Gold 614820c), each containing 40 cores. The hybrid compact FD-WENO code was developed by Wang et al. 35 and was parallelized using one-dimensional domain decomposition. Due to the requirement to implement the finite difference schemes in the code, the maximum number of the cores that can be used in this case is 64. The total central processing unit (CPU) time is 506 469 s, and the wall clock time is 7 888 s. Our DUGKS code, which was parallelized using the two-dimensional domain decomposition strategy, 52 can be run easily with 256 cores with 16 divisions in each of the two directions (more processors can be used if needed). The total CPU time is 16 658 809 s, and the wall clock time is 65 890 s. This comparison indicates that currently the DUGKS code is inefficient when compared to the hybrid compact FD-WENO code. Two main reasons are given here. First, the compressible DUGKS code has not been optimized, and many implementation details could be improved. Second, the hybrid compact FD-WENO code only solves five hydrodynamic variables including density (ρ), temperature (T), and three velocity components (ux, uy, uz). While in the DUGKS code, we track the particle distribution functionsg and h. The hydrodynamic variables are obtained by integratingg andh over the particle velocity space. Using the E3Q77A9 discrete particle velocity model, 154 discrete distribution functions were solved, namely,gα andhα, α = 1, 2, . . ., 77. The relative ratio of the number of variables between the two methods is about 1:30, which reflects roughly the CPU time ratio. We are in the process of optimizing the DUGKS code including decreasing the number of the discrete particle velocities, but such efforts are beyond the scope of the current paper, and they will be reported in the future. The results from the current paper will serve as a benchmark for such future optimization efforts.

VII. CONCLUSIONS AND DISCUSSIONS
We have proposed a new kinetic model with both an adjustable Prandtl number and a tunable ratio of bulk to shear viscosity and applied the model to simulate CDHIT under the DUGKS approach. Different fifth-order WENO schemes are used to reconstruct the distribution function at the interface, and different treatments of the interfacial derivatives are implemented. Simulation results from various implementations are compared with each other and with the literature data. A newly developed 3D lattice velocity model with 77 discrete particle velocities and a Gauss-Hermite quadrature order at 9 has been used to evaluate the velocity moments of the particle distribution functions. This new lattice velocity model can ensure an accurate evaluation of the heat flux at the mesoscopic level.
We apply the re-designed DUGKS approach to simulate two CDHIT cases with different initial turbulent Mach numbers and Taylor microscale Reynolds numbers. The simulation results agree well with the results from the literature. For the local flow region with high compressible effects, the NSF system may no longer be accurate in describing the local flow state, while the Boltzmann equation can handle both the continuum flow regime and the flows in the non-continuum regimes. Due to the coupled treatment of the particle transport and collision processes, the DUGKS can successfully simulate the flows with multiscale physics, with a second-order accuracy in both space and time. The hydrodynamic variables represent only the low-order moments of the distribution functions. Although the use of more discrete velocities increases the computational cost, more information can be obtained from the particle distribution functions at a given point in the DUGKS.
As a direct application of our DNS data, we study the effect of bulk viscosity on turbulence statistics and structures. The DUGKS

ARTICLE
scitation.org/journal/phf results show that the compressibility effect in CDHIT could be damped with the increase in the bulk to shear viscosity ratio, similar to what has been found in compressible stationary HIT. 37,50 Our results show that the DUGKS approach can be a reliable tool for simulating CDHIT at low and moderate turbulent Mach numbers. In a companion paper, we first report results in applying the DUGKS approach to simulate compressible stationary HIT. 53 More parametric studies and the code optimization are needed in the future. The nth order Hermite polynomial 54 is defined by where the weighting function ω(ξ, T 0 ) is ω(ξ, T 0 ) = exp(−ξ 2 /(2RT 0 ))/(2πRT 0 ) D/2 and T 0 is the reference temperature. The Hermite polynomials of different orders are orthogonal to each other in the following sense: where i represents the abbreviation of i 1 i 2 ⋯in and δ n ij is equal to 1 if and only if i is the permutation of j.
From the Chapman-Enskog expansion, we observe that the fourth-order Hermite expansion of the equilibrium g eq and the second-order Hermite expansion of the equilibrium h eq are sufficient for accurate calculation of the heat flux. They are given by where A 0 -A 4 are given by Based on the orthogonal property of Hermite polynomials, Shan et al. 55,56 developed a general solution of the discrete velocity model by solving linear equations for weights of velocities with predefined abscissas on the lattice points. The three-dimensional discrete velocity model E 9 3,103 with 103 dicrete particle velocities and a Gauss-Hermite quadrature order at 9 has been proposed. However, the DUGKS does not require the abscissas of the discrete velocities to be coincided with the Cartesian grids. This provides us with an opportunity to design a discrete velocity model with a less number of discrete velocities and the same Gauss-Hermite quadrature order of 9. Recently, a new model E 9 3,77 with a less number of discrete velocities is proposed by Wen et al. 57 They began with a construction based on E 9 3,125 = E 9 1,5 × E 9 1,5 × E 9 1,5 and then examined the true number of constraints in the model for achieving a Gauss-Hermite quadrature order of 9 in 3D. They found that the one lattice velocity group in E 9 3,125 can be removed, leading to E 9 3,77 . The relevant symmetry groups, discrete velocity vectors, number of velocities, and the corresponding weights of the resulting E 9 3,77 model are summarized in Table I,   To improve the spatial accuracy, we employ the fifth-order WENO-JS 40 scheme to reconstruct the particle distribution functions at the cell interface. Due to the nature of linear convection of the Boltzmann equation, the upwind rule should be satisfied by considering the directions of particle velocities in order to enhance the numerical stability.
(B2c) Therefore, the nonlinear weights can be expressed as ,ω JS q = γq (ϵ 0 + βq) 2 , j, q = 0, 1, 2, where γ 0 = 1/10, γ 1 = 3/5, and γ 2 = 3/10 are the linear weights. The small parameter ε 0 is originally used to avoid divisions by zero in the nonlinear weights. As suggested, ε 0 = 10 −6 is often used in the simulation. 38,40 However, it also has a collateral effect on the order of accuracy of the scheme. [58][59][60] Therefore, according to the above discussion, ε 0 should be the order of O((δx/L) 2 ) (with the same dimension as the smoothness factors), where δx is the grid spacing and L is the side length of the cubic box.
Finally, the distribution functionsφ + (x j+1/2 , ξ α , tn) can be a weighted sum, given as ϕ + (x j+1/2 , ξ α , tn) = ω 0φ +(0) (x j+1/2 , ξ α , tn) + ω 1φ +(1) (x j+1/2 , ξ α , tn) If the x component of the particle velocity ξαx < 0 at the cell interface xj +1/2 , the similar procedures can be repeated with three cells in the upstream direction and two cells in the downstream direction. If ξαx = 0, we take the average of the results obtained by two cases mentioned above. For the cell interfaces in the y direction and z direction, the logic is similar to what we have done in the x direction. The reconstruction process is displayed in Fig. 25.
The WENO-JS scheme can lose its optimal order of accuracy at the critical points (where df /dx = 0, but d 2 f /dx 2 ≠ 0) of a smooth function. 41,59 Therefore, Borges et al. 41 and Castro et al. 42 proposed the WENO-Z scheme by redesigning the nonlinear weights that satisfy the necessary and sufficient conditions for an optimal order at the critical points. The WENO-Z scheme is less dissipative than WENO-JS. 60 In addition to WENO-JS and WENO-Z schemes, Henrick et al. 43 proposed the WENO-M scheme by introducing a mapping function to the weights. The above conclusions are drawn from 1D or 2D simulations in the existing literature. For the present 3D implementation using the DUGKS approach, the effect of different WENO schemes on turbulence statistics should also be further studied. In this appendix, we provide the details of the Chapman-Enskog analysis for the redesigned mesoscopic model. From Eqs. (8a) and (8b), the Chapman-Enskog expansion for two reduced distribution functions g and h can be written as From the definitions of the equilibriums and source terms in Eqs. (9) and (10), the following moments can be evaluated directly: ∫ g eq dξ = ρ, ∫ ξg eq dξ = ρu, ∫ (ξ 2 g eq + h eq )dξ = 2ρE, (C2) The remaining task is to find the explicit expression for the heat flux q by using the Chapman-Enskog expansion of two reduced distribution functions given in Eqs. (C1a) and (C1b). First, we note that by using Eqs. (C3) and (C4), we obtain ∫ cc 2 Sgdξ = (D + 2)(χ − ∂g eq ∂t + ξ ⋅ ∇g eq )dξ + ∫ c( ∂h eq ∂t + ξ ⋅ ∇h eq )dξ] Therefore, the Fourier's law can be recovered, as shown in Eq. (13). Finally, it is worth pointing out that the source terms can be designed by applying the Hermite expansion with five constraints in order to recover the NSF system. 47 Obviously, the simplest choice is to take a truncated Hermite polynomial without increasing the Gauss-Hermite quadrature order.