Thermodynamically admissible 13-moment equations

Closure is an important issue in developing simplified physical theories from more fundamental levels of description. Closure is often perceived as a problem but, more positively, it can also be seen as the crucial first step toward describing and understanding the essence of certain phenomena of interest by identifying appropriate variables.¹ Closure is expected to be possible whenever a set of variables allows for an autonomous description of a problem of interest.

In the context of fluid dynamics, the densities of the locally conserved mass, momentum, and energy provide an extremely successful level of description (the five fields of hydrodynamics). However, when we are interested in time and length scales comparable to molecular scales, for example in rarefied gas dynamics, then we need more detailed levels of description.² The Boltzmann equation for the single-particle velocity distribution provides an excellent, very detailed level of description for rarefied gas dynamics. However, for many applications, the Boltzmann equation is unnecessarily detailed and too complicated to solve. As standard hydrodynamics can be recovered from the five moments of the single-particle velocity distribution that correspond to the mass, momentum, and energy densities, it is tempting to use larger sets of moments to construct theories of fluid dynamics going beyond standard hydrodynamics.³ Then, the issue of moment closure arises from the fact that the Boltzmann equation is equivalent to an infinite hierarchy of coupled moment equations.^2,3

How can we recognize suitable sets of moments and valid closure procedures? An obvious criterion, which we refer to as empirical adequacy, requires that all relevant experimental results are reproduced by autonomous moment equations. Instead of experimental data, one might use the well-established Boltzmann equation or computer simulations of rarefied gases as a reference for empirical adequacy. However, there are further requirements from fundamental principles, such as the validity of a suitable version of the second law of thermodynamics for nonequilibrium systems. We refer to this type of criteria as thermodynamic admissibility, where we here rely on general equation for the nonequilibrium reversible–irreversible coupling (GENERIC)^4–6 as a thermodynamic framework.

We consider both empirical adequacy and thermodynamic admissibility as necessary conditions for successful closure. Whereas none of these criteria by itself provides a sufficient condition, the combination of these two criteria offers a promising indicator for identifying successful closures. Of course, thermodynamic admissibility as a further constraint restricts the options for fitting experimental data. Yet, empirical adequacy should never be inferred from a perfect fit of a particular set of experimental data; it should be recognized in the overall performance of a model for all kinds of possible experiments and, in particular, in a complete coverage of all relevant qualitative phenomena.

In producing simplified theories of rarefied gas dynamics from the Boltzmann equation, it is natural to supplement the five hydrodynamic fields associated with the local conservation laws by five-second moments characterizing the momentum flux, and three third moments characterizing the energy flux. Whereas a lot is known about the empirical adequacy of the resulting 13-moment equations, starting with the work of Grad,^2,7,8 very little is known about their thermodynamic admissibility. The construction of 13-moment equations that are consistent with the GENERIC framework of nonequilibrium thermodynamics is the topic of the present paper.

One key issue for thermodynamic admissibility is the existence of an entropy function. In systematic expansion techniques for the Boltzmann equation, such as Grad's moment method,^7,8 one typically loses the positivity of the probability density for the velocity distribution and hence also Boltzmann's entropy involving the logarithm of that probability density. By using anisotropic Gaussian trial functions, this problem can be solved for 10-moment equations.⁹ Also for the linearized 13-moment equations an entropy has been found.¹⁰ While the framework of Rational Extended Thermodynamics³ aims at a thermodynamic structure for transport equations, approximations are used that are akin to the Grad method, so that the second law is met only approximately. The successful construction of an entropy for the nonlinear 13-moment equations requires not only suitable trial functions but also proper treatment of the general tensor structure of all the moments (see Ref. 11, refined through the discussion in Refs. 12 and 13).

