Entropy Conditions Involved in the Nonlinear Coupled Constitutive Method for Solving Continuum and Rarefied Gas Flows

The numerical study of continuum-rarefied gas flows is of considerable interest because it can provide fundamental knowledge regarding flow physics. Recently, the nonlinear coupled constitutive method (NCCM) has been derived from the Boltzmann equation and implemented to investigate continuum-rarefied gas flows. In this study, we first report the important and detailed issues in the use of the H theorem and positive entropy generation in the NCCM. Importantly, the unified nonlinear dissipation model and its relationships to the Rayleigh–Onsager function were demonstrated in the treatment of the collision term of the Boltzmann equation. In addition, we compare the Grad moment method, the Burnett equation, and the NCCM. Next, differences between the NCCM equations and the Navier–Stokes equations are explained in detail. For validation, numerical studies of rarefied and continuum gas flows were conducted. These studies include rarefied and/or continuum gas flows around a two-dimensional (2D) cavity, a 2D airfoil, a 2D cylinder, and a three-dimensional space shuttle. It was observed that the present results of the NCCM are in good agreement with those of the Direct Simulation Monte Carlo (DSMC) method in rarefied cases and are in good agreement with those of the Navier–Stokes equations in continuum cases. Finally, this study can be regarded as a theoretical basis of the NCCM for the development of a unified framework for solving continuum-rarefied gas flows.


Introduction
The numerical study of continuum-rarefied gas flows is of great interest because it can provide fundamental knowledge regarding flow physics and provide a theoretical tool to precisely predict the aerodynamic or aerothermodynamic performance of hypersonic vehicles and/or Micro-Electro-Mechanical Systems [1][2][3][4].For more than 150 years, the Navier-Stokes-Fourier (NSF) equations have been accepted as the model that describes flow dynamics.Unfortunately, the NSF equations have serious limitations in capturing the correct flow physics under highly nonequilibrium conditions, such as for rarefied gas [5].For the investigation of rarefied and/or micro-gas flows, much effort has been put into the development of computational models beyond the NSF framework.The models can be classified into two categories: the full kinetic model and the fluid dynamics model.The most successful methods in the former category are the direct simulation Monte Carlo (DSMC) [6][7][8][9] and Lattice Boltzmann Methods (LBM) [10][11][12][13][14], and they have been extensively used for the computation of rarefied and/or micro-scale gas flows.Although the Direct Simulation Monte Carlo (DSMC) method is commonly used in the investigation of rarefied gas flows, it has unacceptable memory demands in near continuum states because it uses simulated molecules to model the movement of molecules and the collision of gas flows.Great progress has been made in the development of LBM to remove its capability limit in complex far-from-equilibrium flows, of which the regularized LBM work done by Montessori and Succi is typically representative [12].Importantly, LBM has been successfully extended to investigate turbulent gas flows by Chen and Succi [15].On the other hand, the fluid dynamics models are based on the hyperbolic conservation laws and nonconserved variables appearing in the conservation laws, the latter being determined by their evolution equations which can be derived with the help of the Boltzmann equation.Up to present, several models have been developed for the latter category.This category may be subdivided into a few classes: in one of them, the Burnett-type equations [16,17] are used for devising computational models; in another, Grad moment equations [18] and 13 moment equations [19,20] are employed in conjunction with extended thermodynamics; and in still another, the Eu moment equations that serve as the basis of the generalization of thermodynamics [21,22] have been used in a manner consistent with the laws of thermodynamics at every order of approximation employed.Recently, Eu moment equations have been well-developed by Myong to investigate stiff hypersonic nonequilibrium gas flows [22][23][24][25][26].
Recently, we reported and validated the nonlinear coupled constitutive method (NCCM) for the development of a unified scheme for modelling continuum-rarefied gas flows [4].It can be regarded as another option to treat EU equations.In addition, we tried to extend the NCCM to investigate turbulence and have made great progress recently.In the NCCM, the conservation laws of mass, momentum, and energy are derived from the Boltzmann equation via the EU method [21,22] and are similar to those of the NSF equations.However, nonconserved variables associated with thermal nonequilibrium, such as the shear stress tensor and the heat flux vector, are described by the evolution equations.These variables are quite different from those in the NSF equations in conjunction only with the linear constitutive relations of the gradients of velocity and temperature.In the present study, we first provide the theoretical basis of the NCCM, particularly for positive entropy generation and the unified nonlinear dissipation model according to Eu method.This study can be regarded as an extension of the previous study of the NCCM and Eu method, in which we focus on the numerical scheme for solving NCCM equations and its validation in typical continuum-rarefied gas flows.

