3.1. Definition of moments

Moments are weighted averages of the distribution function, some of which have direct physical meaning, such as the mass density $\rho$, or the bulk velocity $v_{i}$. Following Struchtrup (Reference Struchtrup2005), we consider the usual central moments of the distribution function defined as

(3.1)\begin{equation} u_{\langle i_{1}\cdots i_{n}\rangle}^{\alpha}=m\int C^{2\alpha} C_{\langle i_{1}}\cdots C_{i_{n}\rangle}f\,{\rm d}\boldsymbol{c}, \end{equation}

where $C_{i}=c_{i}-v_{i}$ is the peculiar velocity and $m$ is particle mass; indices between angular brackets indicate trace-free and symmetric tensors. Experience from the Boltzmann equation (Struchtrup Reference Struchtrup2005; Torrilhon Reference Torrilhon2016) shows that the following 26 variables serve well to describe gas flows at not too large Knudsen numbers:

(3.2)\begin{equation} \left.\begin{gathered} \rho =mn=u^{0}, \\ \rho v_{i} =m\int c_{i}f\,{\rm d}\boldsymbol{c},\\ \rho\epsilon =\tfrac{3}{2}\rho RT=\tfrac{3}{2}\rho\theta=\tfrac{1}{2} u^{2}, \\ \sigma_{ij} =u_{\langle ij\rangle}^{0},\\ q_{i} =\tfrac{1}{2}u_{i}^{1}, \\ m_{ijk} =u_{\langle ijk\rangle}^{0}, \\ \varDelta =u^{4}-15\rho\theta^{2}, \\ R_{ij} =u_{\langle ij\rangle}^{1}-7\theta\sigma_{ij}. \end{gathered}\right\} \end{equation}

Above, for the first five moments, $\rho$, $\rho v_{i}$ and $\rho \epsilon$ are mass density, momentum density, and density of thermal energy, where $\theta =RT$ is the temperature in units of specific energy, $R=k_{B}/m$ is the gas constant, $k_{B}$ is the Boltzmann constant, and $T$ is thermodynamic temperature. The other variables are all non-equilibrium moments, that is, they vanish in equilibrium states. Stress tensor $\sigma _{ij}$ and heat flux $q_{i}$ are the kinetic contributions to overall stress and energy flux; there will be additional contributions from the Enskog and Vlasov terms that will emerge further below. The order of magnitude method for classical rarefied gas dynamics for Maxwell molecules shows that stress and heat flux are of first order in the Knudsen number, while the variables $m_{ijk}$, $\varDelta$ and $R_{ij}$ are constructed to be of second order (Struchtrup Reference Struchtrup2004, Reference Struchtrup2005). Since we expect a close relationship between moment equations from the EV theory and classical kinetic theory, we prefer to use the same variables.

3.2. Grad distribution function

Moment equations result from multiplying the kinetic equation with appropriate tensor combinations of microscopic velocity, and integration. The moment equations are coupled, since higher moments always appear under space derivatives. When one has chosen a certain set of moments as variables, the corresponding set of moment equations contains higher moments that must be related somehow to the chosen variables. In the Grad method, this closure problem is solved by constructing a distribution function that reproduces the variables, and allows us to determine all other moments through the variables (Grad Reference Grad1949, Reference Grad1958; Struchtrup Reference Struchtrup2005; Kremer Reference Kremer2010).

Just as for the Boltzmann equation, the Grad distribution is constructed as a deviation from a local Maxwellian, whose macroscopic parameters coincide with the corresponding moments of $f$. With the 26 moments listed above, closure is performed by the 26 moment distribution function (Struchtrup Reference Struchtrup2005), which reads

(3.3)\begin{equation} f_{\mid26}=\hat{f}_{M}\left( \rho+\varPhi\right), \end{equation}

with the reduced Maxwellian

(3.4)\begin{equation} \hat{f}_{M}=\frac{1}{m\sqrt{2{\rm \pi}\theta}^{3}}\,{\rm e}^{-({C^{2}}/{2\theta})}, \end{equation}

and the non-equilibrium deviation $\varPhi$ expressed as a polynomial in the peculiar velocity,

(3.5)\begin{align} \varPhi &= \frac{\sigma_{ij}}{2\theta^{2}}C_{\langle i}C_{j\rangle}+\frac{2}{5}\, \frac{q_{k}}{\theta^{2}}C_{k}\left( \frac{C^{2}}{2\theta}-\frac{5}{2}\right) \nonumber\\ &\quad +\frac{m_{ijk}}{6\theta^{3}}\,C_{\langle i}C_{j}C_{k\rangle}+ \frac{R_{ij}}{14\theta^{3}}\,C_{\langle i}C_{j\rangle}\left( \frac{C^{2}}{2\theta}- \frac{7}{2}\right) +\frac{\varDelta}{8\theta^{2}}\left( 1-\frac{2}{3}\, \frac{C^{2}}{\theta}+\frac{1}{15}\,\frac{C^{4}}{\theta^{2}}\right). \end{align}

Here, the first two terms correspond to Grad's classical $13$ moments closure. In equilibrium, all higher moments vanish, $0=\varDelta =q_{k}=\sigma _{ij} =R_{ij}=m_{ijk}$, hence the equilibrium solution is the Maxwellian $f_{{M}}=\rho \hat {f}_{{M}}$. More precisely, it is worth observing that the equilibrium solution of the EV equation is a Maxwell distribution function with uniform $\theta$ and $\boldsymbol {v}$. Depending on the equilibrium value of $\theta$, the equilibrium density $\rho$ is not generally uniform, obeying an integral equation expressing the local balance of mean field and hard-sphere pressure tensor contributions. Also, for later reference, we emphasize that the non-equilibrium correction $\varPhi$ is independent of mass density.

This distribution (3.3) reproduces the 26 moments of the theory (see (3.2)), and the higher moments required for closure are determined as (Struchtrup Reference Struchtrup2005)

(3.6a–c)\begin{equation} u_{i\mid26}^{2}=28\theta q_{i},\quad u_{ijk\mid26}^{1}=9\theta m_{ijk},\quad u_{ijkl\mid26}^{0}=0. \end{equation}

This distribution will also be used for determining the Enskog production terms, which result from taking appropriate moments of the collision term.

A distribution function should be positive. Due to its polynomial structure, for the Grad approximation this is not the case, since $\varPhi$ becomes negative for large values of $C_{i}$. For not too large values of the moments in (3.5), the multiplying Maxwellian in (3.3) suppresses these negative values, and the resulting equations remain meaningful.

The entropy inequality – also known as the H-theorem – results from integration of the kinetic equation against $\ln f$, which does not exist for the Grad distribution, due to the non-positivity. Evaluation of the entropy therefore must rely on Taylor expansion of $\ln \varPhi$ in terms of the moments. This implies that for Grad-type moment equations, an H-theorem can be only approximated, but not guaranteed for all values of the moments. One can show that fully linearized Grad equations as derived from the Boltzmann equation are accompanied by a quadratic entropy; see Struchtrup & Torrilhon (Reference Struchtrup and Torrilhon2007). While an H-theorem for nonlinear moment equations cannot be proved, Grad-type moment equations remain stable and give meaningful results also outside the linear regime; e.g. see the review by Torrilhon (Reference Torrilhon2016). The same behaviour is expected from the moment equations derived below. An examination of the accompanying H-theorem, at least for the fully linearized equations, is planned for the future.

3.3. Moment equations

3.3.1. General form

Moment equations are formed by multiplying the kinetic equation with the appropriate polynomial in velocity, and subsequent integration. The general moment equation for arbitrary moments (3.1) was presented in Struchtrup (Reference Struchtrup2005) for the Boltzmann equation. Since EV differs only in the Vlasov force and the details of the collision term, we can take the general moment equation from Struchtrup (Reference Struchtrup2005), add the Vlasov force $F_{k}$, and keep the general notation for production terms, so that

(3.7)\begin{align} & \frac{{\rm D} u_{i_{1}\cdots i_{n}}^{\alpha}}{{\rm D} t}+2{\alpha} u_{i_{1}\cdots i_{n}k}^{\alpha -1}\left[\frac{{\rm D} v_{k}}{{\rm D} t}-F_{k}-G_{k}\right] \nonumber\\ &\quad +\frac{n\left(2\alpha+2n+1\right) }{2n+1}\,u_{\langle i_{1}\cdots i_{n-1}}^{\alpha} \left[\frac{{\rm D} v_{i_{n}\rangle}}{{\rm D} t}-F_{i_{n}\rangle}-G_{i_{n}\rangle}\right] \nonumber\\ &\quad +\frac{\partial u_{i_{1}\cdots i_{n}k}^{\alpha}}{\partial x_{k}}+ \frac{n}{2n+1}\,\frac{\partial u_{\langle i_{1}\cdots i_{n-1}}^{\alpha+1}}{\partial x_{i_{n}\rangle}}+ 2\alpha u_{i_{1}\cdots i_{n}kl}^{\alpha-1}\,\frac{\partial v_{k}}{\partial x_{l}} \nonumber\\ &\quad +2\alpha\,\frac{n+1}{2n+3}\,u_{\langle i_{1}\cdots i_{n}}^{\alpha}\, \frac{\partial v_{k\rangle}}{\partial x_{k}}+2\alpha\,\frac{n}{2n+1}\,u_{k\langle i_{1}\cdots i_{n-1}}^{\alpha}\, \frac{\partial v_{k}}{\partial x_{i_{n}\rangle}}+nu_{k\langle i_{1}\cdots i_{n-1}}^{\alpha}\, \frac{\partial v_{i_{n}\rangle}}{\partial x_{k}} \nonumber\\ &\quad +u_{i_{1}\cdots i_{n}}^{\alpha}\,\frac{\partial v_{k}}{\partial x_{k}} +\frac{n( n-1) }{4n^{2}-1}\left( 2\alpha+2n+1\right) u_{\langle i_{1}\cdots i_{n-2}}^{\alpha+1}\, \frac{\partial v_{i_{n-1}}}{\partial x_{i_{n}\rangle}}=\mathcal{P}_{i_{1}\cdots i_{n}}^{\alpha}. \end{align}

Here, ${\textrm {D}}/{\textrm {D} t}={\partial }/{\partial t}+v_{k}({\partial }/{\partial x_{k}})$ is the convective time derivative, $G_{i}$ is the external body force, e.g. from gravity, and $F_{i}$ is the Vlasov force from above. The so-called production terms on the right-hand side are integrals over the Enskog collision term,

(3.8)\begin{equation} \mathcal{P}_{i_{1}\cdots i_{n}}^{\alpha}=m\int C^{2\alpha}C_{\langle i_{1}} \cdots C_{i_{n}\rangle}\mathcal{S}_{{En}}\,{\rm d}\boldsymbol{c}. \end{equation}

The general moment equation contains time derivatives of velocity that enter through derivatives of the peculiar velocity $C_{i}=c_{i}-v_{i}$; the balance of momentum will be used to replace these contributions. Equations for particular moments follow by choice of particular values of the velocity power $\alpha$ and the tensor rank $n$, and replacing moments $u_{\langle i_{1}\cdots i_{n}\rangle }^{\alpha }$ through the variables listed in (3.2). Details of the Vlasov force and production terms will be discussed in later sections.

3.3.2. Conservation laws

We begin with the conservation laws, in which we have to write production terms for energy and momentum, since there will be flux contributions (in divergence form for the full conservation laws) from the Enskog terms for momentum and energy. This accounts for the fact that in a hard-sphere collision, energy and momentum are transferred over the distance of the colliding particles’ centre of mass. Since mass is not transferred between the particles, so that there is no production of mass, $m\int \mathcal {S}_{{En}}\,\textrm {d}\boldsymbol {c}=0$, and the production of momentum is $m\int c_{i}\mathcal {S}_{{En}}\,\textrm {d}\boldsymbol {c}=m\int C_{i}\mathcal {S}_{{En}}\,\textrm {d}\boldsymbol {c}=\mathcal {P}_{i}^{0}$.

The conservation laws for mass and momentum, and the balance of thermal energy, read

(3.9)\begin{equation} \left.\begin{gathered} \frac{{\rm D}\rho}{{\rm D} t}+\rho\,\frac{\partial v_{k}}{\partial x_{k}} =0, \\ \rho\,\frac{{\rm D} v_{i}}{{\rm D} t}+\frac{\partial\rho\theta}{\partial x_{i}}+ \frac{\partial\sigma_{ik}}{\partial x_{k}} =\rho\left( F_{i}+G_{i}\right) +\mathcal{P}_{i}^{0},\\ \frac{3}{2}\rho\,\frac{{\rm D}\theta}{{\rm D} t}+\frac{\partial q_{k}}{\partial x_{k}} ={-}\rho\theta\,\frac{\partial v_{k}}{\partial x_{k}}-\sigma_{kl}\, \frac{\partial v_{k}}{\partial x_{l}}+\frac{1}{2}\mathcal{P}^{1}. \end{gathered}\right\} \end{equation}

Here, and further below, we removed the velocity time derivative from (3.7) by means of the momentum balance, which also removes the contributions of external body force and Vlasov force. Hence the force terms will appear only in the momentum balance, but not in the other equations for central moments.

3.3.3. Balance for $\varDelta$

The conservation laws are used to simplify the equations for other variables, by eliminating the time derivative of velocity, and use of the energy equation and mass balance where required. We proceed through several tensor orders, beginning here with the scalar equation for $\varDelta =u^{2}-15\rho \theta ^{2}$, which requires not only the moment equation for $u^{2}$, but also mass and energy balances to eliminate time derivatives of $15\rho \theta ^{2}$. Omitting all detail of computation, we simply state the final result, based on the 26 moments closure (3.3). We obtain the same equation as in Struchtrup (Reference Struchtrup2005), only with additional terms in the production, which come from the Enskog contributions to energy and momentum and will be discussed later:

(3.10)\begin{align} & \frac{{\rm D}\varDelta}{{\rm D} t}+8\theta\,\frac{\partial q_{k}}{\partial x_{k}}+8\theta \sigma_{kl}\,\frac{\partial v_{k}}{\partial x_{l}}+20q_{k}\, \frac{\partial\theta}{\partial x_{k}}-8q_{k}\theta\,\frac{\partial\ln\rho}{\partial x_{k}} -\frac{8q_{k}}{\rho}\,\frac{\partial\sigma_{kl}}{\partial x_{l}}+4R_{kl}\, \frac{\partial v_{k}}{\partial x_{l}}+\frac{7}{3}\varDelta\,\frac{\partial v_{k}}{\partial x_{k}} \nonumber\\ &\quad =\mathcal{P}^{2}-20\theta\,\frac{\mathcal{P}^{1}}{2}-\frac{8q_{k}}{\rho}\,\mathcal{P}_{k}^{0}. \end{align}

3.3.4. Balance for $q_{i}$

The equation for the kinetic contribution to energy flux, $q_{i}$, is the same as in standard kinetic theory (Struchtrup Reference Struchtrup2005), except that the Enskog term yields additional production terms:

(3.11)\begin{align} & \frac{{\rm D} q_{i}}{{\rm D} t}+\frac{5}{2}\rho\theta\,\frac{\partial\theta}{\partial x_{i}}+ \frac{5}{2}\sigma_{ik}\,\frac{\partial\theta}{\partial x_{k}}+\theta\, \frac{\partial\sigma_{ik}}{\partial x_{k}}-\sigma_{ik}\theta\,\frac{\partial \ln\rho}{\partial x_{k}}+\frac{7}{5}q_{k}\,\frac{\partial v_{i}}{\partial x_{k}}+ \frac{2}{5}q_{k}\,\frac{\partial v_{k}}{\partial x_{i}}+\frac{7}{5}q_{i}\, \frac{\partial v_{k}}{\partial x_{k}} \nonumber\\ &\quad +\frac{1}{2}\,\frac{\partial R_{ik}}{\partial x_{k}}+\frac{1}{6}\, \frac{\partial\varDelta}{\partial x_{i}}+m_{ikl}\,\frac{\partial v_{k}}{\partial x_{l}}- \frac{\sigma_{ik}}{\rho}\,\frac{\partial\sigma_{kl}}{\partial x_{l}}= \frac{1}{2}\mathcal{P}_{i}^{1}-\frac{5}{2}\theta\mathcal{P}_{i}^{0}- \frac{\sigma_{ik}}{\rho}\,\mathcal{P}_{k}^{0}. \end{align}

3.3.5. Balances for $\sigma _{ij}$ and $R_{ij}$