Another, more subtle key issue is the Hamiltonian structure of the convection mechanism. Assuming the existence of a Poisson bracket for the 10-moment case, its necessary form has been derived on p. 309 of Ref. 6. However, the Jacobi identity for the resulting Poisson bracket had not been tested at the time. Later tests showed that, contrary to what was expected in Ref. 6, one actually does not arrive at a valid Poisson bracket. For the 10-moment case, this problem has recently been solved in Ref. 14 by introducing an irreversible contribution with Casimir symmetry to correct for velocity slip. In the present paper, the same strategy is applied to the 13-moment equations. The structure of the resulting thermodynamically consistent 13-moment equations is then compared to that of previous empirically adequate 13-moment equations lacking thermodynamic admissibility. Specifically, the GENERIC building blocks are designed such that they reproduce most elements of Grad's 13-moment equations within second-order accuracy, but contain additional higher order terms that give thermodynamic structure, which the Grad equations do not possess.

The structure of this paper is as follows: In Sec. II, we review Grad's original 13-moment equations, including a brief discussion of entropy and Knudsen scaling. Section III begins with a summary of the essentials of the GENERIC framework and then provides all the building blocks for a new, thermodynamically admissible version of the 13-moment equations. A detailed discussion of the resulting equations is offered in Sec. IV, including the entropy balance, the proper hydrodynamic limit, and second-order matching in the Knudsen scaling with Grad's original equations, which appears to not be completely possible. A short discussion and outlook concludes the paper. A class of positive-semidefinite matrices that is useful for the modeling of relaxation is introduced in  Appendix A. Finally,  Appendix B offers a brief discussion of Onsager–Casimir symmetry.

We shall not describe the derivation of the Grad 13-moment equations, but give only a summary of their features. Detailed accounts of derivation and properties of the equations are available in the literature.^2,3,7,8 We will use tensor index notation with the Einstein summation throughout the paper. The 13 variables are moments of the distribution function $f(r,t,c)$⁠, where r and t are the space and time variables, respectively, and c denotes the microscopic velocities of particles. Specifically, we consider mass density ρ, momentum density $Mi=ρvi$⁠, where v_i is center of mass velocity of the gas, and specific internal energy ε, which are given as

Here, m is the mass of a gas particle, $θ=kBmT$ is the temperature in specific energy units with the Boltzmann constant $kB$ and the thermodynamic temperature T, v_k is center of mass velocity, and C_k = c_k − v_k is the so-called peculiar velocity that is the particle velocity in the co-moving frame. The trace-free and symmetric part of the second moment (indicated by angular brackets) defines the kinetic stress tensor as

and the heat flux is the trace of the third moment

The 13 variables of the original Grad-13 theory are

Within GENERIC, we shall use the equivalent set of variables

where the vector w_i will be related to the heat flux q_i, and the temperature tensor Θ_ij is defined through the pressure tensor $Pij$ as

which is conveniently split into its trace and trace-free parts. Here, $p=ρθ$ is the ideal gas pressure; we note that temperature is the trace of the temperature tensor

and that Θ_ij and P_ij are positive definite by definition.

The total energy density is the sum of internal and kinetic contributions,

The closure of the moment equations relies on an approximation of the distribution function that is entirely specified through the chosen variables (⁠$u¯A$ or u_A). Grad constructed his approximation such that it reproduces the 13 variables $u¯A$⁠, as the product of the equilibrium phase density—the Maxwellian—and a polynomial in C_i that accounts for deviation from equilibrium

With this, all other moments become constitutive functions of the chosen variables.

By definition, a distribution function is positive. Due to its polynomial structure, the Grad approximation is not strictly positive, since the bracket in (9) becomes negative for large values of C_i. This behavior is compensated somewhat by the Maxwellian, which suppresses the polynomial for large C_i. With this, the Grad approximation remains meaningful as long as stress σ_ik and heat flux q_i remain sufficiently low.³ 

The Grad-13 moment equations arise by multiplying the Boltzmann equation with the appropriate velocity polynomials and subsequent integration, using the Grad-13 moment distribution to determine all moments that appear.^2,3,7,8,14 This yields the conservation laws for mass and momentum, and the balance law for internal energy (which together with momentum conservation yields conservation of total energy), in their well-known forms²² 

and full balance laws for stress and heat flux, as

