\documentclass[12pt]{article}
\usepackage{graphicx}
\usepackage{color}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{marvosym}
\newcommand{\envelope}{(\raisebox{-.5pt}{\scalebox{1.45}{\Letter}}\kern-1.7pt)}
\usepackage{times}
\usepackage[LY1]{fontenc}

\begin{document}

\title{Grad's 14-moments for dense gases of spheres subjected to inelastic collisions}

\author{Elvira Barbera and Annamaria Pollino\\
Department of Mathematical, Computer, \\ 
Physical and Earth Sciences,\\
University of Messina,\\ 
ebarbera@unime.it, anpollino@unime.it}

\maketitle

\begin{abstract} 
A 14-moments model for dense granular gas is developed in the context of Grad's theory. 
The gas consists of spheres subjected to inelastic collisions. 
This condition implies that the energy is not conserved, but a part of it is transformed into heat 
and this determines a dissipation of temperature.
Furthermore, dense gases are characterized by two phenomena in inelastic collisions: the transfer and the transport of
particle properties.
The balance laws for mass, momentum, energy, stress tensor, heat flux and the fourth scalar moment are derived. 
The fluxes and the source terms are computed for each balance equation following 
the procedures elaborated by Jenkins and Richman.
The results are compared with the dilute elastic and inelastic gases.
This new set of field equations can be useful in different applications.
\end{abstract}

\section{Introduction}

The simplest case of gases studied in kinetic theory are dilute and elastic gases, 
that are gases whose particles have positions and velocities independent of each other and they are subject to elastic collisions. 
In dilute gases, when two particles collide, the contribution of the other particles is not taken into account because they are very spaced between them. 
In addition, the centers of the two colliding particles are assumed to coincide. 
Normally, two phenomena in binary collision can be identified: the transfer of particle properties, such as momentum and energy, and the transport of these properties. In the dilute case, the transport predominates and the transfer is neglected.
In this work, we will focus on a more complex case, represented by dense and granular gases: 
the particles have positions and velocities that depend on each other and the collisions are inelastic. 
This last condition implies that the energy is not conserved, but a part of it is transformed into heat and this determines a temperature dissipation. 
In dense gases, the positions of the centers of two colliding particles are distinguished, 
and in collision the position of the two particles is affected by the presence of neighboring particles. Ultimately, the phenomenon of transfer of particle properties and transport are both considered. Mathematically, this means that there are terms in the balance equations, related to transport and transfer, that are considered. 

The starting point for the present research is the paper of Jenkins and Richman \cite{Jenkins}, 
in which a model for 13-moments dense granular gases is presented in the context of Grad's theory \cite{Grad1,Grad2}. 
Using the techniques of kinetic theory fluxes associated with transport, 
fluxes associated with transfer and source terms are calculated for each balance equation. 
We amplify the model in \cite{Jenkins} by adding a 14th balance law. 
Then we construct a model of differential equations for 14 moments and we determine all fluxes and source terms. 

In literature, the 14-moments equations are already studied by Kremer and Marques \cite{Kremer1,Kremer2} that 
also amplified the 13-moments model by adding a 14th variable. In addition, they studied the decays of
temperature, stress tensor, heat flux and this new 14th moment together with 
the stability of linear perturbation of a time-dependent solution.

Also other authors \cite{ted} studied the cases with more moments with the stability of the homogeneous 
time-dependent solution. 
On the basis of the solution obtained by Kremer and other authors \cite{Kremer1,Kremer2,ted,Brillantov}
an analysis of the propagation of acceleration waves propagating into the spatially
homogeneous case was conducted in \cite{Acc1D}.

In this paper, we follow the Jenkins and Richman method \cite{Jenkins} without any truncation of the field equations 
and we determine a particular solution for the obtained set of field equations.
More in detail the article is organized as follows: 
In Sect.2 we define the first 14 moments of the kinetic theory, that are the unknown fields: 
the density, the velocity, the temperature, the stress tensor, the heat flux and the 14th scalar moment $\rho_{llkk}$. 
In Sect.3 we derive the 14 balance laws from the Boltzmann equation.  
In Sections 4 and 5, we present the setting of the kinetic theory. 
In particular, in Sect.4 we define the inelastic collisions analytically and in Sect.5 we discuss the approximations of the distribution function 
necessary for calculating various production terms. 
Sect.6 is devoted to the computations of all fluxes and source terms generalizing for the 14-moment model the formulas, obtained from Jenkins and Richman, which are linear approximations in the non-equilibrium fields. 
Finally, in Sect.7 a particular solution is obtained and discussed.

\section{Definition of moments}

Kinetic theory \cite{CC} deals with non-equilibrium problems concerning gases or gas mixtures. 
The starting point is the distribution function that is defined in the phase space
consisting of macroscopic variables, such as time and space, and microscopic variables, such as velocity. 
Precisely, we say that the function $f\left( t,\boldsymbol{x},\boldsymbol{c}\right) $ is the single
particle distribution function, if
\begin{equation}
f\left( t,\boldsymbol{x},\boldsymbol{c}\right) \text{d}\boldsymbol{c}
\end{equation}%
represents the number of particles that, at the time $t$, 
are in the position $\boldsymbol{x}$ with velocity 
in the range $\boldsymbol{c}$ and $\boldsymbol{c}+$d$\boldsymbol{c}$.
In terms of this function, we define the so-called moments of
the distribution function:
\begin{equation}
F_{i_{1}i_{2}...i_{N}}\left( t,\boldsymbol{x}\right) =m\int
c_{i_{1}}c_{i_{2}}...c_{i_{N}}f\left( t,\boldsymbol{x},\boldsymbol{c}\right) 
\text{d}\boldsymbol{c},
\end{equation}
with $m$ the molecular mass. 
These moments are related to macroscopic thermodynamic quantities \cite{CC}. 
In fact, the moment of zeroth and first order are related respectively to the density 
$\rho \left( t,\boldsymbol{x}\right)$ 
and the macroscopic velocity $\boldsymbol{v}\left( t,\boldsymbol{x}\right) $ that are
\begin{equation}
\begin{tabular}{c}
$\rho \left( t,\boldsymbol{x}\right) =m\int f\left( t,\boldsymbol{x},
\boldsymbol{c}\right) \text{d}\boldsymbol{c},\medskip$  \\ 
$\rho \left( t,\boldsymbol{x}\right) v_{i}\left( t,\boldsymbol{x}\right)
=m\int c_{i}$ $f\left( t,\boldsymbol{x},\boldsymbol{c}\right) \text{d}
\boldsymbol{c}.$
\end{tabular}
\label{def mom1}
\end{equation}%
Then, defining the relative velocity $\boldsymbol{C}=\boldsymbol{c}-\boldsymbol{v}$, 
it is possible to introduce the internal moments as
\begin{equation}
\rho _{i_{1}i_{2}...i_{N}}\left( t,\boldsymbol{x}\right) =m\int
C_{i_{1}}C_{i_{2}}...C_{i_{N}}f\left( t,\boldsymbol{x},\boldsymbol{c}\right)
d\boldsymbol{C}.
\label{def int_mom}
\end{equation}

The first internal moments are related to well-known macroscopic quantities,
that are the granular temperature, $\theta \left( t,\boldsymbol{x}\right) $,
the stress tensor $\rho _{ij}$ and the heat flux $q_{i}$:
\begin{equation}
\begin{tabular}{l}
$\theta \left( t,\boldsymbol{x}\right) =\frac{1}{3}\frac{m}{\rho }\int
C^{2}f\left( t,\boldsymbol{x},\boldsymbol{c}\right) \text{d}\boldsymbol{c},$
\medskip \\ 
$\rho _{ij}\left( t,\boldsymbol{x}\right) =m\int C_{i}C_{j}f\left( t,
\boldsymbol{x},\boldsymbol{c}\right) \text{d}\boldsymbol{c},\medskip $ \\ 
$q_{i}\left( t,\boldsymbol{x}\right) =\frac{1}{2}m\int C_{i}C^{2}f\left( t,
\boldsymbol{x},\boldsymbol{c}\right) \text{d}\boldsymbol{c}.$
\end{tabular}
\label{def mom2}
\end{equation}

Another quantity, that will be useful later, is the non-equilibrium part of
the so-called fourth moment, that is 
$\Delta \left( t,\boldsymbol{x}\right)=\rho _{llss}-15\rho \theta ^{2}$ where,
following the definition (\ref{def int_mom}), holds
\begin{equation}
\rho _{llss}=m\int C^{4}f\left( t,\boldsymbol{x},\boldsymbol{c}\right) \text{d}\boldsymbol{c}.  
\label{def mom3}
\end{equation}

Using the relation between the molecular and peculiar velocities, $c_{i}=C_{i}+v_{i}$, 
it is possible to derive the relation between the moments $F$'s and the 
corresponding internal moments $\rho$'s, 
so we split the moments in velocity-dependent parts and internal moments:
\begin{equation}
\begin{tabular}{l}
$F_{ij}=\rho _{ij}+\rho v_{i}v_{j},\medskip $ \\ 
$F_{ijk}=\rho _{ijk}+3\rho _{(ij}v_{k)}+\rho v_{i}v_{j}v_{k},\medskip $ \\ 
$F_{ikll}=\rho _{ikll}+4\rho _{(ikl}v_{l)}+6\rho _{(ik}v_{l}v_{l)}+\rho
v^{2}v_{i}v_{k},\medskip $ \\ 
$F_{kllss}=\rho _{kllss}+5\rho _{(klls}v_{s)}+10\rho
_{(kll}v_{s}v_{s)}+10\rho _{(kl}v_{l}v_{s}v_{s)}+\rho v^{4}v_{k}.$
\end{tabular}
\label{decomposition v}
\end{equation}

By this decomposition, 
it is equivalent to consider the $F$'s as field variables or equivalently the 
corresponding internal moments $\rho$'s.
In this paper, we consider the moments  $\rho$'s, as it is usually done inside the Grad's theory \cite{Grad1,Grad2}.

\section{Balance equations}

The distribution function $f\left( t,\boldsymbol{x},\boldsymbol{c}\right) $
obeys the Boltzman equation  
\begin{equation}
\frac{\partial f}{\partial t}+c_{i}\frac{\partial f}{\partial x_{i}}+f_{i}%
\frac{\partial f}{\partial c_{i}}=\mathbb{C}\left( f\right),
\label{Boltzman}
\end{equation}
where
$f_{i}$ represents the specific body force acting on the particles. 
The term $\mathbb{C}\left( f\right) $ is the collisional operator, 
that is regarded as a measure of the change of the distribution function due to collisions between particles. 
The Boltzman equation governs the evolution in time and space of the distribution function and allows us to predict the behavior of the gas under investigation.

The Boltzmann equation, that describes the gas at the microscopic level, 
can be also used to recover macroscopic equations.
If we multiply it by the velocity $c$ and integrate in the space of the velocities, 
it implies a set of balance equations for the moments. 
The purpose of this article is in fact to use the tools of kinetic theory to determine the fluxes and source terms of the balance laws, 
provided that the results here obtained are consistent with the continuum thermodynamic theories.

So, multiplication of the Boltzman equation (\ref{Boltzman}) by $mc_{i_{1}}c_{i_{2}}...c_{i_{N}}$
and integration over the whole range of $\mathbf{c}$ provides the balance
equations for the macroscopic fields 
$F_{i_{1}i_{2}...i_{N}}\left( t,\boldsymbol{x}\right) $, which assume the compact form
\begin{equation}
\frac{\partial F_{i_{1}i_{2}...i_{N}}}{\partial t}+\frac{\partial
F_{ki_{1}i_{2}...i_{N}}}{\partial x_{k}}%
-NF_{(i_{1}i_{2}...i_{N-1}}f_{i_{N})}=P_{i_{1}i_{2}...i_{N}}.
\end{equation}
Here the round brackets indicate the symmetric part of a tensor. The terms
in the right hand side represent the productions that will be discussed
later.

We assume that the gas can be described by the first fourteen moments, that
are density $\rho \left( t,\boldsymbol{x}\right) $, velocity $v_{i}\left( t,%
\boldsymbol{x}\right) $, temperature $\theta \left( t,\boldsymbol{x}\right) $%
, stress tensor $\rho _{ij}\left( t,\boldsymbol{x}\right) $, heat flux $%
q_{i}\left( t,\boldsymbol{x}\right) $ and the non-equilibrium part of the
double trace of the moment of rank four $\Delta \left( t,\boldsymbol{x}%
\right) $. Then, the balance equations for these 14 macroscopic fields are 
\begin{equation}
\begin{tabular}{l}
$\frac{\partial \rho }{\partial t}+\frac{\partial \rho v_{k}}{\partial x_{k}}%
=0,\medskip $ \\ 
$\frac{\partial \rho v_{i}}{\partial t}+\frac{\partial F_{ik}}{\partial x_{k}%
}=\rho f_{i}+P_{i},\medskip $ \\ 
$\frac{\partial F_{ij}}{\partial t}+\frac{\partial F_{ijk}}{\partial x_{k}}%
=2\rho v_{(i}f_{j)}+P_{ij},\medskip $ \\ 
$\frac{\partial F_{ill}}{\partial t}+\frac{\partial F_{ikll}}{\partial x_{k}}%
=3F_{(il}f_{l)}+P_{ill},\medskip $ \\ 
$\frac{\partial F_{llss}}{\partial t}+\frac{\partial F_{kllss}}{\partial
x_{k}}=4F_{(lss}f_{l)}+P_{llss}.$%
\end{tabular}
\label{balance F}
\end{equation}

Equations (\ref{balance F}) represent the classical hierarchy of 14
equations for classical monatomic gases.
This hierarchy together with the hierarchies corresponding to 13 and more moments
for classical monoatomic gasses are described in detail in \cite{giallo}.
 
\section{Inelastic collisions}

In order to evaluate the production terms, following the computations elaborated by Jenkins and Richman 
\cite{Jenkins}, we assume that the microscopic particles are identical spheres of diameter $\sigma $, 
the collisions between them are binary but inelastic. 
In this way, indicating with $\boldsymbol{c}^{1}$ and $\boldsymbol{c}^{1\prime}$ 
the velocities of the particle "1" before and after the collisions and
with $\boldsymbol{c}^{2}$ and $\boldsymbol{c}^{2\prime }$ those of the
particle "2", one has
\begin{equation}
\begin{tabular}{l}
$mc_{i}^{1\prime }=mc_{i}^{1}-J_{i},\medskip $ \\ 
$mc_{i}^{2\prime }=mc_{i}^{2}+J_{i}.$%
\end{tabular}
\label{J}
\end{equation}

If $g_{i}=c_{i}^{1}-c_{i}^{2}$ and $g_{i}^{\prime }=c_{i}^{1\prime}-c_{i}^{2\prime }$ 
are the relative velocities of the centers of the
spheres immediately before and after the collisions and $\boldsymbol{k}$ is
the unit vector directed from the center of particle "1" to the center of
"2" at contact, their components normal to the plane of contact are characterized by
\begin{equation}
\left( \boldsymbol{g}^{\prime }\cdot \boldsymbol{k}\right) =-e\left( 
\boldsymbol{g}\cdot \boldsymbol{k}\right),  
\label{inel}
\end{equation}
where $e$ is the so-called restitution coefficient.
It is a parameter that verifies $0\leq e\leq 1$.
The particular case $e=1$ corresponds to perfectly elastic collisions,
while, when $e<<1$, a strong inelastic collision
is considered.

By combination of (\ref{J}) and (\ref{inel}), it is possible to get 
the velocities of the particles "1" and "2"
before the collision in terms of the velocity after this collision:
\begin{equation}
\begin{tabular}{l}
$c_{i}^{1\prime }=c_{i}^{1}-\frac{1}{2}\left( 1+e\right) \left( \boldsymbol{g
}\cdot \boldsymbol{k}\right) k_{i},\medskip $ \\ 
$c_{i}^{2\prime }=c_{i}^{2}+\frac{1}{2}\left( 1+e\right) \left( \boldsymbol{g
}\cdot \boldsymbol{k}\right) k_{i}.$
\end{tabular}
\label{vel pre-post}
\end{equation}

From (\ref{vel pre-post}), it is easy to obtain some relations
that will be useful later in the determination of the productions, that are%
\begin{equation}
\begin{tabular}{l}
$c_{i}^{1\prime }c_{j}^{1\prime }-c_{i}^{1}c_{j}^{1}$ $\ \ \ \ \ \ \ =\frac{1%
}{2}\left( 1+e\right) \left( \boldsymbol{g}\cdot \boldsymbol{k}\right) \left[
\frac{1}{2}\left( 1+e\right) \left( \boldsymbol{g}\cdot \boldsymbol{k}%
\right) k_{i}k_{j}-\left( k_{i}c_{j}^{1}+k_{j}c_{i}^{1}\right) \right]
,\medskip $ \\ 
$c_{i}^{1\prime }c_{j}^{1\prime }c_{p}^{1\prime
}-c_{i}^{1}c_{j}^{1}c_{p}^{1}=-\frac{1}{2}\left( 1+e\right) \left( 
\boldsymbol{g}\cdot \boldsymbol{k}\right) \left[ 3k_{(i}c_{j}^{1}c_{p)}^{1}-%
\frac{3}{2}\left( 1+e\right) \left( \boldsymbol{g}\cdot \boldsymbol{k}%
\right) k_{(i}k_{j}c_{p)}^{1}+\right. \medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. +\frac{1}{4}\left(
1+e\right) ^{2}\left( \boldsymbol{g}\cdot \boldsymbol{k}\right)
^{2}k_{i}k_{j}k_{p}\right] ,\medskip $ \\ 
$c_{l}^{1\prime }c_{l}^{1\prime }c_{s}^{1\prime }c_{s}^{1\prime
}-c_{l}^{1}c_{l}^{1}c_{s}^{1}c_{s}^{1}=-2\left( 1+e\right) \left( 
\boldsymbol{g}\cdot \boldsymbol{k}\right) \left[
k_{l}c_{l}^{1}c_{s}^{1}c_{s}^{1}-\frac{1}{2}\left( 1+e\right) \left( 
\boldsymbol{g}\cdot \boldsymbol{k}\right) \left( k_{l}c_{l}^{1}\right)
^{2}+\right. \medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \left. -\frac{1}{4}\left( 1+e\right) \left( \boldsymbol{g}%
\cdot \boldsymbol{k}\right) c_{s}^{1}c_{s}^{1}+\frac{1}{4}\left( 1+e\right)
^{2}\left( \boldsymbol{g}\cdot \boldsymbol{k}\right) ^{2}k_{l}c_{l}^{1}-%
\frac{1}{32}\left( 1+e\right) ^{3}\left( \boldsymbol{g}\cdot \boldsymbol{k}%
\right) ^{3}\right] .$
\end{tabular}
\end{equation}
The first two relations are obtained by Jenkins and Richman \cite{Jenkins},
the third one is instead appropriate to the 14th moment.

Similarly, it is possible to obtain other important relations. 
In fact, setting $\Pi\left( \psi \right) =\psi ^{1\prime }+\psi ^{2\prime }-\psi ^{1}-\psi ^{2}$, 
one has
\begin{equation}
\begin{tabular}{l}
$\Pi (c_{i}c_{j})=\frac{1}{2}\left( 1+e\right) \left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) \left[ \left( 1+e\right) \left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) k_{i}k_{j}-\left( k_{i}g_{j}+k_{j}g_{i}\right) \right],\medskip $ \\ 
$\Pi (c_{i}c_{j}c_{k})=3Q_{(i}\Delta (c_{j}c_{k)}),\medskip $ \\ 
$\Pi (c^{4})=\frac{1}{8}\left( 1+e\right) ^{4}\left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) ^{4}-\frac{1}{2}\left( 1+e\right) ^{3}\left( 
\boldsymbol{g}\cdot \boldsymbol{k}\right) ^{4}+\medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ +\left( 1+e\right) ^{2}\left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) ^{2}\left[ Q^{2}+\frac{1}{4}g^{2}+2\left( \boldsymbol{Q%
}\cdot \boldsymbol{k}\right) ^{2}+\frac{1}{2}\left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) ^{2}\right] +\medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ -\left( 1+e\right) \left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) \left[ \frac{1}{2}\left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) g^{2}+2\left( \boldsymbol{g}\cdot \boldsymbol{k}%
\right) Q^{2}+4\left( \boldsymbol{g}\cdot \boldsymbol{Q}\right) \left( 
\boldsymbol{Q}\cdot \boldsymbol{k}\right) \right] $%
\end{tabular}%
\label{ccc}
\end{equation}%
with $Q_{i}=\frac{1}{2}\left( c_{i}^{1}+c_{i}^{2}\right) $  the velocity of the center of mass of the two particles. 
Also here, the first two conditions are evaluated by Jenkins and Richman \cite{Jenkins},
while the third is appropriate to the balance equation for the 14th moment.