Also, the balances for stress $\sigma _{ij}$ and the variable $R_{ij}$ agree with results in (Struchtrup Reference Struchtrup2005), except that additional production terms appear. We find

(3.12)\begin{equation} \frac{{\rm D}\sigma_{ij}}{{\rm D} t}+2\rho\theta\,\frac{\partial v_{\langle i}}{\partial x_{j\rangle}}+\sigma_{ij}\, \frac{\partial v_{k}}{\partial x_{k}}+\frac{4}{5}\, \frac{\partial q_{\langle i}}{\partial x_{j\rangle}}+2\sigma_{k\langle i}\, \frac{\partial v_{j\rangle}}{\partial x_{k}}+\frac{\partial m_{ijk} }{\partial x_{k}}=\mathcal{P}_{ij}^{0}, \end{equation}

and with the trace-free velocity gradient $S_{ij}={\partial v_{\langle i}}/{\partial x_{j\rangle }}$,

(3.13)\begin{align} & \frac{{\rm D} R_{ij}}{{\rm D} t}+\frac{28}{5}\theta\,\frac{\partial q_{\langle i}}{\partial x_{j\rangle}}+ \frac{28}{5}q_{\langle i}\,\frac{\partial\theta}{\partial x_{j\rangle}}-\frac{28}{5}\theta q_{\langle i}\, \frac{\partial\ln\rho}{\partial x_{j\rangle}}-\frac{28}{5}\,\frac{q_{\langle i}}{\rho}\,\frac{\partial \sigma_{j\rangle l}}{\partial x_{l}}+8\sigma_{k\langle i}S_{j\rangle k} \nonumber\\ &\qquad -\frac{14}{3}\,\frac{\sigma_{ij}}{\rho}\,\frac{\partial q_{k}}{\partial x_{k}}- \frac{14}{3}\,\frac{\sigma_{ij}\sigma_{kl}}{\rho}\,\frac{\partial v_{k}}{\partial x_{l}}+ \frac{14}{15}\varDelta\,\frac{\partial v_{\langle i}}{\partial x_{j\rangle}}- \frac{2m_{ijk}}{\rho}\,\frac{\partial\sigma_{kl}}{\partial x_{l}}-2m_{ijk}\theta\, \frac{\partial\ln\rho}{\partial x_{k}} \nonumber\\ &\qquad +7m_{ijk}\,\frac{\partial\theta}{\partial x_{k}}+2\theta\,\frac{\partial m_{ijk}}{\partial x_{k}}+ \frac{8}{7}R_{k\langle i}S_{j\rangle k}+2R_{k\langle i}\, \frac{\partial v_{k}}{\partial x_{j\rangle}}+\frac{5}{3}R_{ij}\,\frac{\partial v_{k}}{\partial x_{k}} \nonumber\\ &\quad =\mathcal{P}_{ij}^{1}-7\theta\mathcal{P}_{ij}^{0}-\frac{14}{3}\, \frac{\sigma_{ij}}{\rho}\,\frac{\mathcal{P}^{1}}{2}-2\,\frac{m_{ijk}}{\rho}\, \mathcal{P}_{k}^{0}-\frac{28}{5}\,\frac{q_{\langle i}}{\rho}\,\mathcal{P}_{j\rangle}^{0}. \end{align}

3.3.6. Balance for $m_{ijk}$

Finally, for the third moment we find

(3.14)\begin{align} & \frac{{\rm D} m_{ijk}}{{\rm D} t}+3\theta\,\frac{\partial\sigma_{\langle ij}}{\partial x_{k\rangle}}-3 \theta\sigma_{\langle ij}\,\frac{\partial\ln\rho}{\partial x_{k\rangle}}+\frac{12}{5}q_{\langle i}\, \frac{\partial v_{j}}{\partial x_{k\rangle}}-\frac{3\sigma_{\langle ij}}{\rho}\, \frac{\partial\sigma_{k\rangle l}}{\partial x_{l}}+\frac{3}{7}\, \frac{\partial R_{\langle ij}}{\partial x_{k\rangle}} \nonumber\\ &\quad +3m_{l\langle ij}\,\frac{\partial v_{k\rangle}}{\partial x_{l}}+m_{ijk}\, \frac{\partial v_{l}}{\partial x_{l}}=\mathcal{P}_{ijk}^{0}- \frac{3\sigma_{\langle ij}}{\rho}\,\mathcal{P}_{k\rangle}^{0}. \end{align}

4.1. Smooth density approximation

As shown in (2.7), the Vlasov force $F_{k}$ depends solely on the density distribution of the particles, independent of their velocities. Therefore, the full distribution function $f$ is not required for its computation. Assuming that the number density $n$ is sufficiently smooth anywhere, including at the phase interface, we can expand the number density in the integral, so that

(4.1)\begin{equation} F_{k}( x_{l}) =\frac{\partial n}{\partial x_{s}}\int_{r>a} \frac{{\rm d}\phi}{{\rm d} r}\,\frac{r_{k}r_{s}}{r}\,{\rm d}\boldsymbol{r}+\frac{1}{6}\, \frac{\partial ^{3}n}{\partial x_{r}\,\partial x_{s}\,\partial x_{t}}\int_{r>a} \frac{{\rm d}\phi}{{\rm d} r}\,\frac{r_{k}r_{r}r_{s}r_{t}}{r}\,{\rm d}\boldsymbol{r}, \end{equation}

where $r_{k}=x_{k}^{1}-x_{k}$ is the distance vector between the interacting particles, and we have switched integration variables (with $x_{k}^{1} =x_{k}+r_{k}$, we find $\textrm {d}\boldsymbol {x}\,\textrm {d}\boldsymbol {x}^{1}=\textrm {d}\boldsymbol {x}\,\textrm {d}\boldsymbol {r}$); anisotropic contributions of the integrals vanish, since we integrate over a sphere.

Integration yields the Vlasov force in the form

(4.2)\begin{equation} F_{k}( x_{l}) =\frac{4{\rm \pi} a^{3}}{3}\,\frac{\phi_{a}}{m}\,\chi_{1}\, \frac{\partial\rho}{\partial x_{k}}+\frac{2{\rm \pi} a^{5}}{15}\,\frac{\phi_{a} }{m}\,\chi_{3}\, \frac{\partial^{3}\rho}{\partial x_{k}\,\partial x_{s}\,\partial x_{s}}, \end{equation}

where for power potentials

(4.3a,b)\begin{equation} \chi_{1}=\frac{\gamma}{\gamma-3},\quad \chi_{3}=\frac{\gamma}{\gamma-5}, \end{equation}

and for the exponential tail

(4.4a,b)\begin{equation} \chi_{1}=\frac{6}{\varphi^{3}}+\frac{6}{\varphi^{2}}+\frac{3}{\varphi }+1,\quad \chi_{3}=\frac{120}{\varphi^{5}}+\frac{120}{\varphi^{4}} +\frac{60}{\varphi^{3}}+\frac{20}{\varphi^{2}}+\frac{5}{\varphi}+1. \end{equation}

4.2. Korteweg stress tensor

In the balance of momentum, the Vlasov force $F_{k}$ is multiplied by mass density $\rho$. It is a simple exercise to rewrite the force term as the divergence of the potential tail contribution to the Korteweg stress tensor, i.e.

(4.5)\begin{equation} \rho F_{k}={-}\frac{\partial\mathcal{T}_{ks}^{K}}{\partial x_{s}}, \end{equation}

with

(4.6)\begin{align} \mathcal{T}_{ks}^{K}={-}\frac{2{\rm \pi} a^{3}}{3}\,\frac{\phi_{a}}{m}\,\chi_{1}\rho^{2}\delta_{ks}- \frac{2{\rm \pi} a^{5}}{15}\,\frac{\phi_{a}}{m}\,\chi_{3} \left[\left(\rho\,\frac{\partial^{2}\rho}{\partial x_{r}\,\partial x_{r}}+\frac{1}{2}\, \frac{\partial\rho}{\partial x_{r}}\,\frac{\partial\rho}{\partial x_{r}}\right) \delta_{ks}-\frac{\partial\rho}{\partial x_{s}}\,\frac{\partial\rho}{\partial x_{k}}\right]. \end{align}

The first term is a van-der-Waals-like contribution describing pressure reduction through attractive forces, and the second term describes extra stress due to density gradients, including surface tension effects in phase interfaces (Anderson et al. Reference Anderson, McFadden and Wheeler1998).

4.3. Korteweg/Vlasov energy

The potential energy due to the attractive forces between particles is a contribution to the overall energy, and appears in the overall energy balance. Energy for the Korteweg stresses has a van-der-Waals-like contribution and a contribution due to density gradients (Anderson et al. Reference Anderson, McFadden and Wheeler1998)

(4.7)\begin{equation} \epsilon_{K}={-}\frac{2{\rm \pi} a^{3}}{3}\,\frac{\phi_{a}}{m}\,\chi_{1}\rho+\frac{{\rm \pi} a^{5}}{15}\, \frac{\phi_{a}}{m}\,\chi_{3}\,\frac{1}{\rho}\,\frac{\partial\rho}{\partial x_{s}}\, \frac{\partial\rho}{\partial x_{s}}, \end{equation}

which obeys the non-conservation law

(4.8)\begin{equation} \rho\,\frac{{\rm D}\epsilon_{K}}{{\rm D} t}+\frac{\partial q_{k}^{K}}{\partial x_{k} }={-}\mathcal{T}_{rs}^{K}\, \frac{\partial v_{r}}{\partial x_{s}},\end{equation}

where

(4.9)\begin{equation} q_{i}^{K}=\frac{2{\rm \pi} a^{5}}{15}\,\frac{\phi_{a}}{m}\,\chi_{3}\rho\, \frac{\partial\rho}{\partial x_{i}}\,\frac{\partial v_{r}}{\partial x_{r}} \end{equation}

is an extra contribution to energy flux (Giovangigli Reference Giovangigli2020).

The balance of kinetic energy results from multiplying the momentum balance with velocity:

(4.10)\begin{align} \rho\,\frac{{\rm D}\tfrac{1}{2}v^{2}}{{\rm D} t}+ \frac{\partial(\rho\theta\delta _{ik}+\sigma_{ik}+\mathcal{T}_{ik}^{K}) v_{i}}{\partial x_{k}}= (\rho\theta\delta_{ik}+\sigma_{ik}+\mathcal{T}_{ik}^{K}) \frac{\partial v_{i}}{\partial x_{k}}+\rho v_{i}G_{i}+v_{i}\mathcal{P}_{i}^{0}. \end{align}

Adding balances for thermal (third equation of the set (3.9)), potential (4.8) and kinetic (4.10) energies yields the balance of total energy in intermediate form:

(4.11)\begin{align} & \rho\,\frac{{\rm D}\left(\dfrac{3}{2}\theta+\dfrac{1}{2}v^{2}+\epsilon_{K}\right)}{{\rm D} t}+ \frac{\partial\left[ q_{k}+q_{k}^{K}+( \rho\theta\delta _{ik}+\sigma_{ik}+ \mathcal{T}_{ik}^{K}) v_{i}\right] }{\partial x_{k} } \nonumber\\ &\quad =\rho v_{i}G_{i}+ \frac{1}{2}\mathcal{P}^{1}+v_{i}\mathcal{P}_{i}^{0}. \end{align}

The energy balance still contains the energy and momentum contributions from the Enskog collision term, which will turn out to be of divergence form (Resibois & De Leener Reference Resibois and De Leener1977), that is, these are additional contributions to stress and heat flux.

5.1. Moment productions

We write in compact notation

(5.2)\begin{equation} P_{\psi}=\int\psi( \boldsymbol{c})\,\mathcal{S}_{{En}}\,{\rm d}\boldsymbol{c}. \end{equation}

Due to symmetry of collision details, we can use the same standard arguments on exchange and renaming of integration variables as for the Boltzmann equation (Struchtrup Reference Struchtrup2005; Kremer Reference Kremer2010) to compact the expression for the production terms into

(5.3)\begin{align} \mathcal{P}_{\psi} &= a^{2}\int\int\int_{0}^{2{\rm \pi}}\int_{0}^{{\rm \pi}/2} \left\{\psi^{\prime}-\psi\right\} Y\left[ n\left( x_{r}-\frac{ak_{r}}{2}\right)\right] f\left( x_{s}-ak_{s},c_{s}^{1}\right) f( x_{t},c_{t}) \nonumber\\ &\quad \times g\cos\theta\sin\theta \,{\rm d}\theta \,{\rm d}\varepsilon \,{\rm d}\boldsymbol{c}^{1}\,{\rm d}\boldsymbol{c}, \end{align}

where $\psi ^{\prime }=\psi ( \boldsymbol {c}^{\prime })$. Apart from the non-localities in the arguments of $n$ and $f$, this agrees with (3.18) in Struchtrup (Reference Struchtrup2005) for the Boltzmann equation. The production term for the mass balance is obtained for $\psi ( \boldsymbol {c}) =m$, so that $\psi ^{\prime }-\psi =0$, hence we see immediately that mass is conserved, $\mathcal {P}_{m}=\mathcal {P}^{0}=0$, as anticipated above.

The last step is insertion of the Grad distribution (3.3) into (5.3), and evaluation of the integral. We need to integrate over the collision vector $\boldsymbol {k}=\{ \cos \varepsilon \sin \theta,\sin \varepsilon \sin \theta,\cos \theta \}$, which occurs in the non-local argument of $n$ and $f$. To make this integration accessible, we perform expansions of the functions $Y[ n( x_{r}-{ak_{r}}/{2}) ]$ and $f( x_{s}-ak_{s},c_{s}^{1})$ around the point $\boldsymbol {x}$, as outlined in the following subsections. As noted by Frezzotti et al. (Reference Frezzotti, Gibelli and Lorenzani2005), it is necessary to retain terms up to $a^{5}$ in the Vlasov force, in order to obtain the correct soft tail contribution to the Korteweg capillary stress tensor. In the production terms, we have already the factor $a^{2}$ in front, hence, for consistency, we expand the integrand up to third-order terms in $a$.

5.2. Expansion of correlation function

We expand the function $Y[ n( x_{r}-{ak_{r}}/{2}) ]$. The Taylor series around $n( x_{r})$ reads

(5.4)\begin{equation} Y\left[ n\left( x_{r}-\frac{ak_{r}}{2}\right) \right] \simeq Y-Y_{r} ak_{r}+Y_{rs}a^{2}k_{r}k_{s}-Y_{rst}a^{3}k_{r}k_{s}k_{t}, \end{equation}

with

(5.5a–d)\begin{equation} Y=Y\left[ n(x_{r})\right],\quad Y_{r}=\frac{1}{2}\,\frac{\partial Y}{\partial x_{r}},\quad Y_{rs}=\frac{1}{8}\,\frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{s}},\quad Y_{rst}=\frac{1}{48}\,\frac{\partial^{3}Y}{\partial x_{r}\,\partial x_{s}\,\partial x_{t}}, \end{equation}

where $Y_{rs}$ and $Y_{rst}$ are fully symmetric.

5.3. Expansion of the distribution function

Also, the expansion of the distribution function is, in principle, straightforward:

(5.6)\begin{equation} f\left( x_{s}-ak_{s},c_{s}\right) =f-\frac{\partial f}{\partial x_{u}}\,ak_{u}+ \frac{1}{2}\,\frac{\partial^{2}f}{\partial x_{u}\,\partial x_{v}}\,a^{2} k_{u}k_{v}-\frac{1}{6}\, \frac{\partial^{3}f}{\partial x_{u}\,\partial x_{v}\,\partial x_{w}}\,a^{3}k_{u}k_{v}k_{w}. \end{equation}

To proceed, we have to insert $f_{|26}$ from (3.3) and consider the derivatives carefully. Since $f_{|26}$ depends explicitly on the 26 variables $\{\rho,v_{i},\theta,\varDelta,\sigma _{ij},q_{i},R_{ij},m_{ijk}\}$, the derivatives of $f_{|26}$ give rise to a large number of linear and nonlinear terms. We will aim for a mainly linear theory, and ignore all nonlinear terms in the non-equilibrium moments and the derivatives of all variables, except density. Density changes strongly over the equilibrium phase interface, and its gradient terms appear in the Vlasov force. For consistency, we have to have the same level of density gradient terms in the Enskog contributions. However, from (3.3), where the non-equilibrium contribution $\phi$ is independent of density, it is clear that $f_{|26}$ is linear in mass density.

From (3.3), for the first derivative of $f_{|26}$ we have

