Coupled constitutive relations: a second law based higher-order closure for hydrodynamics

In the classical framework, the Navier–Stokes–Fourier equations are obtained through the linear uncoupled thermodynamic force-flux relations which guarantee the non-negativity of the entropy production. However, the conventional thermodynamic descrip- tion is only valid when the Knudsen number is sufficiently small. Here, it is shown that the range of validity of the Navier–Stokes–Fourier equations can be extended by incorporating the nonlinear coupling among the thermodynamic forces and fluxes. The resulting system of conservation laws closed with the coupled constitutive relations is able to describe many interesting rarefaction effects, such as Knudsen paradox, transpiration flows, thermal stress, heat flux without temperature gradients, etc., which cannot be predicted by the classical Navier–Stokes–Fourier equations. For this system of equations, a set of phenomenological boundary conditions, which respect the second law of thermodynamics, is also derived. Some of the benchmark problems in fluid mechanics are studied to show the applicability of the derived equations and boundary conditions.

In the classical framework, the Navier-Stokes-Fourier equations are obtained through the linear uncoupled thermodynamic force-flux relations which guarantee the non-negativity of the entropy production. However, the conventional thermodynamic description is only valid when the Knudsen number is sufficiently small. Here, it is shown that the range of validity of the Navier-Stokes-Fourier equations can be extended by incorporating the nonlinear coupling among the thermodynamic forces and fluxes. The resulting system of conservation laws closed with the coupled constitutive relations is able to describe many interesting rarefaction effects, such as Knudsen paradox, transpiration flows, thermal stress, heat flux without temperature gradients, etc., which cannot be predicted by the classical Navier-Stokes-Fourier equations. For this system of equations, a set of phenomenological boundary conditions, which respect the second law of thermodynamics, is also derived. Some of the benchmark problems in fluid mechanics are studied to show the applicability of the derived equations and boundary conditions. expansion method [18], the Grad's moment method [19], regularized moment method [20,21], and entropy maximization [22,23]. Another approximation method, which derives a system of algebraic nonlinear coupled constitutive relations (NCCR) using the so-called balanced closure for the conservation laws was proposed by Myong [24].
The CE expansion method relies on the smallness of the Knudsen number. At zeroth-and first-order approximations the method leads to the Euler and NSF equations, and is fully equivalent to the results of LIT and RT (local equilibrium, etc.). At the second-and thirdorder approximations the method yields the Burnett and super-Burnett equations, respectively, which exhibit instabilities for time-dependent problems [25], and are thermodynamically inconsistent [26,27]. During the last decade, several modified forms of the Burnett equations have been suggested in the literature (e.g. [28][29][30][31]) that are stable; however, at present no proper boundary conditions are available for any of these sets of equations, hence their applicability is limited. On the other hand, the linearized Grad's 13 (also the R13) moment equations comply with the second law of thermodynamics, which also leads to thermodynamically consistent boundary conditions for these set of equations [32]. However, no entropy law has been established for the nonlinear moment equations.
In the following, we propose a phenomenological procedure in which the entropy flux contains a nonlinear contribution in stress and heat flux, which is motivated by results from rational extended thermodynamics (RET) [15]. In contrast to RET, where, similar to the moment method, the set of variables is enlarged to contain non-equilibrium variables, our approach still uses the basic equilibrium variables and fulfills the main conditions of local thermodynamic equilibrium (classical thermal and caloric equations of state, Gibbs equation). However, the additional term in the entropy flux yields additional terms in the entropy generation, which, as will be seen, can still be written as a sum of products of forces and fluxes, but now the fluxes appear explicitly in the forces, i.e. the fluxes are coupled through the additional terms in the forces.
We shall show that the conservation laws together with the resulting coupled constitutive relations (CCR) can capture many interesting non-equilibrium effects, such as Knudsen paradox, transpiration flows, thermal stress, heat flux without temperature gradients, etc., in good agreement with experiments and with kinetic theory, e.g. the solution of the Boltzmann equation.
The remainder of the paper is structured as follows. The derivation of the CCR for the conservation laws is detailed in §2. A thermodynamically consistent set of boundary conditions complementing the system of the conservation laws closed with the CCR (referred to as the CCR model hereafter) is presented in §3. The linear stability of the CCR model is analysed in §4 to show that the CCR model is stable to small perturbations. Classical flow problems of Knudsen minimum, heat transfer in an isothermal lid-driven cavity and normal shock structure are investigated in §5. The paper ends with conclusion in §6.