It can be easily seen that, by contracting $i$ and $j$ in $\eqref{ccc}_1$, we get
\begin{equation}
\Pi\left( E\right)=-\frac{1}{4}m(1-e^2)(\boldsymbol{g}\cdot \boldsymbol{k})^2.
\end{equation}
This term represents the energy production that is due to the collisions between the molecules.
For inelastic collision this quantity is different from zero,
while it vanishes in the elastic case, that is $e=1$.

\section{Determination of the Production Terms}

The statistics of binary collisions is characterized by the complete pair
distribution function $f^{\left( 2\right) }\left( t,\boldsymbol{x}^{1},%
\boldsymbol{c}^{1},\boldsymbol{x}^{2},\boldsymbol{c}^{2}\right) $, so%
\begin{equation}
f^{\left( 2\right) }\left( t,\boldsymbol{x}^{1},\boldsymbol{c}^{1},
\boldsymbol{x}^{2},\boldsymbol{c}^{2}\right) \text{d}\boldsymbol{x}^{1}\text{
d}\boldsymbol{c}^{1}\text{d}\boldsymbol{x}^{2}\text{d}\boldsymbol{c}^{2}
\end{equation}%
is the probable number of pairs of particles that at the time $t$ are
located between $\boldsymbol{x}^{1}$ and $\boldsymbol{x}^{1}+$d$\boldsymbol{x
}^{1}$ and $\boldsymbol{x}^{2}$ and $\boldsymbol{x}^{2}+$d$\boldsymbol{x}
^{2} $ with the velocities between $\boldsymbol{c}^{1}$ and $\boldsymbol{c}
^{1}+$d$\boldsymbol{c}^{1}$ and $\boldsymbol{c}^{2}$ and $\boldsymbol{c}
^{2}+ $d$\boldsymbol{c}^{2}$, respectively.