Boltzmann Equations and the Feature of Irreversibility
The distribution function f i denotes the probability of finding a particle of species i in the range of v i + dv i and r i + dr i at time t.The Boltzmann equation yields the evolution equation for f i , In this equation, R( f i ) is the collision term, F i is the external force, and v i is the velocity of a particle of species i.We can simply rewrite R( f i ) in the following form: where b is the impact parameter and it is the distance between the two parallel lines passing through the center of the collision molecule, r is the number of species, and g ij = v i − v j is the relative speed.C f i , f j is the Boltzmann collision integral, and it accounts for the decrease and increase in probability caused by collisions and reverse collisions, respectively [21].
If we reverse time t to −t, then the velocities must also be reversed, and the Boltzmann equation takes the following form: In this equation, we have Because the time is reversed, the positions of the post-collision and pre-collision terms are exchanged in Equation ( 4).This reversion transforms the Boltzmann collision integral into the negative of the original one before the time reversal.Thus, the time reversal invariance is broken because of the presence of the Boltzmann collision integral.As a result, in the sense that the time reversal invariance is broken, the Boltzmann equation is irreversible and is consequently expected to describe the irreversible macroscopic process towards equilibrium.This feature expression is very important and is represented in a more useful form via the H theorem and the Entropy balance equation [21].

H Theorem and Entropy Balance Equation
The H function was introduced by Boltzmann to show that the solution of the Boltzmann equation generally approaches a unique equilibrium solution.In addition, this function's time derivative can be proved to decrease to zero as the system approaches equilibrium.Moreover, this H function was found to be closely related to the entropy introduced by Clausius.Thus, the nonequilibrium entropy of the system is defined by [21]: The H theorem states that dS dt ≥ 0.
If we use the local equilibrium temperature T e to replace the temperature T in the Maxwell-Boltzmann distribution function, then the local equilibrium distribution function can be obtained: Thus, the local equilibrium entropy can also be expressed by Furthermore, the nonequilibrium entropy can be defined by By the well-known inequality, y ln(y/x) − y + x ≥ 0 (10) We conclude that In other words, the entropy increases towards the equilibrium entropy, as dictated by the second law of thermodynamics.By differentiating S n (t) as a function of time, we obtain Substituting the Boltzmann equation and proceeding in the same manner as that for the proof of the H theorem, we find This finding indicates that the nonequilibrium entropy increases as the system approaches equilibrium.
In addition, this feature can be expressed in the latter entropy balance equation.We define the entropy density ζ by the formula where ζ is the entropy per unit mass of the gas, and the total entropy can be obtained by integrating ρζ over the volume.By differentiating ρζ over time and using the Boltzmann equation on the right hand, the entropy balance equation can be obtained: where J s (r, t) is the entropy flux, and σ ent (r, t) is the entropy production.Next, the same procedure as that used to prove the H theorem can be applied to show that the entropy production is positive: Therefore, the entropy balance equation has a positive source for entropy, i.e., the source creates entropy.The statistical expression for σ ent (r, t) suggests that entropy is created through numerous, random collisions of the particles in the gas because entropy production is intimately associated with energy dissipation.In addition, the energy dissipation may be regarded as a result of particle collisions that occur incessantly and randomly in the gas flow.From the discussion above, one can conclude that the Boltzmann equation is irreversible in the nonequilibrium state.In addition, this feature can be expressed by the H theorem and entropy production.Meanwhile, the entropy production σ ent (r, t) is related to the collision term of the Boltzmann equation [21].Moreover, this production is also associated with energy dissipation.The treatment of the collision term of the Boltzmann equation requires the identification of a model of energy dissipation.