Derivation of coupled constitutive relations
The conservation laws, which are evolution equations for mass density ρ, macroscopic velocity v i and temperature T, read Here, D/Dt denotes the convective time derivative, p is the pressure, Π ik is the viscous stress tensor, q k is the heat flux and F i is the external force per unit mass. Throughout the paper, Einstein summation is assumed over the repeated indices unless stated otherwise. We shall consider monatomic ideal gases only, for which the pressure is p = ρRT with R being the gas constant and specific internal energy is u = c v T with the specific heat c v = 3R/2. Moreover, for monatomic ideal gases, the stress tensor is symmetric and tracefree, i.e. Π kk = 0 [2].
It should be noted that conservation laws (2.1)-(2.3) contain the stress tensor Π ik and heat flux q k as unknowns, hence constitutive equations are required, which link the these quantities to the variables ρ, v i and T.
The second law of thermodynamics states that the total entropy of an isolated system can never decrease over time [10,13], that is where s denotes the specific entropy, Ψ k is the non-convective entropy flux and Σ is the nonnegative entropy generation rate. An important aspect of a constitutive theory is to determine the appropriate relations among the properties s, Ψ k , Σ and the variables ρ, v i , T and their gradients so that the closed conservation laws guarantee the second law of thermodynamics.
(a) Uncoupled constitutive relations: the NSF equations In LIT, the constitutive relations for closing the system of conservation laws (2.1)-(2.3), are obtained such that the second law of thermodynamics is satisfied for all thermodynamic processes. For this, in LIT, local thermodynamic equilibrium is assumed [10,12], which implies the local validity of the Gibbs equation For convenience, we shall write temperature T in energy units as θ = RT, and dimensionless entropy as η = s/R, so that the Gibbs equation ( Comparison of equation (2.7) with equation (2.4) gives the LIT expression for non-convective entropy flux as Ψ k = q k /θ and that for the entropy production term as The angular brackets around indices represent the symmetric and traceless part of a tensor, for example, where δ ij is the Kronecker delta tensor. The entropy production (2.8) assumes a canonical form, i.e.
with the thermodynamic fluxes J α = {Π kl , q k } α and the thermodynamic forces heat flux must be Galilean invariant tensors [13], and it follows that only forces and fluxes of the same tensor type (scalars, vectors, 2-tensors, etc.) can be linked (Curie Principle [10]). Accordingly, where μ and κ are the coefficients of shear viscosity and thermal conductivity, respectively, where factors with θ are absorbed in the coefficients.
It is worth pointing out that for a monatomic ideal gas interacting with power potentials, the viscosity depends on temperature alone as where μ 0 is the viscosity at a reference temperature θ 0 and w is referred to as the viscosity exponent [2,16]. Furthermore, the heat conductivity is proportional to viscosity, κ = 5μ/(2 Pr), where Pr 2/3 denotes the Prandtl number [2]. Relations (2.10) 1 and (2.10) 2 are the Navier-Stokes law and Fourier's law, respectively, and we refer to them as the linear uncoupled constitutive relations-emanating from LIT. When relations (2.10) are substituted in conservation laws (2.1)-(2.3), they yield the well-known compressible NSF equations of hydrodynamics.

(b) Coupled constitutive relations
As mentioned in the Introduction, LIT and RT yield identical results for simple fluids-such as ideal gases-but differ in their assumptions. Indeed, RT assumes the entropy flux Ψ k = q k /θ that is an outcome in LIT. In a recent paper [33], Paolucci & Paolucci considered entropy flux as function of the hydrodynamic variables and their gradients in a complete nonlinear representation, but found that only the classical term Ψ k = q k /θ is compatible with the second law of thermodynamics. Nevertheless, allowing higher-order contributions for stress and heat flux, they found additions which are fully nonlinear in the gradients. Extended Irreversible Thermodynamics [15,34] is similar to LIT, only that non-equilibrium variables, in particular stress and heat flux, are added, with additional contributions to the Gibbs equation, in an attempt to go beyond local thermodynamic equilibrium. RET [15] proceeds differently, but also adds nonequilibrium variables, and yields a non-equilibrium Gibbs equation. The consideration of local non-equilibrium gives additions to the LIT entropy flux, which appears as explicit function of the extended set of variables.
In theories of (ET) for 13 moments, stress tensor Π kl and heat flux q i are field variables [15,34], and the entropy flux is expressed through these. Using representation theorems for isotropic tensor functions [13] and dimensional analysis, we find the most general entropy flux expression for these variables as where the dimensionless coefficients γ , α, β depend on the dimensionless invariants (Π 2 ) ll /p 2 , (Π 3 ) ll /p 3 , q 2 /(p 2 θ ) (note that the invariant Π kk = 0) [13]. Presently, we are only interested in the leading correction to the classical entropy flux q k /θ . Considering Π kl and q k as small and of the same order, Taylor expansion to second order yields where the coefficient of the first term on the right-hand side was chosen such that the classical result is reproduced, and α 0 is a constant. The form (2.13) of the entropy flux appears as a result in RET of 13 moments (with α 0 = 2/5) [15]. Interestingly, the gradient of the additional term α 0 (Π kl q l /pθ) can be expanded such that it yields higher-order additions to the thermodynamic forces. Indeed, introduction of this additional contribution in equation (2.7) yields an extended form of the second law that reads 14) The underlined terms on the left-hand side of equation (2.14) are equal to the underlined terms on the right-hand side of equation (2.14), where α 1 and α 2 are arbitrary numbers, and α * r = 1 − α r , r ∈ {1, 2}. The coefficients α r and α * r distribute the contributions to entropy generation between different force-flux pairs; their values will be obtained from comparison to results from kinetic theory. For α 0 = 0, the underlined terms in equation (2.14) vanish, and equation (2.7) is recovered.
The right-hand side of equation (2.14) is the entropy generation rate, which can-again-be recognized as a sum of products of the fluxes Π kl and q k with generalized thermodynamic forces (in square brackets). The corresponding phenomenological equations that guarantee positivity of the entropy production read Here, the coefficients were chosen such that in the classical limit (i.e. when α 0 = 0), the NSF equations are obtained. Since the fluxes Π kl and q k appear explicitly in the forces, we refer to relations ( The extended hydrodynamic models, e.g. the Grad 13-moment equations, contain additional contribution to the entropy density as well as to the fluxes. On the other hand, in the CCR model a second-order contribution to the entropy flux is considered, while the Gibbs equation, and hence the entropy density, remain unchanged. We emphasize that unlike the Grad 13-moment equations, the full Burnett equations cannot be obtained from the CCR model. Hence the CCR model stands somewhere between the NSF equations (first order in Kn) and the Burnett equations (second order in Kn).
With the extended entropy flux (2.13) and coupling force-flux relations, we find additional linear contributions to stress and heat flux, in contrast to [33] where all additional contributions are strictly nonlinear. The linear contributions are well known in kinetic theory, and are required to describe processes. For instance, the linear contribution ∂q k /∂x l in (2.15) describes thermal stresses, where temperature gradients can induce flow [4]. Similarly, the linear contribution ∂Π kl /∂x l in (2.16) describes stress-induced heat flux, as discussed in §5b.
Moment theories with 13 and more moments, as well as theories of ET, describe these effects as well, by providing full balance laws with time derivatives for higher moments. By contrast, the CCR model does not contain time derivatives. From this, one will in particular expect differences between the CCR model and ET in the description of transient processes, but only small differences for slow and steady-state processes. The benefit of the CCR model compared to available equations from ET and moment method is its full thermodynamic structure without restriction.
(c) Evaluation of the phenomenological coefficients  higher-order effects in gases with some accuracy. Hence, we determine the coefficients in the CCR from comparison with the Burnett equations. The Burnett equations are obtained from the CE expansion in the Knudsen number, which is proportional to the viscosity. The procedure can be easily applied to the constitutive relations (2.15) and (2.16) as follows: stress and heat flux are expanded in terms of the viscosity, so that Inserting this ansatz into equations (2.15) and (2.16), we find at first order the NSF laws and at second order where w is the exponent in the expression of viscosity for power potentials, see equation ( where i and θ i are the Burnett coefficients. Relations (2.21) yield the unknown coefficients as We note that the coefficients α 0 , α 1 , α 2 in the CCR can be fitted to the Burnett coefficients in agreement with the well-known relations between Burnett coefficients, θ 4 = 3 and 3 The Burnett coefficients depend upon the choice of intermolecular potential function appearing in the Boltzmann collision operator, values of these coefficients for inverse-power law potentials can be found, for example, in [2,36]. The values of the phenomenological coefficients α 0 , α 1 and α 2 for the hard-sphere (HS) and Maxwell molecule (MM) gases are given in table 1; for other power potentials, they can be computed from equations (2.21). We emphasize that we have performed the expansion (2.17) only to determine the coefficients α 0 , α 1 , α 2 , but will use the full CCR as given in (2.15) and (2.16).
We also note that, at least for Maxwell molecules, the Burnett equations can be obtained by CE-like expansion of the 13 moment system (based on ET or Grad, higher moment sets yield the same). Hence, in the sense of the CE expansion, the CCR and ET models agree in the transport coefficients for those terms that the CCR model contain, but the CCR model offers a reduced model, which has less (Burnett) contributions than the full 13 moment system. In the following sections, we show that the CCR model-just as Grad method and ET-provides linearly stable equations ( §4), and describe some classical flow problems in sufficient agreement to detailed kinetic models ( §5).

Wall boundary conditions
Just as the process of finding constitutive relations in the bulk, the development of wall boundary conditions is based on the second law. Specifically, one determines the entropy generation at the interface, and finds the boundary conditions as phenomenological laws that guarantee positivity of the entropy generation.
The entropy production rate at the boundary Σ w is given by the difference between the entropy fluxes into and out of the surface [10], i.e.
Here, n k is unit normal pointing from the boundary into the gas , q w k denotes the heat flux in the wall at the interface, and θ w denotes the temperature of the wall at the interface. Here, the wall is assumed to be a rigid Fourier heat conductor, with the entropy flux q w k /θ w and Π w ik = 0. At the interface, the total fluxes of mass, momentum and energy are continuous, due to conservation of these quantities [10,32,37,38], where all quantities with superscript w refer to wall properties, and the others refer to the gas properties. To proceed, we combine entropy generation and continuity conditions by eliminating the heat flux in the wall q w k and the pressure p w , and find, after insertion of the entropy flux (2.13), Here, V i = v i − v w i is the slip velocity, with V i n i = 0, and T = θ − θ w is the temperature jump. To write the entropy generation properly as sum of products of forces and fluxes, it is necessary to decompose the stress tensor and heat flux into their components in the normal and tangential directions as [32] Π ij = Π nn ( 3 2 n i n j − 1 2 δ ij ) +Π ni n j +Π nj n i +Π ij , and q i = q n n i +q i , (3.4) where q n = q k n k , Π nn = Π kl n k n l , and such thatq k n k =Π nk n k =Π kk =Π ij n j =Π ij n i = 0, see appendix A for more details. Substituting equations (3.4) and (3.5) into equation (3.3), the entropy generation can be written as a sum of two contributions:  where P = p − α 0 Π nn . As in LIT of the bulk, the positivity of entropy production is ensured by phenomenological equations Here, ς 1 and ς 2 are non-negative coefficients, which can be obtained either from experiments or from kinetic theory models. We determine ς 1 and ς 2 from the Maxwell accommodation model [39]. This model employs the accommodation coefficient χ , which is defined such that a fraction χ of the molecules hitting the wall returns to the gas in equilibrium with the wall properties (wall Maxwellian), whereas the remaining fraction (1 − χ ) is specularly reflected. Comparison with boundary conditions from this model shows that these coefficients are related to the velocity-slip coefficient η VS and the temperature-jump coefficient η TJ as [32] The values of the velocity-slip coefficient η VS and the temperature-jump coefficient η TJ as obtained from the linearized Boltzmann equation [40][41][42][43] are given in table 2. Conditions (3.6) are the required, and thermodynamically consistent, boundary conditions for the CCR model. The first term in the brackets of boundary condition (3.6) 1 relates the shear stress (Π ni ) to the tangential velocity slip (V i ) while the second term describes thermal transpiration-a flow induced by the tangential heat flux. It is evident from boundary condition (3.6) 1 that the NSF equations, for which α 0 = 0, do not allow thermal transpiration within the framework of thermodynamically consistent boundary conditions. Boundary condition (3.6) 2 relates temperature jump (T ) and viscous heating (Π ni V i ) to the normal heat flux q n [44].
Several authors have suggested second-order boundary conditions for the Navier-Stokes equations, which are in a form, possibly with other values of coefficients, identical to (3.6) [45]. Our derivation suggests that in order for the boundary conditions to be compatible with the second law of thermodynamics, the NSF equations do not suffice, while they arise quite naturally from the CCR.
We note that a full theory of boundary conditions is still missing in the theory of ET. Indeed, since the transport equations in ET are hyperbolic, one has to study the characteristics of the system in order to determine how many, and which, boundary conditions are required. We are not aware of complete sets of boundary conditions for nonlinear ET at present. For moment equations and linearized version of ET, one can use ideas from kinetic theory to determine boundary conditions. The approach suggested by Grad [19] violates Onsager symmetry conditions, and ideas as those above can be used to adjust some coefficients in order to guarantee the proper symmetries [32].

Linear stability analysis
In this section, we analyse the stability of the CCR model (

(a) Linear dimensionless equations
For linear stability analysis, the CCR model is linearized by introducing small perturbations in the field variables from their values in an equilibrium rest state. We write where ρ 0 and θ 0 are the values of ρ and θ in the equilibrium while the remaining variables vanish in the equilibrium rest state; hats denote the dimensionless perturbations in the field variables from their values in the equilibrium rest state; and √ θ 0 is the scale for making the velocity dimensionless. These dimensionless perturbations are assumed to be much smaller than unity so that the linear analysis remains valid. Consequently, the pressure is linearized as p = p 0 (1 +p) with p 0 = ρ 0 θ 0 andp ≈ρ +θ . Furthermore, a relevant length scale L is introduced so that the dimensionless space variable isx i = x i /L and the dimensionless time ist = t √ θ 0 /L. Also, the external force is assumed to be small (of the order of small perturbations) in a linear analysis, and the dimensionless external force is given byF i = F i L/θ 0 . Substituting the values of field variables from equations (4.1) in the CCR model, introducing the dimensionless space and time variables, and retaining only the linear terms in the perturbation variables, one obtains with Π kl = −2 Kn ∂v k ∂x l + α 0 ∂q k ∂x l and q k = − 5 2

Kn Pr
In equations (4.2)-(4.5) and henceforth, the hats are dropped for better readability, and is the Knudsen number with μ 0 being the viscosity of the gas in equilibrium. Thus, all quantities in equations (4.2)-(4.5) and henceforth are dimensionless unless otherwise mentioned.

(b) Plane harmonic waves
Now we consider a one-dimensional process (in the x-direction) without any external forces (i.e. F i = 0) and assume a plane wave solution of the form for equations (4.2)-(4.5), where ψ = {ρ, v x , θ, Π xx , q x } , • ψ is a vector containing the complex amplitudes of the respective waves, ω is the (dimensionless) frequency and k is the (dimensionless) wavenumber. Substitution of the plane wave solution (4.7) into equations (4.2)-(4.5) gives a system of algebraic equations Aψ = 0, where  For non-trivial solutions, the determinant of matrix A must vanish, i.e. det (A) = 0. This leads to the dispersion relation, which is the relationship between ω and k: For temporal stability, a disturbance is considered in space; consequently, the wavenumber k is assumed to be real while the frequency ω can be complex. From the plane wave solution (4.7), it is clear that temporal stability requires the imaginary part of the frequency to be non-negative, i.e. Im(ω) ≥ 0, for all wavenumbers. In other words, if Im(ω) < 0, then a small disturbance in space will blow up in time. Figure 1 illustrates the stability diagram for the CCR model, NSF equations and Burnett equations shown by blue, red and green colours, respectively. The results for the NSF equations and Burnett equations are included for comparison. Assuming Maxwell molecules, we have Pr = 2 3 and α 0 = 2 5 ; the Knudsen number is set to Kn = 1, that is the mean free path is used as relevant length scale. The figure exhibits ω(k) = Re(ω) + i Im(ω)-obtained from the dispersion relations for the CCR model, NSF equations and Burnett equations-in the complex plane with k as parameter. The grey region in the figure is the region, in which the stability conditions are not fulfilled, hence it depicts the unstable region, whereas the white region represents the stable one.
It is evident from figure 1 that the NSF equations (red) as well as the CCR model (blue) are linearly stable in time for all wavenumbers since the roots of their dispersion relations always have non-negative imaginary parts; on the other hand, the Burnett equations (green) are unstable in time, with negative roots for large wavenumbers. Grad-type moment systems and ET systems are linearly stable (e.g. [2,46]).

Classical flow problems (a) Knudsen minimum
The Knudsen minimum is observed in a force-driven Poiseuille flow of rarefied gases, in which, for given force, the mass flow rate of a gas first decreases with the Knudsen number, attains a minimum value around a critical Knudsen number and then increases with the Knudsen number. We shall investigate this problem through the CCR model.
Let us consider the steady state (∂(·)/∂t = 0) of a gas confined between two isothermal, fully diffusive walls of a channel. Let the walls be located at (dimensionless positions) y = ∓1/2 and be kept at a (dimensionless) temperature θ w = 1. The flow is assumed to be fully developed and driven by a uniform (dimensionless) external force F in the positive x-direction parallel to the walls; all field variables are independent of x; and the velocity component in the y-direction is zero, i.e. v y = 0. The problem can be described through the (linear-dimensionless) CCR model (equations (4.2)-(4.5)). For the problem under consideration, the mass balance equation (4.2) is identically satisfied and the rest of the equations simplify to Kn α 0 ∂q y ∂y , and Kn Pr Boundary conditions (3.6) at the walls (i.e. at y = ∓ 1 2 ) in the linear-dimensionless form reduce to ± Π xy = −ς 1 (v x + α 0 q x ) and ± q y = −ς 2 α 0 Π yy . 1 Consequently, the mass flow rate of the gas is Here, the additional factor of 1/ √ 2 in the mass flow rate is included in order to compare the present results with those obtained from the linearized Boltzmann equation (LBE) in [5] where the authors scale velocity by √ 2θ . As the walls of the channel are fully diffusive, the accommodation coefficient χ is unity, and hence from equation (3.7) 1 , ς 1 = ( √ 2/π )/η VS . The mass flow rate from the NSF equations with slip boundary conditions is obtained by setting α 0 = 0, which gives (1/2 √ 2)F((1/6)(1/Kn) + (1/ς 1 )), while for NSF with no-slip boundary conditions one finds (1/12 √ 2)(1/Kn)F. Figure 2 shows the mass flow rate of a hard-sphere gas plotted over the Knudsen number Kn = 4 √ 2 Kn/5 for F = 1, as obtained from the CCR model (blue solid line), from the NSF equations with (red dashed line) and without (magenta solid line) slip boundary conditions, and from the LBE (symbols). Again, the Knudsen number Kn is multiplied with a factor of 4 √ 2/5 in order to compare the results with those from [5]. The parameters Pr and η VS for hard-spheres are taken from tables 1 and 2, respectively. It is clear from the figure that the CCR model predicts the Knudsen minimum and closely follow the mass flow rate profile from the LBE up to Kn 1. On the other hand, the NSF equations with or without slip do not predict the Knudsen minimum at all; in particular, the mass flow rate from the NSF equations with no-slip boundary conditions is not even close to that from the LBE equations. The mass flow rate from the CCR model

(b) Lid-driven cavity flow
The lid-driven cavity is a well-known test problem in rarefied gas dynamics, in which a gasunder no external force-is confined to a square enclosure of (dimensionless) length 1. The boundaries at x = 0, x = 1 and y = 0 are stationary while the upper boundary at y = 1 is moving in the x-direction with a velocity v lid . The boundaries are kept at a constant temperature, which is equal to the initial temperature of the gas inside the cavity so that the dimensionless temperature of the walls θ w = 1. The flow in the cavity is assumed to be in steady state and independent of the z-direction, i.e. ∂(·)/∂t = 0 and ∂(·)/∂z = 0. The problem is solved through the CCR model ((2.1)-(2.3) along with (2.15) and (2.16) in dimensionless form) numerically using a finite difference scheme whose details can be found in [9]. Figure 3 shows the vertical and horizontal components of the velocity along the horizontal and vertical centrelines of the cavity, respectively, computed through the CCR model and the NSF equations for Kn = 0.1/ √ 2 and v lid = 0.21 (50 m s −1 in dimensional units). The results are compared to DSMC simulations for hard spheres, hence we chose the phenomenological coefficients for hard spheres from table 1. The velocity profiles from the CCR model as well as from the NSF equations are in good agreement with the DSMC simulations [8].
It is well known that the classical NSF laws cannot describe heat transfer phenomena in a liddriven cavity (see e.g. [9,47,48] and references therein). Therefore, either extended models, such as Grad-type moment equations or the R13 equations [2,17], or other advanced constitutive laws ought to be used for the closure in order to describe the heat transfer characteristics.
We study in particular the temperature field and heat flux induced in the lid-driven cavity and compare the results from the CCR model to those from DSMC presented in [8]. Figure 4 illustrates the heat flux lines plotted over the temperature contours, and compares the predictions of (from left to right) the CCR model, DSMC and the NSF equations. The NSF equations predict that the heat flows from hot (top-right corner) to cold (top-left corner) in an orthogonal direction to the temperature contours. However, both the CCR model and DSMC predict heat flux from cold region to hot region, which is non-Fourier effect.  The non-Fourier heat flux can easily be understood from the linear terms itself, which dominate the nonlinear terms. From the (linearized) momentum balance equation (4.3) (ignoring the time derivative and external force terms), the divergence of the stress is given by Substitution of this in the (linearized) constitutive relation for the heat flux (4.5) 2 yields Kn Pr The first term on the right-hand side of equation (5.6) is the Fourier's contribution to the heat flux while the second term is the non-Fourier contribution to the heat flux due to the pressure gradient, which is responsible for heat transfer from the cold region to the hot one. For one-dimensional normal shock structure, one considers a reference frame moving with the shock and hence all the quantities depend only on a single variable ξ = x − v x t [49,50]. The dimensionless hydrodynamic variables (denoted with bars) in upstream (ξ → −∞) areρ u = 1, v u = √ 5/3 Ma,θ u = 1, wherev 0 is the dimensionless velocity in the x-direction and Ma is the Mach number. Furthermore, the dimensionless hydrodynamic variables in downstream (ξ → ∞),ρ d ,v d andθ d , follow from the Rankine-Hugoniot conditions [46,49]: The shock profiles are obtained by solving the different models (CCR, NSF and Grad) numerically using a central difference scheme employed in [46]. Figure 6 illustrates the profiles for (a) density and (b) entropy density for Ma = 2, w = 1 and Kn = 1 compared with the DSMC data given in [46]. In order to make ξ dimensionless, the mean free path of the upstream region λ 0 = 0.0014 m [2,46] has been used. Both the CCR (blue solid line) and NSF (red dashed line) models give smooth shock profiles while the Grad  not observed in the DSMC simulations (symbols). The sub-shock in Grad 13-moment equations is due to the hyperbolicity of the equations and occurs above the characteristic speed [46,49]. Clearly, both the CCR and NSF models mismatch the DSMC results. Nevertheless, the density profile from the CCR model agrees well with the DSMC results in the upstream part of the shock while being overpredicted in the downstream of the shock whereas that from the NSF equations has clear mismatch from the DSMC simulations both in the upstream and downstream parts of the shock. Notwithstanding, it is worthwhile to note that one of the main assumption made in the CCR model is the validity of local equilibrium, which implies that the entropy density η depends only on the equilibrium variables and is given by the Gibbs equation (2.5), i.e.
This assumption may not be valid in strong non-equilibrium processes, such as problems involving shocks. Figure 6b exhibits the entropy density, which turns out to be non-monotonous for both the CCR (blue solid line) and NSF (red dashed line) models. Here we want to emphasize that while entropy is not monotonous, it certainly has positive production. Although, we could not find the Boltzmann/DSMC solution for the entropy density in the literature, it has been reported, for instance in [50], that the entropy density computed from the Boltzmann equation is monotonous. This again asserts that the local equilibrium assumption for entropy may not be valid for strongly non-equilibrium processes, and non-equilibrium variables must also be included in the expression of entropy (e.g. [33,50]) along with the entropy flux.

Conclusion
Combining ideas of different approaches to Irreversible Thermodynamics, in particular LIT, RT and RET, we derived an improved set of constitutive relations for stress tensor and heat fluxthe CCR. The model describes processes in mildly rarefied gases in sufficient approximation, and reproduces important rarefaction effects, such as the Knudsen minimum, and non-Fourier heat transfer, which cannot be described by classical hydrodynamics (NSF). By construction, the resulting transport equations are accompanied by a proper entropy inequality with non-negative entropy generation for all processes and are linearly stable. This is a clear distinction to other models for rarefied gases, such as the Burnett equations, which are unstable due to lack of a proper entropy, or Grad-type moment equations, which are accompanied by proper entropy inequalities, and are stable, only in the linear case [15].
Thermodynamically consistent boundary conditions for the CCR model have been developed as well, which describe velocity slip, temperature jump and transpiration flow at the boundaries.
The CCR add several higher order terms to the NSF system, in bulk and at the boundary. The CCR model is accompanied by an entropy inequality, in which the entropy remains the equilibrium entropy as integrated from the equilibrium Gibbs equation, but entropy flux and entropy generation exhibit higher order correction terms. The model gives a good description of some important rarefaction effects, but does not provide as fine resolution as the full Boltzmann equation, or higher order moment equations (e.g. the R13 and R26 equations [17,21]). In particular, Knudsen layers are not resolved, and only appear indirectly in the corrected jump and slip coefficients.
The development of macroscopic transport equations for rarefied gases at larger orders in the Knudsen number with a full formulation of the second law of thermodynamics is an important project within the field of non-equilibrium thermodynamics. Our results for the normal shock problem show that the assumption of local equilibrium may not be legitimate since problems involving shocks are intrinsically in strong non-equilibrium. A possible way to resolve this issue is to consider an extended form of entropy-involving non-equilibrium variables-along with non-equilibrium contribution to the entropy flux, which is planned for future work. On the other hand, one will accept an approximation of the second law, if the system of equations provides sufficient accuracy for the description of processes. For instance, the R13 equations provide good accuracy, but have a proven second law only in the linearized case. While this guarantees linear stability, little can be said about the nonlinear behaviour. On the other hand, cases where the full second law-linear and nonlinear-is enforced, are either not amenable to analytic closure [23], or not sufficiently accurate [51].
In steady state, the linearized CCR model reduces to the linearized Grad's 13-moment equations, which have been studied extensively in the literature. In particular, Green's functions solutions were obtained for the steady-state linearized Grad equations by Lockerby & Collyer [52]. The numerical framework based upon these fundamental solutions can readily be implemented for the linearized CCR model allowing for three-dimensional steady computation at remarkably low computational cost. The appropriate numerical recipes for the nonlinear CCR model would be challenging-especially due to the non-local coupling of stress and heat flux-which will be the subject of future research.
The successful development of thermodynamically consistent transport equations for rarefied gases at higher orders will only be possible by using the best of several approaches to irreversible thermodynamics and new ideas. We hope that the development of the CCR system based on ideas of LIT, RT and RET will be a useful step towards this end. We also envisage that these and similar ideas will prove useful for further development of recent moment models for monatomic gas mixtures [53,54] and granular gases [55].
Data accessibility. The research material (source code and data) can be accessed in the electronic supplementary material.