The production terms on the right-hand side are for Maxwell molecules, and μ is the shear viscosity of the gas.

Within kinetic gas theory, the Boltzmann entropy density, the non-convective entropy flux of the gas, and the entropy production rate are given by^2,15,16

with a scaling factor $y=e(mh)3$⁠, where e and h are Euler's and Planck's constants, respectively, and the Boltzmann collision term $S(f,f)$⁠. The corresponding balance of entropy, also known as H-theorem, reads

Evaluation with the Maxwellian gives the well-known equilibrium entropy density of the ideal gas

with the constant $η0=[52+ln (2π3m4h3)]$⁠, while flux and production vanish (strictly speaking, the logarithmic terms should be combined because the argument of the ln-function must be dimensionless).

The entropy (13) is defined for arbitrary nonequilibrium states, which have positive distribution functions. With its loss of positivity, the logarithm of the Grad distribution $f|13$ in (9) does not exist for all C_i and the entropy density can only be determined by approximation for small stress and heat flux through Taylor expansion of $ln f|13$⁠. To second order, entropy, entropy flux, and entropy production are found as³ 

Due to the approximations, insertion of these expressions into the entropy balance (14) does not produce a proper second law for the Grad-13 equations. Indeed, the above quadratic approximations in the nonequilibrium variables (σ_kl, q_k) yield a suitable second law for the linearized Grad-13 moment equations,^10,17 but not for the full nonlinear Grad-13 system (10)–(12).

We are not aware of a successful attempt to find a complete entropy inequality to accompany the Grad-13 equations (or any other Grad-type system). The loss of positivity of the distribution function, and the corresponding loss of entropy, removes proper physics when the equations are overstretched. Indeed, while close to equilibrium the Grad-13 equations are symmetric hyperbolic, far from equilibrium, that is, for large stresses and heat fluxes, they lose hyperbolicity and instead have complex eigenvalues,³ which leads to unphysical solutions. This implies that the Grad equations do not exclude solutions, which might violate the H-theorem. This does not render the equations meaningless from the perspective of empirical adequacy: The Grad system aims at approximating the Boltzmann equation through moment equations, and with this, also the H-theorem is only approximated. The equations describe gas rarefaction effects such as thermal stresses, or non-gradient heat fluxes, within limits of sufficiently small stress σ_kl and heat flux $qk$⁠.

For meaningful solutions, the Grad equations should remain hyperbolic. Müller and Ruggeri³ demonstrate that hyperbolicity is lost for large (σ_kl, q_k) and determine a hyperbolicity radius that estimates limits on (σ_kl, q_k). Proper thermodynamic structure—a valid H-theorem—and mathematical structure—hyperbolicity—are related; hence, the breakdown of hyperbolicity gives a strong hint that far from equilibrium Grad-13 will violate the second law.

Scaling arguments are helpful for a deeper understanding of the limitations of Grad-13. Indeed, considering scaling with the Knudsen number $Kn$ (the ratio of mean free path and desired length scale of resolution) as smallness parameter, it can be shown that Grad-13 is the proper model when one is interested in second-order accuracy $O(Kn2)$⁠, while higher order accuracy must account for additional moments and their equations.^2,18,19 We here present elements of the arguments as required to proceed.

The Grad-13 equations describe transport in mildly rarefied gases. They reduce to the classical Navier–Stokes equations by means of the Chapman–Enskog (CE) expansion, which relies on small mean free path of the gas molecules. Since the mean free path is proportional to the viscosity $μ$⁠, we can rescale the production terms for stress and heat flux as $−1ϖpμσij,−231ϖpμqi$⁠, with a formal smallness parameter ϖ that plays the role of the Knudsen number; the parameter can be set to unity after the expansion. For the CE expansion, stress and heat flux are written as a series

Next, the expansions are inserted into the respective transport equations, and matching of powers of ϖ yields the leading contributions for stress and heat flux as

When considering stress and heat flux only to first order, we obtain the well-known equations of Navier–Stokes and Fourier with shear viscosity μ and heat conductivity $154μ$⁠. Note that both transport coefficients have the same dimension due to our use of $θ=kBmT$ as measure for temperature.