In the case of dilute gases, the positions and the velocities of the
particles are assumed to be independent, therefore for rarefied gases the pair distribution function is assumed as the product
$f^{\left( 2\right) }\left( t,\boldsymbol{x}^{1},\boldsymbol{c}^{1},%
\boldsymbol{x}^{2},\boldsymbol{c}^{2}\right) =f\left( t,\boldsymbol{x}^{1},%
\boldsymbol{c}^{1}\right) f\left( t,\boldsymbol{x}^{2},\boldsymbol{c}%
^{2}\right) $. 
For such gases the transfer of particles proprieties, like momentum and energy in
collisions is neglected in comparison to the transport of particles
proprieties.

Instead, for dense gases the transfer of particles proprieties in collisions
is so important as the transport of particles proprieties. 
The probable position of two colliding particles is
influenced by the presence of neighboord particles.
To account for the interactions of neighboring particles when two particles collide, the equilibrium value of a radial distribution function $g_0$ is introduced into the relation that binds $f^{(2)}$ and the velocity distribution function of each particle $f$:
\begin{equation}
f^{\left( 2\right) }\left( t,\boldsymbol{x}^{1},\boldsymbol{c}^{1},%
\boldsymbol{x}^{2},\boldsymbol{c}^{2}\right) =g_{0}\left( x\right) f\left( t,%
\boldsymbol{x}^{1},\boldsymbol{c}^{1}\right) f\left( t,\boldsymbol{x}^{2},%
\boldsymbol{c}^{2}\right).
\end{equation}
\ \ \ \ \ \ \ 
There are different evaluations of the radial distribution function $g_0$, for example Charnahan
and Starling \cite{CS} evaluated the probability of collision in terms of the solid
volume fraction $\nu =n\pi \sigma ^{3}/6$ as 
\begin{equation}
g_{0}\left( \nu \right) =\frac{1}{1-\nu }+\frac{3\nu }{2\left( 1-\nu \right)
^{2}}+\frac{\nu ^{2}}{2\left( 1-\nu \right) ^{3}}.
\end{equation}

