Transition to the Fluid Dynamic Limit: Mathematical Models and Simulation Results

1. Introduction

Numerical simulations of flows in the transition from the gas kinetic (mesoscopic) to the fluid dynamic regime (macroscopic) represent a demanding task. On the mesoscopic side, the collision term operates on the 6-dimensional phase space and is very complex. On the macroscopic side, including a collision term in the fluid dynamic equations makes the problem become numerically stiff and hard to resolve in time.

In this situation simplyfying models are en vogue. The most popular ones are systems of BGK type which treat mesoscopic effects applying simple one- or more-parameter relaxation rules. Theoretical models try to find formulas from asymptotic analysis, taking into account the small mesoscopic perturbations as terms of an expansion with respect to some small parameter $ϵ$ (Hilbert, Chapman-Enskog expansion). From this we find the Euler and the Navier Stokes equations as zeroth and first order approximations. Improved models, e.g. for a more detailed refinement of the pressure tensor are expected from higher order terms.

In the course of the present paper we demonstrate that all these models fail in certain practically relevant situations. E.g., the $ϵ$-expansion ansatz fails in some condensation-evaporation problem for a gas mixture (see section 6.2), leading to some unphysical ghost effect. The BGK ansatz yields systematic errors even in simple heat layer problems (section 6.1). In the following we propose a new ansatz (trace solutions) which works on the tangent space to the set of all Maxwellians (rather than the classical ansatz “solution=Maxwellian plus orthogonal perturbation”, supplemented with some exponential relaxation rule). As a motivation we discuss the relaxation of a simple local perturbation problem (see section 5.2). The crucial feature of the new approach is that (close to the fluid dynamic limit) the local macroscopic moments are not reflected by the local Maxwellian but by a specific moment perturbation.

The idea to construct intermediate states between arbitrary density functions and Maxwellians is not new. E.g. Shakhov [12] proposed such a system which is at present used for numerical simulations [14]. A comparable intention lies behind the idea of the extended BGK model with an additional relaxation parameter to match the correct Prandtl number. However, a better understanding of the transfer to fluid dynamics is not a question of parameter matching, and the above approaches do not yield a theoretical basis and a safe mathematical ground. Proposing a modified structure of kinetic solutions, the present paper claims to provide a new attempt to better understand the passage from the rarefied gas description to the macroscopic limit.

We derive closure relations which take the form of nonlinear first order (rather than second order) differential terms and thus are completely different from the parabolic second order terms of the Navier-Stokes system. The results allow to interpret the differences of various kinetic models (here: nonlinear and linearized Boltzmann collision operators and the BGK relaxation operator) in the fluid dynamical limit. As particular results we point at a purely nonlinear effect of the Boltzmann collision operator which is not reflected in the Navier-Stokes approach (section 5.3), and we demonstrate that the BGK system produces systematic artificial effects which are non-local and which do not vanish in the fluid dynamical limit (section 6.1).

The derivation of the trace ansatz requires some tedious mathematical arguments, which are necessary for a sufficient understanding of the mathematical framework. The technical part has been worked out in [5] (with Lea Bold as coauthor who has contributed the numerical results of section 5.3 which are crucial for the justification of the trace ansatz). In the present paper we have tried to shift the focus from the technical part to the applications (in particular section 6.2), while keeping the paper self-contained.

Having derived the trace ansatz, we have performed numerical simulation studies for testing and comparing the various models mentioned above. As a numerical tool we have used discrete velocity models (in the version as described in [3]) on $5\times 5\times 5$- resp. $7\times 7\times 7$-velocity grids which are sufficient in our cases. These are applicable in identical form for all of the discussed models. (To implement collision terms on grids in an efficient manner with valid flow parameters, the application of computer algebra tools is recommended, as demonstrated in [6] for the 2D case and [7] for the 3D case.) These systems are well-investigated, with systematic errors in numerical simulations being well-understood and under control [11]. Grids of the above size are preferrable mainly for small to moderate Mach numbers, creep flows, etc. The numerical results are in full agreement with the theoretical results.

The paper is organized as follows. Section 2 gives a short review over the abovementioned kinetic models. Section 3 introduces the general (non-closed) moment system derived from gas kinetics and derives the Euler system. Section 4 reinterprets the steps of the first order Chapman-Enskog approach for the Navier-Stokes system. In this way the central point of procedure can be generalized and becomes applicable to the full nonlinear collision operator without recourse to a formal series expansion. A new concept of balanced states is introduced, which are elements of the tangent space to the manifold of the Maxwellians and which replace (resp. supplement) the first order terms of the density function as provided by the Chapman-Enskog procedure. In Section 5 we introduce the concept of traces of kinetic solutions. Traces are comparable to projections of kinetic solutions onto the tangent spaces. The differential structure of the underlying dynamics provides a powerful tool to describe distributions in the neighborhood of Maxwellians. In order to keep the underlying ideas as concise and understandable as possible, we restrict theory and numerical examples one-dimensional geometries. Section 6 applies the results to the heat layer problem and discusses an evaporation condensation problem for gas mixtures.

2. Kinetic Models

The most fundamental kinetic model equation for a density function $f=f(t,x,v)$ (x space-, $v={\left(v_{\xi }\right)}_{\xi =1}^d$ velocity-parameter) of a single-species gas in phase space $\mathrm{l}{\mathrm{R}}^d\times \mathrm{l}{\mathrm{R}}^d$ ($d=2,3$) is the classical Boltzmann equation

$$\begin{array}{c}\hfill {\displaystyle ({\partial }_t+v{\nabla }_x)f=\frac{1}{ϵ}·C_{J}f}\end{array}$$

(1)

The Boltzmann collision operator $C_{J}f=J[f,f]$ is a bilinear operator integrating over all possible conservative two-particle collisions. It satisfies the conservation laws

$$\begin{array}{ccc}\hfill {\displaystyle \langle \mathrm{l}1,C_{J}f\rangle =0}& & \left(\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right)\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill \langle v,C_{J}f\rangle =0& & \left(\mathrm{m}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{u}\mathrm{m}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right)\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill \langle v^2,C_{J}f\rangle =0& & \left(\mathrm{e}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{y}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{e}\mathrm{r}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\right)\hfill \end{array}$$

(4)

The Boltzmann collision operator is ruled by the H-Theorem stating that in a space homogeneous environment all densities converge for $t\nearrow \infty$ to equilibrium functions e (these are functions satisfying $C_{J}e=0$); all equilibrium functions are Maxwellians, i.e. functions of the form

$$\begin{array}{c}\hfill {\displaystyle e\left(v\right)=e[\rho ,\overline{v},T]\left(v\right)=\frac{\rho }{{\sqrt{2\pi T}}^d}·exp\left(-|v-\overline{v}{|}^2/2T\right)}\end{array}$$

(5)

uniquely determined by its macroscopic moments

$$\begin{array}{ccc}\hfill {\displaystyle \rho =\langle 1,e\rangle }& & \left(\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{t}\mathrm{y}\right)\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill \rho \overline{v}=\langle v,e\rangle & & \left(\mathrm{m}\mathrm{o}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{u}\mathrm{m}\right)\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill T=\frac{1}{d·\rho }·\langle {(v-\overline{v})}^2,e\rangle & & \left(\mathrm{t}\mathrm{e}\mathrm{m}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{e}\right)\hfill \end{array}$$

(8)

In the following we denote by

$$\begin{array}{ccc}\hfill {\displaystyle \mathcal{M}}& =& span(1,v_{\xi }(\xi =1,\dots ,d),v^2)\hfill \end{array}$$

the space of collision invariants, and by ${\mathcal{M}}^{\perp }$ its orthogonal complement. We make use of the

(2.1) Decomposition Lemma: 

Given a nonnegative density f and a Maxwellian e, there exists a unique decomposition

$$\begin{array}{ccc}\hfill {\displaystyle f}& =& eM·r+f_{\perp },r={\left(r_i\right)}_{i=0}^{d+1},f_{\perp }\in {\mathcal{M}}^{\perp }\hfill \end{array}$$

(9)

with M defined as

$$\begin{array}{ccc}\hfill {\displaystyle M}& =& (1,v_{\xi }(\xi =1,\dots ,d),v^2)\hfill \end{array}$$