Conservation Laws in the NCCM
The Boltzmann equation fulfils the requirements of the mass, momentum, and energy conservation laws.By differentiating the densities of mass, momentum, and energy, the conservation equations of gas flows can be obtained: In this equation, P is the stress tensor and is defined by P = pI + ∆I + P. F i is the external force per unit mass.In the gas flows, the external force is usually omitted.Thus, the conservation laws are simplified into [22][23][24][25][26]: In this equation, P and Q are nonconservation variables with molecular expressions that do not yield a collisional invariant.

Evolution Equations of the Nonconservation Variables in NCCM
To derive the evolution equations for the nonconservation variables, we must define the velocity moment of order l R (abc...l) i Next, differencing this moment with time, we have The last collision term of Equation ( 20) is given by: In this equation, Finally, the derivative of the velocity moment can be obtained: Substituting the stress tensor and the heat flux term into Equation ( 22), we have Next, we also obtain the transport equation of the stress tensor and the heat flux: In this equation, In addition, the external force is usually omitted.Moreover, it can be seen that there are four terms in these transport equations: the derivative with time, the gradient of high-order terms, and the flux and collision parts.Usually, we use the viscous stress and the excess normal stress to replace the stress tensor.For this purpose, the following formulations should be defined in the molecular expressions for the moments: when the traceless part of the stress tensor is averaged over velocity, we denote them by the symbol where a = 1, 2, 3, and i = 1, 2, 3, . . ., n.The meaning of the term Φ (a) i is as follows: Thus, the transport equation of nonconservation variables can be obtained: Here, Λ k denotes the collision term and can be expressed as

Treatment of the Distribution Function in the Boltzmann Equation
At this stage, the most important work remaining is determining how to treat the collision term.The first step is the definition of the distribution function.The following primary concepts were implemented in this definition.
(1) The definition of the distribution function from nonequilibrium to equilibrium is weak in form and dynamic.In other words, the distribution function has no exact and specific expression in the process of nonequilibrium.
Entropy 2017, 19, 683 7 of 19 (2) The distribution function should fulfil not only the conservations of mass, momentum, and energy, but also the H theorem and positive entropy generation.
To fulfill the abovementioned rules and essential construction of the thermodynamic theory of irreversible processes, the definition of the distribution function is expressed via Eu Methods [17] as follows: In this equation, i = v i − u denotes the thermal motion velocity.Ψ i (r) denotes the potential energy.
is the molecular expressions for the moments of the nonequilibrium variable, and X α denotes the weighting coefficient.
In addition, in Equation ( 30), H i can be regarded as the local equilibrium, H denotes the nonequilibrium condition, and µ i connects the local equilibrium and nonequilibrium conditions to fulfill the conservation laws and positive entropy generation.
The advantage of this exponential form of the distribution function is obvious; in the physical sense, it is the only form that satisfies the additive property of the entropy and entropy production as well as the calortropy and calortropy production, all of which are in the logarithmic form of the distribution function; in the mathematical sense, it assures the non-negativity of the distribution function regardless of the level of approximations.
If we substitute this definition of the distribution function into the entropy balance equation, the entropy generation can be obtained: This form renders a clear, physical interpretation of entropy production: the entropy production is a direct measure of the energy dissipation arising from molecular collisions in the system, which in turn gives rise to a dissipative evolution of nonconserved macroscopic fluxes (moments).
In addition, the local equilibrium distribution function can be given in a form similar to the exponential form as follows: By applying the following dimensionless processes, The local equilibrium distribution function can be further written as Moreover, the molecular expressions of the viscous stress, excess normal stress, and heat flux can be rewritten again.
The unknown functions of X (α) i now must be calculated in terms of the macroscopic variables obeying the evolution equations of mass, momentum, and energy (conservation laws).We substitute the nonequilibrium canonical form of the distribution function into the Boltzmann equation and multiply by h (α) i to obtain the energy dissipation: In this equation, Λ (γ) i denotes the energy dissipation and represents the collision term of the transportation equation for nonconservation variables.
Consider the transport equation of Comparing this equation with the transportation equation for nonconservation variables, it can be rewritten as Equation ( 38) is equivalent to Equation (39), leading to the conclusion that Here, b is a constant.Equations ( 36) and (39) are both the transport equations of nonconservation variables coupled with the Boltzmann equation, with the distinction being (38); thus, b = 0, and we have In the lowest-order approximation of the dissipation term, by applying the model of small perturbation, we can set The moments of h (α) i are treated as orthogonal tensor Hermite polynomials; as a result, we have Finally, we obtain