(5.7)\begin{equation} \frac{\partial f_{|26}}{\partial x_{u}}=\frac{\partial\hat{f}_{M}}{\partial x_{u}} \left( \rho+\varPhi\right) +\hat{f}_{M}\left( \frac{\partial\rho }{\partial x_{u}}+ \frac{\partial\varPhi}{\partial x_{u}}\right), \end{equation}

where

(5.8)\begin{equation} \frac{\partial\hat{f}_{M}}{\partial x_{u}}=\hat{f}_{M} \left[\left(\frac{C^{2}}{2\theta}-\frac{3}{2}\right) \frac{1}{\theta}\,\frac{\partial \theta}{\partial x_{u}}+ \frac{C_{i}}{\theta}\,\frac{\partial v_{i}}{\partial x_{u}}\right]. \end{equation}

The term $({\partial \hat {f}_{M}}/{\partial x_{u}})\,\varPhi$ is nonlinear in the non-equilibrium moments and gradients of temperature and velocity, hence will be dropped in linearization. Within the linearization, the derivative of $\varPhi$ is simply the same expression for $\varPhi$ as in (3.5), except that the moments $\varDelta$, $\sigma _{ij}$, $q_{i}$, $R_{ij}$, $m_{ijk}$ are replaced by their gradients.

Explicitly, we find for the first derivative (semi-linearized)

(5.9)\begin{equation} \frac{\partial f_{|26}}{\partial x_{u}}=\hat{f}_{M}\left[\begin{array}[c]{@{}c@{}} \dfrac{\partial\rho}{\partial x_{u}}+\dfrac{\rho}{\theta}\,\dfrac{\partial v_{i}}{\partial x_{u}}\,C_{i}+ \dfrac{\rho}{\theta}\,\dfrac{\partial\theta}{\partial x_{u}}\left(\dfrac{C^{2}}{2\theta}-\dfrac{3}{2}\right) \\ {}+\dfrac{1}{8\theta^{2}}\,\dfrac{\partial\varDelta}{\partial x_{u}} \left( 1-\dfrac{2}{3}\,\dfrac{C^{2}}{\theta}+\dfrac{1}{15}\,\dfrac{C^{4}}{\theta^{2}}\right) \\ {}+\dfrac{2}{5}\,\dfrac{1}{\theta^{2}}\,\dfrac{\partial q_{k}}{\partial x_{u}}\, C_{k}\left(\dfrac{C^{2}}{2\theta}-\dfrac{5}{2}\right) +\dfrac{1}{2\theta^{2} }\, \dfrac{\partial\sigma_{ij}}{\partial x_{u}}\,C_{\langle i}C_{j\rangle}\\ {}+\dfrac{1}{14\theta^{3}}\,\dfrac{\partial R_{ij}}{\partial x_{u}}\,C_{\langle i}C_{j\rangle} \left(\dfrac{C^{2}}{2\theta}-\dfrac{7}{2}\right) +\dfrac{1}{6\theta^{3}}\, \dfrac{\partial m_{ijk}}{\partial x_{u}}\,C_{\langle i}C_{j}C_{k\rangle} \end{array}\right]. \end{equation}

Generalization to higher derivatives yields (we require this for $n=1,2,3$)

(5.10)\begin{align} &\frac{\partial^{n}f_{|26}}{\partial x_{u_{1}}\cdots\partial x_{u_{n}}}\nonumber\\ &\quad = \hat{f}_{M}\left[\begin{array}[c]{@{}c@{}} \dfrac{\partial^{n}\rho}{\partial x_{u_{1}}\cdots\partial x_{u_{n}} }+ \mathcal{A}_{\left[ u_{1}\cdots u_{n}\right] i}\,\dfrac{C_{i}}{\theta} +\dfrac{1}{\theta}\,\mathcal{B}_{\left[ u_{1}\cdots u_{n}\right] }\left(\dfrac{C^{2}}{2\theta}-\dfrac{3}{2}\right) \\ {}+\dfrac{1}{8\theta^{2}}\,\dfrac{\partial^{n}\varDelta}{\partial x_{u_{1}} \cdots\partial x_{u_{n}}}\left( 1-\dfrac{2}{3}\,\dfrac{C^{2}}{\theta}+ \dfrac{1}{15}\,\dfrac{C^{4}}{\theta^{2}}\right) \\ {}+\dfrac{2}{5}\,\dfrac{1}{\theta^{2}}\,\dfrac{\partial^{n}q_{k}}{\partial x_{u_{1} } \cdots\partial x_{u_{n}}}\,C_{k}\left( \dfrac{C^{2}}{2\theta}-\dfrac{5}{2}\right) + \dfrac{1}{2\theta^{2}}\, \dfrac{\partial^{n}\sigma_{ij}}{\partial x_{u_{1}} \cdots\partial x_{u_{n}}}\,C_{\langle i}C_{j\rangle}\\ {}+\dfrac{1}{14\theta^{3}}\,\dfrac{\partial^{n}R_{ij}}{\partial x_{u}\, \partial x_{v}\,\partial x_{w}}\,C_{\langle i}C_{j\rangle} \left(\dfrac{C^{2}}{2\theta }-\dfrac{7}{2}\right) + \dfrac{1}{6\theta^{3}}\, \dfrac{\partial^{n}m_{ijk} }{\partial x_{u_{1}} \cdots\partial x_{u_{n}}}\,C_{\langle i}C_{j}C_{k\rangle} \end{array}\right] \end{align}

where

(5.11)\begin{equation} \left.\begin{gathered} \mathcal{A}_{\left[ u\right] i} =\left( \frac{\partial\rho v_{i} }{\partial x_{u}}-v_{i}\, \frac{\partial\rho}{\partial x_{u}}\right) =\left(\rho\,\frac{\partial v_{i}}{\partial x_{u}}\right), \\ \mathcal{A}_{\left[ uv\right] i} =\left( \frac{\partial^{2}\rho v_{i} }{\partial x_{u}\,\partial x_{v}}-v_{i}\, \frac{\partial^{2}\rho}{\partial x_{u}\,\partial x_{v}}\right), \\ \mathcal{A}_{\left[ uvw\right] i} =\left( \frac{\partial^{3}\rho v_{i} }{\partial x_{u}\, \partial x_{v}\,\partial x_{w}}-v_{i}\,\frac{\partial^{3}\rho }{\partial x_{u}\,\partial x_{v}\,\partial x_{w}}\right), \end{gathered}\right\} \end{equation}

and

(5.12)\begin{equation} \left.\begin{gathered} \mathcal{B}_{\left[ u\right] } =\left( \frac{\partial\rho\theta }{\partial x_{u}}-\theta\, \frac{\partial\rho}{\partial x_{u}}\right) =\left(\rho\,\frac{\partial\theta}{\partial x_{u}}\right), \\ \mathcal{B}_{\left[ uv\right] } =\left( \frac{\partial^{2}\rho\theta }{\partial x_{u}\,\partial x_{v}}-\theta\, \frac{\partial^{2}\rho}{\partial x_{u}\,\partial x_{v}}\right), \\ \mathcal{B}_{\left[ uvw\right] } =\left( \frac{\partial^{3}\rho\theta }{\partial x_{u}\, \partial x_{v}\,\partial x_{w}}-\theta\,\frac{\partial^{3}\rho }{\partial x_{u}\, \partial x_{v}\,\partial x_{w}}\right). \end{gathered}\right\} \end{equation}

Here, square brackets around indices are just a notational convenience, used for identification of different index groups.

5.4. Expansion of the production terms

We now insert the above expansions into the expression (5.3) for the moment productions, and consider contributions up to fifth-order terms in particle diameter $a$. We sort by order in $a$, to write the moment production (where $f$ stands for $f_{|26}$) as

(5.13)\begin{align} & \mathcal{P}_{\psi}=a^{2}Y\int\int f^{1}f\left[ g\int_{0}^{2{\rm \pi}} \int_{0}^{{\rm \pi}/2}\left( \psi^{\prime}-\psi\right) \cos\theta \sin\theta \,{\rm d}\theta \,{\rm d}\varepsilon\right] {\rm d}\boldsymbol{c}^{1}\,{\rm d}\boldsymbol{c} \nonumber\\ &\quad -a^{3}\int\int\left( Y_{r}f^{1}f+Y\,\frac{\partial f^{1}}{\partial x_{r}}\,f\right) \left[ g\int_{0}^{2{\rm \pi}}\int_{0}^{{\rm \pi}/2}\left( \psi^{\prime}-\psi\right) k_{r} \cos\theta\sin\theta \,{\rm d}\theta \,{\rm d}\varepsilon\right] {\rm d}\boldsymbol{c}^{1}\,{\rm d}\boldsymbol{c} \nonumber\\ &\quad +a^{4}\int\int\left( Y_{rs}f^{1}f+Y_{r}\,\frac{\partial f^{1}}{\partial x_{s}}\,f+ \frac{1}{2}Y\,\frac{\partial^{2}f^{1}}{\partial x_{r}\,\partial x_{s}}\,f\right) \nonumber\\ &\quad \times\left[ g\int_{0}^{2{\rm \pi}}\int_{0}^{{\rm \pi}/2}\left( \psi^{\prime} -\psi\right) k_{r}k_{s}\cos\theta\sin\theta \,{\rm d}\theta \,{\rm d}\varepsilon\right] {\rm d}\boldsymbol{c}^{1}\,{\rm d}\boldsymbol{c} \nonumber\\ &\quad -a^{5}\int\int\left(\begin{array}[c]{c} Y_{rst}f^{1}f+Y_{rs}\,\dfrac{\partial f^{1}}{\partial x_{t}}\,f \\ {}+\dfrac{1}{2}Y_{r}\,\dfrac{\partial^{2}f^{1}}{\partial x_{s}\,\partial x_{t}} f+ \dfrac{1}{6}Y\,\dfrac{\partial^{3}f^{1}}{\partial x_{r}\,\partial x_{s}\,\partial x_{t}}\,f \end{array}\right) \nonumber\\ &\quad \times\left[ g\int_{0}^{2{\rm \pi}}\int_{0}^{{\rm \pi}/2}\left( \psi^{\prime}-\psi\right) k_{r}k_{s}k_{t}\cos\theta\sin\theta \,{\rm d}\theta \,{\rm d}\varepsilon\right] {\rm d}\boldsymbol{c}^{1}\,{\rm d}\boldsymbol{c}. \end{align}

There will be further linearization in the terms quadratic in the distribution functions and their derivatives, where we will ignore nonlinear terms in non-equilibrium contributions. The first line above is the standard production for hard-sphere molecules; the others are correction terms from the Enskog term.

The Vlasov force has fifth-order terms that remain in equilibrium, and are required for the computation of the interface. For our first approach to moment equations for EV, in our linearization of $\mathcal {P}_{\psi }$ we will consider only those fifth-order terms that remain in equilibrium, while non-equilibrium contributions will be considered only up to fourth order in $a$.

The evaluation of these integrals is rather involved, and it is impossible to discuss all the details. In the following sections, we outline the main steps for the computation of the $\mathcal {P}_{\psi }$.

5.5. Velocity transformation

Computation of the production integrals is a standard procedure in kinetic theory; here, it is just more involved, due to the larger number of contributions to the productions. Solution of all integrals implies change of variables from $\boldsymbol {C},\boldsymbol {C}^{1}$ to

(5.14a,b)\begin{equation} g_{i}=C_{i}^{1}-C_{i},\quad G_{i}=\tfrac{1}{2}( C_{i}+C_{i}^{1}), \end{equation}

or the inverse,

(5.15a,b)\begin{equation} C_{i}=G_{i}-\tfrac{1}{2}g_{i},\quad C_{i}^{1}=G_{i}+\tfrac{1}{2}g_{i}. \end{equation}

Use of the peculiar velocity implies that we look at the collision terms from the rest frame of the gas, which is natural when we use central moments.

We also have

(5.16a–c)\begin{equation} C^{2}=G^{2}-G_{i}g_{i}+\frac{g^{2}}{4},\quad ( C^{1})^{2}=G^{2}+G_{i}g_{i}+\frac{g^{2}}{4},\quad C^{2}+( C^{1})^{2}=2G^{2}+\frac{g^{2}}{2} \end{equation}

and

(5.17)\begin{equation} {\rm d}\boldsymbol{c}\,{\rm d}\boldsymbol{c}^{1}={\rm d}\boldsymbol{C}\,{\rm d}\boldsymbol{C}^{1}= \left\Vert\frac{\partial( \boldsymbol{C,C}^{1}) }{\partial ( \boldsymbol{G,g}) }\right\Vert {\rm d}\boldsymbol{G}\,{\rm d}\boldsymbol{g}=\left\Vert \begin{array}[c]{cc} \delta_{ij} & -\tfrac{1}{2}\delta_{ij}\\[2pt] \delta_{ij} & +\tfrac{1}{2}\delta_{ij} \end{array}\right\Vert {\rm d}\boldsymbol{G}\,{\rm d}\boldsymbol{g}={\rm d}\boldsymbol{G}\,{\rm d}\boldsymbol{g}. \end{equation}

5.6. $k$-integration

We have $\psi ( C_{k}) =mC_{i_{1}}\cdots C_{i_{n}}$ and $C_{i}^{\prime }=C_{i}+k_{i}g\cos \theta$. In (5.13), the direction vector $k_{i}$ appears only in $\psi ^{\prime }$ and in the integrals in square brackets, but not in the distribution function terms. Therefore, we can deal with the $k$-integrals first.

The explicit dependence of $\psi$ on $k$ is as follows. For $\psi (C_{k}) =m$,

(5.18)\begin{equation} \frac{1}{m}\left( \psi^{\prime}-\psi\right) =0; \end{equation}

for $\psi ( C_{k}) =mC_{i}$,

(5.19)\begin{equation} \frac{1}{m}\left( \psi^{\prime}-\psi\right) =g\cos\theta k_{i}; \end{equation}

for $\psi ( C_{k}) =mC_{i}C_{j}$,

(5.20)\begin{equation} \frac{1}{m}\left( \psi^{\prime}-\psi\right) =2C_{(i}k_{j)}g\cos\theta +k_{i}k_{j}g^{2}\cos^{2}\theta; \end{equation}

for $\psi ( C_{k}) =mC_{i}C_{j}C_{k}$,

(5.21)\begin{equation} \frac{1}{m}\left( \psi^{\prime}-\psi\right) =3C_{(i}C_{j}k_{k)}g\cos \theta+3C_{(i}k_{j}k_{k)}g^{2}\cos^{2}\theta+k_{i}k_{j}k_{k}g^{3}\cos^{3}\theta; \end{equation}

for $\psi ( C_{k}) =mC^{2}C_{i}C_{j}$,

(5.22)\begin{align} \frac{1}{m}\left( \psi^{\prime}-\psi\right) &= 2C_{i}C_{j}C_{k}k_{k} g\cos \theta+2C^{2}C_{(i}k_{j)}g\cos\theta+C_{i}C_{j}g^{2}\cos^{2}\theta \nonumber\\ &\quad +4C_{k}C_{(i}k_{j)}k_{k}g^{2}\cos^{2}\theta+C^{2}k_{i}k_{j}g^{2}\cos^{2} \theta+2C_{(i}k_{j)}g^{3}\cos^{3}\theta \nonumber\\ &\quad +2C_{k}k_{i}k_{j}k_{k}g^{3}\cos^{3}\theta+k_{i}k_{j}g^{4}\cos^{4}\theta. \end{align}

Insertion of the above into the production terms shows that we are facing integrals, for $n\leq 3$ and $m\leq 3$, of the type

(5.23)\begin{equation} \left[ I_{m,n}^{\alpha}\right] _{i_{1}\cdots i_{n+m}}=g^{n+\alpha+1} \int_{0}^{2{\rm \pi}}\int_{0}^{{\rm \pi}/2}k_{i_{1}}\cdots k_{i_{n+m}}\cos^{n+\alpha+1} \theta\sin\theta \,{\rm d}\theta\,{\rm d}\varepsilon, \end{equation}

where $k_{k}=\{\cos \varepsilon \sin \theta,\sin \varepsilon \sin \theta,\cos \theta \}_{k}$ and $k_{k}g_{k}=g\cos \theta$. Explicit expressions for the integrals required are given in Appendix A.

5.7. Explicit productions

The next step is to insert the integrals $[I_{m,n}^{\alpha }]_{i_{1}\cdots i_{n+m}}$ into (5.13), which gives explicit integrals of $f^{1}f$ and $({\partial ^{n}f_{|26}^{1}}/{\partial x_{u_{1}}\cdots \partial x_{u_{n}}})\,f$ over polynomials in $g_{i}$ and $C_{i}$; for these integrations, $C_{i}$ must be expressed through $g_{i}$ and $G_{i}$. These integrals were evaluated on the computer, using Wolfram Mathematica^TM. Results are inserted and simplified to produce the final expressions for the production terms.