(10)

The decomposition takes the form

$$\begin{array}{ccc}\hfill {\displaystyle f}& =& e+f_{\perp }\hfill \end{array}$$

(11)

iff e has the same macroscopic moments as f.

• ($eM·r$ is a short hand notation for an element in $e\mathcal{M}$, fully written as $e·(r_0+{\sum }_{\xi =1}^d{r}_{\xi }{v}_{\xi }+r_{d+1}{v}^2)$.)

Proof: 

r can be uniquely determined by calculating mass, momenta and temperature on both sides. □

• The set of all Maxwellians is denoted by $\mathcal{E}$. Given $e=e[\rho ,\overline{v},T]\in \mathcal{E}$, the tangent space in e is defined by

  $$\begin{array}{c}\hfill {\displaystyle \mathcal{T}\left[e\right]=\left\{e(1+\sum _{\xi =1}^d{r}_{\xi }{(v-\overline{v})}_{\xi }+r_{d+1}{(v-\overline{v})}^2),{\left(r_{\xi }\right)}_{\xi =1}^{d+1}\in \mathrm{l}{\mathrm{R}}^{d+1}\right\}\subset e\mathcal{M}}\end{array}$$

  (12)

The full tangent space $\mathcal{TE}$ is the union

$$\begin{array}{ccc}\hfill {\displaystyle \mathcal{TE}}& =& \bigcup _{e\in \mathcal{E}}\mathcal{T}\left[e\right]\hfill \end{array}$$

(13)

Two simplified alternatives to the nonlinear Boltzmann collision operator are the linearized collision operator $C_L$ and the BGL relaxation model$C_{BGK}$. Both are based on the decomposition (3). Given $f=e+f_{\perp }$, the Boltzmann collision operator is

$$\begin{array}{c}\hfill {\displaystyle J[f,f]=J[e,e]+J[e,f_{\perp }]+J[f_{\perp },e]+J[f_{\perp },f_{\perp }]=J[e,f_{\perp }]+J[f_{\perp },e]+J[f_{\perp },f_{\perp }]}\end{array}$$

Dropping the term which is quadratic in $f_{\perp }$, we end up with the linearization of $C_J$ around e,

$$\begin{array}{ccc}\hfill {\displaystyle C_{L}f=L_{e}f=L_e{f}_{\perp }}& =& J[e,f_{\perp }]+J[f_{\perp },e]\hfill \end{array}$$

(14)

$L_e$ transfers f exponentially in time to its equilibrium e. The well-known properties of $L_e$ are as follows.

(2.2) The linearized collision operator: 

Write $L_{e}f=:K\left(e^{-1}f\right)$. Then

• (a) $ker\left(K\right)={\left(im\left(K\right)\right)}^{\perp }=\mathcal{M}$
• (b) K is self-adjoint (e.g. in a weighted $L^2$ space) and negative semidefinite.
• (c) The restriction $K:{\mathcal{M}}^{\perp }\to {\mathcal{M}}^{\perp }$ has an inverse $K^{-1}$.
• (d) The linear system

  $$\begin{array}{c}\hfill {\displaystyle L_e\psi =\varphi ·e}\end{array}$$

  (15)

  is solvable iff $\varphi ·e\in {\mathcal{M}}^{\perp }$. In this case, the solution is

  $$\begin{array}{c}\hfill {\displaystyle \psi =L_e^{†}\varphi =e{K}^{-1}(\varphi ·e)}\end{array}$$

  (16)

  It satisfies

  $$\begin{array}{c}\hfill {\displaystyle \langle \varphi ,\psi \rangle =\langle \varphi ·e,K^{-1}(\varphi ·e)\rangle <0}\end{array}$$

  (17)
• For a more thorough treatment of $L_e$, see e.g. [9,10]. For corresponding results in the case of Discrete Velocity Models see [2,3].
• The simplest kind of exponential decay of f to its equilibrium e is given by the one-parameter BGK model

  $$\begin{array}{ccc}\hfill {\displaystyle C_{BGK}f}& =& \lambda ·(e-f)\hfill \end{array}$$

  (18)
• The solution of the linear system corresponding to (4) follows immediately from calculation.

(2.3) Linear system for BGK: 

Suppose $\varphi ·e\in {\mathcal{M}}^{\perp }$. The unique solution of

$$\begin{array}{c}\hfill {\displaystyle C_{BGK}\psi =\varphi ·e}\end{array}$$

(19)

is

$$\begin{array}{c}\hfill {\displaystyle \psi =e·\left(1-\frac{1}{\lambda }·\varphi \right),\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\langle \varphi ,\psi \rangle =-\frac{1}{\lambda }·\langle \varphi ,\varphi ·e\rangle <0}\end{array}$$

(20)

3. The moment hierarchy and the Euler system

Taking moments $\langle m,.\rangle$ ($m\in \mathcal{M}$) of the Boltzmann equation or one of its models, and taking into account that $\langle m,Cf\rangle =0$ for any of the above collision operators we end up with the following moment hierarchy.

• (3.4) General moment system:

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_t\rho +\nabla ·\left(\rho \overline{v}\right)}& =& 0\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\partial }_t\left(\rho {\overline{v}}_{\xi }\right)+\sum _{\eta }{\partial }_{\eta }(\rho {\overline{v}}_{\xi }{\overline{v}}_{\eta }+p_{\xi \eta })& =& 0\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\partial }_t\left(\rho ·|\overline{v}{|}^2+d·p\right)+\sum _{\xi }{\partial }_{\xi }\left[\rho {\overline{v}}_{\xi }{\left|\overline{v}\right|}^2+q_{\xi }+2{\left(P\overline{v}\right)}_{\xi }+d·p{\overline{v}}_{\xi }\right]& =& 0\hfill \end{array}$$

(23)

with the macroscopic moments density $\rho$, moment vector$\rho ·\overline{v}$, and (kinetic) energy E defined by

$$\begin{array}{ccc}\hfill {\displaystyle E}& =& {0.5\langle \left|v\right|}^2,f\rangle =0.5·(\rho {\overline{v}}^2+d\rho T)\hfill \end{array}$$

(24)

(9) to (11) is not a closed system, since it contains the closure moments

$$\begin{array}{ccc}\hfill {\displaystyle P={\left(p_{\xi ,\nu }\right)}_{1\le \xi ,\nu \le d}}& =& {\left(\langle {\pi }_{\xi ,\nu },f\rangle \right)}_{1\le \xi ,\nu \le d}\left(\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{o}\mathrm{r}\right)\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill q={\left(q_{\xi }\right)}_{\xi =1}^d& =& {\left(\langle {\sigma }_{\xi },f\rangle \right)}_{\xi =1}^d\left(\mathrm{h}\mathrm{e}\mathrm{a}\mathrm{t}\mathrm{f}\mathrm{l}\mathrm{u}\mathrm{x}\right)\hfill \end{array}$$

(26)

described by the centered moments

$$\begin{array}{ccc}\hfill {\displaystyle {\pi }_{\xi ,\nu }={\pi }_{\xi ,\nu }\left(\overline{v}\right)}& =& (v_{\xi }-{\overline{v}}_{\xi })(v_{\nu }-{\overline{v}}_{\nu })l\le \xi ,\nu \le d\hfill \end{array}$$

(27)

$$\begin{array}{ccc}\hfill {\sigma }_{\xi }={\sigma }_{\xi }\left(\overline{v}\right)& =& \frac{1}{2}(v_{\xi }-{\overline{v}}_{\xi }){|v-\overline{v}|}^2\hfill \end{array}$$

(28)

For given $e=e[\rho ,\overline{v},T]$, $\xi \ne \nu$, the function ${\pi }_{\xi ,\nu }e\in {\mathcal{M}}^{\perp }$. For the other centered moments, the projections onto ${\mathcal{M}}^{\perp }$ are given by

$$\begin{array}{ccc}\hfill {\displaystyle {\widehat{\pi }}_{\xi }}& =& {\pi }_{\xi \xi }-\frac{1}{d}\sum _{\nu }{\pi }_{\nu \nu }\hfill \end{array}$$

(29)