Treatment of the Collision Term in the Boltzmann Equation
The distribution function can be rewritten in the following form: Thus, the entropy production can be expressed as In this instance, x i,j = x i + x j , y i,j = x * i + x * j .Next, we normalize the entropy production, In addition, we define The entropy production can be further expressed as If we neglect the second-or higher-order terms, then we have when the gas flows approach the near-equilibrium state from the nonequilibrium state, κ approaches zero.Therefore, we take the limit: In other words, when the gas flows approach the near-equilibrium state, In addition, the Rayleigh-Onsager dissipation function denotes the dissipation energy in near-equilibrium states; thus, κ 2 is the Rayleigh-Onsager dissipation function.
In this equation, From the above discussion, we conclude that entropy production is caused by nonlinear energy dissipation: When this dissipation approaches the near-equilibrium state, q(κ) = 1.In addition, the nonlinear energy dissipation is limited to the linear Rayleigh-Onsager dissipation function of κ 2 .In the far-from-equilibrium state, nonlinear factors q(κ) gradually increase.Of course, we can prove that q(κ) = 1 in the Grad equation; thus, the Grad moment method cannot solve the problem of the far-equilibrium state of gas flow.

Nonlinear Coupled Constitutive Method
The transport equation of nonconservation variable Φ (k) is: This new hydrodynamic equation can be rewritten as follows: The Discontinuous Galerkin (DG) scheme is used to solve the above equations; the specific numerical process is given as follows: (1) compute Π 0 , ∆ 0 , Q 0 (which are Newton's law of shear and bulk viscosity and Fourier's law of heat flux, respectively) based on conserved variables U; (2) update Π, ∆, Q based on Π 0 , ∆ 0 , Q 0 in Equation (58) by the Runge-Kutta method; (3) according to the updated Π, ∆, Q, update U in Equation (57) by the Runge-Kutta method; and (4) return to step (1) until the convergent error is satisfied.

Comparisons among the Viscous Stress in the NCCM, the Grad Moment Method, and the Burnett Equation
Figure 1 shows comparisons of the one-dimensional dimensionless normal viscous stress in the NCCM at steady state for the Grad moment method, the second-order Burnett equation, and the third-order Burnett equation.The following observations are made:

Results
In the section above, we can see that the NCCM equations approach the NSF equation for near-equilibrium states.In addition, this equation can be validated by the DSMC in far-equilibrium states.For validation, we first provide the verification and validation of the NCCM based on experimental results.Next, we present the test cases of hypersonic gas flows for both a low and a high Kn number.The comparison between the NCCM and the NSF equation at a low Kn number and the comparison between the NCCM and the DSMC method at a high Kn number were both performed in hypersonic gas flows.Furthermore, a comparison among the NCCM results, the DSMC, and the experiment involving three-dimensional (3D) complicated hypersonic vehicles in moment was also conducted.All convergent error is set to be 10 −6 .

Verification and Validation
For verification and validation of the present NCCM equations, the air flows around a NACA 0012 airfoil were first considered.The chord length is regarded as the characteristic length to define Re and the Kn number.To verify grid independence, as shown in Table 1, four sets of meshes are used