5.7.1. Density, $\psi (C_{k}) =m$

As we have seen already, there is no contribution to mass balance:

(5.24)\begin{equation} \mathcal{P}_{m}=\mathcal{P}^{0}=0. \end{equation}

5.7.2. Momentum, $\psi (C_{k}) =mC_{i}$

The production terms for momentum and energy do not vanish, since in the collision momentum and energy are transferred across the diameter, $a$, that is, the distance between the centres of mass of the colliding hard spheres. However, for collision invariants, non-local production terms can be recast in divergence form (Karkheck & Stell Reference Karkheck and Stell1981). This property is inherited by the local expansion at any order, of course.

We take a closer look at the computation of the momentum production term, but for the other productions we will just present the results. After insertion of the integrals $[ I_{m,n}^{\alpha }] _{i_{1}\cdots i_{n+m}}$, the momentum production term reads

(5.25)\begin{align} \mathcal{P}_{mC_{i}} &= \mathcal{P}_{i}^{0}=\frac{\rm \pi}{2}ma^{2}Y\int \int f^{1}f\left\{ gg_{i}\right\} {\rm d}\boldsymbol{G}\,{\rm d}\boldsymbol{g} \nonumber\\ &\quad -\frac{2{\rm \pi}}{15}ma^{3}\int\int\left( Y_{r}f^{1}f+Y\, \frac{\partial f^{1}}{\partial x_{r}}\,f\right) \left\{ g^{2}\delta_{ir}+2g_{i}g_{r}\right\} {\rm d}\boldsymbol{G}\,{\rm d}\boldsymbol{g} \nonumber\\ &\quad +\frac{\rm \pi}{12}ma^{4}\int\int\left( Y_{rs}f^{1}f+Y_{r}\, \frac{\partial f^{1}}{\partial x_{s}}\,f+\frac{1}{2}Y\, \frac{\partial^{2}f^{1}}{\partial x_{r}\,\partial x_{s}}f\right) \left\{3g\delta_{(ir}g_{s)}+\frac{g_{i}g_{r}g_{s}}{g}\right\} {\rm d}\boldsymbol{G}\,{\rm d}\boldsymbol{g} \nonumber\\ &\quad -\frac{2{\rm \pi}}{35}ma^{5}\int\int\left(\begin{array}[c]{c} Y_{rst}f^{1}f+Y_{rs}\,\dfrac{\partial f^{1}}{\partial x_{t}}\,f\\ {}+\dfrac{1}{2}Y_{r}\,\dfrac{\partial^{2}f^{1}}{\partial x_{s}\,\partial x_{t}} f+ \dfrac{1}{6}Y\,\dfrac{\partial^{3}f^{1}}{\partial x_{r}\,\partial x_{s}\,\partial x_{t}}\,f \end{array}\right) \nonumber\\ &\quad \times\left\{ g^{2}\delta_{(ir}\delta_{st)}+4\delta_{(ir}g_{s} g_{t)}\right\} {\rm d}\boldsymbol{G}\,{\rm d}\boldsymbol{g}. \end{align}

Obviously, there is a larger number of integrals over $f^{1}f$ and $({\partial ^{n}f^{1}}/{\partial x^{n}})\,f$ that must be determined and inserted; see Appendix B for a brief discussion and some results. Inserting the integrals, and making use of (5.5a–d), we finally find the momentum collision term, as expected in divergence form:

(5.26)\begin{align} \mathcal{P}_{mC_{i}}={-}\frac{\partial}{\partial x_{r}}\left\{\begin{array}[c]{c} \left[\dfrac{2{\rm \pi}}{3}\,\dfrac{\rho a^{3}}{m}Y\left( \rho\theta\delta_{ir} +\dfrac{2}{5}\sigma_{ir}\right) \right] \\ {}-\left[\dfrac{2{\rm \pi}}{3}\dfrac{a^{4}}{m\sqrt{{\rm \pi}\theta}}\,\rho^{2}Y\dfrac{6}{5} \left(\theta\dfrac{\partial\delta_{(ir}v_{s)}}{\partial x_{s}}+ \dfrac{1}{10}\,\dfrac{\partial\dfrac{\delta_{(ir}q_{s)}}{\rho}}{\partial x_{s}}+ \dfrac{1}{42}\,\dfrac{\partial\dfrac{m_{irs}}{\rho}}{\partial x_{s}}\right) \right] \\ {}+\left[\dfrac{\rm \pi}{60}\,\dfrac{a^{5}}{m}\left(\dfrac{\partial^{2}Y\rho ^{2}\theta}{\partial x_{s}\,\partial x_{t}}+3Y\, \dfrac{\partial^{2}\rho^{2}\theta }{\partial x_{s}\,\partial x_{t}}-12Y\dfrac{\partial\rho}{\partial x_{s}} \dfrac{\partial\rho\theta}{\partial x_{t}}\right) \delta_{(ir}\delta_{st)}\right] \end{array}\right\}. \end{align}

We emphasize that to reach the above, we have ignored nonlinear terms between non-equilibrium contributions. Moreover, as stated before, we have kept only those fifth-order terms that contribute in equilibrium, where the density gradient is non-zero in the interface, while we removed fifth-order terms that contribute only in non-equilibrium, assuming that their influence is decidedly smaller.

5.7.3. 2-tensor, $\psi ( C_{k}) =mC_{i}C_{j}$

The Enskog contribution to energy is half of the trace, $\psi (C_{k}) =\frac {1}{2}mC^{2}$, and the contribution to the stress tensor is the trace-free part, $\psi ( C_{k}) =mC_{\langle i}C_{j\rangle }$; we list both contributions separately. Moreover, for more transparent presentation of the results, we show contributions of different orders in the expansion.

Order $a^{2}$

The leading terms correspond to results from the Boltzmann collision term for hard spheres, except that the correlation function $Y$ appears as a correction factor. There is no contribution to energy, only to stress:

(5.27)\begin{gather} \mathcal{P}_{({1}/{2})mC^{2}}^{\left( 2\right) } =0, \end{gather}

(5.28)\begin{gather}\mathcal{P}_{mC_{\langle i}C_{j\rangle}}^{\left( 2\right)} ={-}\frac {16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Y\left( \sigma_{ij}+\frac{1}{28}\,\frac{R_{ij}}{\theta}\right). \end{gather}

With the relaxation time $\tau$ defined through

(5.29)\begin{equation} \frac{1}{\tau}=\frac{16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Y, \end{equation}

we can write the last expression as $\mathcal {P}_{mC_{\langle i}C_{j\rangle } }^{( 2) }=-({1}/{\tau })\left (\sigma _{ij}+\frac {1}{28} ({R_{ij}}/{\theta })\right )$. The same relaxation time will appear in other second-order expressions below.

Order $a^{3}$

To third order, there is a contribution to thermal energy, which is related to energy transfer over the distance $a$ of the colliding particles:

(5.30)\begin{gather} P_{({1}/{2})mC^{2}}^{\left( 3\right)} ={-}\frac{2{\rm \pi}}{3}\,\frac{a^{3}} {m}\left\{ \rho^{2}\theta Y\,\frac{\partial v_{r}}{\partial x_{r}}+\frac{3} {5}\,\frac{\partial\rho Yq_{r}}{\partial x_{r}}\right\}, \end{gather}

(5.31)\begin{gather}P_{mC_{\langle i}C_{j\rangle}}^{\left( 3\right)} ={-}\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\left\{\begin{array}[c]{c} \dfrac{4}{5}\rho^{2}\theta Y\,\dfrac{\partial v_{\langle i}}{\partial x_{j\rangle }}+ \dfrac{12}{25}\,\dfrac{\partial\rho Yq_{\langle i}}{\partial x_{j\rangle}}+\dfrac{6}{25}\, \dfrac{q_{\langle i}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{j\rangle}}\\ {}+\dfrac{6}{35}\,\dfrac{\partial\rho Ym_{ijr}}{\partial x_{r}}+\dfrac{8}{35}\, \dfrac{m_{ijr}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{r}}. \end{array}\right\}. \end{gather}

The contribution to total energy results from adding the corresponding term for kinetic energy; see (4.11). Within the linearizations performed, the contribution to total energy can be written in divergence form as

(5.32)\begin{equation} \mathcal{P}_{({1}/{2})mC^{2}}^{\left( 3\right) }+v_{i}\mathcal{P}_{mC_{i} }^{\left( 3\right) }={-} \frac{\partial}{\partial x_{r}}\left\{ \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\left[ \rho\theta v_{r}+ \frac{3}{5}q_{r}+\frac{2}{5}v_{i}\sigma_{ir}\right]\right\}. \end{equation}

This form is shown to give evidence that the Enskog collision term leads not to an actual energy production, but rather to energy transfer, as described by the flux under the divergence.

Order $a^{4}$

The contributions to energy and stress equations to fourth order read

(5.33)\begin{align} & \mathcal{P}_{({1}/{2})mC^{2}}^{\left( 4\right)} \nonumber\\ &\quad =\frac{\partial}{\partial x_{r}} \left\{\frac{2{\rm \pi}}{15}\,\frac{a^{4}}{m\sqrt{\rm \pi}}\,\rho^{2}Y\left[5\sqrt{\theta}\, \frac{\partial\theta}{\partial x_{r}}+\frac{1}{12}\,\frac{1}{\sqrt{\theta}}\, \frac{\partial}{\partial x_{r}}\left(\dfrac{\varDelta}{\rho}\right) +2\sqrt{\theta}\,\frac{\partial}{\partial x_{s}}\left(\dfrac{\sigma_{rs}}{\rho}\right) +\frac{1}{7}\,\frac{1}{\sqrt{\theta}}\,\frac{\partial}{\partial x_{s}}\left(\dfrac{R_{rs}}{\rho}\right)\right]\right\}, \end{align}

(5.34)\begin{gather} \mathcal{P}_{mC_{\langle i}C_{j\rangle}}^{\left(4\right)}=\frac{4}{105}\,\frac{\sqrt{\rm \pi}a^{4}}{m}\,\frac{1}{\sqrt{\theta}} \left(\begin{array}[c]{c} 14\theta\,\dfrac{\partial}{\partial x_{\langle i}}\left( \rho^{2}Y\,\dfrac{\partial\theta}{\partial x_{j\rangle}}\right) + \dfrac{7}{30}\,\dfrac{\partial }{\partial x_{\langle i}} \left(\rho^{2}Y\,\dfrac{\partial}{\partial x_{j\rangle}}\left(\dfrac{\varDelta}{\rho}\right)\right) \\ {}+2\theta\rho\,\dfrac{\partial}{\partial x_{r}}\left(Y\,\dfrac{\partial\sigma_{ij} }{\partial x_{r}}\right) -12\theta\sigma_{ij}\,\dfrac{\partial}{\partial x_{r} }\left(Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) -\dfrac{5}{2}\rho \theta\sigma_{ij}\,\dfrac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\\ {}+8\theta\rho\,\dfrac{\partial}{\partial x_{(r}}\left( Y\,\dfrac{\partial \sigma_{r\langle i}\delta_{j\rangle s}}{\partial x_{s)}}\right) -20\theta\sigma_{r\langle i}\delta_{j\rangle s}\, \dfrac{\partial}{\partial x_{(r}}\left( Y\,\dfrac{\partial\rho}{\partial x_{s)}}\right) \\ {}-3\rho\theta\sigma_{r\langle i}\delta_{j\rangle s}\,\dfrac{\partial^{2} Y}{\partial x_{r}\,\partial x_{s}}\\ {}+\dfrac{1}{7}\rho\,\dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial R_{ij} }{\partial x_{r}}\right) -\dfrac{1}{2}R_{ij}\,\dfrac{\partial}{\partial x_{r} }\left(Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) -\dfrac{5}{56}\rho R_{ij}\,\dfrac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\\ {}+\dfrac{4}{7}\rho\,\dfrac{\partial}{\partial x_{(r}}\left( Y\,\dfrac{\partial R_{r\langle i}\delta_{j\rangle s}}{\partial x_{s)}}\right) -R_{r\langle i}\delta_{j\rangle s}\,\dfrac{\partial}{\partial x_{(r}}\left( Y\,\dfrac {\partial\rho}{\partial x_{s)}}\right) \\ {}-\dfrac{3}{28}\rho R_{r\langle i}\delta_{j\rangle s}\,\dfrac{\partial^{2} Y}{\partial x_{r}\,\partial x_{s}} \end{array}\right). \end{gather}

Within the linearizations performed, we find the fourth-order contribution to energy as $\mathcal {P}_{({1}/{2})mC^{2}}^{( 4) }+v_{i} \mathcal {P}_{mC_{i}}^{( 4) }\simeq \mathcal {P}_{({1}/{2})mC^{2} }^{( 4) }$, which is of divergence form. Obviously, the linearization removed terms that bring the full expression to divergence form.

All fifth-order terms vanish in equilibrium, and are ignored.

5.7.4. 3-tensor, $\psi ( C_{k}) =mC_{i}C_{j}C_{k}$

For the Enskog contribution to energy flux, we need half of the trace, $\psi ( C_{k}) =\frac {1}{2}mC^{2}C_{i}$, and for the contribution to the third moment we need the trace-free part, $\psi ( C_{k}) =mC_{\langle i}C_{j}C_{k\rangle }$; we list both contributions separately. The overall contribution to the heat flux equation (3.11) is the combination

(5.35)\begin{equation} \mathcal{P}_{i}^{HF}=\frac{1}{2}\mathcal{P}_{i}^{1}-\frac{5}{2}\theta \mathcal{P}_{i}^{0}-\frac{\sigma_{ik}}{\rho}\mathcal{P}_{k}^{0}= \mathcal{P}_{({1}/{2})mC^{2}C_{i}}-\frac{5}{2}\theta \mathcal{P}_{mC_{i}}-\frac{\sigma_{ik}}{\rho}\mathcal{P}_{mC_{k}}, \end{equation}

and for the contribution to the equations for $m_{ijk}$ we need the combination $\mathcal {P}_{ijk}^{0}-{3\sigma _{\langle ij}}/{\rho }\mathcal {P}_{k\rangle }^{0}$. These will be given as well.

Again, for more transparent presentation of the results, we show contributions of different orders in the expansion.

Order $a^{2}$

The leading terms correspond to results from the Boltzmann collision term for hard spheres, except that the correlation $Y$ function appears as a correction factor:

(5.36)\begin{gather} \mathcal{P}_{({1}/{2})mC^{2}C_{i}}^{\left( 2\right) } ={-}\frac{2} {3}\,\frac{16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Yq_{i}, \end{gather}

(5.37)\begin{gather}\mathcal{P}_{mC_{\langle i}C_{j}C_{k\rangle}}^{\left( 2\right)} ={-}\frac{3}{2}\,\frac{16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Ym_{ijk}. \end{gather}

With the relaxation time $\tau$ introduced for the stress contribution, we can write these as $\mathcal {P}_{({1}/{2})mC^{2}C_{i}}^{( 2) }=-\frac {2}{3}({q_{i}}/{\tau })$ and $\mathcal {P}_{mC_{\langle i} C_{j}C_{k\rangle }}^{( 2) }=-\frac {3}{2}({1}/{\tau })m_{ijk}$; these have the same factors $\frac {2}{3}$ and $\frac {3}{2}$ that one finds for hard spheres and Maxwell molecules; see Gupta & Torrilhon (Reference Gupta and Torrilhon2012).

Order $a^{3}$

We find

(5.38)\begin{align} & \mathcal{P}_{({1}/{2})mC^{2}C_{i}}^{\left( 3\right) } \nonumber\\ &\quad ={-}\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\left\{\begin{array}[c]{c} 4\theta\,\dfrac{\partial Y\rho^{2}\theta}{\partial x_{k}}-\dfrac{3}{2}\theta^{2}\, \dfrac{\partial\rho^{2}Y}{\partial x_{k}}+\dfrac{1}{10}\, \dfrac{\partial\rho Y\varDelta}{\partial x_{k}}+\dfrac{1}{30}\, \dfrac{\varDelta}{\rho}\,\dfrac{\partial\rho ^{2}Y}{\partial x_{k}}\\ {}+\dfrac{8}{5}\,\dfrac{\partial\rho Y\theta\sigma_{kr}}{\partial x_{r}}+ \dfrac{1}{10}\,\dfrac{\theta\sigma_{kr}}{\rho}\, \dfrac{\partial\rho^{2}Y}{\partial x_{r} }+\dfrac{6}{35}\, \dfrac{\partial\rho YR_{kr}}{\partial x_{r}}+\dfrac{11}{140}\, \dfrac{R_{kr}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{r}} \end{array}\right\}, \end{align}

(5.39)\begin{gather} \mathcal{P}_{mC_{\langle i}C_{j}C_{k\rangle}}^{\left( 3\right) }={-}\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\left[\begin{array}[c]{c} \dfrac{18}{35}\,\dfrac{\partial\rho Y\theta\sigma_{\langle ij}}{\partial x_{k\rangle}}+ \dfrac{9}{5}\,\dfrac{\theta\sigma_{\langle ij}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{k\rangle}}\\ {}+\dfrac{36}{245}\,\dfrac{\partial\rho YR_{\langle ij}}{\partial x_{k\rangle} }+ \dfrac{93}{490}\,\dfrac{R_{\langle ij}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{k\rangle}} \end{array}\right]. \end{gather}

The corresponding production term in the heat flux equation (3.11) reads, after appropriate linearization,

(5.40)\begin{equation} \mathcal{P}_{i}^{HF\left( 3\right) }={-}\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\left[\begin{array}[c]{c} \dfrac{3}{2}\rho^{2}Y\theta\,\dfrac{\partial\theta}{\partial x_{i}}+ \dfrac{1}{10}\,\dfrac{\partial\rho Y\varDelta}{\partial x_{i}}+\dfrac{1}{30}\, \dfrac{\varDelta}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{i}}\\ {}+\dfrac{3}{5}\,\dfrac{\partial\rho Y\theta\sigma_{ir}}{\partial x_{r}}-\dfrac{9}{10}\, \dfrac{\sigma_{ik}}{\rho}\,\dfrac{\partial\rho^{2}\theta Y}{\partial x_{k} }\\ {}+\dfrac{6}{35}\,\dfrac{\partial\rho YR_{ir}}{\partial x_{r}}+\dfrac{11}{140}\, \dfrac{R_{ir}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{r}} \end{array}\right]. \end{equation}

The contribution for the $m_{ijk}$ equation (3.14) becomes

(5.41)\begin{equation} \mathcal{P}_{ijk}^{0}-\frac{3\sigma_{\langle ij}}{\rho}\,\mathcal{P}_{k\rangle }^{0}={-}\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\left[\begin{array}[c]{c} \dfrac{18}{35}\theta\,\dfrac{\partial\rho Y\sigma_{\langle ij}}{\partial x_{k\rangle}}-\dfrac{6}{5}\, \dfrac{\theta\sigma_{\langle ij}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{k\rangle}}+ \dfrac{36}{245}\,\dfrac{\partial\rho YR_{\langle ij}}{\partial x_{k\rangle}}\\ {}+\dfrac{93}{490}\,\dfrac{R_{\langle ij}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{k\rangle}} \end{array}\right]. \end{equation}

Order $a^{4}$

For the fourth-order contribution to heat flux, we find the full production as

(5.42)\begin{equation} \mathcal{P}_{k}^{HF\left( 4\right) }=\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta }\left[\begin{array}[c]{c} -\dfrac{19}{15}\theta\left[\dfrac{\partial}{\partial x_{r}}\left( \rho^{2}Y\, \dfrac{\partial v_{k}}{\partial x_{r}}\right) +2\dfrac{\partial}{\partial x_{(k}} \left( \rho^{2}Y\,\dfrac{\partial v_{r}}{\partial x_{r)}}\right) \right]\\ {}-\dfrac{53}{70}\rho\,\dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial q_{k}}{\partial x_{r}}\right) -\dfrac{1}{2}q_{k}\,\dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) -\dfrac{2}{25}\rho q_{k}\dfrac{\partial^{2}Y}{\partial x_{r}\partial x_{r}}\\ {}-\dfrac{7}{15}q_{r}\,\dfrac{\partial}{\partial x_{(k}}\left( Y\,\dfrac{\partial\rho }{\partial x_{r)}}\right) +\dfrac{169}{105}\,\rho\,\dfrac{\partial}{\partial x_{(k}}\left( Y\,\dfrac{\partial q_{r}}{\partial x_{r)}}\right) -\dfrac{2}{75}\rho q_{r}\,\dfrac{\partial^{2}Y}{\partial x_{r}\,\partial x_{k}}\\ {}+\dfrac{9}{70}\,\rho\,\dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial m_{krs}}{\partial x_{s}}\right) -\dfrac{1}{6}m_{irs}\,\dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial\rho}{\partial x_{s}}\right) -\rho\,\dfrac{1}{105}\,m_{krs}\,\dfrac{\partial^{2}Y}{\partial x_{r}\partial x_{s}} \end{array}\right]. \end{equation}

The total contribution to the fourth-order terms for the $m_{ijk}$ equation (3.14) reads

(5.43)\begin{equation} \mathcal{P}_{m_{ijk}}=\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[\begin{array}[c]{c} \dfrac{4}{35}\theta\,\dfrac{\partial}{\partial x_{\langle i}} \left( \rho^{2}Y\,\dfrac{\partial v_{j}}{\partial x_{k\rangle}}\right) \\ {}-\dfrac{22}{35}\,\rho\,\dfrac{\partial}{\partial x_{\langle i}} \left(Y\,\dfrac{\partial q_{j}}{\partial x_{k\rangle}}\right) -\dfrac{2}{5}q_{\langle i}\,\dfrac{\partial}{\partial x_{j}} \left( Y\,\dfrac{\partial\rho}{\partial x_{k\rangle}}\right) -\dfrac{4}{175}\,\rho q_{\langle i}\,\dfrac{\partial^{2} Y}{\partial x_{j}\,\partial x_{k\rangle}}\\ {}+\dfrac{2}{70}\rho\,\dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial m_{ijk}}{\partial x_{r}}\right) -\dfrac{3}{5}m_{ijk}\,\dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) -\dfrac{1}{7}\rho m_{ijk}\,\dfrac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\\ {}+\dfrac{6}{70}\rho\,\dfrac{\partial}{\partial x_{r}} \left( Y\,\dfrac{\partial m_{r\langle ij}}{\partial x_{k\rangle}}\right) +\dfrac{6}{70}\rho\, \dfrac{\partial}{\partial x_{\langle i}}\left( Y\,\dfrac{\partial m_{jk\rangle r}}{\partial x_{r}}\right) -\dfrac{3}{7}m_{r\langle ij}\,\dfrac{\partial }{\partial x_{k\rangle}} \left( Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) \\ {}-\dfrac{3}{7}m_{r\langle ij}\,\dfrac{\partial}{\partial x_{r}} \left(Y\,\dfrac{\partial\rho}{\partial x_{k\rangle}}\right) -\dfrac{6}{35}\rho m_{r\langle ij}\,\dfrac{\partial^{2}Y}{\partial x_{k\rangle}\,\partial x_{r}} \end{array}\right]. \end{equation}

5.7.5. Trace of 4-tensor, $\psi ( C_{k}) =mC^{2}C_{i} C_{j}$

Again, we distinguish the trace and trace-free parts, at different orders. These are the contributions to equations (3.10) and (3.13) for higher moments $\varDelta$ and $R_{ij}$.

Order $a^{2}$

The leading terms correspond to results from the Boltzmann collision term for hard spheres, except that the correlation function appears as a correction factor:

(5.44)\begin{gather} \mathcal{P}_{mC^{4}}^{\left( 2\right) } ={-}\frac{2}{3}\,\frac{16}{5}\, \frac{a^{2}\rho\sqrt{{\rm \pi}\theta}}{m}\,Y\varDelta, \end{gather}

(5.45)\begin{gather}\mathcal{P}_{mC^{2}C_{\langle i}C_{j\rangle}}^{\left( 2\right) } ={-}\frac{16}{5}\,\frac{a^{2}\rho\sqrt{{\rm \pi}\theta}}{m}\,Y\left[ \frac{15}{2} \theta\sigma_{ij}+\frac{247}{168}R_{ij}\right]. \end{gather}

Once more, we find the hard-sphere relaxation time ${1}/{\tau }=\frac {16}{5}({\rho \sqrt {{\rm \pi} \theta }a^{2}}/{m})Y$. The coefficients ($\frac {2}{3},\frac {15}{2},\frac {247}{168}$) agree with the results for hard spheres by Gupta & Torrilhon (Reference Gupta and Torrilhon2012).

Order $a^{3}$

The third-order contributions are

(5.46)\begin{gather} \mathcal{P}_{mC^{4}}^{\left( 3\right) }={-}\frac{2}{3}\,\frac{{\rm \pi}\theta a^{3} }{m} \left\{\frac{84}{5}\,\frac{\partial Y\rho q_{r}}{\partial x_{r}}+ \frac{8}{5}\,\frac{q_{r}}{\rho}\,\frac{\partial Y\rho^{2}}{\partial x_{r}} +20Y\theta\rho^{2}\,\frac{\partial v_{r}}{\partial x_{r}}\right\}, \end{gather}

(5.47)\begin{gather}\mathcal{P}_{mC^{2}C_{\langle i}C_{j\rangle}}^{\left( 3\right) }={-}\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\theta\left\{\begin{array}[c]{@{}c@{}} \dfrac{132}{25}\,\dfrac{\partial\rho Yq_{\langle i}}{\partial x_{j\rangle}} +\dfrac{112}{25}\,\dfrac{q_{\langle i}}{\rho}\,\dfrac{\partial\rho^{2}Y}{\partial x_{j\rangle}}+ \dfrac{28}{5}Y\theta\rho^{2}\,\dfrac{\partial v_{\langle i} }{\partial x_{j\rangle}}\\ {}+\dfrac{66}{35}\,\dfrac{\partial Y\rho m_{ijr}}{\partial x_{r}}+\dfrac{71}{35}\, \dfrac{m_{ijr}}{\rho}\,\dfrac{\partial Y\rho^{2}}{\partial x_{r}} \end{array}\right\}. \end{gather}

For the balance laws (3.10), (3.13) for $\varDelta$ and $R_{ij}$, we find the required combinations up to third order, after linearization:

(5.48)\begin{align} \mathcal{P}^{2}-20\theta\,\frac{\mathcal{P}^{1}}{2}-\frac{8q_{k}}{\rho }\,\mathcal{P}_{k}^{0} &={-}\frac{2}{3}\left[ \frac{16}{5}\,\frac{a^{2}\rho \sqrt{{\rm \pi}\theta}}{m}\,Y\right]\varDelta \nonumber\\ &\quad -\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\, \theta\left[ \frac{24}{5}\,\frac{\partial\rho Yq_{r}}{\partial x_{r}}-\frac{32}{5}\, \frac{q_{k}}{\rho}\,\frac{\partial Y\rho^{2}}{\partial x_{k}}\right], \end{align}

(5.49)\begin{align} &\mathcal{P}_{ij}^{1}-7\theta\mathcal{P}_{ij}^{0}-\frac{14}{3}\, \frac{\sigma_{ij}}{\rho}\,\frac{\mathcal{P}^{1}}{2}-2\,\frac{m_{ijk}}{\rho}\,\mathcal{P}_{k}^{0}-\frac{28}{5}\, \frac{q_{\langle i}}{\rho}\,\mathcal{P}_{j\rangle}^{0}\nonumber\\ &\quad ={-}\frac{16}{5}\,\frac{a^{2}\rho\sqrt{{\rm \pi}\theta}}{m}\,Y \left[\frac{1}{2}\theta\sigma_{ij}+\frac{205}{168}R_{ij}\right] -\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\theta\left[ \frac{48}{25}\, \frac{\partial\rho Yq_{\langle i}}{\partial x_{j\rangle}}\right.\nonumber\\ &\qquad \, -\left.\frac{14}{5}\,\frac{q_{\langle i}}{\rho}\, \frac{\partial\rho^{2}Y}{\partial x_{j\rangle}}+\frac{24}{35}\, \frac{\partial\rho Ym_{ijr}}{\partial x_{r}}-\frac{11}{7}\,\frac{m_{ijk}}{\rho}\, \frac{\partial\rho^{2}Y}{\partial x_{k}}\right]. \end{align}

The fourth-order contributions are not shown separately, but are included in the full equations shown below.

6.1. Mass balance

The mass balance is not affected by the Enskog collision term; it has the usual form,

(6.1)\begin{equation} \frac{\partial\rho}{\partial t}+\frac{\partial\rho v_{k}}{\partial x_{k}}=0. \end{equation}

6.2. Momentum balance

The momentum balance is best written in conservative form, with the total momentum flux in the divergence. For order $a^{5}$, we consider only equilibrium contributions:

(6.2)\begin{equation} \rho\frac{Dv_{i}}{Dt}+\frac{\partial}{\partial x_{k}}\left(\begin{array}[c]{c} \left[ \rho\theta\left( 1+\dfrac{2{\rm \pi}}{3}\,\dfrac{a^{3}}{m}\,\rho Y\right) -\dfrac{2{\rm \pi} a^{3}}{3}\,\dfrac{\phi_{a}}{m}\,\chi_{1}\rho^{2}\right] \delta_{ik}\\ {}+\left[ 1+\dfrac{2}{5}\,\dfrac{2{\rm \pi}}{3}\,\dfrac{a^{3}}{m}\,\rho Y\right] \sigma_{ik}\\ {}-\dfrac{4}{5}\,\dfrac{a^{4}\sqrt{\rm \pi}}{m}\,\dfrac{\rho^{2}Y}{\sqrt{\theta}} \left(\theta\,\dfrac{\partial\delta_{(ik}v_{s)}}{\partial x_{s}}+\dfrac{1}{10}\, \dfrac{\partial\dfrac{\delta_{(ik}q_{s)}}{\rho}}{\partial x_{s}}+ \dfrac{1}{42}\,\dfrac{\partial\dfrac{m_{iks}}{\rho}}{\partial x_{s}}\right) \\ {}+\dfrac{\rm \pi}{60}\,\dfrac{a^{5}}{m}\left(\dfrac{\partial^{2}Y\rho^{2}\theta }{\partial x_{s}\,\partial x_{t}}+3Y\, \dfrac{\partial^{2}\rho^{2}\theta}{\partial x_{s}\,\partial x_{t}}-12Y\, \dfrac{\partial\rho}{\partial x_{s}}\,\dfrac{\partial \rho\theta}{\partial x_{t}}\right) \delta_{(ik}\delta_{st)}\\ {}-\dfrac{2{\rm \pi}}{15}\,\dfrac{a^{5}}{m}\,\phi_{a}\chi_{3}\left[\left( \rho\, \dfrac{\partial^{2}\rho}{\partial x_{r}\,\partial x_{r}}+\dfrac{1}{2}\, \dfrac{\partial\rho}{\partial x_{r}}\,\dfrac{\partial\rho}{\partial x_{r}}\right) \delta_{ik}-\dfrac{\partial\rho}{\partial x_{i}}\,\dfrac{\partial\rho}{\partial x_{k}}\right] \end{array}\right) =\rho G_{i}. \end{equation}

This equation has the usual form of the momentum balance,

(6.3)\begin{equation} \rho\,\frac{{\rm D} v_{i}}{{\rm D} t}+\frac{\partial\left( p\delta_{ik}+\varPi_{ik}\right) }{\partial x_{k}}= \rho G_{i}, \end{equation}

with a generalized stress tensor $p\delta _{ik}+\varPi _{ik}$ that has many contributions that can be read from (6.2). The first line describes the van-der-Waals-like pressure

(6.4)\begin{equation} p=\rho\theta\left( 1+\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) -\frac{2{\rm \pi} a^{3}}{3}\,\frac{\phi_{a}}{m}\,\chi_{1}\rho^{2}, \end{equation}

with a volume correction stemming from the Enskog term, and the pressure-reducing Vlasov contribution. The remainder is the non-equilibrium stress $\varPi _{ik}$, where the second line is the contribution from kinetic stress $\sigma _{ij}$ with a correction due to the Enskog term. The last line is the fifth-order Korteweg stress, which is corrected by the fifth-order Enskog terms in the preceding line. The third line describes higher contributions to stress, including a Navier–Stokes-like contribution from the velocity gradient, a Burnett-like term with the gradient of the kinetic heat flux (i.e. thermal stress), and a higher contribution with $m_{ijk}$.