$$\begin{array}{ccc}\hfill {\widehat{\sigma }}_{\xi }={\widehat{\sigma }}_{\xi }(\overline{v},T)& =& {\sigma }_{\xi }-\frac{(d+2)·T}{2}·(v_{\xi }-{\overline{v}}_{\xi })\hfill \end{array}$$

(30)

Define pressure p by

$$\begin{array}{ccc}\hfill {\displaystyle p}& =& d^{-1}\langle |v-\overline{v}{|}^2,f\rangle =\rho T\hfill \end{array}$$

(31)

The easiest way to close system (9) to (11) is to assume that f is in its local equilibrium which gives

$$\begin{array}{ccc}\hfill {\displaystyle P}& =& p·I=\rho T·I,q=0\hfill \end{array}$$

(32)

with the resulting moment system

• (3.5) Euler equations:

  $$\begin{array}{ccc}\hfill {\displaystyle {\partial }_t\rho +\nabla ·\left(\rho \overline{v}\right)}& =& 0\hfill \end{array}$$

  (33)

  $$\begin{array}{ccc}\hfill {\partial }_t\left(\rho \overline{v}\right)+\nabla ·(\rho \overline{v}\otimes \overline{v})+\nabla \left(\rho T\right)& =& 0\hfill \end{array}$$

  (34)

  $$\begin{array}{ccc}\hfill {\partial }_t\left(\rho \right|\overline{v}{|}^2+d·\rho T)+\nabla ·(\rho \overline{v}|\overline{v}{|}^2)+(d+2)\nabla ·(\rho T·\overline{v})& =& 0\hfill \end{array}$$

  (35)

4. Navier-Stokes equations and balanced states

A prominent role in the refinement of the Euler system play the solutions ${\psi }_{\xi ,\nu }$ ($\xi \ne \nu$), ${\widehat{\psi }}_{\xi }$ and ${\psi }_{\xi }^{\sigma }$ in ${\mathcal{M}}^{\perp }$ of the linear equations

$$\begin{array}{ccc}\hfill {\displaystyle L{\psi }_{\xi ,\nu }}& =& {\pi }_{\xi ,\nu }e\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill L{\widehat{\psi }}_{\xi }& =& {\widehat{\pi }}_{\xi ,\nu }e\hfill \end{array}$$

(37)

$$\begin{array}{ccc}\hfill L{\psi }_{\xi }^{\sigma }& =& {\widehat{\sigma }}_{\xi }e\hfill \end{array}$$

(38)

(For the solvability of the equations see (2.2)(c).)

• A convenient way to derive the Navier-Stokes system is the Chapman Enskog expansion of first order. Plugging the ansatz

  $$\begin{array}{c}\hfill {\displaystyle f=f_0+ϵf_1}\end{array}$$

  into the Boltzmann equation and setting equal the terms of equal power of $ϵ$. For ${ϵ}^{-1}$ we find

  $$\begin{array}{c}\hfill {\displaystyle J[f_0,f_0]=0}\end{array}$$

  and thus $f_0=e\in \mathcal{E}$. For ${ϵ}^1$ we find

  $$\begin{array}{c}\hfill {\displaystyle ({\partial }_t+v·{\nabla }_x)e=L\left[e\right]f_1}\end{array}$$

  (39)
• For its solution we perform the decomposition (2) on the whole system. The part acting on ${\mathcal{M}}^{\perp }$ reads

  $$\begin{array}{c}\hfill {\displaystyle {\left(({\partial }_t+v·{\nabla }_x)e\right)}_{\perp }=L\left[e\right]f_1}\end{array}$$

  (40)

  with

  $$\begin{array}{ccc}\hfill {\displaystyle {\left(({\partial }_t+v·{\nabla }_x)e\right)}_{\perp }}& =& {(v·{\nabla }_{x}e)}_{\perp }\hfill \\ & =& \rho ·\sum _{\xi }\sum _{\nu }{\partial }_{\xi }{\overline{v}}_{\nu }{\left(v_{\xi }{\partial }_{{\overline{v}}_{\nu }}e[1,\overline{v},T]\right)}_{\perp }+\rho ·\sum _{\xi }{\left(v_{\xi }{\partial }_{T}e[1,\overline{v},T]·{\partial }_{\xi }T\right)}_{\perp }\hfill \\ & =& \frac{\rho }{T}·\sum _{\xi }\sum _{\nu }{\partial }_{\xi }{\overline{v}}_{\nu }·{\pi }_{\xi ,\nu }e+\frac{\rho }{T^2}·\sum _{\xi }{\partial }_{\xi }T·{\widehat{\sigma }}_{\xi }e\hfill \end{array}$$

  (41)
• The solution of (19) is

  $$\begin{array}{ccc}\hfill {\displaystyle f_1^{NS}}& =& \frac{1}{T}·\sum _{\xi }\left[\sum _{\nu \ne \xi }{\partial }_{\xi }{\overline{v}}_{\nu }·{\psi }_{\xi \nu }+{\partial }_{\xi }{\overline{v}}_{\xi }·{\widehat{\psi }}_{\xi }\right]·e+\frac{ϵ}{2T^2}·\sum _{\xi }{\partial }_{\xi }T·{\psi }_{\xi }^{\sigma }·e\hfill \end{array}$$

  (42)
• The Navier-Stokes correction to the Euler equations we now get by calculating the moment system (3.4) of the function $f=e+ϵf^{\left(1\right)}$. Taking into account the inequalities (6), we find the correction terms to the closure moments of the Euler system,

  $$\begin{array}{ccc}\hfill {\displaystyle p_{\xi ,\xi }^{NS}}& =& -ϵ·{\mu }^{\left(1\right)}·\frac{d}{d-1}·\frac{1}{T}·\left({\partial }_{\xi }{\overline{v}}_{\xi }-\frac{1}{d}\nabla ·\overline{v}\right)\hfill \end{array}$$

  (43)

  $$\begin{array}{ccc}\hfill p_{\xi ,\nu }^{NS}& =& -ϵ{\mu }^{\left(2\right)}·\frac{1}{T}·\left({\partial }_{\xi }{\overline{v}}_{\nu }+{\partial }_{\nu }{\overline{v}}_{\xi }\right)(\xi \ne \nu )\hfill \end{array}$$

  (44)

  $$\begin{array}{ccc}\hfill q_{\xi }^{NS}& =& -ϵ·\lambda ·\frac{1}{2T^2}·{\partial }_{\xi }T\hfill \end{array}$$

  (45)

  with the positive viscosity resp. heat coefficients

  $$\begin{array}{ccc}\hfill {\displaystyle {\mu }^{\left(1\right)}}& =& -\langle {\widehat{\pi }}_{\xi },{\widehat{\psi }}_{\xi }\rangle \hfill \end{array}$$

  (46)

  $$\begin{array}{ccc}\hfill {\mu }^{\left(2\right)}& =& -\langle {\pi }_{\xi ,\nu },{\psi }_{\xi ,\nu }\rangle \hfill \end{array}$$

  (47)

  $$\begin{array}{ccc}\hfill \lambda & =& -\langle {\widehat{\sigma }}_{\xi },{\psi }_{\xi }^{\sigma }\rangle \hfill \end{array}$$

  (48)
• All the above results are classical and need no further explanation. Here, we find a way to reinterpret the Chapman Enskog results and to generalize them, without making use of the series expansion.

(4.6) Lemma: Let $e=e(t_0,x_0)$ be a fixed Maxwellian and $s=s(t_0,x_0)=v·{\nabla }_{x}e$ the corresponding flow term. The expansion of first order $f^{\left(0\right)}+ϵ·f^{\left(1\right)}$ of Chapman-Enskog is the unique asymptotic solution ($\tau \nearrow \infty$) of the time-homogeneous initial value problem (IVP)

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_{\tau }f+ϵ·s_{\perp }}& =& L_{e}f,f\left(t_0\right)=e\hfill \end{array}$$

Proof: Decompose f uniquely into the form

$$\begin{array}{ccc}\hfill {\displaystyle f\left(\tau \right)}& =& f_{\parallel }\left(\tau \right)+f_{\perp }\left(\tau \right),f_{\parallel }\left(\tau \right)\in e·\mathcal{M},f_{\perp }\in {\mathcal{M}}^{\perp }\hfill \end{array}$$