Results
In the section above, we can see that the NCCM equations approach the NSF equation for near-equilibrium states.In addition, this equation can be validated by the DSMC in far-equilibrium states.For validation, we first provide the verification and validation of the NCCM based on experimental results.Next, we present the test cases of hypersonic gas flows for both a low and a high Kn number.The comparison between the NCCM and the NSF equation at a low Kn number and the comparison between the NCCM and the DSMC method at a high Kn number were both performed in hypersonic gas flows.Furthermore, a comparison among the NCCM results, the DSMC, and the experiment involving three-dimensional (3D) complicated hypersonic vehicles in moment was also conducted.All convergent error is set to be 10 −6 .

Verification and Validation
For verification and validation of the present NCCM equations, the air flows around a NACA 0012 airfoil were first considered.The chord length is regarded as the characteristic length to define Re and the Kn number.To verify grid independence, as shown in Table 1, four sets of meshes are used to calculate the drag coefficient.Next, a set of unstructured triangle meshes, with 100 cells in the airfoil-conforming direction and 100 cells in the wall-normal direction, is implemented in the computational domain.Then, the smallest cell size is less than one percent of chord.Free stream conditions are applied in the outside circle, which has a radius that is 10 times the length of the chord.
Figure 2 presents the dimensionless velocity contours of the NSF equation and the NCCM together with the experimental [27] and DSMC results [8].All of these results capture the shock wave with a large thickness.The velocity in the after stream of the shockwave is under-predicted by the NSF equation.Both the NCCM and DSMC-IP (Information Preservation) results show good agreement with the experiment, especially in capturing the value of the 0.9 contour at the trailing edge.

Cavity Flow
Lid-driven cavity flows are studied at high Knudsen numbers.The cavity has a fixed wall temperature Tw = 290 K.The Knudsen number is defined as the ratio of the mean free path to the length of the cavity's side wall, and the mean free path of gas is evaluated for a hard sphere model.Similar to the previous work [28], the wall velocity is 300 m/s U   for cases at Kn = 0.671 and 6.712.

Cavity Flow
Lid-driven cavity flows are studied at high Knudsen numbers.The cavity has a fixed wall temperature T w = 290 K.The Knudsen number is defined as the ratio of the mean free path to the length of the cavity's side wall, and the mean free path of gas is evaluated for a hard sphere model.Similar to the previous work [28], the wall velocity is U ∞ = 300m/s for cases at Kn = 0.671 and 6.712.The solutions are compared with the DSMC ones.The computational domain is discretized using a mesh of 100 × 100 cells in physical space.The results of the present study and the DSMC at Kn = 0.671 and 6.712 are presented in Figures 3 and 4, where the temperature contours are shown.Excellent agreements are obtained.

Hypersonic Gas Flow in a Low Kn Number
In this section, two-dimensional (2D) NCCM is validated for hypersonic gas flows past a circular cylinder.The typical triangular mesh is used in 2D cases, and the outer circular radius that corresponds to a good computational domain is chosen as Router = 30 R. The cell size in unstructured mesh is 200 × 120.In other words, there are 200 points on the cylinder surface and 120 points in the radial direction of the computational domain.The linear element (P 1 ) and the quadratic element (P 2 ) are initially tested for the 2D case; because their numerical results are not distinguishable, the linear element is chosen for all 2D simulations.
The input parameters for the first case are Ma = 5.48, Kn = 0.001, and Pr = 0.607, and the working gas is argon.The contours of Ma number, density, and heat flux are shown in Figure 5; typical distributions around the cylinder are found.Good agreement between the NCCM and the NSF Tw=290 K.The Knudsen number is defined as the ratio of mean free path to the length of cavity side wall, 1 and the mean free path of gas is evaluated for a hard sphere model.Similar to the previous work [29], the wall 2 velocity is for cases at Kn=0.671 and 6.712.The solutions are compared with the DSMC ones. 3 The computational domain is discretized using a mesh of cells in physical space.Results of the 4 present study and the DSMC at Kn = 0.671 and 6.712 are presented in Figs.3-4, where the temperature contours