In equilibrium, we have $v_{i}=\sigma _{ik}=q_{k}=m_{ijk}=\varDelta =R_{ik}=0$. Then the equilibrium density field follows from the first, fourth and fifth lines of (6.2). The solution allows for smooth phase interfaces, where the Korteweg stress stabilizes the non-monotonic van-der-Waals-like pressure. Details of one-dimensional interfaces are presented by Frezzotti et al. (Reference Frezzotti, Gibelli and Lorenzani2005) .

6.3. Balance of thermal energy

The balance of thermal energy includes terms up to fourth order, with frictional heating:

(6.5)\begin{align} & \frac{3}{2}\rho\,\frac{D\theta}{Dt}+\frac{\partial}{\partial x_{k}}\left[\begin{array}[c]{c} \left(1+\dfrac{3}{5}\,\dfrac{2{\rm \pi}}{3}\,\dfrac{a^{3}}{m}\,\rho Y\right) q_{k}\\ {}-\dfrac{a^{4}\sqrt{\rm \pi}}{m}\,\dfrac{\rho^{2}Y}{\sqrt{\theta}}\left(\dfrac{2}{3}\theta\, \dfrac{\partial\theta}{\partial x_{k}}+\dfrac{4}{15}\theta\, \dfrac{\partial\dfrac{\sigma_{ks}}{\rho}}{\partial x_{s}}+\dfrac{1}{90}\, \dfrac{\partial\dfrac{\varDelta}{\rho}}{\partial x_{k}}+\dfrac{2}{105}\, \dfrac{\partial\dfrac{R_{ks}}{\rho}}{\partial x_{s}}\right) \end{array}\right] \nonumber\\ &\quad ={-}\left( \varPi_{kl}-\mathcal{T}_{rs}^{K}\right) \frac{\partial v_{k} }{\partial x_{l}}. \end{align}

This has the usual form of the balance of internal energy,

(6.6)\begin{equation} \frac{3}{2}\,\rho\,\frac{D\theta}{Dt}+\frac{\partial Q_{k}}{\partial x_{k}} ={-}\hat{\varPi}_{kl}\,\frac{\partial v_{k}}{\partial x_{l}}, \end{equation}

with a generalized heat flux $Q_{k}$ that has several contributions that can be read from (6.5). The first line is the contribution from kinetic heat flux $q_{i}$ with a correction due to the Enskog term. The second line describes higher contributions to heat flux, including a Fourier-like contribution from the temperature gradient, a Burnett-like term with the kinetic stress tensor $\sigma _{ij}$ (non-Fourier heat flux), and higher-order contributions related to $\varDelta$ and $R_{ij}$.

The energy exchange term on the right-hand side contains only the non-Vlasov stress tensor contribution $\hat {\varPi }_{kl}=( \varPi _{kl}-\mathcal {T} _{rs}^{K})$ and not the full stress tensor $\varPi _{kl}$ that appears in the balance of momentum. As discussed earlier, the Vlasov stresses appear in the balance of Korteweg energy (see (4.8)), but do not contribute to frictional heating.

As written in (6.5), the right-hand side contains nonlinear terms related to friction, but in the linear approximation used throughout, these are ignored, and the energy exchange reduces to the work for volume change:

(6.7)\begin{equation} -\hat{\varPi}_{kl}\,\frac{\partial v_{k}}{\partial x_{l}}={-}\left( \varPi_{kl}-\mathcal{T}_{rs}^{K}\right) \frac{\partial v_{k}}{\partial x_{l}} \simeq{-}\rho\theta\left( 1+\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \frac{\partial v_{k}}{\partial x_{k}}. \end{equation}

The balance of total energy has balance law form, with the total energy flux under the divergence. Within the EV model, one has to account for thermal, kinetic and Vlasov energy, and the equation reads, in compact form,

(6.8)\begin{equation} \rho\frac{{\rm D}\left( \tfrac{3}{2}\theta+\tfrac{1}{2}v^{2}+\epsilon_{K}\right) }{{\rm D} t}+ \frac{\partial( Q_{k}+q_{k}^{K}+\varPi_{ik}v_{i}) }{\partial x_{k}}=\rho v_{i}G_{i}. \end{equation}

6.4. Stress balance

The linearized balance for the kinetic stress tensor includes a larger number of terms up to fourth order:

(6.9)\begin{align} & \frac{\partial\sigma_{ij}}{\partial t}+2\left[ \rho\theta \left( 1+\frac{2}{5}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \right] \frac{\partial v_{\langle i}}{\partial x_{j\rangle}}+\frac{4}{5}\, \frac{\partial}{\partial x_{\langle i}}\left[ \left( 1+\frac{3}{5}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) q_{j\rangle}\right] \nonumber\\ &\qquad +\frac{\partial}{\partial x_{k}}\left[ \left( 1+\frac{6}{35}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) m_{ijk}\right] +\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m} \left[ \frac{6}{25}\,q_{\langle i}\delta_{j\rangle r}+\frac{8}{35}\,m_{ijr}\right] \frac{1}{\rho}\,\frac{\partial\rho^{2}Y}{\partial x_{r}} \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{8}{15}\rho^{2}\, \frac{\partial}{\partial x_{\langle i}}\left( Y\,\frac{\partial\theta}{\partial x_{j\rangle}}\right) +\frac{16}{15}\rho Y\,\frac{\partial\rho}{\partial x_{\langle i}}\, \frac{\partial\theta}{\partial x_{j\rangle}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m\sqrt{\theta}}\left[ \frac{2}{225}\rho\, \frac{\partial}{\partial x_{\langle i}}\left( Y\,\frac{\partial\varDelta}{\partial x_{j\rangle}}\right) -\frac{2}{225}\varDelta\,\frac{\partial}{\partial x_{\langle i}} \left( Y\,\frac{\partial\rho}{\partial x_{j\rangle}}\right) \right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{8}{105}\rho\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\sigma_{ij}}{\partial x_{r}}\right) -\frac{16}{35}\sigma_{ij}\,\frac{\partial}{\partial x_{r}}\left(Y\,\frac{\partial\rho}{\partial x_{r}}\right) -\frac{2}{21}\rho\sigma_{ij}\,\frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[\begin{array}[c]{c} \dfrac{16}{105}\rho\,\dfrac{\partial}{\partial x_{\langle i}} \left(Y\,\dfrac{\partial\sigma_{j\rangle r}}{\partial x_{r}}\right) +\dfrac{16}{105}\rho\, \dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial\sigma_{r\langle i}}{\partial x_{j\rangle}}\right) \\ {}-\dfrac{8}{21}\sigma_{r\langle i}\,\dfrac{\partial}{\partial x_{j\rangle}} \left(Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) -\dfrac{8}{21}\sigma_{r\langle i}\,\dfrac{\partial}{\partial x_{r}} \left(Y\,\dfrac{\partial\rho}{\partial x_{j\rangle}}\right) -\dfrac{4}{35}\rho \sigma_{r\langle i}\,\dfrac{\partial^{2}Y}{\partial x_{j\rangle}\,\partial x_{r}} \end{array}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m\sqrt{\theta}}\left[\frac{4}{735}\rho\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial R_{ij}}{\partial x_{r}}\right) -\frac{2}{105}R_{ij}\,\frac{\partial}{\partial x_{r}} \left(Y\,\frac{\partial\rho}{\partial x_{r}}\right) -\frac{1}{294}\rho R_{ij}\, \frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m\sqrt{\theta}}\left[\begin{array}[c]{c} \dfrac{8}{735}\rho\,\dfrac{\partial}{\partial x_{\langle i}} \left(Y\,\dfrac{\partial R_{j\rangle r}}{\partial x_{r}}\right) +\dfrac{8}{735}\rho\, \dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial R_{r\langle i}}{\partial x_{j\rangle}}\right) \\ {}-\dfrac{2}{105}R_{r\langle i}\,\dfrac{\partial}{\partial x_{r}} \left(Y\,\dfrac{\partial\rho}{\partial x_{j\rangle}}\right) -\dfrac{2}{105}R_{r\langle i}\,\dfrac{\partial}{\partial x_{j\rangle}} \left(Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) -\dfrac{1}{245}\rho R_{r\langle i}\, \dfrac{\partial^{2}Y}{\partial x_{j\rangle}\,\partial x_{r}} \end{array}\right] \nonumber\\ &\quad ={-}\left[\frac{16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Y\right] \left(\sigma_{ij}+\frac{1}{28}\,\frac{R_{ij}}{\theta}\right). \end{align}

6.5. Heat flux balance

The linearized balance for kinetic energy flux also includes many terms up to fourth order:

(6.10)\begin{align} & \frac{\partial q_{i}}{\partial t}+\frac{5}{2}\rho\theta\left( 1+\frac{3}{5}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \frac{\partial\theta}{\partial x_{i}}+ \theta\,\frac{\partial}{\partial x_{k}}\left[ \left( 1+\frac{3}{5}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \sigma_{ik}\right] \nonumber\\ &\qquad +\frac{1}{2}\,\frac{\partial}{\partial x_{k}}\left[ \left( 1+\frac{12}{35}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) R_{ik}\right] +\frac{1}{6}\, \frac{\partial}{\partial x_{i}}\left[ \left( 1+\frac{3}{5}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) \varDelta\right] \nonumber\\ &\qquad -\frac{\sigma_{ik}}{\rho}\,\theta\,\frac{\partial}{\partial x_{k}} \left[ \left(1+\frac{9}{10}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \rho\right] +\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\left[ \frac{1}{30}\,\frac{\varDelta}{\rho}\, \delta_{ik}+\frac{11}{140}\,\frac{R_{ik}}{\rho}\right] \frac{\partial\rho^{2}Y}{\partial x_{k}} \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ -\frac{19}{15}\,\theta\, \frac{\partial}{\partial x_{r}}\left( \rho^{2}Y\,\frac{\partial v_{i}}{\partial x_{r}}\right) -\frac{38}{15}\theta\,\frac{\partial}{\partial x_{(i}}\left(\rho^{2}Y\, \frac{\partial v_{r}}{\partial x_{r)}}\right) \right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{9}{50}\,\rho\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial q_{i}}{\partial x_{r}}\right) -\frac{1}{2}q_{i}\,\frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\rho}{\partial x_{r}}\right) -\frac{2}{25}\rho q_{i}\,\frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{9}{25}\,\rho\, \frac{\partial}{\partial x_{(i}}\left( Y\,\frac{\partial q_{r}}{\partial x_{r)}}\right) -\frac{11}{15}q_{r}\,\frac{\partial}{\partial x_{(i}}\left( Y\,\frac{\partial\rho}{\partial x_{r)}}\right) -\frac{2}{75}\rho q_{r}\,\frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{i}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{9}{70}\,\rho\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial m_{irs}}{\partial x_{s}}\right) -\frac{1}{6}m_{irs}\,\frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\rho}{\partial x_{s}}\right) -\frac{1}{105}\rho m_{irs}\,\frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{s}}\right] \nonumber\\ &\quad ={-}\frac{2}{3}\left[ \frac{16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Y\right] q_{i}. \end{align}

6.6. $m_{ijk}$ balance

The linearized transport equation up to fourth order reads

(6.11)\begin{align} & \frac{\partial m_{ijk}}{\partial t}+3\theta\,\frac{\partial}{\partial x_{\langle i}} \left[\left( 1+\frac{6}{35}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \sigma_{jk\rangle}\right] -3\,\frac{\theta\sigma_{\langle ij}}{\rho}\, \frac{\partial}{\partial x_{k\rangle}}\left[ \rho+\frac{2}{5}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho^{2}Y\right] \nonumber\\ &\qquad +\frac{3}{7}\,\frac{\partial}{\partial x_{\langle i}}\left[ \left( 1+\frac{12}{35}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) R_{jk\rangle}\right] +\frac{93}{490}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\frac{R_{\langle ij}}{\rho}\, \frac{\partial\rho^{2}Y}{\partial x_{k\rangle}} \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{4}{35}\,\theta\, \frac{\partial}{\partial x_{\langle i}}\left( \rho^{2}Y\, \frac{\partial v_{j}}{\partial x_{k\rangle}}\right) \right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{54}{175}\,\rho\, \frac{\partial}{\partial x_{\langle i}}\left( Y\,\frac{\partial q_{j}}{\partial x_{k\rangle}}\right) -\frac{2}{5}q_{\langle i}\,\frac{\partial}{\partial x_{j}} \left( Y\,\frac{\partial\rho}{\partial x_{k\rangle}}\right) -\frac{4}{175}\rho q_{\langle i}\,\frac{\partial^{2}Y}{\partial x_{j}\,\partial x_{k\rangle}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{2}{70}\,\rho\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial m_{ijk}}{\partial x_{r}}\right) -\frac{3}{5}m_{ijk}\,\frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\rho}{\partial x_{r}}\right) -\frac{1}{7}\rho m_{ijk}\,\frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{6}{70}\,\rho\,\frac{\partial }{\partial x_{\langle i}}\left( Y\,\frac{\partial m_{jk\rangle r}}{\partial x_{r}}\right) -\frac{3}{7}m_{r\langle ij}\,\frac{\partial}{\partial x_{k\rangle}} \left( Y\,\frac{\partial\rho}{\partial x_{r}}\right) -\frac{3}{35}\rho m_{r\langle ij}\, \frac{\partial^{2}Y}{\partial x_{k\rangle}\,\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{6}{70}\,\rho\,\frac{\partial}{\partial x_{r}} \left( Y\,\frac{\partial m_{r\langle ij}}{\partial x_{k\rangle}}\right) -\frac{3}{7}m_{r\langle ij}\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\rho}{\partial x_{k\rangle}}\right) -\frac{3}{35}\rho m_{r\langle ij}\,\frac{\partial^{2}Y}{\partial x_{k\rangle}\,\partial x_{r}}\right] \nonumber\\ &\quad ={-}\frac{3}{2}\left[ \frac{16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Y\right] m_{ijk}. \end{align}

6.7. $\varDelta$ balance

The linearized transport equation up to fourth order reads

(6.12)\begin{align} & \frac{{\rm D}\varDelta}{{\rm D} t}+8\theta\,\frac{\partial}{\partial x_{k}} \left[ \left(1+\frac{3}{5}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) q_{k}\right] -8\theta\,\frac{q_{k}}{\rho}\,\frac{\partial}{\partial x_{k}} \left[ \left(1+\frac{4}{5}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \rho\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[ \frac{4}{3}\rho^{2}\theta\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\theta}{\partial x_{r}}\right) + \frac{8}{3}\theta\rho Y\,\frac{\partial\rho}{\partial x_{r}}\,\frac{\partial\theta}{\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[ \frac{8}{15}\rho^{2}\theta\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\frac{\sigma_{rs}}{\rho}}{\partial x_{s}}\right) -\frac{32}{3}\rho\theta Y\,\frac{\partial\rho}{\partial x_{r}}\, \frac{\partial\dfrac{\sigma_{rs}}{\rho}}{\partial x_{s}}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[ \frac{1}{3}\rho\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\varDelta}{\partial x_{r}}\right) -\frac{31}{45}\varDelta\,\frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\rho}{\partial x_{r}}\right) -\frac{4}{45}\varDelta\rho\,\frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[ \frac{4}{7}\rho\,\frac{\partial}{\partial x_{r}} \left( Y\,\frac{\partial R_{rs}}{\partial x_{s}}\right) -\frac{76}{105}R_{rs}\,\frac{\partial}{\partial x_{r}} \left(Y\,\frac{\partial\rho}{\partial x_{s}}\right) -\frac{4}{105}\rho R_{rs}\, \frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{s}}\right] \nonumber\\ &\quad ={-}\frac{2}{3}\left[ \frac{16}{5}\,\frac{a^{2}\rho\sqrt{{\rm \pi}\theta}}{m}\,Y\right]\varDelta. \end{align}

6.8. $R_{ij}$ balance

The linearized transport equation up to fourth order reads