In presence of dense gases Jenkins and Richman in \cite{Jenkins} showed that
that the production terms can be expressed as%
\begin{equation}
P_{i_{1}i_{2}...i_{N}}=\Psi _{i_{1}i_{2}...i_{N}}-\frac{\partial \Theta
_{ki_{1}i_{2}...i_{N}}}{\partial x_{k}}.  \label{productions}
\end{equation}%
The quantities $\Psi _{i_{1}i_{2}...i_{N}}=\Psi \left(
mc_{i_{1}}c_{i_{2}}...c_{i_{N}}\right) $\ are defined as%
\begin{equation}
\begin{tabular}{l}
$\Psi \left( \psi \right) =\frac{1}{2}\int \int \int \left( \psi ^{1\prime
}+\psi ^{2\prime }-\psi ^{1}-\psi ^{2}\right) f^{\left( 2\right) }\left( 
\boldsymbol{x}-\sigma \boldsymbol{k},\boldsymbol{c}^{1},\boldsymbol{x},%
\boldsymbol{c}^{2}\right) \medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\sigma ^{2}\left( \boldsymbol{g}\cdot \boldsymbol{k}\right) \text{d}%
\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d}\boldsymbol{c}^{2}.$%
\end{tabular}
\label{psi}
\end{equation}
They represent the typical source terms that are already present in
monatomic gases. The other terms, that are not present in kinetic theory of
dilute gases, depend on the dense nature of the material under
consideration. One has $\Theta _{si_{1}i_{2}...i_{N}}=\Theta _{s}\left(
mc_{i_{1}}c_{i_{2}}...c_{i_{N}}\right) $ with 
\begin{equation}
\begin{tabular}{l}
$\Theta _{s}\left( \psi \right) =-\frac{\sigma }{2}\int \int \int \left(
\psi ^{1\prime }-\psi ^{1}\right) k_{s}\left( 1-\frac{\sigma }{2!}k_{j}\frac{%
\partial }{\partial x_{j}}+\frac{\sigma ^{2}}{3!}k_{j}k_{m}\frac{\partial }{%
\partial x_{j}\partial x_{m}}+...\right) \medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ f^{\left(
2\right) }\left( \boldsymbol{x},\boldsymbol{c}^{1},\boldsymbol{x}+\sigma 
\boldsymbol{k},\boldsymbol{c}^{2}\right) \sigma ^{2}\left( \boldsymbol{g}%
\cdot \boldsymbol{k}\right) \text{d}\boldsymbol{k}\text{d}\boldsymbol{c}^{1}%
\text{d}\boldsymbol{c}^{2}.$%
\end{tabular}
\label{teta}
\end{equation}

By substitution of (\ref{productions}-\ref{teta}) into the balance equations for the moments (\ref{balance F}), we get 
\begin{equation}
\begin{tabular}{l}
$\frac{\partial \rho }{\partial t}+\frac{\partial \rho v_{k}}{\partial x_{k}}%
=0,\medskip $ \\ 
$\frac{\partial \rho v_{i}}{\partial t}+\frac{\partial }{\partial x_{k}}%
\left[ F_{ik}+\Theta _{ik}\right] =\rho f_{i},\medskip $ \\ 
$\frac{\partial F_{ij}}{\partial t}+\frac{\partial }{\partial x_{k}}\left[
F_{ijk}+\Theta _{ijk}\right] =2\rho v_{(i}f_{j)}+\Psi _{ij},\medskip $ \\ 
$\frac{\partial F_{ill}}{\partial t}+\frac{\partial }{\partial x_{k}}\left[
F_{ikll}+\Theta _{ikll}\right] =3F_{(il}f_{l)}+\Psi _{ill},\medskip $ \\ 
$\frac{\partial F_{llss}}{\partial t}+\frac{\partial }{\partial x_{k}}\left[
F_{kllss}+\Theta _{kllss}\right] =4F_{(lss}f_{l)}+\Psi _{llss}.$%
\end{tabular}
\label{balance F produz}
\end{equation}
In the particular case of dilute gases, the $\Theta$s vanish.


We define the internal quantities $\theta $ and $\psi $ exactly 
like in the definitions (\ref{psi}) and (\ref{teta}) for the quantities $\Psi $ and $\Theta $,
but using the peculiar velocity $\boldsymbol{C}$ instead of $\boldsymbol{c}$.
So, taking again into account the relation $C_i=c_i-v_i$,
the following decomposition for the production terms holds 
\begin{equation}
\begin{tabular}{l}
$\Theta _{ij}=\theta _{ij},\medskip $ \\ 
$\Theta _{ijk}=\theta _{ijk}+2\theta _{k(i}v_{j)},\medskip $ \\ 
$\Theta _{ikll}=\theta _{ikll}+3\theta _{k(il}v_{l)}+3\theta
_{k(i}v_{l}v_{l)},\medskip $ \\ 
$\Theta _{kllss}=\theta _{kllss}+4\theta _{k(lls}v_{s)}+6\theta
_{k(ll}v_{s}v_{s)}+4\theta _{kl}v_{l}v_{s}v_{s},\medskip $ \\ 
$\Psi _{ij}=\psi _{ij},\medskip $ \\ 
$\Psi _{ill}=\psi _{ill}+3\psi _{(il}v_{l)}\medskip $ \\ 
$\Psi _{ikll}=\psi _{ikll}+4\psi _{(ikl}v_{l)}+\psi _{(ik}v_{l}v_{l)}.$%
\end{tabular}
\label{dec 2}
\end{equation}