Hypersonic Gas Flow in a Low Kn Number
In this section, two-dimensional (2D) NCCM is validated for hypersonic gas flows past a circular cylinder.The typical triangular mesh is used in 2D cases, and the outer circular radius that corresponds to a good computational domain is chosen as Router = 30 R. The cell size in unstructured mesh is 200 × 120.In other words, there are 200 points on the cylinder surface and 120 points in the radial direction of the computational domain.The linear element (P 1 ) and the quadratic element (P 2 ) are initially tested for the 2D case; because their numerical results are not distinguishable, the linear element is chosen for all 2D simulations.
The input parameters for the first case are Ma = 5.48, Kn = 0.001, and Pr = 0.607, and the working gas is argon.The contours of Ma number, density, and heat flux are shown in Figure 5; typical T w =290 K.The Knudsen number is defined as the ratio of mean free path to the length of cavity side wall, 1 and the mean free path of gas is evaluated for a hard sphere model.Similar to the previous work [29], the wall 2 velocity is for cases at Kn=0.671 and 6.712.The solutions are compared with the DSMC ones. 3 The computational domain is discretized using a mesh of cells in physical space.Results of the 4 present study and the DSMC at Kn = 0.671 and 6.712 are presented in Figs.3-4, where the temperature contours

Hypersonic Gas Flow in a Low Kn Number
In this section, two-dimensional (2D) NCCM is validated for hypersonic gas flows past a circular cylinder.The typical triangular mesh is used in 2D cases, and the outer circular radius that corresponds to a good computational domain is chosen as R outer = 30 R. The cell size in unstructured mesh is 200 × 120.In other words, there are 200 points on the cylinder surface and 120 points in the radial Entropy 2017, 19, 683 14 of 19 direction of the computational domain.The linear element (P 1 ) and the quadratic element (P 2 ) are initially tested for the 2D case; because their numerical results are not distinguishable, the linear element is chosen for all 2D simulations.
The input parameters for the first case are Ma = 5.48, Kn = 0.001, and Pr = 0.607, and the working gas is argon.The contours of Ma number, density, and heat flux are shown in Figure 5; typical distributions around the cylinder are found.Good agreement between the NCCM and the NSF equation was observed.These factors include the position of the shock wave and the contour lines (both in the forebody and the afterbody).This case is for a low Kn number and in the region of near-equilibrium states.

Hypersonic Gas Flow at a High Kn Number
The input parameters for these two cases are Ma = 20, Kn = 0.01 and Ma = 20, Kn = 1.0.The working gas is argon.The heat fluxes around the solid surface are given in Figure 6.Good agreement of the results of the NCCM and the DSMC and UGKS (Unified gas-kinetic scheme) methods [29] was observed.These two cases are of a high Kn number and in the region of far-equilibrium states.Furthermore, a 3D NCCM is validated for hypersonic gas flows past 3D complicated hypersonic vehicles as shown in Figures 7 and 8.The typical tetrahedron mesh is used in 3D cases, and the outer circular radius for a good computational domain is chosen as Router = 30 R. In this equation, R is the characteristic size of the vehicles.The input parameters for the first case are Ma = 20.0,Kn = 2.531, and Pr = 0.707, and the working gas is air.Good agreement between the results of the NCCM and the DSMC method was observed.This case is of a high Kn number and in the region of far-equilibrium states.

Hypersonic Gas Flow at a High Kn Number
The input parameters for these two cases are Ma = 20, Kn = 0.01 and Ma = 20, Kn = 1.0.The working gas is argon.The heat fluxes around the solid surface are given in Figure 6.Good agreement of the results of the NCCM and the DSMC and UGKS (Unified gas-kinetic scheme) methods [29] was observed.These two cases are of a high Kn number and in the region of far-equilibrium states.