(6.13)\begin{align} & \frac{\partial R_{ij}}{\partial t}+\frac{28}{5}\theta\, \frac{\partial}{\partial x_{\langle i}}\left[ \left( 1+\frac{12}{35}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) q_{j\rangle}\right] -\frac{28}{5}\theta\,\frac{q_{\langle i}}{\rho}\, \frac{\partial}{\partial x_{j\rangle}}\left[ \left( 1+\frac{1}{2}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \rho\right] \nonumber\\ &\qquad +2\theta\,\frac{\partial}{\partial x_{k}}\left[ \left( 1+\frac{12}{35}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) m_{ijk}\right] -2\,\frac{m_{ijk}\theta}{\rho}\, \frac{\partial}{\partial x_{k}}\left[ \left(1+\frac{11}{14}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) \rho\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{8}{15}\rho^{2}\theta\, \frac{\partial}{\partial x_{\langle i}}\left( Y\,\frac{\partial\theta}{\partial x_{j\rangle}}\right) + \frac{16}{15}\rho\theta Y\,\frac{\partial\rho}{\partial x_{\langle i}}\, \frac{\partial\theta}{\partial x_{j\rangle}}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[ \frac{2}{15}\rho\, \frac{\partial}{\partial x_{\langle i}}\left( Y\,\frac{\partial\varDelta}{\partial x_{j\rangle}}\right) -\frac{38}{225}\varDelta\,\frac{\partial}{\partial x_{\langle i}} \left(Y\,\frac{\partial\rho}{\partial x_{j\rangle}}\right) -\frac{2}{225}\rho \varDelta\,\frac{\partial^{2}Y}{\partial x_{\langle i}\,\partial x_{j\rangle}}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[ \frac{8}{105}\rho\theta\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial\sigma_{ij}}{\partial x_{r}}\right) -\frac{4}{15}\theta\sigma_{ij}\,\frac{\partial}{\partial x_{r}} \left( Y\,\frac{\partial\rho}{\partial x_{r}}\right) -\frac{1}{21}\rho \theta\sigma_{ij}\,\frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[\begin{array}[c]{c} \dfrac{16}{105}\rho\theta\,\dfrac{\partial}{\partial x_{\langle i}} \left(Y\,\dfrac{\partial\sigma_{j\rangle r}}{\partial x_{r}}\right) +\dfrac{16}{105}\rho\theta\, \dfrac{\partial}{\partial x_{r}}\left( Y\,\dfrac{\partial \sigma_{r\langle i}}{\partial x_{j\rangle}}\right) \\ {}-\dfrac{4}{15}\theta\sigma_{r\langle i}\,\dfrac{\partial}{\partial x_{j\rangle}} \left(Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) -\dfrac{4}{15}\theta\sigma_{r\langle i}\,\dfrac{\partial}{\partial x_{r}} \left(Y\,\dfrac{\partial\rho}{\partial x_{j\rangle}}\right) \\ {}-\dfrac{2}{35}\rho\theta\sigma_{r\langle i}\,\dfrac{\partial^{2}Y}{\partial x_{j\rangle}\,\partial x_{r}} \end{array}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[\frac{4}{49}\rho\, \frac{\partial}{\partial x_{r}}\left( Y\,\frac{\partial R_{ij}}{\partial x_{r}}\right) -\frac{149}{245}R_{ij}\,\frac{\partial}{\partial x_{r}} \left(Y\,\frac{\partial\rho}{\partial x_{r}}\right) -\frac{129}{980}\rho R_{ij}\, \frac{\partial^{2}Y}{\partial x_{r}\,\partial x_{r}}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[\begin{array}[c]{c} \dfrac{8}{49}\rho\,\dfrac{\partial}{\partial x_{\langle i}} \left(Y\,\dfrac{\partial R_{j\rangle r}}{\partial x_{r}}\right) + \dfrac{8}{49}\rho\,\dfrac{\partial}{\partial x_{r}} \left(Y\,\dfrac{\partial R_{r\langle i}}{\partial x_{j\rangle}}\right) \\ {}-\dfrac{257}{735}R_{r\langle i}\,\dfrac{\partial}{\partial x_{j\rangle}} \left(Y\,\dfrac{\partial\rho}{\partial x_{r}}\right) -\dfrac{257}{735}R_{r\langle i}\,\dfrac{\partial}{\partial x_{r}} \left(Y\,\dfrac{\partial\rho}{\partial x_{j\rangle}}\right) \\ {}-\dfrac{137}{1470}\rho R_{r\langle i}\,\dfrac{\partial^{2}Y}{\partial x_{j\rangle}\partial x_{r}} \end{array}\right] \nonumber\\ &\quad ={-}\left[ \frac{16}{5}\,\frac{a^{2}\rho\sqrt{{\rm \pi}\theta}}{m}\,Y\right] \left(\frac{1}{2}\theta\sigma_{ij}+\frac{205}{168}R_{ij}\right). \end{align}

6.9. Ideal gas limit

The above equations are valid, within the limitations of the underlying theory, for all states from compressed liquid over two-phase flows and real gas, to ideal gas. For the case of the latter, they agree with the Grad 26 moment equations as presented in Struchtrup (Reference Struchtrup2005) within the limits of the performed linearizations. Indeed, the ideal gas case of the above is obtained from (6.1)–(6.13) by removing all terms with powers $a^{3}$ and higher, and setting $Y=1$ (with $p=\rho \theta$, ${1}/{\tau }=\frac {16}{5}({\rho \sqrt {{\rm \pi} \theta }a^{2}}/{m})$), as

(6.14)\begin{gather} \left.\begin{gathered} \frac{\partial\rho}{\partial t}+\frac{\partial\rho v_{k}}{\partial x_{k}}=0, \\ \rho\,\frac{{\rm D} v_{i}}{{\rm D} t}+\frac{\partial}{\partial x_{k}}\left( \rho\theta \delta_{ik}+\sigma_{ik}\right) =\rho G_{i}, \\ \frac{3}{2}\rho\,\frac{{\rm D}\theta}{{\rm D} t}+\frac{\partial q_{k}}{\partial x_{k}} ={-}\rho\theta\,\frac{\partial v_{k}}{\partial x_{k}}, \\ \frac{\partial\sigma_{ij}}{\partial t}+2\rho\theta\,\frac{\partial v_{\langle i}}{\partial x_{j\rangle}}+ \frac{4}{5}\,\frac{\partial q_{\langle i}}{\partial x_{j\rangle}}+ \frac{\partial m_{ijk}}{\partial x_{k}} ={-}\frac{1}{\tau }\left( \sigma_{ij}+\frac{1}{28}\, \frac{R_{ij}}{\theta}\right), \\ \frac{\partial q_{i}}{\partial t}+\frac{5}{2}\rho\theta\,\frac{\partial\theta }{\partial x_{i}}+ \theta\,\frac{\partial\sigma_{ik}}{\partial x_{k} }-\underline{\theta\sigma_{ik}\, \frac{\partial\ln\rho}{\partial x_{k}}}+\frac{1}{2}\,\frac{\partial R_{ik}}{\partial x_{k}}+ \frac{1}{6}\,\frac{\partial\varDelta}{\partial x_{i}} ={-}\frac{2}{3}\,\frac{q_{i}}{\tau}, \\ \frac{\partial m_{ijk}}{\partial t}+3\theta\frac{\partial\sigma_{\langle ij} }{\partial x_{k\rangle}}+ \frac{3}{7}\frac{\partial R_{\langle ij}}{\partial x_{k\rangle}}-\underline{3\theta\sigma_{\langle ij}\, \frac{\partial\ln\rho }{\partial x_{k\rangle}}} ={-}\frac{3}{2}\,\frac{m_{ijk}}{\tau}, \\ \frac{\partial\varDelta}{\partial t}+8\theta\,\frac{\partial q_{k}}{\partial x_{k} }-\underline{8\theta q_{k}\, \frac{\partial\ln\rho}{\partial x_{k}}} ={-}\frac{2}{3}\,\frac{\varDelta}{\tau}, \\ \frac{\partial R_{ij}}{\partial t}+\frac{28}{5}\theta\, \frac{\partial q_{\langle i}}{\partial x_{j\rangle}}-\underline{\frac{28}{5}\theta q_{\langle i}\, \frac{\partial\ln\rho}{\partial x_{j\rangle}}}-\underline{2m_{ijk}\theta\, \frac{\partial\ln\rho}{\partial x_{k}}}+2\theta\,\frac{\partial m_{ijk} }{\partial x_{k}} ={-}\frac{1}{\tau}\left( \frac{1}{2}\theta\sigma_{ij}+\frac{205}{168}R_{ij}\right). \end{gathered}\right\} \end{gather}

Comparison with the equations in Chapter 9 of Struchtrup (Reference Struchtrup2005) shows that the above without the underlined terms are the linearized version of the Grad 26 equations for the ideal gas. The underlined terms are nonlinear with the density gradient and appear in the nonlinear Grad 26 equations, which, however, contain other nonlinear terms that were removed during the above derivation. The collision terms here are for a hard-sphere gas, while in Struchtrup (Reference Struchtrup2005) they are given for Maxwell molecules.

6.10. Wall boundary conditions

The EV26 equations presented above ((6.1), (6.2), (6.5), (6.9)–(6.13)) contain many terms with higher derivatives. For the computation of boundary-value problems, they must be furnished with appropriate boundary conditions.

Development of boundary conditons for the ideal gas moment equations is well understood; e.g. see Torrilhon & Struchtrup (Reference Torrilhon and Struchtrup2008) and Torrilhon (Reference Torrilhon2016), with moment boundary conditions developed from boundary conditions for the Boltzmann equation and the Grad distribution function.

A similar approach is expected to find boundary conditions for EV26. However, due to the higher-order derivatives appearing in the equations, as a result of the Taylor expansions (essentially the terms with $a^4$), it is expected that EV26 requires a larger number of boundary conditions than Grad26 (or R13).

Sufficiently far from the interface and walls, the bulk liquid and vapour behaviour are governed by the classical NSF laws. For studying non-equilibrium phase interfaces far from walls, it suffices to impose the NSF equations at the system boundary. This approach is used in results presented below.

Classical kinetic theory boundary conditions are effectively based on impulsive wall scattering models (Cercignani Reference Cercignani1975), and do not consider Vlasov like long-distance interaction between wall and a vapour or liquid particle. Hence both the kinetic boundary conditions for the EV equation and the development of a complete set of boundary conditions for its 26 moments approximation will require further thought.

7.1. Equilibrium

For additional context, we will have a brief look at equilibrium rest states, where all time derivatives vanish, and we expect homogeneous temperature, and vanishing non-equilibrium variables, so that

(7.2)\begin{equation} \frac{\partial\theta}{\partial x_{i}}=v_{i}=\sigma_{ij}=q_{i}=m_{ijk}=\varDelta=R_{ij}=0. \end{equation}

Then apart from the momentum balance, all equations are identically fulfilled. In a one-dimensional setting, where $\rho =\rho ( x)$, the momentum balance reduces to an equation for density (no body force $G_{i}$):

(7.3)\begin{gather} \frac{\partial}{\partial x}\left\{\begin{array}[c]{@{}c@{}} \rho\theta\left( 1+\dfrac{2{\rm \pi}}{3}\rho Y\right) -\dfrac{2{\rm \pi} a^{3}}{3}\chi_{1}\rho^{2}\\ {}+\dfrac{\rm \pi}{60}\theta\left[\dfrac{\partial^{2}Y\rho^{2}}{\partial x\,\partial x}+3Y \left(\dfrac{\partial^{2}\rho^{2}}{\partial x\partial x}-4\,\dfrac{\partial\rho}{\partial x}\, \dfrac{\partial\rho}{\partial x}\right) \right]\,{-}\,\dfrac{2{\rm \pi}}{15}\chi_{3} \left(\rho\,\dfrac{\partial^{2}\rho}{\partial x\,\partial x}-\dfrac{1}{2}\, \dfrac{\partial\rho}{\partial x}\,\dfrac{\partial\rho}{\partial x}\right) \end{array}\right\}{=}\,0. \end{gather}

That is, the divergence of the overall normal stress $\{ \bullet \}$ (the expression between the braces in (7.3)) vanishes. This equation allows phase equilibrium between saturated liquid and vapour phases, with continuous density variation across the interface.

For the case without interface, that is, for a single phase, the gradients of density vanish, and the pressure

(7.4)\begin{equation} p=\rho\theta\left( 1+\frac{2{\rm \pi}}{3}\rho Y\right) -\frac{2{\rm \pi}}{3}\chi_{1} \rho^{2} \end{equation}

is homogeneous.

It is well known that van-der-Waals-like equations of state such as this are non-monotonic, and allow for phase equilibrium (Struchtrup Reference Struchtrup2014). We denote the specific volume as $v=1/\rho$. The critical point is found from the condition for a horizontal inflection point in the $p$–$v$ diagram,

(7.5)\begin{equation} \left(\frac{\partial p}{\partial v}\right) _{\theta,crit}= \left(\frac{\partial^{2}p}{\partial v^{2}}\right) _{\theta,crit}=0, \end{equation}

as

(7.6a–c)\begin{equation} \theta_{crit}=0.7546,\quad p_{crit}=0.06748,\quad v_{crit}=4.014. \end{equation}

For subcritical temperatures, $\theta <\theta _{crit}$, equilibrium solutions exist, where a continuous interface connects two homogeneous regions at different mass densities – these are the states of saturated liquid and vapour, which both are at the unique saturation pressure $p_{sat}(\theta )$. As an example, figure 1 shows the liquid–vapour equilibrium for a system of length $L=40$ at temperature $\theta =0.65$, obtained from numerical solution of (7.3). The location of the interface depends on the average density $\bar {\rho }=({1}/{L})\int _{0}^{L}\rho \,{\textrm {d}x}$ of the system, which must lie between the saturated liquid and vapour densities, $\rho _{V}<\bar {\rho }<\rho _{L}$.

Figure 1. Density profile in liquid–vapour equilibrium for $\theta =0.65$.

According to (7.3), the overall normal stress is a constant. The pressure $p( \theta,\rho )$ of (7.4), however, varies throughout the interface, as shown in figure 2, where it is balanced by the Korteweg stresses. The pressure in the bulk phases, and hence the overall stress, equals the saturation pressure $p_{sat}$, as indicated in the figure.

Figure 2. Pressure profile $p( \theta,\rho )$ over the interface. The density curve is underlaid to identify the location of the liquid–vapour interface.

Solving (7.3) at other subcritical temperatures gives a complete set of data for saturation pressure $p_{sat}( \theta )$ and the corresponding saturated liquid and vapour volumes $v_{L}( \theta ) =1/\rho _{L}( \theta )$ and $v_{V}( \theta ) =1/\rho _{V}( \theta )$. These are collected in the $p$–$v$ diagram of figure 3, which shows the saturated liquid and vapour lines.

Figure 3. Saturation curve (blue) and three isothermal lines $p(\theta,v)$ (red) in the $p$–$v$ diagram. Saturated liquid and vapour volumes $v_{L}$ and $v_{V}$, and saturation pressure $p_{sat}$, are indicated for $\theta =0.65$. The black dot is the critical point.

The $p$–$v$ diagram also shows the isotherms of pressure $p( \theta,v)$ of (7.4) with the well-known non-monotonic shape of the subcritical isotherm (e.g. $\theta =0.65$). Thermodynamic stability requires $({\partial p}/{\partial v})_{\theta }<0$ (Struchtrup Reference Struchtrup2014), hence the pressure (7.4) alone is not stable. The interface is stabilized through the Korteweg stresses in (7.3), which balance the non-monotonic parts of the pressure, resulting in the continuous interface structure shown in figure 1.

We note that the determination of the saturation state for the EV gas relies on the evaluation of the balance of momentum (7.3). The typical thermodynamic argument employs Maxwell's equal area rule (Struchtrup Reference Struchtrup2014), where the saturation state follows from the condition

(7.7)\begin{equation} \int_{L}^{V}p( \theta,v)\,{\rm d} v=p_{sat}\left( v_{V}-v_{L}\right). \end{equation}

For $\theta =0.65$, the equal area rule gives $p_{sat}=0.03098$, $\rho _{L}=0.5205$, $\rho _{V}=0.06552$; in comparison, the EV results have minute relative errors of 0.09 %, 0.13 %, 0.005 %.

Thermodynamic equilibrium requires that the Gibbs free energies of both phases agree, and this implies the equal area rule (Struchtrup Reference Struchtrup2014). Hence the equations as derived stand in (slight) contradiction to thermodynamics. Aifantis & Serrin (Reference Aifantis and Serrin1983) presented a condition on capillary stress terms to agree with the equal area rule. Evaluation shows that the capillary stresses from the Vlasov term fulfil the condition, but the fifth-order Enskog contributions do not. A small ad hoc correction to the Enskog stress term of the form

(7.8)\begin{equation} P_{ij,{correction}}=\frac{\rm \pi}{60}\,\frac{a^{5}}{m}\left( -\frac{1}{2}\theta\rho^{2}\, \frac{{\rm d}^{2}Y}{{\rm d}\rho^{2}}\,\frac{\partial\rho}{\partial x_{s} }\, \frac{\partial\rho}{\partial x_{t}}\right) \delta_{(ik}\delta_{st)} \end{equation}

suffices to retrieve the equal area rule; see Appendix C.

7.2. Isothermal Couette flow

As a first test for the set of 26 moment equations for the EV equation, we show a numerical solution for the case of a steady two-phase and isothermal Couette flow (Rah & Eu Reference Rah and Eu2001; Frezzotti & Rossi Reference Frezzotti and Rossi2012), along a liquid–vapour interface, where we compare results to a direct simulation Monte Carlo (DSMC). In the considered test problem, an infinite and planar liquid film, of nominal thickness $x_{l}$, is in contact with its vapour, through the liquid–vapour interface. A shear flow is induced in the two-phase system in such a way that, everywhere in the flow field, the velocity has a constant direction, parallel to the vapour liquid interface. Specifically, we use a coordinate system such that the interface normal points in the $x_{1}$ direction and the liquid film occupies the region $\mathcal {L}=\{ (x_{1},x_{2},x_{3})\in R^{3}:0\leq x_{1}< x_{l}\}$. The velocity, which only has component $v_{2}( x_{1})$, is assumed to be zero at $x_{1}=0$, and to have constant gradient ${\partial v_{2}}/{\partial x_{1}}=\alpha$ far from the interface. Frictional heating is proportional to $\alpha ^{2}$ and can be ignored in the linearized approximation, assuming infinitely small shear rate. Hence the transport equations can be solved for the isothermal case, where specifically,

(7.9)\begin{equation} \frac{\partial\theta}{\partial x}=\varDelta=0. \end{equation}

Then the only non-vanishing variables are

(7.10)\begin{equation} \left\{ \rho( x),\ v_{2}( x), \sigma_{12}( x),\ q_{2}( x), m_{112}( x),\ R_{12}( x) \right\}. \end{equation}

Within the approximations used, the density $\rho ( x)$ is undisturbed and remains at the equilibrium solution, as given by (7.3), solely determined by the assumed temperature and liquid film thickness. The other variables assume non-equilibrium values, governed by the following set of linear differential equations, in which $x_{1}$ has been renamed $x$.

Momentum conservation, direction $x_{2}$:

(7.11)\begin{align} & \frac{\partial}{\partial x}\left[ \left( 1+\frac{2}{5}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) \sigma_{12} \right.\nonumber\\ &\quad \left.{}-\frac{a^{4}\sqrt{\rm \pi}}{m}\, \frac{\rho^{2}Y}{\sqrt{\theta}}\left( \frac{4}{15}\,\theta\,\frac{\partial v_{2} }{\partial x}+ \frac{2}{75}\,\frac{\partial}{\partial x}\left( \frac{q_{2}}{\rho}\right) +\frac{2}{105}\, \frac{\partial}{\partial x}\left( \frac{m_{112} }{\rho}\right) \right) \right] =0. \end{align}

Balance equation for kinetic shear stress $\sigma _{12}$:

(7.12)\begin{align} & \left[\rho\theta\left( 1+\frac{2}{5}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \right] \frac{\partial v_{2}}{\partial x}+\frac{2}{5}\,\frac{\partial}{\partial x} \left[\left( 1+\frac{3}{5}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) q_{2}\right] \nonumber\\ &\qquad +\frac{\partial}{\partial x}\left[ \left( 1+\frac{6}{35}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) m_{112}\right] +\frac{2{\rm \pi}}{3}\,\frac{a^{3}} {m}\left[ \frac{3}{25}q_{2}+\frac{8}{35}m_{112}\right] \frac{1}{\rho}\, \frac{\partial\rho^{2}Y}{\partial x} \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{8}{35}\,\rho\, \frac{\partial}{\partial x}\left( Y\,\frac{\partial\sigma_{12}}{\partial x}\right) -\frac{88}{105}\,\sigma_{12}\,\frac{\partial}{\partial x}\left( Y\,\frac{\partial\rho}{\partial x}\right) -\frac{16}{105}\rho\sigma_{12}\,\frac{\partial^{2}Y}{\partial x^{2}}\right] \nonumber\\ &\qquad-\frac{\sqrt{\rm \pi}a^{4}}{m\sqrt{\theta}}\left[ \frac{4}{245}\,\rho\, \frac{\partial}{\partial x}\left( Y\,\frac{\partial R_{12}}{\partial x}\right) -\frac{4}{105}R_{12}\,\frac{\partial}{\partial x}\left( Y\,\frac{\partial\rho}{\partial x}\right) -\frac{4}{735}\rho R_{12}\,\frac{\partial^{2}Y}{\partial x^{2}}\right] \nonumber\\ &\quad ={-}\left[ \frac{16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Y\right] \left(\sigma_{12}+\frac{1}{28}\,\frac{R_{12}}{\theta}\right). \end{align}

Balance equation for kinetic heat flux component $q_{2}$:

(7.13)\begin{align} & \theta\,\frac{\partial}{\partial x}\left[ \left( 1+\frac{3}{5}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) \sigma_{12}\right] -\frac{\sigma_{12}}{\rho}\,\theta\, \frac{\partial}{\partial x}\left[ \left( 1+\frac{9}{10}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) \rho\right] \nonumber\\ &\qquad +\frac{1}{2}\,\frac{\partial}{\partial x}\left[ \left( 1+\frac{12}{35}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) R_{12}\right] +\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\frac{11}{140}\,\frac{R_{12}}{\rho}\,\frac{\partial\rho^{2}Y}{\partial x} \nonumber\\ &\qquad +\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\,\frac{19}{15}\theta \left[\frac{\partial}{\partial x}\left( \rho^{2}Y\,\frac{\partial v_{2}}{\partial x}\right) \right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{9}{50}\,\rho\, \frac{\partial}{\partial x}\left( Y\,\frac{\partial q_{2}}{\partial x}\right) -\frac{1}{2}q_{2}\,\frac{\partial}{\partial x}\left( Y\,\frac{\partial\rho}{\partial x}\right) -\frac{2}{25}\rho q_{2}\,\frac{\partial^{2}Y}{\partial x^{2}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{9}{70}\,\rho\, \frac{\partial}{\partial x}\left( Y\,\frac{\partial m_{112}}{\partial x}\right) -\frac{1}{6}m_{122}\,\frac{\partial}{\partial x}\left( Y\,\frac{\partial\rho}{\partial x}\right) -\frac{1}{105}\rho m_{112}\,\frac{\partial^{2}Y}{\partial x^{2}}\right] \nonumber\\ &\quad ={-}\frac{2}{3}\left[ \frac{16}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Y\right] q_{2}. \end{align}

Balance equation for kinetic third-order moment $m_{112}$:

(7.14)\begin{align} & \frac{8}{5}\theta\,\frac{\partial}{\partial x}\left[ \left( 1+\frac{6}{35}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \sigma_{12}\right] -\frac{8}{5}\,\frac{\theta\sigma_{12}}{\rho}\,\frac{\partial}{\partial x} \left[\rho+\frac{2}{5}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho^{2}Y\right] \nonumber\\ &\qquad +\frac{8}{35}\,\frac{\partial}{\partial x}\left[ \left( 1+\frac{12}{35}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) R_{12}\right] +\frac{124}{1225}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\frac{R_{12}}{\rho}\,\frac{\partial\rho^{2}Y}{\partial x} \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{32}{525}\rho\theta Y\, \frac{\partial\rho}{\partial x}\,\frac{\partial v_{2}}{\partial x}+ \frac{16}{525}\rho^{2}\theta\,\frac{\partial}{\partial x} \left(Y\,\frac{\partial v_{2}}{\partial x}\right) \right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{72}{875}\rho\, \frac{\partial}{\partial x}\left( Y\,\frac{\partial q_{2}}{\partial x}\right) -\frac{8}{75}q_{2}\,\frac{\partial}{\partial x}\left( Y\,\frac{\partial\rho}{\partial x}\right) -\frac{16}{2625}\rho q_{2}\,\frac{\partial^{2}Y}{\partial x^{2}}\right] \nonumber\\ &\qquad -\frac{\sqrt{\rm \pi}a^{4}}{m}\sqrt{\theta}\left[ \frac{3}{25}\rho\, \frac{\partial}{\partial x}\left( Y\,\frac{\partial m_{112}}{\partial x}\right) -\frac{37}{35}m_{112}\,\frac{\partial}{\partial x}\left( Y\,\frac{\partial\rho}{\partial x}\right) -\frac{41}{175}\rho m_{112}\,\frac{\partial^{2}Y}{\partial x^{2}}\right] \nonumber\\ &\quad ={-}\frac{24}{5}\,\frac{\rho\sqrt{{\rm \pi}\theta}a^{2}}{m}\,Ym_{112}. \end{align}

Balance equation for kinetic fourth-order moment $R_{12}$:

(7.15)\begin{align} & \frac{14}{5}\theta\,\frac{\partial}{\partial x}\left[ \left( 1+\frac{12} {35}\,\frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) q_{2}\right] -\frac{14}{5}\, \theta\,\frac{q_{2}}{\rho}\,\frac{\partial}{\partial x}\left[ \left(1+\frac{1}{2}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) \rho\right] \nonumber\\ &\qquad +2\theta\,\frac{\partial}{\partial x}\left[ \left( 1+\frac{12}{35}\, \frac{2{\rm \pi}}{3}\,\frac{a^{3}}{m}\,\rho Y\right) m_{112}\right] -2\,\frac{m_{112}\theta}{\rho}\, \frac{\partial}{\partial x}\left[ \left( 1+\frac{11}{14}\,\frac{2{\rm \pi}}{3}\, \frac{a^{3}}{m}\,\rho Y\right) \rho\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[ \frac{8}{35}\rho\theta\, \frac{\partial}{\partial x}\left( Y\,\frac{\partial\sigma_{12}}{\partial x}\right) -\frac{8}{15}\theta\sigma_{12}\,\frac{\partial}{\partial x} \left(Y\,\frac{\partial\rho}{\partial x}\right) -\frac{8}{105}\rho\theta\sigma_{12}\, \frac{\partial^{2}Y}{\partial x^{2}}\right] \nonumber\\ &\qquad -\frac{a^{4}\sqrt{\rm \pi}}{m}\sqrt{\theta}\left[\frac{12}{49}\rho\, \frac{\partial}{\partial x}\left( Y\,\frac{\partial R_{12}}{\partial x}\right) -\frac{704}{735}R_{12}\,\frac{\partial}{\partial x}\left( Y\,\frac{\partial\rho}{\partial x}\right) -\frac{131}{735}\rho R_{12}\,\frac{\partial^{2}Y}{\partial x^{2}}\right] \nonumber\\ &\quad ={-}\left[ \frac{16}{5}\,\frac{a^{2}\rho\sqrt{{\rm \pi}\theta}}{m}\,Y\right] \left(\frac{1}{2}\theta\sigma_{12}+\frac{205}{168}R_{12}\right). \end{align}

The system of five linear second-order ordinary differential equations (7.11)–(7.15) defines a boundary-value problem, once it is prescribed that at $x=0$ and far from the interface the Navier–Stokes solution for a uniform shear flow holds.

The linearity allows us to consider approximations of intermediate complexity, with respect to the above full set of five equations. The simplest one, named 2-modes approximation, consists in forcing $q_{2}=m_{112}=R_{12}=0$ and describing the flow by (7.11) and (7.12). It is interesting to note that the 2-modes approximation provides the correct Navier–Stokes solution in the regions of uniform density but it does not coincide with the isothermal DIM description of the interface flow.

Grad's 13-moments approximation is obtained by describing the flow by the fields $v_{2}(x)$, $\sigma _{12}(x)$, $q_{2}(x)$, governed by (7.11)–(7.13) with $m_{112}=R_{12}=0$.

The equations associated with the approximations described above have been solved by turning them into linear algebraic systems, through quite standard discretization methods. The coefficients of the unknowns and their derivatives have been computed by approximating the unperturbed density profile by a hyperbolic tangent profile, matched to the solution of (7.3) (Frezzotti et al. Reference Frezzotti, Gibelli and Lorenzani2005).

Moment equations solutions have been compared to DSMC (stochastic) solutions of the EV equation, tailored to the considered Couette flow configuration.

In DSMC, two infinite planar and parallel walls are located at positions $x=\pm x_{w}$, where $x$ is the coordinate spanning the direction normal to the walls. The central part of the gap between the walls is occupied by a liquid slab of width $2x_{l}$, surrounded by two symmetrically located gas regions of width $x_{w}-x_{l}$. A simple steady shear flow is generated in the two-phase fluid by assigning velocities $\pm U_{w} \hat {\boldsymbol {y}}$ to the walls at $x=\pm x_{w}$, where $\hat {\boldsymbol {y}}$ is a unit vector parallel to the walls. It is further assumed that gas molecules hitting walls are re-emitted with velocities distributed according Maxwellians characterized by wall velocities $\pm U_{w}\hat {\boldsymbol {y}}$ and temperature $T_{w}$. A thermostat keeps the central part of the liquid slab at temperature $T_{w}$, thus absorbing the energy dissipated into the system by shear stresses.

The results presented below have been obtained by time-averaging the velocities of 24 million particles over about $10^{5}$ time steps, after the onset of a steady solution. The large sample size is necessary to resolve the small deviation from equilibrium required by the comparison of the DSMC solution of the nonlinear EV equation with solutions of the linearized moment equations. In this respect, it is worth mentioning that although the DSMC shear rate, tuned by $U_{w}$, has been set as small as permitted by the noise-to-signal ratio, the DSMC solution is not exactly isothermal. In the preliminary example shown below, the normalized equilibrium temperature has been set equal to $0.6$. The DSMC shows a slight temperature maximum of $0.605$ next to the liquid–vapour interface. Although further numerical work is required, we believe that the present accuracy level is sufficient to judge the quality of the moment formulation described above.

We proceed with comparing results from EV26 and DSMC. Due to the large sample size required, the computational time for DSMC is relatively large, about 24 hours, on a workstation equipped with a 10-core CPU. We note that DSMC requires unsteady three-dimensional simulations (with periodic boundaries in the $y$, $z$ directions). In contrast, this one-dimensional problem can be solved in about one minute with the moment equations, which benefit from a lower number of variables, the dimensional reduction, steadiness and linearity of the problem formulation.

Figure 4 shows velocity $v_{2}$, kinetic shear stress $\sigma _{12}$, kinetic heat flux $q_{2}$, and moments $m_{112}$, $R_{12}$ for dimensionless temperature $\theta =0.6\simeq 0.8\times \theta _{crit}$, as computed from the 26 moment system, the reduced system with 13 moments (variables $\rho,\ v_{2}, \sigma _{12},\ q_{2}$), and the further reduced 2-mode system (variables $\rho,\ v_{2},\ \sigma _{12}$). The results are compared with the DSMC solution of the EV equation itself (Frezzotti Reference Frezzotti1997) that serves as a reference.

Figure 4. Two-phase Couette flow: comparison of EV DSMC and linearized moment equations, in the case $\theta =0.6$. (a) Shear velocity $v_{2}$; (b) kinetic shear stress $\sigma _{12}$; (c) kinetic heat flux component $q_{2}$; (d) third- and fourth-order moments $m_{112}$ and $R_{12}$.

The main features of the flow field can be seen in figure 4(a), where the velocity field exhibits the structure of two uniform shear flows in the liquid and vapour regions, respectively, connected by the interface flow region. The macroscopically sharp change in velocity over the interface is a resolved slip effect, which would be dealt with as a slip boundary condition in a traditional hydrodynamic context. As can be seen from figures 4(a) and  4(b), moment theories with 26 and 13 moments give very good agreement to DSMC for velocity $v_{2}$ and kinetic stress $\sigma _{12}$. The 2-mode model does not succeed in approximating the velocity change across the interface, predicting a too-weak velocity slip.

The kinetic heat flux $q_{2}$ shown in figure 4(c) is due to forced gas motion, and not driven by a temperature gradient. Such a contribution to heat flux does not occur in classical hydrodynamics (Navier–Stokes–Fourier), but is well known in rarefied gas dynamics (Struchtrup Reference Struchtrup2005). The heat flux computed from the 26 moment system is in reasonably good agreement with the DSMC result, while the 13 moments theory clearly overpredicts the heat flux. This indicates the importance of including higher-order moments ($m_{112},R_{12}$), which have clear influence on process detail around the interface. As can be seen in figure 4(d), the 26 moment system gives a good approximation for $m_{112}$, while somewhat overpredicting the value of $R_{12}$ in the vapour in front of the interface. While velocity and stress variations are confined to the interfacial region as defined through the density variation, the higher moments ($q_{2} ,m_{112},R_{12}$) exhibit variations further into the vapour – these are Knudsen layer contributions.