Then, taking into account these decompositions and (\ref{decomposition v}), 
it is possible to obtain the 14 balance equations for the 14 internal moments:
\begin{equation}
\begin{tabular}{l}
$\frac{\text{d}\rho }{\text{d}t}+\rho \frac{\partial v_{k}}{\partial x_{k}}%
=0,\medskip $ \\ 
$\rho \frac{\text{d}v_{i}}{\text{d}t}+\frac{\partial }{\partial x_{k}}\left[
\rho _{ik}+\theta _{ik}\right] =\rho f_{i},\medskip $ \\ 
$\frac{\text{d}\rho _{ij}}{\text{d}t}+\frac{\partial }{\partial x_{k}}\left[
\rho _{ijk}+\theta _{ijk}\right] +\rho _{ij}\frac{\partial v_{k}}{\partial
x_{k}}+2\left[ \rho _{k(i}+\theta _{k(i}\right] \frac{\partial v_{j)}}{%
\partial x_{k}}=\psi _{ij},\medskip $ \\ 
$\frac{\text{d}\rho _{ill}}{\text{d}t}+\frac{\partial }{\partial x_{k}}\left[
\rho _{ikll}+\theta _{ikll}\right] +\rho _{ill}\frac{\partial v_{k}}{%
\partial x_{k}}+3\left[ \rho _{k(il}+\theta _{k(il}\right] \frac{\partial
v_{l)}}{\partial x_{k}}+\medskip $ \\ 
$\hspace{1.8in}-3\frac{\rho _{(il}}{\rho }\frac{\partial }{\partial x_{k}}%
\left[ \rho _{l)k}+\theta _{l)k}\right] =\psi _{ill},\medskip $ \\ 
$\frac{\text{d}\rho _{llss}}{\text{d}t}+\frac{\partial }{\partial x_{k}}%
\left[ \rho _{kllss}+\theta _{kllss}\right] +\rho _{llss}\frac{\partial v_{k}%
}{\partial x_{k}}+4\left[ \rho _{ksll}+\theta _{ksll}\right] \frac{\partial
v_{s}}{\partial x_{k}}+\medskip $ \\ 
$\hspace{1.8in}-8\frac{q_{l}}{\rho }\frac{\partial }{\partial x_{k}}\left[
\rho _{lk}+\theta _{lk}\right] =\psi _{llss}.$
\end{tabular}
\label{balance ins}
\end{equation}

The first equation represents the conservation law of mass, the second is
the balance law of momentum. The quantity $\rho _{ik}+\theta _{ik}$
is the total pressure tensor: the sum of the quantity due to the transport
of momentum between collisions and those related to the transfer of momentum. 

The trace of the third equation represents
the balance law of energy
\begin{equation}
\frac{3}{2}\rho \frac{\text{d}\theta }{\text{d}t}+\frac{1}{2}\frac{\partial %
\left[ \rho _{kll}+\theta _{kll}\right] }{\partial x_{k}}+\left[ \rho
_{kl}+\theta _{kl}\right] \frac{\partial v_{l}}{\partial x_{k}}=\frac{1}{2}%
\psi _{ll},
\end{equation}%
where $\theta $ is the so-called granular temperature, the sum $\frac{1}{2}%
\left[ \rho _{kll}+\theta _{kll}\right] $ is the heat flux, with the
transport and collisional parts, the last term $\psi _{ll}$ represents the
dissipation due to the inelastic nature of the collision.

The traceless part of equation $(\ref{balance ins})_3$ is instead 
the balance law for the stress tensor and it assumes the form
\begin{equation}
\frac{\text{d}\rho _{<ij>}}{\text{d}t}+\frac{\partial \left[ \rho
_{<ij>k}+\theta _{<ij>k}\right] }{\partial x_{k}}+\rho _{<ij>}\frac{\partial
v_{k}}{\partial x_{k}}+2\left[ \rho _{k<i}+\theta _{k<i}\right] \frac{%
\partial v_{j>}}{\partial x_{k}}=\psi _{<ij>}.
\end{equation}%
Square brackets in the indexes indicate traceless part of a tensor.

The remaining equations are instead the balance laws for the heat flux and the 14th moment,
since $q_i=\frac{1}{2}\rho _{ill}$ and $\Delta=\rho _{ssll}-15\rho\theta^2$.

In next section, we will close the set of balance equations (\ref{balance ins}) by the method of Grad's.

\section{Closure by Grad}

Equations (\ref{balance ins}) represent a system of 14 partial differential equations for the 14
fields $\rho $, $v_{i}$, $\theta $, $\rho _{<ij>}$, $q_{k}=1/2\rho _{kll}$ and $\Delta $. 
Unfortunately, these equations are not closed for the occurrence of the quantities $\rho _{<ijk>}$, $\rho _{<ik>ll}$, $\rho _{kllss}$, $\theta _{ik}$, 
$\theta _{ijk}$, $\theta _{ikll}$, $\theta _{kllss}$, $\psi _{ij}$, $\psi _{ill}$ and $\psi _{llss}$. 
In this section, we will obtain constitutive relations for these additional quantities 
using the methods of Grad's  \cite{Grad1}.
First of all, the single distribution function is written as
\begin{equation}
f\left( t,\boldsymbol{x},\boldsymbol{c}\right) =\left( 1-a_{i}\frac{\partial }{\partial c_{i}}+
a_{ij}\frac{\partial ^{2}}{\partial c_{i}\partial c_{j}}-a_{ijk}\frac{\partial ^{3}}{\partial c_{i}\partial c_{j}\partial c_{k}}
+...\right) f_{0}\left( t,\boldsymbol{x},\boldsymbol{c}\right),
\label{svil grad}
\end{equation}
where $f_{0}\left( t,\boldsymbol{x},\boldsymbol{c}\right) $ is the Maxwellian distribution function given by
\begin{equation}
f_{0}\left( t,\boldsymbol{x},\boldsymbol{c}\right) =\frac{n}{\left( 2\pi
\theta \right) ^{\frac{3}{2}}}e^{-\frac{C^{2}}{2\theta }},
\label{maxwellian}
\end{equation}
where $n$ is the number density defined as $n=\rho/m$.
Insertion of (\ref{maxwellian}) into (\ref{svil grad}) yields, after simple calculations, the single
distribution function in terms of the Hermite polynomials:
\begin{equation}
\begin{tabular}{l}
$f\left( t,\boldsymbol{x},\boldsymbol{c}\right) =f_{0}\left\{ 1+\frac{a_{i}}{%
\theta }C_{i}+\frac{a_{ij}}{\theta ^{2}}\left[ C_{i}C_{j}-\theta \delta _{ij}%
\right] +\frac{a_{ijk}}{\theta ^{3}}\left[ C_{i}C_{j}C_{k}-3\theta
C_{(i}\delta _{jk)}\right] +\right. \medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. +\frac{%
a_{ijks}}{\theta ^{4}}\left[ C_{i}C_{j}C_{k}C_{s}-6\theta \delta
_{(ij}C_{k}C_{s)}+3\theta ^{2}\delta _{(ij}\delta _{ks)}\right] +...\right\}
.$%
\end{tabular}
\label{svil grad2}
\end{equation}
The coefficients $\boldsymbol{a}$ can be evaluated by insertion of (\ref{svil grad2}), 
truncated to the fourth order, into the definitions of the moments (\ref{def mom1},\ref{def mom2},\ref{def mom3}), so it follows 
\begin{equation}
\begin{tabular}{llll}
$a_{i}=0,$ & $a_{ij}=\frac{\rho _{<ij>}}{2\rho },$ & $a_{ijk}=\frac{q_{i}}{5\rho }\delta _{jk},$ & 
$a_{ijks}=\frac{\Delta }{120\rho }\delta _{ij}\delta_{ks}.$
\end{tabular}
\end{equation}

With these quantities, the single velocity distribution function, appropriate to 14 moments, takes the form
\begin{equation}
\begin{tabular}{l}
$f\left( t,\boldsymbol{x},\boldsymbol{c}\right) =f_{0}\left\{ 1+\frac{\rho
_{<ij>}}{2\rho \theta ^{2}}C_{i}C_{j}+\frac{q_{i}C_{i}}{5\rho \theta ^{3}}%
\left[ C^{2}-5\theta \right] +\right. \medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. \frac{%
\Delta }{120\rho \theta ^{4}}\left[ C^{4}-10\theta C^{2}+15\theta ^{2}\right]
\right\} .$%
\end{tabular}
\label{grad 14}
\end{equation}

Once this approximate distribution function is obtained, its insertion into the definition of the moments 
$\rho _{ijk}$, $\rho _{<ij>ll}$ and $\rho _{illss}$ yields the determination of the first constitutive relations: 
\begin{equation}
\begin{tabular}{l}
$\rho _{ijk}=\frac{2}{5}\left( q_{i}\delta _{jk}+q_{j}\delta
_{ik}+q_{k}\delta _{ij}\right) ,\medskip $ \\ 
$\rho _{<ij>ll}=7\theta \rho _{<ij>},\medskip $ \\ 
$\rho _{illss}=28\theta q_{i}.$
\end{tabular}
\label{constitutive}
\end{equation}

\section{Collision integrals}

In \cite{Jenkins} Jenkins and Richman showed that the production terms can
be evaluated from (\ref{psi}) and (\ref{teta}) obtaining 
\begin{equation}
\begin{tabular}{l}
$\psi \left( \varphi \right) =\frac{g_{0}\left( \boldsymbol{x}\right) }{2}%
\int \int \int \Pi \left( \varphi \right) f^{\left( 1\right) }\left( 
\boldsymbol{x},\boldsymbol{c}^{1}\right) f^{\left( 1\right) }\left( 
\boldsymbol{x},\boldsymbol{c}^{2}\right) \times \medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \times \left[ 1+\frac{\sigma }{2}k_{i}\frac{\partial }{%
\partial x_{i}}\ln \frac{f^{\left( 1\right) }\left( \boldsymbol{x},%
\boldsymbol{c}^{2}\right) }{f^{\left( 1\right) }\left( \boldsymbol{x},%
\boldsymbol{c}^{1}\right) }\right] \sigma ^{2}\left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) \text{d}\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d%
}\boldsymbol{c}^{2}$%
\end{tabular}%
\end{equation}%
with $\Pi \left( \varphi \right) =\varphi ^{1\prime }+\varphi ^{2\prime
}-\varphi ^{1}-\varphi ^{2}$ for the source terms and 
\begin{equation}
\begin{tabular}{l}
$\theta _{i}\left( \varphi \right) =-\sigma \frac{g_{0}\left( \boldsymbol{x}%
\right) }{2}\int \int \int \left( \varphi ^{1\prime }-\varphi ^{1}\right)
k_{i}f^{\left( 1\right) }\left( \boldsymbol{x},\boldsymbol{c}^{1}\right)
f^{\left( 1\right) }\left( \boldsymbol{x},\boldsymbol{c}^{2}\right) \times
\medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \times \left[ 1+\frac{\sigma }{2}k_{j}\frac{\partial 
}{\partial x_{j}}\left( \ln \frac{f^{\left( 1\right) }\left( \boldsymbol{x},%
\boldsymbol{c}^{2}\right) }{f^{\left( 1\right) }\left( \boldsymbol{x},%
\boldsymbol{c}^{1}\right) }\right) \right] \sigma ^{2}\left( \boldsymbol{g}%
\cdot \boldsymbol{k}\right) \text{d}\boldsymbol{k}\text{d}\boldsymbol{c}^{1}%
\text{d}\boldsymbol{c}^{2}.$%
\end{tabular}%
\end{equation}

We provide approximate expression for $\theta _{i}(\psi)$ and $\psi(\varphi)$ that are linear in the perturbations 
and the spatial gradients of velocity and temperature. 
So we obtain for the 14-moment case
\begin{equation}
\begin{tabular}{l}
$\theta _{i}\left( \varphi \right) =A_{i}\left( \varphi \right) +B_{i}\left(
\varphi \right) +\rho _{jk}B_{ijk}\left( \varphi \right)
+q_{j}B_{ijll}\left( \varphi \right) +\Delta \hat{B}_{i}\left( \varphi
\right), \medskip $ \\ 
$\psi \left( \varphi \right) =E\left( \varphi \right) +F\left( \varphi
\right) +\rho _{jk}F_{jk}\left( \varphi \right) +q_{j}F_{j}\left( \varphi
\right) +\Delta \hat{F}\left( \varphi \right), $%
\end{tabular}
\label{espresslineari}
\end{equation}
where the integrals are given by
\begin{equation}
\begin{tabular}{l}
$A_{i} =-\frac{\sigma^3 }{2}g_{0}\int \int \int \left(
\varphi ^{1\prime }-\varphi ^{1}\right) k_{i}f_{01}f_{02} 
(\boldsymbol{g}\cdot \boldsymbol{k}) \text{d}\boldsymbol{k}\text{d}
\boldsymbol{c}^{1}\text{d}\boldsymbol{c}^{2},\medskip $ \\ 