Hypersonic Gas Flow at a High Kn Number
The input parameters for these two cases are Ma = 20, Kn = 0.01 and Ma = 20, Kn = 1.0.The working gas is argon.The heat fluxes around the solid surface are given in Figure 6.Good agreement of the results of the NCCM and the DSMC and UGKS (Unified gas-kinetic scheme) methods [29] was observed.These two cases are of a high Kn number and in the region of far-equilibrium states.Furthermore, a 3D NCCM is validated for hypersonic gas flows past 3D complicated hypersonic vehicles as shown in Figures 7 and 8.The typical tetrahedron mesh is used in 3D cases, and the outer circular radius for a good computational domain is chosen as Router = 30 R. In this equation, R is the characteristic size of the vehicles.The input parameters for the first case are Ma = 20.0,Kn = 2.531, and Pr = 0.707, and the working gas is air.Good agreement between the results of the NCCM and the DSMC Furthermore, a 3D NCCM is validated for hypersonic gas flows past 3D complicated hypersonic vehicles as shown in Figures 7 and 8.The typical tetrahedron mesh is used in 3D cases, and the outer circular radius for a good computational domain is chosen as R outer = 30 R. In this equation, R is the characteristic size of the vehicles.The input parameters for the first case are Ma = 20.0,Kn = 2.531, and Pr = 0.707, and the working gas is air.Good agreement between the results of the NCCM and the DSMC method was observed.This case is of a high Kn number and in the region of far-equilibrium states.The present authors also validated hypersonic gas flows around a 3D Apollo 6 command module with Ma = 5.0, Kn = 0.5 in the earlier study [4], and good agreement between the results of the NCCM and the DSMC method was observed in this case.Here, we just give the results as can be seen in Figure 9 to prove that the NCCM can solve the problems of far equilibrium.More details can be seen in reference [4].The present authors also validated hypersonic gas flows around a 3D Apollo 6 command module with Ma = 5.0, Kn = 0.5 in the earlier study [4], and good agreement between the results of the NCCM and the DSMC method was observed in this case.Here, we just give the results as can be seen in Figure 9 to prove that the NCCM can solve the problems of far equilibrium.More details can be seen in reference [4].The present authors also validated hypersonic gas flows around a 3D Apollo 6 command module with Ma = 5.0, Kn = 0.5 in the earlier study [4], and good agreement between the results of the NCCM and the DSMC method was observed in this case.Here, we just give the results as can be seen in Figure 9 to prove that the NCCM can solve the problems of far equilibrium.More details can be seen in reference [4].
The present authors also validated hypersonic gas flows around a 3D Apollo 6 command module with Ma = 5.0, Kn = 0.5 in the earlier study [4], and good agreement between the results of the NCCM and the DSMC method was observed in this case.Here, we just give the results as can be seen in Figure 9 to prove that the NCCM can solve the problems of far equilibrium.More details can be seen in reference [4].

Discussion
The processes for the derivation of the NCCM equations from the Boltzmann equation were given in detail in the present study.The two features implemented in this derivation are discussed in detail.First, the H theorem and positive entropy generation are considered in the irreversible process.In addition, the definition of the distribution function is formal and dynamic.Second, based on the constraint of positive entropy generation, the unified highly nonlinear dissipation model from the nonequilibrium state to the equilibrium state is derived; this model turns into the Rayleigh-Onsager dissipation function near the equilibrium state.Next, the nonlinear coupled constitutive relations of viscous stress and heat flux, which strictly fulfil the second law of thermodynamics, were derived.In addition, coupled with the conservative laws of mass, moment, and energy, a unified scheme for continuum and rarefied gas flow was proposed.A comparison among the Grad moment method, the Burnett equation, and the NCCM was presented.In addition, the differences and relationships between the NCCM equations and the Navier-Stokes-Fourier equations were explained in detail.
For validation, numerical studies of rarefied and continuum gas flows were conducted.These studies include rarefied and/or continuum gas flows around a cylinder and a 3D model of a space shuttle.It was observed that the present results of the NCCM equations were in good agreement with those of the DSMC and UGKS methods in rarefied cases and were in good agreement with those of the Navier-Stokes equations in continuum cases.All these results indicate that the present NCCM equations provide a new approach for developing a unified scheme for solving continuum-rarefied gas flows.