For the following, it will be important to note that both, σ_ij and q_i, have no zeroth order—that is, equilibrium—contributions, but non-vanishing first-order contributions. Hence, we can say, both are of first order in the Knudsen number, that is, $O(ϖ)$⁠. This understanding will now be employed to discuss the full Grad's 13-moment equations (11) and (12) for stress and heat flux.

We recognize the first-order characteristics formally by writing $σij=ϖσ¯ij, qi=ϖq¯i, and μ=ϖμ¯,$ where $σ¯ij, q¯i, μ¯$ are considered to be of order $O(1)$⁠, and the factor $ϖ$ explicitly accounts for the actual order. Inserting this into (11) and (12), we obtain scaled transport equations for (⁠$σ¯ij, q¯k$⁠) as

The Navier–Stokes and Fourier laws result from keeping only the first-order terms with factor $ϖ1$ on the rhs, that is, from ignoring the terms with factors $ϖ2$ and $ϖ3$ on the lhs. The terms with power $ϖ2$ on the lhs account for the second-order corrections to $(σij,qk)$⁠. These are complete in the sense that consideration of additional moments will not add any terms at this order.

The underlined term in the heat flux balance is the only third-order term appearing in Grad-13. However, additional third-order terms appear in the equations when higher moments are added to the theory.^18,19 With this, adding more moments to the list of variables results in higher order accuracy for stress and heat flux, and the predictions. We note that the conservation laws remain unchanged in their structure, and their accuracy is directly related to accuracy of stress and heat flux.

In the present contribution, we focus on the Grad-13 equations, which from the viewpoint of accuracy in terms of Knudsen number powers are limited to second order. Taking this at face value, the (underlined) third-order term should not play a role, so that it can be removed. On the other hand, second-order accuracy is maintained when third-order terms are added to the equations. This opens a strategy for reconfiguration of the Grad-13 equations by adding higher order terms such that the revised equations are accompanied by a properly formulated second law. Specifically, we shall aim to add terms that embed the revised Grad-13 equations into GENERIC.

In the framework of GENERIC, irreversible dynamics is generated by the entropy, while reversible dynamics is generated by the energy, as in Hamiltonian dynamics. We give a summary of GENERIC elements that suits our desire to reformulate the Grad-13 equations such that they fit into this general framework.

In GENERIC, when the variables are denoted by u_A, the balance laws in local form read

where $L$ is the antisymmetric Poisson matrix, and $M$ is the friction matrix, and we assume a strictly local form of the total energy and entropy functionals

where ϵ and η are the local energy and entropy densities appearing as generators in the evolution equation (23).

Poisson matrix and friction matrix must fulfill a number of requirements:⁶ 

The Poisson matrix must obey the degeneracy requirement

to ensure entropy conservation under reversible dynamics.

The antisymmetry of $LAB$ is defined through the Poisson bracket ${A,B}$⁠,

where $A(uA)$ and $B(uB)$ are arbitrary functionals of the variables u_A. Finally, $LAB$ has to obey the Jacobi identity

The friction matrix $MAB$ appears in the transport equations linearly. Each dissipative mechanism contributes a friction matrix $MABk$ so that the overall matrix $MAB$ is the sum

All friction matrices must obey the degeneracy requirement

which ensures energy conservation in irreversible processes.

The individual $MABk$ must obey Onsager–Casimir relations; that is, they are either symmetric or antisymmetric. The dissipative bracket for $MABk$ is defined as

For a symmetric friction matrix $MABk,sym$⁠, symmetry implies

where non-negative entropy production requires that

For an antisymmetric building block $MABk,anti$⁠, we have the non-dissipative bracket

with

According to the discussion in Sec. 3.2 of Ref. 6, an antisymmetric contribution $MABk,anti$ corresponds to the irreversible evolution without entropy production associated with the well-known Casimir symmetry in linear irreversible thermodynamics.²⁰ Note that a total friction matrix consisting of both symmetric and antisymmetric contributions associated with well-defined irreversible processes has no well-defined symmetry.