$B_{i}=-\frac{\sigma ^{4}}{4}g_{0}\int \int \int
\left( \varphi ^{1\prime }-\varphi ^{1}\right) k_{i}k_{m}f_{01}f_{02}\frac{%
\partial }{\partial x_{m}}\left( \ln \frac{f_{02}}{f_{01}}\right) 
\left( \boldsymbol{g}\cdot \boldsymbol{k}\right) \text{d}\boldsymbol{k}%
\text{d}\boldsymbol{c}^{1}\text{d}\boldsymbol{c}^{2},\medskip $ \\ 

$B_{ijk}=-\frac{\sigma^3 }{4 \rho}g_{0}\int \int \int \left(
\varphi ^{1\prime }-\varphi ^{1}\right) k_{i}\left( f_{01}\frac{\partial
^{2}f_{02}}{\partial c_{j}^{2}\partial c_{k}^{2}}+f_{02}\frac{\partial
^{2}f_{01}}{\partial c_{j}^{1}\partial c_{k}^{1}}\right) \left( 
\boldsymbol{g}\cdot \boldsymbol{k}\right) \text{d}\boldsymbol{k}\text{d}
\boldsymbol{c}^{1}\text{d}\boldsymbol{c}^{2},\medskip $ \\ 

$B_{ijll}=\frac{\sigma^3 g_{0}}{10\rho}\int \int \int
\left( \varphi ^{1\prime }-\varphi ^{1}\right) k_{i}\left( \ f_{01}\frac{%
\partial ^{3}f_{02}}{\partial c_{j}^{2}\partial c_{l}^{2}\partial c_{l}^{2}}%
+f_{02}\frac{\partial ^{3}f_{01}}{\partial c_{j}^{1}\partial
c_{l}^{1}\partial c_{l}^{1}}\right) \left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) \text{d}\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d%
}\boldsymbol{c}^{2},\medskip $ \\ 

$\hat{B}_{i}=-\frac{\sigma^3 g_{0}}{240\rho}\int \int \int
\left( \varphi ^{1\prime }-\varphi ^{1}\right) k_{i}\left( \ f_{01}\frac{%
\partial ^{4}f_{02}}{\partial c_{l}^{2}\partial c_{l}^{2}\partial
c_{s}^{2}\partial c_{s}^{2}}+f_{02}\frac{\partial ^{4}f_{01}}{\partial
c_{l}^{1}\partial c_{l}^{1}\partial c_{s}^{1}\partial c_{s}^{1}}\right)
(\boldsymbol{g}\cdot \boldsymbol{k}) \text{d}%
\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d}\boldsymbol{c}^{2}$%
\end{tabular}
\label{AB}
\end{equation}
and
\begin{equation}
\begin{tabular}{l}
$E =\frac{\sigma ^{2}g_{0}}{2}\int \int \int \Pi \left( \varphi
\right) f_{01}f_{02}\left( \boldsymbol{g}\cdot \boldsymbol{k}%
\right) \text{d}\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d}\boldsymbol{c%
}^{2},\medskip $ \\ 

$F =\frac{\sigma ^{3} g_{0}}{4}\int \int \int \Pi \left(
\varphi \right) k_{m}f_{01}f_{02}\frac{\partial }{\partial x_{m}}\left( \ln 
\frac{f_{02}}{f_{01}}\right) \sigma ^{2}\left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) \text{d}\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d%
}\boldsymbol{c}^{2},\medskip $ \\ 

$F_{ij} =\frac{\sigma ^{2}g_{0}}{4\rho}\int \int \int \Pi \left(
\varphi \right) \left( f_{01}\frac{\partial ^{2}f_{02}}{\partial
c_{i}^{2}\partial c_{j}^{2}}+f_{02}\frac{\partial ^{2}f_{01}}{\partial
c_{i}^{1}\partial c_{j}^{1}}\right) \left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) \text{d}\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d%
}\boldsymbol{c}^{2},\medskip $ \\ 

$F_{ijk} =-\frac{\sigma ^{2}g_{0}}{10\rho}\int \int \int \Pi \left(
\varphi \right) \left( f_{01}\frac{\partial ^{3}f_{02}}{\partial
c_{i}^{2}\partial c_{j}^{2}\partial c_{k}^{2}}+f_{02}\frac{\partial
^{3}f_{01}}{\partial c_{i}^{1}\partial c_{j}^{1}\partial c_{k}^{1}}\right)
\left( \boldsymbol{g}\cdot \boldsymbol{k}\right) \text{d}%
\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d}\boldsymbol{c}^{2},\medskip $\\
 
$\hat{F}=\frac{\sigma ^{2}g_{0}}{240 \rho}\int \int \int \Pi \left(
\varphi \right) \left( f_{01}\frac{\partial ^{4}f_{02}}{\partial
c_{l}^{2}\partial c_{l}^{2}\partial c_{s}^{2}\partial c_{s}^{2}}+f_{02}\frac{%
\partial ^{4}f_{01}}{\partial c_{l}^{1}\partial c_{l}^{1}\partial
c_{s}^{1}\partial c_{s}^{1}}\right) \left( \boldsymbol{g}\cdot 
\boldsymbol{k}\right) \text{d}\boldsymbol{k}\text{d}\boldsymbol{c}^{1}\text{d%
}\boldsymbol{c}^{2}$%
\end{tabular}
\label{EF}
\end{equation}%
with $f_{01}=f_{0}\left( \boldsymbol{x},\boldsymbol{c}^{1}\right) $ and $%
f_{02}=f_{0}\left( \boldsymbol{x},\boldsymbol{c}^{2}\right) $.

The integrals (\ref{AB})$_{1-4}$ and (\ref{EF})$_{1-4}$ are already been evaluated
in \cite{Jenkins}. 
For completeness, we write here their values together with the new terms due to the presence of the 14th field $\Delta $.

By insertion of the Grad's distribution function (\ref{grad 14}) and by evaluation of well-known Gaussian integrals \cite{CC},
we get, after long computations, the following values for the flux and the source terms of each balance equation,
\begin{equation}
\begin{tabular}{l}

$\theta _{ik}=2(1+e)\upsilon g_{0}\rho \theta \delta _{ik}-\frac{2(1+e)g_{0}\sqrt{\pi \theta }\rho
^{2}d_{p}^{4}}{15m}\left( 2\frac{\partial v_{(i}}{\partial x_{k)}}+\frac{%
\partial v_{l}}{\partial x_{l}}\delta _{ik}\right)+\frac{2(1+e)g_{0}\pi \rho d_{p}^{3}}{15m}\rho _{<ik>},\medskip $ \\ 

$\theta _{ijk}=-\frac{2(1+e)g_{0}\sqrt{\pi \theta }
\rho ^{2}d_{p}^{4}}{5m}\frac{\partial \theta }{\partial x_{(i}}\delta
_{jk)}+\frac{4(1+e)g_{0}\pi \rho d_{p}^{3}
}{75m}\left[ q_{i}\delta _{jk}+\frac{9}{2}q_{(j}\delta _{k)i}\right], \medskip $ \\ 

$\theta _{ikll}=(1+e)\left( 10-3e+3e^{2}\right) \upsilon g_{0}\rho \theta^{2}\delta _{ik}+\medskip $ \\ 
$-\frac{(1+e)g_{0}\sqrt{\pi \theta }\theta
\rho ^{2}d_{p}^{4}}{15m}\left[ \left( 23-3e+8e^{2}\right) \frac{\partial
v_{(i}}{\partial x_{k)}}+\left( 19-4e+4e^{2}\right) \frac{\partial v_{l}}{%
\partial x_{l}}\delta _{ik}\right]\medskip $ \\ 
$+\frac{(1+e)\left(
43-21e+12e^{2}\right) g_{0}\pi \theta \rho d_{p}^{3}}{30m}\rho
_{<ik>}-\frac{(1+e)\left(
16-3e+3e^{2}\right) g_{0}\pi \rho d_{p}^{3}}{180m}\Delta \delta_{ik}, \medskip$\\