Since $f_{\parallel }\left(\tau \right)\in \mathrm{k}\mathrm{e}\mathrm{r}\left(L_e\right)$ and $s_{\perp }\in {\mathcal{M}}^{\perp }$, we find

$$\begin{array}{ccc}\hfill {\displaystyle f_{\parallel }\left(\tau \right)}& =& e\hfill \end{array}$$

and $f_{\perp }\left(\tau \right)$ solution of the IVP in ${\mathcal{M}}^{\perp }$,

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_{\tau }{f}_{\perp }+ϵ·s_{\perp }}& =& L_e{f}_{\perp },f_{\perp }\left(t_0\right)=0\hfill \end{array}$$

Since $L_e:{\mathcal{M}}^{\perp }\to {\mathcal{M}}^{\perp }$ is invertible and negative definite, the asymptotics

$$\begin{array}{c}\hfill {\displaystyle f_{\perp }^{\infty }=\underset{\tau \nearrow \infty }{lim}{f}_{\perp }\left(\tau \right)}\end{array}$$

exists and satisfies

$$\begin{array}{ccc}\hfill {\displaystyle ϵ·s_{\perp }}& =& L_e{f}_{\perp }^{\infty }\hfill \end{array}$$

Thus $f_{\perp }^{\infty }=ϵ·f^{\left(1\right)}.\square$

• Thus the Chapman Enskog procedure as sketched above may be seen as a special case of the following three-step procedure for the calculation of closure relations.

(4.7) Scheme for closure relations:(1) Let $(t_0,x_0)$ and $f=f(t_0,x_0)$ be given and fixed. Calculate a reduced discription $\mathbf{f}$ of f containing all information concerning the relevant moments.

• (2) For an appropriate approximation C of the Boltzmann collision operator solve the IVP

  $$\begin{array}{ccc}\hfill {\displaystyle {\partial }_{\tau }g+ϵ·{\left(v·{\nabla }_x\mathbf{f}\right)}_{\perp }}& =& Cg,g\left(0\right)=\mathbf{f}\hfill \end{array}$$

  (49)

(iii) Determine $g^{\infty }={lim}_{\tau \nearrow \infty }g\left(\tau \right)$ and calculate from this the closure relations for f.

• (In the Chapman-Enskog case, $\mathbf{f}=e$ and $C=L_e$.) We can interpret (22) as a local solution to the kinetic equation: Introduce a microscopic time scale $\tau =(t-t_0)/ϵ$ and solve the kinetic system. In lowest order in $ϵ$, we freeze the source term ${\left(v·{\nabla }_x\mathbf{f}\right)}_{\perp }$ at $\tau =0$ and solve the system for $\tau \nearrow \infty$. In the following we call the distributions $g^{\infty }$ balanced states, since they balance the trend to equilibrium given by the collision operator with the perturbing effects produced by the streaming term.

5. Traces in tangent space

5.1. Traces

The definition of the linearized collision operator and the results of the above Chapman Enskog procedure above as well as the BGK model suggest a decomposition of f of the form (3) rather than (2). In this case, the central element e contains all macroscopic moments, but no contribution to the closure moments. In view of the balance equations described above with source term ${\left(v_x{\partial }_{x}e\right)}_{\perp }$, we propose an alternative decomposition. For simplicity we restrict in the following to spatially one-dimensional problems (with space variable x) and with bulk velocities $\overline{v}={({\overline{v}}_x,0)}^T$ resp. $\overline{v}={({\overline{v}}_x,0,0)}^T$.

• Given a triple $[\overline{\rho },{\overline{v}}_x,\overline{T}]$ with corresponding Maxwellian $\overline{e}=e[\overline{\rho },{\overline{v}}_x,\overline{T}]$, the tangent space$\mathcal{T}\left[\overline{e}\right]$ of $\mathcal{E}$ in $\overline{e}$ is defined as the set of all elements in

$$\begin{array}{c}\hfill {\displaystyle \mathcal{EM}=\bigcup _{e\in \mathcal{E}}e·\mathcal{M}}\end{array}$$

with the same macroscopic moments $[\overline{\rho },{\overline{v}}_x,\overline{T}]$ as $\overline{e}$. It is a two-dimensional manifold in $\mathcal{EM}$ and takes the form

$$\begin{array}{ccc}\hfill {\displaystyle \mathcal{T}\left[\overline{e}\right]}& =& \{\widehat{e}+{\gamma }^{\left(1\right)}·{\widehat{e}}^{\left(1\right)}+{\gamma }^{\left(2\right)}·{\widehat{e}}^{\left(2\right)},\alpha ,\beta \in \mathrm{l}\mathrm{R}\}\hfill \end{array}$$

(50)

$$\begin{array}{ccc}\hfill \widehat{e}& =& e[\overline{\rho },{\widehat{v}}_x,\widehat{T}]\hfill \end{array}$$

(51)

$$\begin{array}{ccc}\hfill {\widehat{e}}^{\left(1\right)}& =& \frac{1}{\widehat{T}}·(v_x-{\widehat{v}}_x)·\widehat{e}\hfill \end{array}$$

(52)

$$\begin{array}{ccc}\hfill {\widehat{e}}^{\left(2\right)}& =& \frac{1}{2{\widehat{T}}^2}·\left(|v-\widehat{v}{|}^2-d·\widehat{T}\right)·\widehat{e}\hfill \end{array}$$

(53)

with $({\gamma }^{\left(1\right)},{\gamma }^{\left(2\right)})$ and $({\widehat{v}}_x,\widehat{T})$ connected to $({\overline{v}}_x,\overline{T})$ via

$$\begin{array}{ccc}\hfill {\displaystyle {\widehat{v}}_x}& =& {\overline{v}}_x-{\gamma }^{\left(1\right)}\hfill \end{array}$$

(54)

$$\begin{array}{ccc}\hfill \widehat{T}& =& \overline{T}+\frac{1}{d}·\left({\left({\gamma }^{\left(1\right)}\right)}^2-{\gamma }^{\left(2\right)}\right)\hfill \end{array}$$

(55)

Elements of $\mathcal{T}\left[\overline{e}\right]$ are addressed as

$$\begin{array}{ccc}\hfill {\displaystyle \mathbf{f}[{\gamma }^{\left(1\right)},{\gamma }^{\left(2\right)}]}& =& \widehat{e}+{\gamma }^{\left(1\right)}·{\widehat{e}}^{\left(1\right)}+{\gamma }^{\left(2\right)}·{\widehat{e}}^{\left(2\right)}\hfill \end{array}$$

The union of tangent spaces is denoted as

$$\begin{array}{ccc}\hfill {\displaystyle \mathcal{TE}}& =& \bigcup _{e\in \mathcal{M}}\mathcal{T}\left[e\right]\hfill \end{array}$$

(56)

In the following, the parameter ${\gamma }^{\left(2\right)}$ is not relevant in the case of the full Boltzmann collision operator. Thus we set ${\gamma }^{\left(2\right)}=0$, ${\gamma }^{\left(1\right)}:=\gamma$ and define as the trace of $\overline{e}$ in $\mathcal{TE}$ the mapping

$$\begin{array}{c}\hfill {\displaystyle \gamma \to \mathbf{f}\left(\gamma \right)=\widehat{e}+\gamma ·{\widehat{e}}^{\left(1\right)}}\end{array}$$

(57)

$\mathbf{f}\left(\gamma \right)$ has the same macroscopic moments as $\overline{e}$, however it contributes to the closure moments as follows.

• (5.8) Closure relations for trace: The closure moments of the trace elements (in relation to the macroscopic moments ${\overline{v}}_x$, $\overline{T}$ of $\mathbf{f}\left(\gamma \right)$) are

$$\begin{array}{ccc}\hfill {\displaystyle \langle {\widehat{\pi }}_x\left(\overline{v}\right),\mathbf{f}\left(\gamma \right)\rangle }& =& -\frac{d-1}{d}·\rho {\gamma }^2\hfill \end{array}$$

(58)

$$\begin{array}{ccc}\hfill \langle {\widehat{\sigma }}_x(\overline{v},\overline{T}),\mathbf{f}\left(\gamma \right)\rangle & =& -\rho {\gamma }^3\hfill \end{array}$$

(59)

Proof: These formulas follow from