While their entropy, entropy flux, and entropy production can be approximated to second order, see Eqs. (16)–(18), the Grad-13 equations do not possess an accompanying second law with positive production. Accordingly, it is not possible to find GENERIC building blocks that just reproduce the full Grad-13 equations (10)–(12).

As already indicated above, instead we aim at finding GENERIC building blocks that reproduce Grad-13 exactly to second order, that is, those with powers $(ϖ,ϖ2)$ in Eqs. (21) and (22), while allowing additional higher order terms, with powers $ϖ3$ or higher. In the spirit of the Knudsen number scaling, these higher terms are expected to contribute very little to the solutions. In other words, we aim at correcting the Grad-13 equations by higher order contributions, such that the resulting equations possess a proper second law, while maintaining their second-order accuracy.

We construct the GENERIC version of the 13-moment equations not for the original variables $(ρ,vi,θ,σij,qi)$ of the Grad model, but instead consider the related set (5) with

Compared to the third moments q_i defined in Eq. (3), the new variables w_i have the major advantage that they possess vector transformation behavior under the general space deformations occurring in fluid dynamics, whereas q_i defines a Cartesian vector only. This behavior is crucial for finding a Poisson bracket with a degenerate entropy, which, as a scalar, can only depend on a scalar combination of $Θjk$ and w_j under general space transformations. Evaluation of (35) with the Grad-13 distribution gives

which to leading order reduces to

The energy density of the monatomic ideal gas is the sum of internal and kinetic energy (8), which for these variables assumes the form

with the gradient

The convex entropy density is postulated to be of the form¹¹ 

where $S¯(φ)$ as a dimensionless scalar under general space transformations is assumed to be a dimensionless function of the scalar quantity

with negative slope

The desired consistency with the Grad-13 equations is achieved by comparing the second-order expansion of the entropy (40) in σ_ij and q_i with Eq. (16) by assuming a linear dependence of $S¯(φ)$ on $φ$⁠. Note that by its definition (6), the temperature tensor is positive definite, and that by definition the flux variable w_i has unit $[θ]$⁠, while the unit of heat flux q_i is $[pθ]$⁠.

The entropy gradient with respect to the variables u_A is

To complete the model, we require a Poisson matrix $LAB$ and a set of friction matrices $MABk$ that reproduce the various terms in the Grad-13 equations. These must be inserted in to the balance laws (23) for the variables u_A. We note that the balance of total energy and entropy balance results as

In Secs. III C–III F, we present matrices $LAB$ and $MABk$ that were constructed to reproduce terms in the Grad-13 equations. We present these matrices one-by-one together with the resulting contributions to the transport equations.

A proper Poisson matrix that fulfills the degeneracy requirement (25) and the Jacobi identity (27) is given by (see Table 2.1 on p. 66 of Ref. 6 for lower convected tensors and vectors)

where indices in round brackets indicate symmetric tensors. It has been argued in Sec. 5(a) of Ref. 10 that lower convected derivatives appear naturally for moments of momentum. The contribution to the transport equations for the variables $uA=(ρ,Mi,Θij,wi)A$ is found as

where for compact notation, we have introduced the velocity gradient as

The contributions to energy and entropy balance are of divergence form

The degeneracy of the entropy with density (40) for the Poisson matrix (46) can be verified by a lengthy but straightforward calculation. The validity of the entropy degeneracy is the motivation for introducing variables with the proper general transformation behavior for scalars, vectors, and tensors.

The convection of moments obtained from the Boltzmann equation does not coincide with the lower convection behavior implemented by the Poisson matrix (46). In our previous contribution on the 10-moment equations,¹⁴ we hence proposed to introduce a “convective correction” term. Inspired by the thermodynamically consistent realization of Schowalter derivatives describing slip effects in rheology (see pp. 119f of Ref. 6), we proposed to implement the convective correction via an antisymmetric or Casimir-type contribution to the friction matrix, which is irreversible but non-dissipative (no entropy production; see  Appendix B for a brief discussion of Onsager–Casimir symmetry). Incorporation of the convective correction as a reversible term is impossible because of a violation of the Jacobi identity, which is widely accepted as a criterion for reversible behavior (in the sense of being under “mechanistic control”). As an extension of the friction matrix $Mcc$ of our previous contribution¹⁴ from 10 to 13 moments by introducing suitable additional contributions, we propose