$\theta _{kllss}=-\frac{2(1+e)\left( 13-6e+4e^{2}\right)g_{0}\sqrt{\pi \theta }\theta \rho^{2}d_{p}^{4}}{3m} \frac{\partial \theta }{%
\partial x_{k}}+\frac{4(1+e)\left( 26-12e+9e^{2}\right)
g_{0}\pi \theta \rho d_{p}^{3}}{15m}q_{k},$%
\end{tabular}
\label{fluxes3}
\end{equation}
and the final value of the source term,
\begin{equation}
\begin{tabular}{l}
$\psi_i=0, \medskip $\\
$\psi _{ij}=-\frac{4}{3}\frac{\left(1-e^{2}\right)g_{0}\rho ^{2}d_{p}^{2}\sqrt{\pi \theta }\theta}{m} \delta _{ij}\medskip $ \\ 
$+\frac{2\left( 1+e\right) \left( 2-e\right) g_{0}\rho ^{2}d_{p}^{3}\pi \theta}{5m}\frac{\partial v_{(i}}{\partial x_{j)}}+
\frac{\left( 1+e\right) \left( 1-3e\right) g_{0}\rho ^{2}d_{p}^{3}\pi\theta}{15m} \frac{\partial v_{k}}{\partial x_{k}}\delta _{ij}\medskip $ \\ 
$ -\frac{4(1+e)\left( 3-e\right)g_{0}\sqrt{\pi \theta }\rho d_{p}^{2}}{5m}\rho _{<ij>}-\frac{\left( 1-e^{2}\right) g_{0}\rho d_{p}^{2}\sqrt{\pi \theta }}
{60m\theta }\Delta \delta _{ij},\medskip $\\

$\psi _{ill}=\frac{\left( 1+e\right) \left(13-9e\right) g_{0}\rho ^{2}d_{p}^{3}\pi \theta}{6m}\frac{\partial \theta }{\partial x_{i}}
-\frac{2(1+e)\left( 49-33e\right)g_{0}\sqrt{\pi \theta }\rho d_{p}^{2}}{15m}q_{i}, \medskip $ \\
 
$\psi _{llss}=-4\frac{\left( 1-e^{2}\right) \left( 9+2e^{2}\right) g_{0}\rho ^{2}d_{p}^{2}\sqrt{\pi \theta }\theta ^{2}}{m}
+\frac{\left( 1-e^{2}\right) \left( 19+5e^{2}\right)g_{0}\rho ^{2}d_{p}^{3}\pi \theta ^{2}}{2m}\frac{\partial v_{k}}{\partial x_{k}} \medskip $ \\

$-\left[ \frac{\left( 1-e^{2}\right)e^{2}}{2m}+\frac{\left( 1+e\right) \left( 271-207e\right) }{60m}\right]
g_{0}\rho d_{p}^{2}\sqrt{\pi \theta }\Delta.$%
\end{tabular}
\label{sources3}
\end{equation}

Equations (\ref{balance ins}) together with the constitutive relations (\ref{constitutive}) and the relations \eqref{fluxes3} and \eqref{sources3}
form a closed set of 14 field equations for the 14 fields.
These equations are appropriate to a dilute gas whose molecules are subjected to inelasic collisions.

We compared the constitutive relations (\ref{fluxes3}) and (\ref{sources3}) with the 14-moments calculations of Kremer \cite{Kremer2}.
He calculated not all the terms in (\ref{fluxes3}) and used a different Grad's approximation. Anyway, the common terms coincides.
Also in \cite{Chen} some similar calculations are carried out as a second-order moment method. The authors evaluated some terms of 
(\ref{fluxes3}) and (\ref{sources3}) and the corresponding terms coincide.
Gupta et all in \cite{ted} studied the dilute case but they considered more moments.
Also in this case, the corresponding terms coincide.

\section{A numerical application}


In order to show a simple application of the previous calculations,
we study system \eqref{balance ins}, together with the relations \eqref{fluxes3} and \eqref{sources3}, in the one-dimensional stationary case, assuming that the fields depend only on $x$ and the velocity vanishes. 
In this case, the field equations become
\begin{equation}
\begin{tabular}{l}
$\frac{\text{d}}{\text{d}x}\left[ \left( 1+\chi \right) p+\left( 1+\frac{2}{5%
}\chi \right) \rho _{<11>}\right] =0,\medskip $ \\ 
$\frac{\text{d}}{\text{d}x}\left[ \left( 1+\frac{3}{5}\chi \right)
q_{1}-\chi \varphi \frac{\text{d}\theta }{\text{d}x}\right] =-6\frac{\left(
1-e\right) \chi \varphi }{d_{p}^{2}}\theta -\frac{3}{40}\frac{\left(
1-e\right) \chi \varphi }{d_{p}^{2}}\frac{\Delta }{\rho \theta },\medskip $
\\ 
$\frac{8}{15}\frac{\text{d}}{\text{d}x}\left[ \left( 1+\frac{9}{10}\chi
\right) q_{1}-\chi \varphi \frac{\text{d}\theta }{\text{d}x}\right] =-\frac{%
12}{5}\frac{\left( 3-e\right) \chi \varphi }{d_{p}^{2}}\frac{\rho _{<11>}}{%
\rho },\medskip $ \\ 
$\frac{\text{d}}{\text{d}x}\left[ \left( 5+\frac{10-3e+3e^{2}}{2}\chi
\right) p\theta +\left( 7+\frac{43-21e+4e^{2}}{10}\chi \right) \theta \rho
_{<11>}\right. +\medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \left. +\frac{1}{3}\left( 1+\frac{16-3e+3e^{2}}{20}%
\chi \right) \Delta \right] =-\frac{2}{5}\frac{\left( 49-33e\right) \chi
\varphi }{d_{p}^{2}}\frac{q_{1}}{\rho }+\frac{\left( 13-9e\right) \chi }{2}%
\rho \theta \frac{\text{d}\theta }{\text{d}x},\medskip $ \\ 
$\frac{\text{d}}{\text{d}x}\left[ \left( 28+4\frac{26-12e+9e^{2}}{5}\chi
\right) \theta q_{1}-2\left( 16-6e+4e^{2}\right) \chi \varphi \theta \frac{%
\text{d}\theta }{\text{d}x}\right] =\medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ \ =-12\frac{\left( 1-e\right) \left(
9+2e^{2}\right) \chi \varphi }{d_{p}^{2}}\theta ^{2}-\left[ \frac{3\left(
1-e\right) e^{2}}{2}+\frac{271-207e}{20}\right] \frac{\chi \varphi \Delta }{%
d_{p}^{2}\rho }$%
\end{tabular}
\label{1D}
\end{equation}%
with 
\begin{equation}
\begin{tabular}{l}
$\chi =\frac{\pi }{3}\left( 1+e\right) \frac{\rho }{m}d_{p}^{3}g_{0},%
\medskip $ \\ 
$\varphi =\rho d_{p}\sqrt{\frac{\theta }{\pi }}.$%
\end{tabular}
\end{equation}
These quantities depend on the gas under consideration. $\chi $ is an dimensionless parameter
that will be discussed later.
For simplicity, we introduce the following dimensionless quantities 
\begin{equation}
\begin{tabular}{ccccc}
$\hat{x}=\frac{x}{L},$ &
$\hat{p}=\frac{p}{P},$ & 
$\hat{\theta}=\frac{\theta }{\theta _{0}},$ & 
$\hat{\sigma}=\frac{\rho_{<11>}}{P}, $ \medskip \\
$\hat{q}=\frac{q}{P\sqrt{\theta _{0}}},$ & 
$\hat{\Delta}=\frac{\Delta}{P\theta _{0}\sqrt{\theta _{0}}},$ & 
$\lambda=\frac{d_{p}}{\sqrt{\pi }L}$ 
\end{tabular}
\end{equation}
where $L$ is the lenght of the domain, $P$ the boundary pressure and $\theta _{0}$ the boundary temperature. 
So equations (\ref{1D}) become
\begin{equation}
\begin{tabular}{l}
$\frac{\text{d}}{\text{d}\hat{x}}\left[ \left( 1+\chi \right) \hat{p}+\left(
1+\frac{2}{5}\chi \right) \hat{\sigma}\right] =0,\medskip $ \\ 
$\frac{\text{d}}{\text{d}\hat{x}}\left[ \left( 1+\frac{3}{5}\chi \right) 
\hat{q}-\chi \lambda\hat{\rho}\sqrt{\hat{\theta}}\frac{\text{d}\hat{%
\theta}}{\text{d}\hat{x}}\right] =-\frac{3\left( 1-e\right) \chi \hat{\rho}%
\sqrt{\hat{\theta}}}{\pi \lambda}\left[ 2\hat{\theta}+\frac{1}{40}\frac{%
\hat{\Delta}}{\hat{\rho}\hat{\theta}}\right] ,\medskip $ \\ 
$\frac{8}{15}\frac{\text{d}}{\text{d}\hat{x}}\left[ \left( 1+\frac{9}{10}%
\chi \right) \hat{q}-\chi \lambda \hat{\rho}\sqrt{\hat{\theta}}\frac{%
\text{d}\hat{\theta}}{\text{d}\hat{x}}\right] =-\frac{12}{5}\frac{\left(
3-e\right) \chi \sqrt{\hat{\theta}}}{\pi \lambda \hat{\sigma}},\medskip $
\\ 
$\frac{\text{d}}{\text{d}\hat{x}}\left[ \left( 5+\frac{10-3e+3e^{2}}{2}\chi
\right) \hat{p}\hat{\theta}+\left( 7+\frac{43-21e+4e^{2}}{10}\chi \right) 
\hat{\theta}\hat{\sigma}\right. +\medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \left. +\frac{1}{3}\left( 1+\frac{16-3e+3e^{2}}{20}\chi
\right) \hat{\Delta}\right] =-\frac{2}{5}\frac{\left( 49-33e\right) \chi }{%
\pi \lambda}\sqrt{\hat{\theta}}\hat{q}+\frac{\left( 13-9e\right) \chi }{%
2}\hat{\rho}\hat{\theta}\frac{\text{d}\hat{\theta}}{\text{d}\hat{x}}%
,\medskip $ \\ 
$\frac{\text{d}}{\text{d}x}\left[ \left( 28+4\frac{26-12e+9e^{2}}{5}\chi
\right) \hat{\theta}\hat{q}-2\left( 16-6e+4e^{2}\right) \chi \lambda
\hat{p}\sqrt{\hat{\theta}}\frac{\text{d}\hat{\theta}}{\text{d}%
\hat{x}}\right] =\medskip $ \\ 
$\ \ \ \ \ \ \ \ \ \ \ \ \ \ =-\frac{12\left( 1-e\right) \left(
9+2e^{2}\right) \chi \hat{p}\hat{\theta}\sqrt{\hat{\theta}}}{\pi \lambda
}-\left[ \frac{3\left( 1-e\right) e^{2}}{2}+\frac{271-207e}{20}\right] \frac{
\chi \sqrt{\hat{\theta}}\hat{\Delta}}{\pi \lambda}.$
\end{tabular}
\label{1Dad}
\end{equation}