$$\begin{array}{ccc}\hfill {\displaystyle \langle {\widehat{\pi }}_x\left(\overline{v}\right),\mathbf{f}\left(\gamma \right)\rangle }& =& \frac{d-1}{d}·\rho ·\left({({\widehat{v}}_x-{\overline{v}}_x)}^2+2\gamma ·({\widehat{v}}_x-{\overline{v}}_x)\right)=-\frac{d-1}{d}·\rho ·{({\widehat{v}}_x-{\overline{v}}_x)}^2\hfill \end{array}$$

and

$$\begin{array}{ccc}\hfill {\displaystyle \langle {\widehat{\sigma }}_x(\overline{v},\overline{T}),\mathbf{f}\left(\gamma \right)\rangle }& =& -\rho ·\frac{{\widehat{v}}_x-{\overline{v}}_x}{2}·\left(\frac{d-2}{d}{({\widehat{v}}_x-{\overline{v}}_x)}^2+(d+2)·(\widehat{T}-\overline{T})\right)\square \hfill \end{array}$$

Given $\mathbf{f}\left(\gamma \right)$, we denote as $e\left(\gamma \right)$, ${\widehat{\pi }}_x\left(\gamma \right)$ and ${\widehat{\sigma }}_x\left(\gamma \right)$ the Maxwellian and the closure moments corresponding to the parameters ${\widehat{v}}_x=-\gamma$ and $\widehat{T}=\overline{T}+{\gamma }^2/d$. The normed vectors ${\mathbf{n}}_{\pi }$ and ${\mathbf{n}}_{\sigma }$ are for given $\gamma$ to be interpreted according to the formulas

$$\begin{array}{ccc}\hfill {\displaystyle {\widehat{\pi }}_x\left(\gamma \right)·\widehat{e}\left(\gamma \right)}& =& \parallel {\widehat{\pi }}_x\left(\gamma \right)\widehat{e}\left(\gamma \right)\parallel ·{\mathbf{n}}_{\pi }\hfill \end{array}$$

(60)

$$\begin{array}{ccc}\hfill {\widehat{\sigma }}_x\left(\gamma \right)·\widehat{e}\left(\gamma \right)& =& \parallel {\widehat{\sigma }}_x\left(\gamma \right)\widehat{e}\left(\gamma \right)\parallel ·{\mathbf{n}}_{\sigma }\hfill \end{array}$$

(61)

Furthermore we define as the local tangent plane the affine linear space

$$\begin{array}{c}\hfill {\displaystyle \mathbf{f}\left(\gamma \right)+\mathrm{s}\mathrm{p}\mathrm{a}\mathrm{n}({\mathbf{n}}_{\pi },{\mathbf{n}}_{\sigma })}\end{array}$$

(62)

Given $\mathbf{z}\in span({\mathbf{n}}_{\pi },{\mathbf{n}}_{\sigma })$, the angle ${\mathrm{\Phi }}_{\mathbf{z}}$ defined by

$$\begin{array}{c}\hfill {\displaystyle \mathbf{z}=cos\left({\mathrm{\Phi }}_{\mathbf{z}}\right)·{\mathbf{n}}_{\pi }+sin\left({\mathrm{\Phi }}_{\mathbf{z}}\right)·{\mathbf{n}}_{\sigma }}\end{array}$$

(63)

describes the position of $\mathbf{f}\left(\gamma \right)+\mathbf{z}$ in the local tangent plane.

We collect some results concerning the local structure of $\mathbf{f}\left(\gamma \right)$ which can be proved by straightforward calculations and by Taylor expansion.

• (5.9) Differential structure of trace : (a) $\mathbf{f}$ is differentiable with derivative

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_{\gamma }\mathbf{f}=\parallel {\partial }_{\gamma }\mathbf{f}\parallel ·{\mathbf{n}}_{\partial \mathbf{f}}}& =& -\frac{\gamma }{{\widehat{T}}^2}·\parallel {\widehat{\pi }}_x(-{\widehat{v}}_x)·\widehat{e}\parallel ·{\mathbf{n}}_{\pi }+\frac{2{\gamma }^2}{d·{\widehat{T}}^3}·\parallel \widehat{\sigma }(-{\widehat{v}}_x,\widehat{T})·\widehat{e}\parallel ·{\mathbf{n}}_{\sigma }\hfill \end{array}$$

(64)

(b) The leading order approximation of $\mathbf{f}\left(\gamma \right)$ for $\left|\gamma \right|\ll 1$ in the local tangent plane $(\gamma =0)$ is

$$\begin{array}{ccc}\hfill {\displaystyle \mathbf{f}\left[\gamma \right]}& =& \overline{e}·\left(1-\frac{{\gamma }^2}{2·{\overline{T}}^2}·{\widehat{\pi }}_x\left(0\right)+\frac{2{\gamma }^3}{3d·{\overline{T}}^3}·{\widehat{\sigma }}_x(0,\overline{T})\right)\hfill \end{array}$$

(65)

(c) For $\left|\gamma \right|\ll 1$,

$$\begin{array}{ccc}\hfill {\displaystyle \mathrm{t}\mathrm{a}\mathrm{n}\left({\varphi }_{\partial \mathbf{f}}\left(\gamma \right)\right)}& =& -\frac{2\gamma }{d·\widehat{T}}·\frac{\parallel {\widehat{\sigma }}_x\widehat{e}\parallel }{\parallel {\widehat{\pi }}_x\widehat{e}\parallel }=:-C_{\partial \mathbf{f}}·\gamma +\mathcal{O}\left({\gamma }^2\right)\hfill \end{array}$$

(66)

$$\begin{array}{ccc}\hfill {\mathrm{\Phi }}_{\partial \mathbf{f}}& =& \left\{\begin{array}{ccc}-C_{\partial \mathbf{f}}·\gamma +\mathcal{O}\left({\gamma }^2\right)& \mathrm{f}\mathrm{o}\mathrm{r}& \gamma <0\\ \pi -C_{\partial \mathbf{f}}·\gamma +\mathcal{O}\left({\gamma }^2\right)& \mathrm{f}\mathrm{o}\mathrm{r}& \gamma >0\end{array}\right.\hfill \end{array}$$

(67)

Thus at $\gamma =0$ the curve $\mathbf{f}\left(\gamma \right)$ is not differentiable and thus not a regular parametrization.

• Proof: (a) follows from straightforward calculations (using the definition (30) of $\mathbf{f}\left(\gamma \right)$ and formulas (27), (28)) which yield

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_{\gamma }\widehat{e}}& =& \left(-\frac{\gamma }{\widehat{T}}-\frac{v_x-{\widehat{v}}_x}{\widehat{T}}+\frac{\gamma }{d·{\widehat{T}}^2}·{(v-\widehat{v})}^2\right)·\widehat{e}\hfill \end{array}$$

(68)

and

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_{\gamma }\left(\gamma ·\frac{v_x-{\widehat{v}}_x}{\widehat{T}}\right)}& =& \frac{\gamma }{\widehat{T}}+\left(\frac{1}{\widehat{T}}-\frac{2{\gamma }^2}{d·{\widehat{T}}^2}\right)·(v_x-{\widehat{v}}_x)\hfill \end{array}$$

(69)

(b) is obtained by integrating (33) and using $\mathbf{f}\left(0\right)=\overline{e}$ and ${\partial }_{\gamma }\mathbf{f}\left(0\right)=0.\square$

• For the derivation of closure moments for the nonlinear collision operator it is useful to follow the trace of $C_J$ along $\mathbf{f}\left(\gamma \right)$ which is defined as the projection of $C_J$ given by

  $$\begin{array}{ccc}\hfill {\displaystyle {\mathbf{C}}_J\mathbf{f}\left(\gamma \right)}& =& \langle C_J\mathbf{f}\left(\gamma \right),{\mathbf{n}}_{\partial \mathbf{f}}\rangle ·{\mathbf{n}}_{\partial \mathbf{f}}\hfill \\ & =& \mathrm{c}\mathrm{o}\mathrm{s}\left({\varphi }_{\partial \mathbf{f}}\right)\langle C_J\mathbf{f}\left(\gamma \right),{\mathbf{n}}_{\pi }\rangle ·{\mathbf{n}}_{\partial \mathbf{f}}+sin\left({\varphi }_{\partial \mathbf{f}}\right)·\langle C_J\mathbf{f}\left(\gamma \right),{\mathbf{n}}_{\sigma }\rangle ·{\mathbf{n}}_{\partial \mathbf{f}}\hfill \end{array}$$

  (70)