with the vorticity tensor

The contribution to the transport equations for the u_A reads

This matrix fulfills the degeneracy condition (29) that it is has no contribution to the energy balance, and also, it has no contribution to entropy balance

The next friction contribution is constructed to give heat flux in the energy balance as well as the cross-coupling between heat flux and stress tensor in their respective moment equations, which we know from the Grad-13 equations. It is of the form

where a natural ansatz for the nonzero matrix elements reads

Here, the dimensionless parameters Z₁, Z_2, and c can be scalar functions of the variable $φ$⁠. In these matrix elements, the expression in square brackets has been chosen such that it provides the most general closure approximation for the third-moment tensor $∫CiCjCkfdc$ [cf. Eq. (15) of Ref. 11]. The straightforward implementation of this general closure approximation through $Mcross$ completes the discussion of 13-moment equations with entropy initiated in Refs. 11 and 13, which are both based on simplified closures. The corresponding general form of the entropy current (3) of Ref. 13 implied by $Mcross$ is given by

or, more explicitly,

with

We recognize that the contributions with Z₂ in (56) and (57) are of higher order; that is, they will not play a role in matching Grad-13 to second order in $Kn$⁠. For simplicity of the following line of arguments, we set $Z2=0$⁠, and we moreover assume that the parameters Z₁ and c are constants.

By construction, $MABcross$ is energy degenerate, $MABcross∂ϵ∂uB=0$⁠, while the contribution to the GENERIC balance laws becomes

The contribution to total energy balance gives the negative divergence of the heat flux as

with the heat flux

To identify the coefficients c, Z₁, and $dS¯(φ)/dφ$⁠, we study the leading order contributions and match these two terms in the Grad equations. To leading order the heat flux reduces to

and by comparison with (37), we identify

The leading contribution to stress balance is the trace-free part times $ρ$⁠,

and comparison with (11) shows that c = 1, which will be used from now on.

The leading contribution to the w_i-balance can be written as

where we used $∂Θkl−1∂xs=−Θkn−1Θml−1∂Θnm∂xs$⁠. With the temperature tensor given by (6), this can be further reduced to give the leading contribution to the heat flux equation (multiplication with $p/2$⁠)

Comparison with the Grad-13 equation (12) shows that setting $Z1=−4$ results in recovery of the first three terms, while the last term $(−σik∂θ∂xk)$ is extra, as will be discussed further below.

Thus, the matching at leading order yields

where the value for c implies symmetry of the third-moment tensor in its three indices and the values for $dS¯(φ)/dφ$ and $Q¯1$ coincide with previous findings.¹¹ We note that $S¯$ itself must vanish in equilibrium, $S¯(φ=0)=0$⁠, so that the proper equilibrium entropy results from (40). We will continue with the simplest choice for $S¯$ that also gives the proper equilibrium value and set

While $MABcross$ is antisymmetric in the Casimir sense, it contributes to entropy flux (but not to entropy production). We find

with the non-convective entropy flux

which is consistent with (59).

For the relaxation terms in the equations for stress and heat flux, we use

with $Aij kl=δik Θjl$⁠. The matrix $MABrelax$ gives energy degeneracy, $MABrelax∂ϵ∂uB=0$⁠, is Onsager symmetric, and its positive semidefiniteness of the friction matrix follows from the results of  Appendix A [see, in particular, the examples given in (A9) and (A10)].

The production terms are found as

where we already used (70). Indeed, the matrix was constructed such, that to second order the production terms reduce to the relaxation terms of the Grad equations

The non-linear second-order contribution $[23pμ(1θΘ⟨ij⟩)wj]$ to the w_i-balance is required for full matching, as will become clear further below.