Appendix A. The Derivation of the Nonlinear Parameter in the NCCM
The nonlinear parameter of R plays an important role in the NCCM equations.Note that the R in Equation ( 58) is different from that in [16].The following gives the modification.
c can be regarded as the parameter related the molecular model.

Appendix B. Positive Entropy Production in Eu Moment Equations
According to Equation (49), if the third order is included, the entropy production σent can be given by the form σent = κ exp 1 2 κ 2 − κ 2 sinh κ − κ 2 /2κ + κ 3 /6κ + κ 3 /3 The higher-order cumulant approximation must be taken in such a way that σent is always positive.Although the second-or third-order cumulant approximation for entropy production has potentially interesting features, rather involved calculations are required for them since tensors of high ranks appear in the integrals to be evaluated.Most of the studies done so far are based on the first-order cumulant approximation for σent , which is already highly nonlinear and able to account for interesting experimental observations as have been shown by the applications of the present theory.
Since first-order cumulant κ is positive, the entropy production given in Equation (50) clearly is positive and vanishes only at equilibrium where X (α) i = 0 for all i and α, as will be shown shortly.In fact, κ 2 turns out to be the Rayleigh-Onsager dissipation function if the 13 moments are taken for macroscopic variables.Since the Rayleigh-Onsager dissipation function is directly related to the entropy production in the linear theory of irreversible processes, Equation ( 50) is already a significant generalization of the nonlinear regime of the entropy production appearing in the linear theory.

( 1 )
The NCCM and the Burnett and Grad moment methods show nonlinear constitutive relations which is quite different from the NSF equation.However, the NCCM and the Burnett and Grad moment methods have the same linear relations as those of Newton's laws near the equilibrium state.In other words, the NSF equation can be regarded as the low-order approximation of NCCM, Burnett, and Grad.Further discussion can be found in[4,22].(2)All three of these nonlinear relationships are consistent in the near-equilibrium region.Because both of the latter methods can solve the problem of near-equilibrium gas flows, the NCCM should be correct in the near-equilibrium region.(3)The NCCM has a different nonlinear trend than those of the Burnett equation and the Grad moment method in the far-from-equilibrium region.Note that both the Burnett equation and the Grad moment method cannot be used in far-from-equilibrium states.In contrast, the NCCM shows the opposite trends as those of the Grad moment and Burnett methods.This finding has been validated by the DSMC method, as shown in Figure1.This feature indicates that the NCCM can be used to describe the gas flow in the far-from-equilibrium state.

Figure 1 .
Figure 1.Comparisons of the normal viscous stress in the nonlinear coupled constitutive method (NCCM) and the Grad, Burnett, and Navier-Stokes-Fourier (NSF) methods.

Figure 1 .
Figure 1.Comparisons of the normal viscous stress in the nonlinear coupled constitutive method (NCCM) and the Grad, Burnett, and Navier-Stokes-Fourier (NSF) methods.

Figure 5 .
Figure 5.The contours of Ma number, density, and Qx (heat flux in the × direction) of High Mach Number gas flow around a Cylinder at Ma = 5.48 and Kn = 0.001.NSF (Top), NCCM (Bottom).

Figure 5 .
Figure 5.The contours of Ma number, density, and Q x (heat flux in the × direction) of High Mach Number gas flow around a Cylinder at Ma = 5.48 and Kn = 0.001.NSF (Top), NCCM (Bottom).

Figure 5 .
Figure 5.The contours of Ma number, density, and Qx (heat flux in the × direction) of High Mach Number gas flow around a Cylinder at Ma = 5.48 and Kn = 0.001.NSF (Top), NCCM (Bottom).

Figure 9 .
Figure 9.The comparisons of temperature and pressure coefficients on Ma = 5.0, Kn = 0.5.

Figure 9 .
Figure 9.The comparisons of temperature and pressure coefficients on Ma = 5.0, Kn = 0.5.