(5.10) Trace of nonlinear operator: (a) With $a_{\pi ,\pi }^J=-\langle J[{\widehat{e}}^{\left(1\right)},{\widehat{e}}^{\left(1\right)}],{\mathbf{n}}_{\pi }\rangle$,

$$\begin{array}{c}\hfill {\displaystyle {\mathbf{C}}_J\mathbf{f}\left(\gamma \right)=sign\left(\gamma \right)·cos\left({\varphi }_{\partial \mathbf{f}}\right)·{\gamma }^2·a_{\pi ,\pi }^J·{\mathbf{n}}_{\partial \mathbf{f}}}\end{array}$$

(71)

(b) Under the approximations of (5.9)(b),(c) for small $\gamma$,

$$\begin{array}{ccc}\hfill {\displaystyle {\mathbf{C}}_J\mathbf{f}\left(\gamma \right)}& =& -{\gamma }^2{a}_{\pi ,\pi }^J·{\mathbf{n}}_{\pi }+\mathcal{O}\left({\gamma }^3\right)\hfill \end{array}$$

and the solution $\mathbf{f}\left(\gamma \right(t\left)\right)$ of the homogeneous trace Boltzmann equation

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_t\mathbf{f}\left(\gamma \left(t\right)\right)}& =& {\mathbf{C}}_J\mathbf{f}\left(\gamma \left(t\right)\right)\hfill \end{array}$$

(72)

is in lowest order given by

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_t\gamma }& =& -\lambda \gamma ,\lambda =a_{\pi ,\pi }^J·{\overline{T}}^2·{\parallel {\widehat{\pi }}_x\overline{e}\parallel }^{-1}\hfill \end{array}$$

(73)

Proof: (a) From the definition of $\mathbf{f}\left(\gamma \right)$ (30) and the fact that ${\widehat{e}}^{\left(1\right)}\in ker\left(L_{\widehat{e}}\right)$ follows

$$\begin{array}{c}\hfill {\displaystyle C_J\mathbf{f}\left(\gamma \right)={\gamma }^2·J[{\widehat{e}}^{\left(1\right)},{\widehat{e}}^{\left(1\right)}]}\end{array}$$

In a coordinate system centered around $\widehat{v}$, $J[{\widehat{e}}^{\left(1\right)},{\widehat{e}}^{\left(1\right)}]$ is an even function with respect to $\widehat{v}$. Thus

$$\begin{array}{c}\hfill {\displaystyle \langle {\widehat{\sigma }}_x(\widehat{v},\widehat{T}),J[{\widehat{e}}^{\left(1\right)},{\widehat{e}}^{\left(1\right)}]\rangle =\langle {\mathbf{n}}_{\sigma },J[{\widehat{e}}^{\left(1\right)},{\widehat{e}}^{\left(1\right)}]\rangle =0}\end{array}$$

(b) follows from inserting the approximations (5.9)(b),(c).$\square$

• In our context, the most significant difference between the nonlinear collision operator on the one side, and the linearized and the BGK operator (shortly termed as “relaxation operators”) on the other side is that

$$\begin{array}{c}\hfill {\displaystyle C_J\mathbf{f}\left(\gamma \right)=\frac{{\gamma }^2}{{\widehat{T}}^2}·J[(v_x-{\widehat{v}}_x)·\widehat{e},(v_x-{\widehat{v}}_x)·\widehat{e}]}\end{array}$$

(74)

is centered around $\widehat{e}$ (with the asymptotics $\overline{e}$ “unvisible”), while the Maxwellian $\overline{e}$ appears in the definition of the relaxation operators and thus marks it as their asymptotic end point. For the same reason we find closed formulas for the solutions of the homogenous equations in the relaxation case,

$$\begin{array}{ccc}\hfill {\displaystyle f\left(t\right)}& =& \overline{e}+exp\left(t{L}_{\overline{e}}\right)·(f_0-\overline{e})\left(linearizedequation\right)\hfill \end{array}$$

(75)

$$\begin{array}{ccc}\hfill f\left(t\right)& =& \overline{e}+exp(-\lambda t)·(f_0-\overline{e})\left(BGKequation\right)\hfill \end{array}$$

(76)

while we only can formulate the solution of the homogeneous nonlinear equation in differential form.

5.2. Relaxation problems

Consider the initial value problem

$$\begin{array}{c}\hfill {\displaystyle {\partial }_{t}f=Cf,f\left(0\right)=f_0=\widehat{e}+ϵ·\varphi }\end{array}$$

with C one of the collision operators discussed above. Depending on $\varphi$, there exists a unique Maxwellian $\overline{e}$ with the same macroscopic moments as $f_0$, which is the asymptotic state of the solution for $t\nearrow \infty$. In the relaxation case we have exponential decay to $\overline{e}$, (33), (34). The dynamics of the nonlinear solution is more detailed, depending on the form of $\varphi$.

• Suppose

$$\begin{array}{ccc}\hfill {\displaystyle \varphi }& =& \widehat{e}·m+{\varphi }_{\perp },m\in \mathcal{M},{\varphi }_{\perp }\in {\mathcal{M}}^{\perp }\hfill \end{array}$$

(77)

Then

$$\begin{array}{ccc}\hfill {\displaystyle J[f_0,f_0]}& =& J[\widehat{e},\widehat{e}]+ϵ·\left(J[\widehat{e},m\widehat{e}]+J[m\widehat{e},\widehat{e}]\right)+ϵ·\left(J[\widehat{e},{\varphi }_{\perp }]+J[{\varphi }_{\perp },\widehat{e}]\right)+\mathcal{O}\left({ϵ}^2\right)\hfill \end{array}$$

(78)

with

$$\begin{array}{ccc}\hfill {\displaystyle J[\widehat{e},\widehat{e}]}& =& 0\hfill \end{array}$$

(79)

$$\begin{array}{ccc}\hfill J[\widehat{e},m\widehat{e}]+J[m\widehat{e},\widehat{e}]& =& L_{\widehat{e}}m\widehat{e}=0\hfill \end{array}$$

(80)

$$\begin{array}{ccc}\hfill J[\widehat{e},{\varphi }_{\perp }]+J[{\varphi }_{\perp },\widehat{e}]& =& L_{\widehat{e}}{\varphi }_{\perp }\hfill \end{array}$$

(81)

For the short time behaviour we may neglect the $\mathcal{O}\left({ϵ}^2\right)$-terms. Thus the initialization phase (initial layer) is dominated by the exponential relaxation of ${\varphi }_{\perp }$ to the “wrong” equilibrium $\widehat{e}$. What follows is a dynamics in $\mathcal{TE}$ which runs on a slower time scale and in which the initial function $\widehat{e}(1+ϵm)$ approaches its new equilibrium $\overline{e}$. The difference between fast exponential decay for the linearized model and of the slower approach to equilibrium for the nonlinear system can be easily verified e.g. running simulations with Discrete Velocity Models (DVM). Figure 1 a, b display the relaxation rates for randomly perturbed initial conditions in the cases 1, $\varphi \in {\mathcal{M}}^{\perp }$, and case 2, $\varphi \in \widehat{e}·\mathcal{M}$. Case 1 exhibits exponential decay, with rates within a narrow band. This situation can be reasonably well imitated with a BGK model (straight line).

Case 2 is more complex and we may recognize three phases. Phase I represents the relaxation of the perturbation $\varphi$ to $(1+ϵ\mathbf{m})\widehat{e}$, phase II (“convex phase”) the transition from $(1+ϵ\mathbf{m})\widehat{e}$ to $\overline{e}$, and phase III the exponential relaxation to $\overline{e}$. In the special case $\widehat{e}·m={\widehat{e}}^{\left(1\right)}$, phase II marks the transition of the trace element $\mathbf{f}\left(ϵ\right)$ to $\mathbf{f}\left(0\right)=\overline{e}$. The slope of phase III is given by the parameter $\lambda$ of (73). As a special feature we notice that the life time of the moment perturbation $ϵ\mathbf{m}\widehat{e}$ is much larger than for a model with relaxation rates based on the mean free path. The BGK model is not capable to cover this situation and we have to decide whether to match to the initial Maxwellian or the final Maxwellian, depending on whether we want to simulate the short-time or the long-time behaviour.