In Fig.1 we show the solution of these field equations obtained with $e=0.9$, $\hat{\theta}_0=1$, 
$\hat{\theta}_1=1.1$,  $\lambda=0.05$ and  $\chi=0.1$. Due to some numerical problems some non-linear effect in the equations are not taken into account. In Fig.1a the solution is shown together with the line connecting the two boundary temperatures.  
The solution shows that the temperature greatly reduces in the center of the domain due to its dissipation in the collision and then it grows due to the boundary temperature. The decrease in temperature  is stronger then in \cite{ETdilute} since here a dense case is considered. The new variable $\Delta$ is proportional to the temperature and it depends on the restitution coefficient $e$. For values of $e \approx 1 $ this field variable tends to vanish, as shown in Fig.1d where the field $\Delta$ is depicted for different labeled values of $e$. 
In Fig.1b the heat flux is shown and its behavior is approximately linear. The stress tensor $\sigma$ is instead approximately constant and, for the chosen values of the constants $\hat{\sigma} \approx 0.0651$ holds. 
This is also an interesting result since in this case the equations herein obtained imply stress tensor components although the gas is at rest. This is a consequence of the inelastic collisions, as also shown in dilute granular gases, but here the effect is stronger since again the gas is more dense.    
\begin{figure}[hhh]
\begin{center}
\includegraphics[width=1 \textwidth]{Fig staz 14 denso.eps}
\end{center}
\caption[]{Solution for the temperature, heat flux and the variable $\hat{\Delta}$ obtained from the field equations 
(\protect\ref{1Dad}) with the values $e=0.9$, $\hat{\theta}_0=1$, 
$\hat{\theta}_1=1.1$,  $\lambda=0.05$ and  $\chi=0.1$. Fig1d: Solution for $\hat{\Delta}$ obtained with different labeled values of the restitution coefficient $e$. }
\label{aaa}
\end{figure}

\section{Concluding remarks}

In this paper, we started from the Grad's 13-moment model for dense and granular gases of Jenkins and Richman \cite{Jenkins} 
and we extended it to a model of 14 moments, adding the 14th fields $\Delta$ and the balance equation for this new variable. 
The $\theta_i$ fluxes, related to the transfer of the particle moments, 
and the $\psi_i$ source terms for each balance equation have been determined analytically, generalizing those obtained by Jenkins and Richman in \cite{Jenkins}.
A particular solution is determined in order to understand, in a particular case, the effect of this new variable $\Delta$. 

A future research perspective would be to treat the model on the basis of continuum thermodynamic theories, 
in particular the Rational Extended Thermodynamics \cite{giallo,Ruggeri Sugyama 1, Ruggeri Sugyama 2}. 
It has been shown that Extended Thermodynamics using universal physical principles is able to recover some results in agreement with the Grad's theory.
It would be interesting to see if this theory is able to obtain the results of this manuscript.

It would also be remarkable to see the applicability of this model to more complex cases. 
It has been shown that the solutions of Grad's theory for rarefied gasses become more interesting 
and different from the Classical Thermodynamic fields when the gas is in curved domains or in particular situations, 
see \cite{helicoidal} and the references therein.   

\begin{thebibliography}{20}

\bibitem{Jenkins} Jenkins, J. T., \& Richman, M. W. (1985). Grad's 13-moment
system for a dense gas of inelastic spheres. Arch. Rat. Mech. Anal, 87(4),
355-377.

\bibitem{Grad1} 
Grad, H.: On the kinetic theory of rarefied gases.
Comm. Pure and Applied Mathematics \textbf{2}, 331-407 (1949).

\bibitem{Grad2} 
Grad, H.: Principles of the Kinetic Theory of Gases. In: Handbuch der Physik \textbf{XII}, pp. 205-294.
Springer, Heidelberg (1958).


\bibitem{Kremer1} Kremer, G.M., Marques Jr, W. (2011). Fourteen moment theory for granular
gases. Kinet. Relat. Models, 4(1), 317-331

\bibitem{Kremer2} Kremer, G.M.  (2019) Fourteen moment method for moderately dense granular
gases, AIP Conference Proceedings 2132, 080002

\bibitem{ted} Gupta V.K., Shukla P., Torrilhon M., (2018). Higher-order moment theories 
for dilute granular gases of smooth hard spheres, J. FluidMech. 836, 451-501.

\bibitem{Brillantov} Brillantov N.V., Pöschel T. 2004, Kinetic Theory of Granular Gases. Oxford
University Press.

\bibitem{Acc1D} Barbera E., Pollino, A. (2023), Mathematical Investigation of 1D Discontinuity
Waves in Dilute Granular Gases, Mathematics 2023, 11, 4935.
https://doi.org/10.3390/math11244935

\bibitem{CC} Chapman, S., Cowling, T. G. (1990). The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases. Cambridge university press.

\bibitem{giallo} M\"{u}ller, I., Ruggeri, T. (2003). Rational Extended
Thermodynamics (Vol. 37). Springer Science \& Business Media.

\bibitem{CS} Carnahan, N. F., Starling, K. E. (1969). Equation of state for nonattracting rigid spheres. The Journal of chemical physics, 51(2), 635-636.

\bibitem{Chen} Chen, J., Wang, S., Sun, D., Lu, H., Gidaspow, D., Yu, H. (2012). 
A second‐order moment method applied to gas–solid risers. AIChE journal, 58(12), 3653-3675.

\bibitem{ETdilute} Barbera E., Pollino A, Rational Extended Thermodynamics of dilute granular gases, submitted (2024).

\bibitem{Ruggeri Sugyama 1} Ruggeri,T., Sugiyama,M.: Rational Extended Thermodynamics Beyond the Monatomic Gas. 
Springer, New York, 2015.

\bibitem{Ruggeri Sugyama 2} Ruggeri,T., Sugiyama,M.:Classical and relativistic rational extended thermodynamics 
of gases. Springer, 2021.

\bibitem{helicoidal} Barbera, E., Brini, F. (2018). Stationary heat transfer in helicoidal flows of a rarefied gas. Europhysics Letters, 120(3), 34001.


\end{thebibliography}

\end{document}