5.3. Relaxation with source, balanced states

In view of the equations to be solved for the closure relations let’s have a look on solutions of the equation with source,

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_{t}f+ϵs}& =& Cf,f\left(0\right)=\overline{e},s={\widehat{\sigma }}_x·\overline{e}\hfill \end{array}$$

(82)

and in particular their asymptotic states $f^{(\infty )}$ for $t\nearrow \infty$ (balanced states). For the relaxation systems we easily find

$$\begin{array}{ccc}\hfill {\displaystyle f_{BGK}\left(t\right)}& =& \overline{e}-ϵ{\lambda }^{-1}·s·(1-exp(-\lambda t))\nearrow f_{BGK}^{\infty }=\overline{e}-ϵ{\lambda }^{-1}·s\hfill \\ \hfill f_L\left(t\right)& =& \overline{e}+exp\left(t{L}_{\overline{e}}\right)f\left(0\right)+ϵ·(1-exp\left(t{L}_{\overline{e}}\right)){\psi }_x^{\sigma }\hfill \\ & & \nearrow f_L^{\infty }=\overline{e}-ϵ{\psi }_x^{\sigma }({\psi }_x^{\sigma }\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{o}\mathrm{f}{L}_{\overline{e}}{\psi }_x^{\sigma }={\widehat{\sigma }}_x·\overline{e},\mathrm{s}\mathrm{e}\mathrm{e}\mathrm{S}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}...)\hfill \end{array}$$

Both solutions satisfy the orthogonality relation

$$\begin{array}{ccc}\hfill {\displaystyle \langle f^{(\infty )},{\widehat{\pi }}_x\rangle }& =& 0\hfill \end{array}$$

This is no surprise, since the source is even with respect to $v_x$, and L maps even into even functions, while ${\widehat{\sigma }}_x$ is odd. So it is remarkable to observe that orthogonality is not given for the nonlinear solution, as can be seen in Figure 2. (This observation was made in [8].) To understand this phenomenon, it helps to look at trace solutions. We easily see that $\gamma =0$ is a balanced solution to the trace version of (36). However, it is not stable in the sense of dynamical systems, and there is another solution. Necessary and sufficient for balancing solutions is the orthogonality relation

$$\begin{array}{c}\hfill {\displaystyle \langle s-\mathbf{C}\mathbf{f}\left(\gamma \right)\rangle \perp {\partial }_{\gamma }\mathbf{f}\left(\gamma \right)}\end{array}$$

(83)

resp.

$$\begin{array}{ccc}\hfill {\displaystyle \langle s,{\mathbf{n}}_{\partial \mathbf{f}}\rangle }& =& \langle {\mathbf{C}}_J\mathbf{f}\left(\gamma \right),{\mathbf{n}}_{\pi }\rangle \hfill \end{array}$$

(84)

which in leading order leads to

$$\begin{array}{ccc}\hfill {\displaystyle ϵ·\parallel {\partial }_{\gamma }{\mathbf{f}\left(\gamma \right)\parallel }^{-1}·\frac{2{\gamma }^2}{d·{\overline{T}}^3}·{\parallel {\widehat{\sigma }}_x\overline{e}\parallel }^2}& =& {\gamma }^2·a_{\pi ,\pi }^J\hfill \end{array}$$

(85)

With $\parallel {\partial }_{\gamma }\mathbf{f}\parallel =\mathcal{O}\left(\gamma \right)$ we finally end up with

(5.11) Balanced state (Relaxation with source): Under the smallness assumptions for $\gamma$ in (5.9) and with an appropriate constant $C_{\sigma }>0$, the leading order term of the stable solution is given by

$$\begin{array}{c}\hfill {\displaystyle {\gamma }_{bs}=ϵ·C_{\sigma }}\end{array}$$

(86)

The solution $\mathbf{f}\left({\gamma }_{bs}\right)$ has an additonal contribution to the pressure tensor of the order $\mathcal{O}\left({ϵ}^2\right)$ as given in (58) .

• This example is an indicator that the concept of traces is relevant for the investigation of the structure of kinetic solutions in the context of closure relations. We can generalize the above principle to other sources and kinetic models.

(5.12) Balancing condition: Suppose given a source $ϵs$ and a collision operator $Cf$ with traces

$$\begin{array}{ccc}\hfill {\displaystyle ϵ·\mathbf{s}\left(\gamma \right)}& =& ϵ·s_{\pi }·{\mathbf{n}}_{\pi }+ϵ·s_{\sigma }·{\mathbf{n}}_{\sigma }\hfill \end{array}$$

(87)

$$\begin{array}{ccc}\hfill \mathbf{Cf}\left(\gamma \right)& =& c_{\pi }·{\mathbf{n}}_{\pi }+c_{\sigma }·{\mathbf{n}}_{\sigma }\hfill \end{array}$$

(88)

$\gamma$ marks a balanced state if and only if the orthogonality condition

$$\begin{array}{ccc}\hfill {\displaystyle ϵ·\mathbf{s}-\mathbf{Cf}}& \perp & {\partial }_{\gamma }\mathbf{f}=d_{\pi }·{\mathbf{n}}_{\pi }+d_{\sigma }·{\mathbf{n}}_{\sigma }\hfill \end{array}$$

(89)

is satisfied which is equivalent to the balancing condition

$$\begin{array}{ccc}\hfill {\displaystyle (ϵ·s_{\pi }-c_{\pi })·d_{\pi }}& =& -(ϵ·s_{\sigma }-c_{\sigma })·d_{\sigma }\hfill \end{array}$$

(90)

6. Applications

6.1. Heat layers with zero flow

We consider the steady flow between two totally reflecting walls at different temperatures $T_0=T(x=-1)>T_1=T(x=1)$ with a temperature difference $\Delta T=T_0-T_1>0$. Since $\overline{v}=0$, the moment equations read

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_x{p}_{xx}}& =& 0\hfill \\ \hfill {\partial }_x{q}_x& =& 0\hfill \end{array}$$

The closure moments for the solution are determined from the asymptotic solution $f_{corr}$ for the model at hand of

$$\begin{array}{ccc}\hfill {\displaystyle {\partial }_x{f}_{corr}+\frac{ϵ}{2T^2}·{\widehat{\sigma }}_x}& =& Cf\hfill \end{array}$$

In the Navier-Stokes description, $f_{corr}$ contains no contribution to ${\widehat{\pi }}_x$, and from equations ... we find the heat layer solution

$$\begin{array}{ccc}\hfill {\displaystyle p_{xx}=p=\rho ·T}& =& const\hfill \\ \hfill \lambda \left(T\right)·{\partial }_x\left(T^{-1}\right)& =& const\hfill \end{array}$$

Comparing the temperature profiles of numerical simulations for the nonlinear and the BGK system (Figure 3) we find that the BGK profile is considerably flatter than that of the nonlinear case. The latter one is close to the constant gradient profile (solid line) linearly connecting the wall temperatures. An explanation for th the difference is given by comparing the closure moments $\langle {\widehat{\pi }}_x,f_{corr}\rangle$ for both models (Figure 4). In Figure 4(a) we find a small almost constant contribution (maybe weakly dependent on the temperature) in the nonlinear case. Zooming into Figure 4(b) confirms this. However, for BGK we find a distinct, almost constant gradient over the whole field of calculation. An explanation for this is again provided by the trace description which allows to draw a rough picture pointing out the differences between nonlinear collision operator and BGK model. A full discussion of the trace solutions for both models in the fluid dynamic limit is provided in [5]. We do not repeat the tedious calculations of [5] but explain by a short argument why the BGK system does not provide the correct temperature gradient in the fluid dynamic limit. Suppose $\mathbf{t}\left[\gamma \right]=\overline{\mathbf{e}}+\mathrm{\Phi }$ is the steady solution of the balance equation for the BGK system,

$$\begin{array}{ccc}\hfill {\displaystyle ϵ·v_x{\partial }_x(\overline{\mathbf{e}}+\mathrm{\Phi })}& =& \lambda ·(\overline{\mathbf{e}}-\mathbf{t}\left[\gamma \right])\hfill \end{array}$$

(91)

The first term on the left hand side represents in the limit a finite heat flux modeled by some term of the form $B·{\sigma }_x·\overline{\mathbf{e}}$. For small $\gamma$ we may use the approximation (65) for $\mathbf{t}\left[\gamma \right]$ shortly denoted as $(a{\gamma }^2·{\pi }_x+b{\gamma }^3·{\sigma }_x)·\overline{\mathbf{e}}$. Furthermore, $v_x{\partial }_x\mathrm{\Phi }$ yields an unknown contribution $A·{\pi }_x\overline{\mathbf{e}}$ which comes out as a solvability condition for the above equation,

$$\begin{array}{ccc}\hfill {\displaystyle ϵ·(A{\pi }_x+B{\sigma }_x)}& =& \lambda (a{\gamma }^2·{\pi }_x+b{\gamma }^3·{\sigma }_x)\hfill \end{array}$$

(92)

Its solution is

$$\begin{array}{c}\hfill {\displaystyle \gamma ={ϵ}^{1/3}·{\left(B/\left(\lambda b\right)\right)}^{1/3}=:c·{ϵ}^{1/3}}\end{array}$$

(93)

and

$$\begin{array}{c}\hfill {\displaystyle A=\lambda ac{ϵ}^{-1/3}}\end{array}$$

(94)

Thus under the above assumptions A is singular in the limit $ϵ\searrow 0$ which is not realistic. The only way out is that c and with this B are $ϵ$-dependent and vanish in the limit. This explains qualitatively the flattened curve for the temperature profile of the BGK model. In the limit it is supposed to yield a constant temperature profile.

Remark: The results show that trace theory does not fit into the classical series expansion ansatz since it contains fractional powers of $ϵ$.

6.2. Gas mixtures: An evaporation condensation problem

The initial motivation for the present work came from the study of a gas mixture problem. The two-component gas consists of vapor (species A) which can evaporate or condensate at the flat confining walls of an infinite strip, and of some noncondensable component (species B) which is totally reflected (with energy exchange) at the walls. Given a pressure difference, the gas is attracted by one of the walls, and some portion of species A condenses, which results in a flow through the wall. Suppose there is only a small (almost neglegible amount of species B in the mixture. Calculating the fluid dynamic limit applying standard expansion techniques (here: Hilbert expansion) leads to a ghost effect: An infinitely thin boundary layer of the noncondensable completely cuts down the flux of species A, although in the limit the flow velocity of the ensemble is zero. Details of this example are described in [1] and the literature cited there. In [4] we have demonstrated that the reason for the confusion lies in some bifurcation structure of the solution space in the case of small velocities. We are now able to present the full solution to the above problem using the trace theory approach.

Deriving a trace theory for a gas mixture could closely follow the lines developed in the preceding sections.We will not do this here but restrict to a subproblem. We have numerically simulated the evaporation condensation problem and ended up with results part of which is given in Figure 5. Here we find some boundary layer close to the left wall (which is large scale and thus not a kinetic boundary layer).Far apart from the wall the whole ensemble seems to tend to some equilibrium state. This cannot be a common Maxwellian, since Species A and B exhibit different velocities and temperatures. Such an equilibrium state does not exist in classical theory. However, trace theory approach results in a common balanced state $({\mathbf{f}}^A,{\mathbf{f}}^B)$ of the mixture far on the right of the perturbing wall. In particular, we find the dependence of the condensation rate for A on the interaction between both species, and on the energy exchange of species B on the left wall. Here, we roughly sketch the corresponding calculations.

Suppose the restrictions ${\mathbf{f}}^A$ and ${\mathbf{f}}^B$ to their species have moments $({\rho }^A,-{\overline{v}}^A,{\overline{T}}^A)$ and $({\rho }^B,0,{\overline{T}}^B)$ with corresponding Maxwellians ${\rho }^A{\overline{\mathbf{e}}}^A$ and ${\rho }^B{\overline{\mathbf{e}}}^B$. Constructing $\mathbf{f}$ we assume that both ${\mathbf{f}}^A$ and ${\mathbf{f}}^B$ are trace elements of their respective Maxwellians, ${\mathbf{f}}^A={\mathbf{t}}^A\left[{\gamma }^A\right]$ and ${\mathbf{f}}^B={\mathbf{t}}^B\left[{\gamma }^B\right]$ with corresponding parameters given by formulas (54) and (55),

$$\begin{array}{c}\hfill {\displaystyle {\widehat{T}}^S={\overline{T}}^S+d^{-1}{\left({\gamma }^S\right)}^2,{\widehat{v}}^S={\overline{v}}^S-{\gamma }^S,S\in \{A,B\}}\end{array}$$

(95)

Compatibility conditions are the equalities ${\tilde{T}}^A={\tilde{T}}^B=:\tilde{T}$ and ${\tilde{v}}_x^A={\tilde{v}}_x^B=:{\tilde{v}}_x$, i.e.

$$\begin{array}{c}\hfill {\displaystyle {\gamma }^B={\gamma }^A-{\overline{v}}_x^A,d·({\overline{T}}^A-{\overline{T}}^B)={\overline{v}}_x^A·({\overline{v}}_x^A-2{\gamma }^A)}\end{array}$$

(96)

These fix ${\gamma }^B$ once ${\gamma }^A$ is given, and the temperature difference in relation to the velocity of A. Both trace elements are now represented by and perturbations of the same Maxwellian $\tilde{\mathbf{e}}=\mathbf{e}[{\tilde{v}}_x,\tilde{T}]$.

Let a nonlinear two-species collision operator be given as

$$\begin{array}{ccc}\hfill {\displaystyle J[\mathbf{f},\mathbf{f}]}& =& J_A[{\mathbf{f}}^A,{\mathbf{f}}^A]+J_{AB}[{\mathbf{f}}^A,{\mathbf{f}}^B]+J_{BA}[{\mathbf{f}}^A,{\mathbf{f}}^B]+J_B[{\mathbf{f}}^B,{\mathbf{f}}^B]\hfill \end{array}$$

(97)

with restrictions to the species being

$$\begin{array}{ccc}\hfill {\displaystyle {\left.J[\mathbf{f},\mathbf{f}]\right|}_A}& =& J^A[{\mathbf{f}}^A,{\mathbf{f}}^A]+J_{BA}[{\mathbf{f}}^B,{\mathbf{f}}^A]\hfill \end{array}$$

(98)

$$\begin{array}{ccc}\hfill {\left.J[\mathbf{f},\mathbf{f}]\right|}_B& =& J_{AB}[{\mathbf{f}}^A,{\mathbf{f}}^B]+J_B[{\mathbf{f}}^B,{\mathbf{f}}^B]\hfill \end{array}$$

(99)

The operators $J_A$ and $J_B$ are the usual one-species collision operators, while $J_{AB}$ and $J_{BA}$ represent the interspecies interactions. The above trace construction is a balanced solution with respect to this operator, if it satisfies the balance conditions

$$\begin{array}{ccc}\hfill {\displaystyle \langle {\left.J[\mathbf{f},\mathbf{f}]\right|}_A,{\pi }_x^A\rangle }& =& 0\hfill \end{array}$$

(100)

$$\begin{array}{ccc}\hfill \langle {\left.J[\mathbf{f},\mathbf{f}]\right|}_B,{\pi }_x^B\rangle & =& 0\hfill \end{array}$$

(101)

where the closure moments ${\pi }_x^S$ have to be taken relative to their corresponding bulk velocities ${\overline{v}}_x^S$. From this we extract conditions of the form

$$\begin{array}{c}\hfill {\displaystyle {\alpha }^A{\rho }^A{\left({\gamma }^A\right)}^2={\rho }^B{\beta }^A,{\alpha }^B{\rho }^B{\left({\gamma }^B\right)}^2={\rho }^A{\beta }^B}\end{array}$$

(102)

with ${\alpha }^S$ and ${\beta }^S$ positive constants to be derived from the collision operator. (For ${\alpha }^S$ compare (71). ${\beta }^S$ can be derived similarly from the interspecies operators.). From the first condition we can calculate ${\gamma }^A$, the second one relates the velocity ${\overline{v}}_x^A$ to the temperature $\tilde{T}$.