Higher-order moment theories for dilute granular gases of smooth hard-spheres

TL;DR

This paper develops and analyzes higher-order moment equations for dilute granular gases using Grad's method, revealing stability properties and constitutive relations that improve understanding of granular flow behavior.

Contribution

It introduces a 26-moment theory for granular gases and compares its predictions with existing models, providing insights into stability and critical system size.

Findings

Homogeneous cooling state follows Haff's law for temperature.

Higher-order moments decay faster than temperature.

Grad's 26-moment theory predicts a smooth critical wavenumber profile.

Abstract

Grad's method of moments is employed to develop higher-order Grad's moment equations---up to first 26-moments---for granular gases within the framework of the (inelastic) Boltzmann equation. The homogeneous cooling state of a freely cooling granular gas is investigated with the Grad's 26-moment equations in a semi-linearized setting and it is shown that the granular temperature in the homogeneous cooling state still decays according to Haff's law while the other higher-order moments decay on a faster time scale. The constitutive relations for stress and heat flux (the Navier--Stokes and Fourier relations) are obtained by performing a Chapman--Enskog-like expansion on the Grad's 26-moment equations and compared with those existing in the literature. The linear stability of the homogeneous cooling state is analyzed through the Grad's 26-moment system and various sub-systems by decomposing them into longitudinal and transverse systems. It is found that one eigenmode in both longitudinal and transverse systems in case of inelastic gases is unstable. By comparing the eigenmodes from various theories, it is established that the 13-moment eigenmode theory predicts that the unstable eigenmode remains unstable for all wavenumbers below a certain coefficient of restitution while any other higher-order moment theory shows that this mode becomes stable above some critical wavenumber for all values of coefficient of restitution. In particular, the Grad's 26-moment theory leads to a smooth profile for the critical wavenumber in contrast to the other considered theories. Furthermore, the critical system size obtained through the Grad 26-moment and existing theories are also in excellent agreement.

Equations

$$\begin{array}{rl}{\mathbf{c}}^{\prime }& =\mathbf{c}-\frac{1+e}{2}(\hat{\mathbf{k}}\cdot \mathbf{g})\hat{\mathbf{k}},\\ {\mathbf{c}}_1^{\prime }& ={\mathbf{c}}_1+\frac{1+e}{2}(\hat{\mathbf{k}}\cdot \mathbf{g})\hat{\mathbf{k}},\end{array}\}$$

$$\frac{\partial f}{\partial t}+c_i\frac{\partial f}{\partial x_i}+F_i\frac{\partial f}{\partial c_i}=d^2\underset{{\mathbb{R}}^3}{\int }\underset{S^2}{\int }\left(\frac{1}{e^2}{f}^{\prime \prime }{f}_1^{\prime \prime }-f{f}_1\right)(\hat{\mathbf{k}}\cdot \mathbf{g})\Theta (\hat{\mathbf{k}}\cdot \mathbf{g})\mathrm{d}\hat{\mathbf{k}}\mathrm{d}{\mathbf{c}}_1,$$

$$\begin{array}{rl}{\mathbf{c}}^{\prime \prime }& =\mathbf{c}-\frac{1+e}{2e}(\hat{\mathbf{k}}\cdot \mathbf{g})\hat{\mathbf{k}},\\ {\mathbf{c}}_1^{\prime \prime }& ={\mathbf{c}}_1+\frac{1+e}{2e}(\hat{\mathbf{k}}\cdot \mathbf{g})\hat{\mathbf{k}}.\end{array}\}$$

$$n(t,\mathbf{x})$$

$$n(t,\mathbf{x})\mathbf{v}(t,\mathbf{x})$$

$$\frac{3}{2}n(t,\mathbf{x})T(t,\mathbf{x})$$

$$⟨\psi ⟩(t,\mathbf{x})=\frac{1}{n(t,\mathbf{x})}\int \psi (t,\mathbf{x},\mathbf{c})f(t,\mathbf{x},\mathbf{c})\mathrm{d}\mathbf{c}.$$

$$1=⟨1⟩,\mathbf{v}=⟨\mathbf{c}⟩\text{and}T=⟨\frac{1}{3}m{C}^2⟩.$$

$$\frac{\mathrm{D}}{\mathrm{D}t}\int \psi f\mathrm{d}\mathbf{c}+\underline{\frac{\partial }{\partial x_i}\int \psi C_{i}f\mathrm{d}\mathbf{c}}+\frac{\partial v_i}{\partial x_i}\int \psi f\mathrm{d}\mathbf{c}-\int \left(\frac{\mathrm{D}\psi }{\mathrm{D}t}+C_i\frac{\partial \psi }{\partial x_i}+F_i\frac{\partial \psi }{\partial c_i}\right)f\mathrm{d}\mathbf{c}=\mathcal{P}(\psi ),$$

$$\mathcal{P}(\psi )$$

$$\displaystyle=d^{2}\iiint\!\big{(}\psi^{\prime}-\psi\big{)}ff_{1}(\hat{\bm{k}}\cdot\bm{g})\,\Theta(\hat{\bm{k}}\cdot\bm{g})\,\mathrm{d}\hat{\bm{k}}\,\mathrm{d}\bm{c}\,\mathrm{d}\bm{c}_{1},$$

$$\frac{\mathrm{D}n}{\mathrm{D}t}+n\frac{\partial v_i}{\partial x_i}=0,$$

$$\frac{\mathrm{D}{v}_i}{\mathrm{D}t}+\frac{1}{mn}\left[\frac{\partial {\sigma }_{ij}}{\partial x_j}+\frac{\partial (nT)}{\partial x_i}\right]-F_i=0,$$

$$\frac{\mathrm{D}T}{\mathrm{D}t}+\frac{2}{3n}\left[\frac{\partial q_i}{\partial x_i}+{\sigma }_{ij}\frac{\partial v_i}{\partial x_j}+nT\frac{\partial v_i}{\partial x_i}\right]=-\zeta T.$$

$${\sigma }_{ij}=m\int C_{⟨i}{C}_{j⟩}f\mathrm{d}\mathbf{c}\text{and}{q}_i=\frac{1}{2}m\int C^2{C}_{i}f\mathrm{d}\mathbf{c},$$

$$\displaystyle\zeta=-\frac{m\,d^{2}}{6\,n\,T}\iiint\!\Big{(}(C^{\prime})^{2}+(C_{1}^{\prime})^{2}-C^{2}-C_{1}^{2}\Big{)}ff_{1}(\hat{\bm{k}}\cdot\bm{g})\,\Theta(\hat{\bm{k}}\cdot\bm{g})\,\mathrm{d}\hat{\bm{k}}\,\mathrm{d}\bm{c}\,\mathrm{d}\bm{c}_{1}.$$

$$u_{i_1{i}_2\dots i_n}^a=m\int C^{2a}{C}_{⟨i_1}{C}_{i_2}\dots C_{i_n⟩}f\mathrm{d}\mathbf{c},a,n\in {\mathbb{N}}_0.$$

$$u^0=mn=\rho ,u_i^0=0,u^1=3nT=3\rho \theta ,u_{ij}^0={\sigma }_{ij},u_i^1=2q_i,$$

$$m_{ijk}=u_{ijk}^0,R_{ij}=u_{ij}^1-7\theta {\sigma }_{ij},w=u^2-15\rho {\theta }^2,$$

$$\frac{\mathrm{D}{\sigma }_{ij}}{\mathrm{D}t}+\frac{\partial m_{ijk}}{\partial x_k}+\frac{4}{5}\frac{\partial q_{⟨i}}{\partial x_{j⟩}}+{\sigma }_{ij}\frac{\partial v_k}{\partial x_k}+2{\sigma }_{k⟨i}\frac{\partial v_{j⟩}}{\partial x_k}+2\rho \theta \frac{\partial v_{⟨i}}{\partial x_{j⟩}}={\mathcal{P}}_{ij}^0,$$

$$\frac{\mathrm{D}{q}_i}{\mathrm{D}t}+\frac{1}{2}\frac{\partial R_{ij}}{\partial x_j}+\frac{1}{6}\frac{\partial w}{\partial x_i}+\theta \frac{\partial {\sigma }_{ij}}{\partial x_j}+\frac{5}{2}{\sigma }_{ij}\frac{\partial \theta }{\partial x_j}+\frac{5}{2}\rho \theta \frac{\partial \theta }{\partial x_i}+m_{ijk}\frac{\partial v_j}{\partial x_k}$$

$$-{\sigma }_{ij}\frac{1}{\rho }\left(\frac{\partial {\sigma }_{jk}}{\partial x_k}+\theta \frac{\partial \rho }{\partial x_j}\right)+\frac{7}{5}{q}_i\frac{\partial v_j}{\partial x_j}+\frac{7}{5}{q}_j\frac{\partial v_i}{\partial x_j}+\frac{2}{5}{q}_j\frac{\partial v_j}{\partial x_i}=\frac{1}{2}{\mathcal{P}}_i^1,$$

$$\frac{\mathrm{D}{m}_{ijk}}{\mathrm{D}t}+\frac{\partial u_{ijkl}^0}{\partial x_l}+\frac{3}{7}\frac{\partial R_{⟨ij}}{\partial x_{k⟩}}+3\theta \frac{\partial {\sigma }_{⟨ij}}{\partial x_{k⟩}}-3\frac{1}{\rho }{\sigma }_{⟨ij}\left(\frac{\partial {\sigma }_{k⟩l}}{\partial x_l}+\theta \frac{\partial \rho }{\partial x_{k⟩}}\right)$$

$$+m_{ijk}\frac{\partial v_l}{\partial x_l}+3m_{l⟨ij}\frac{\partial v_{k⟩}}{\partial x_l}+\frac{12}{5}{q}_{⟨i}\frac{\partial v_j}{\partial x_{k⟩}}={\mathcal{P}}_{ijk}^0,$$

$$\frac{\mathrm{D}{R}_{ij}}{\mathrm{D}t}+2u_{ijkl}^0\frac{\partial v_k}{\partial x_l}+\frac{\partial u_{ijk}^1}{\partial x_k}+\frac{2}{5}\frac{\partial u_{⟨i}^2}{\partial x_{j⟩}}+R_{ij}\frac{\partial v_k}{\partial x_k}-\frac{28}{5}\theta \frac{\partial q_{⟨i}}{\partial x_{j⟩}}-\frac{28}{5}{q}_{⟨i}\frac{\partial \theta }{\partial x_{j⟩}}$$

$$+4\theta {\sigma }_{k⟨i}\frac{\partial v_k}{\partial x_{j⟩}}+4\theta {\sigma }_{k⟨i}\frac{\partial v_{j⟩}}{\partial x_k}-\frac{8}{3}\theta {\sigma }_{ij}\frac{\partial v_k}{\partial x_k}-\frac{14}{3}\frac{1}{\rho }{\sigma }_{ij}\frac{\partial q_k}{\partial x_k}-\frac{14}{3}\frac{1}{\rho }{\sigma }_{ij}{\sigma }_{kl}\frac{\partial v_k}{\partial x_l}$$

$$-7\theta \frac{\partial m_{ijk}}{\partial x_k}+\frac{6}{7}{R}_{⟨ij}\frac{\partial v_{k⟩}}{\partial x_k}+\frac{4}{5}{R}_{k⟨i}\frac{\partial v_k}{\partial x_{j⟩}}+2R_{k⟨i}\frac{\partial v_{j⟩}}{\partial x_k}+\frac{14}{15}w\frac{\partial v_{⟨i}}{\partial x_{j⟩}}$$

$$-2m_{ijk}\frac{\partial \theta }{\partial x_k}-2\frac{1}{\rho }{m}_{ijk}\left(\frac{\partial {\sigma }_{kl}}{\partial x_l}+\theta \frac{\partial \rho }{\partial x_k}\right)-\frac{28}{5}\frac{1}{\rho }{q}_{⟨i}\left(\frac{\partial {\sigma }_{j⟩k}}{\partial x_k}+\theta \frac{\partial \rho }{\partial x_{j⟩}}\right)$$

$$={\mathcal{P}}_{ij}^1-7\theta {\mathcal{P}}_{ij}^0-\frac{14}{3}\frac{1}{\rho }{\sigma }_{ij}\left(\frac{1}{2}{\mathcal{P}}^1\right),$$

$$\frac{\mathrm{D}w}{\mathrm{D}t}+\frac{\partial u_i^2}{\partial x_i}-20\theta \frac{\partial q_i}{\partial x_i}-8q_i\frac{\partial \theta }{\partial x_i}+\frac{7}{3}w\frac{\partial v_i}{\partial x_i}+4R_{ij}\frac{\partial v_i}{\partial x_j}$$

$$+8\theta {\sigma }_{ij}\frac{\partial v_i}{\partial x_j}-8\frac{1}{\rho }{q}_i\left(\frac{\partial {\sigma }_{ij}}{\partial x_j}+\theta \frac{\partial \rho }{\partial x_i}\right)={\mathcal{P}}^2-20\theta \left(\frac{1}{2}{\mathcal{P}}^1\right),$$

$${\mathcal{P}}_{i_1{i}_2\dots i_n}^a=m{d}^2\iiint$$

$$\times f{f}_1(\hat{\mathbf{k}}\cdot \mathbf{g})\Theta (\hat{\mathbf{k}}\cdot \mathbf{g})\mathrm{d}\hat{\mathbf{k}}\mathrm{d}\mathbf{c}\mathrm{d}{\mathbf{c}}_1$$

$$f_M\equiv f_M(t,\mathbf{x},\mathbf{c})=n{\left(\frac{1}{2\pi \theta }\right)}^{3/2}\exp \left(-\frac{C^2}{2\theta }\right),$$

$$f_{|\text{G26}}$$

$$+\frac{1}{28}\frac{R_{ij}}{\rho {\theta }^3}{C}_i{C}_j\left(\frac{C^2}{\theta }-7\right)+\frac{1}{8}\frac{w}{\rho {\theta }^2}\left(1-\frac{2}{3}\frac{C^2}{\theta }+\frac{1}{15}\frac{C^4}{{\theta }^2}\right)].$$

Full text

1 Introduction

A conglomeration of discrete macroscopic particles characterized by dissipative collisions is termed as a granular material. Granular materials are prevalent in various industries—for instance, in chemical, agriculture and food industries—as well as in nature—for instance, in asteroid belt, sand dunes, debris, etc. Under substantially strong driving forces (e.g., vibration, shearing, etc.), granular materials are in rapid flow regime, in which they exhibit fluid-like behavior and often referred to as granular gases (Campbell, 1990; Goldhirsch, 2003). In the rapid flow regime, the particles of a granular material move randomly and, similar to a monatomic gas, it can be assumed that the collisions among them are binary and instantaneous. However, in contrast to a monatomic gas, the collisions between any two particles of a granular material are inherently inelastic, and thereby energy is dissipated during each collision. The inelastic collisions in granular materials lead to many interesting phenomena—for instance, standing wave patterns (Melo et al., 1995; Umbanhowar et al., 1996), clustering (Kudrolli et al., 1997), fingering (Pouliquen et al., 1997), mixing and segregation (Ottino & Khakhar, 2000; Mullin, 2000; Breu et al., 2003), shear banding (Mueth et al., 2000; Shukla & Alam, 2009), jamming (Corwin et al., 2005), density waves (Liss et al., 2002; Alam et al., 2009), etc. Nevertheless, the mechanics of granular material is still not well-understood, although there have been significant developments in the last couple of decades.

The analogy between granular and molecular fluids has motivated several researchers to devise theoretical methods—for studying granular fluids—based on kinetic theory within the framework of the Boltzmann equation, see e.g., Jenkins & Richman (1985a, b); Goldshtein & Shapiro (1995); Brey et al. (1996); Sela & Goldhirsch (1998); Brey et al. (1998a); Garzó & Santos (2003); Brilliantov & Pöschel (2004); Bisi et al. (2004); Kremer et al. (2014). Kinetic theory for granular flows was first introduced in the two seminal papers by Jenkins & Savage (1983) and Lun et al. (1984) where the former deals with the theory for nearly elastic granular flows while the latter with arbitrary inelasticity. The reported works on the kinetic theory for granular fluids also attempted to extend the two well-known approximation methods in kinetic theory for monatomic gases, namely the Chapman–Enskog (CE) expansion (Chapman & Cowling, 1970) and Grad’s method of moments (Grad, 1949b). Despite the fact that the higher-order approximations (Burnett and beyond) resulting from the CE expansion lead to unstable equations while the Grad’s method of moments always yields linearly stable sets of equations in the case of monatomic gases (Bobylev, 1982), the former has extensively been studied even for granular gases. Goldshtein & Shapiro (1995) applied the CE expansion to obtain Euler-like hydrodynamic equations for rough granular flows. The papers by Brey et al. (1998a) and Garzó & Dufty (1999) apply the CE expansion to dilute and dense granular gases, respectively, to derive the first order (Navier–Stokes and Fourier) constitutive relations for the system of mass, momentum and energy balance equations. Both the works assume that the space and time dependence of the distribution function can be expressed completely in terms of the hydrodynamic fields and represent the distribution function in a formal series of a uniformity parameter which measures the strength of spatial gradients of the hydrodynamic fields. Sela et al. (1996) performed a generalized expansion on the distribution function in powers of two small parameters, namely the Knudsen number and the degree of inelasticity $\epsilon=1-e^{2}$ with $e$ being the coefficient of (normal) restitution and determined the constitutive relations up to Burnett order for the system of mass, momentum and energy balance equations specialized to a simply sheared two-dimensional hard-sphere granular gas via the CE expansion. Subsequently, with the same method, Sela & Goldhirsch (1998) (see also Gupta, 2011) determined the constitutive relations up to Burnett order for a smooth hard-sphere granular gas in three dimensions. However, both these works restrict the coefficient of restitution for being very close to $1$ in order to get the zeroth order solution in the expansion as Maxwellian, similar to a monatomic gas. Lutsko (2005) employed the CE expansion to obtain the constitutive relations for dense granular gases with arbitrary energy loss models. His model produced the homogeneous cooling state (HCS) solution at the zeroth order in the expansion while Navier–Stokes and Fourier constitutive relations at first order. Following the approach of expanding the distribution function in powers of uniformity parameter, Khalil et al. (2014) recently obtained the Burnett order constitutive relations for a smooth granular gas described with (inelastic) Maxwell interaction potential via the CE expansion. It is worth to point out that, unlike Sela & Goldhirsch (1998), their approach does not restrict the coefficient of restitution for being close to $1$. Recently, Kremer et al. (2014) employed the CE expansion to obtain the transport coefficients up to Navier–Stokes order for a granular gas of rough hard-spheres.

Despite the success of the Grad’s method of moments over the CE expansion in the case of monatomic gases (see e.g., Torrilhon, 2016, and references therein), the former in the case of granular flows has received much less attention than the latter which has been investigated extensively. Nevertheless, over the last few years, researchers have shown interest in exploring the former as well (see e.g., Kremer & Marques Jr., 2011; Garzó, 2013; Saha & Alam, 2014). The pioneering work on Grad’s method of moments for granular flows is due to Jenkins & Richman who extended the method to granular flows and derived the Grad 13-moment (G13) equations for a dense and smooth granular gas in Jenkins & Richman (1985a) and, subsequently, the Grad 16-moment equations for a dense and rough two-dimensional granular gas in Jenkins & Richman (1985b)—the extra field variables in the latter are due to rotational motion of rough particles on collision. Bisi et al. (2004) further extended the Grad’s method of moments to weakly inelastic granular flows with variable coefficient of restitution, in one-dimension though. As it is well-established that the scalar fourth moment is necessary for the proper description of a granular gas, Risso & Cordero (2002) included the scalar fourth moment and—for a granular gas—obtained Grad moment equations for $9$ hydrodynamic fields in two dimensions (corresponding to $14$ field variables in three dimensions), and applied them to investigate two problems, namely the HCS and a granular system steadily heated by two parallel walls. Kremer & Marques Jr. (2011) presented the Grad 14-moment (G14) theory for dilute granular gases and by employing this theory, they investigated the HCS, found the constitutive relations for five moment equations, and performed eigenmode analysis on thirteen field theory, where they assumed a constant value for the scalar fourth moment. Garzó (2013) demonstrated the G14 theory for moderately dense granular flows, although the work aimed at finding the Navier–Stokes level constitutive relations through the G14 equations. Recently, Saha & Alam (2014) developed the Grad’s method of moments based on anisotropic Gaussian and employed it to investigate the non-Newtonian stress, collisional dissipation and heat flux in a sheared two-dimensional granular flow.

The main reasons—among others—that the Grad’s method of moments for granular flows is receiving less attention could be that $(i)$ the production terms, which emanate through the Boltzmann collision integral, in the moment equations are very difficult to evaluate for a general interaction potential, $(ii)$ how many and which moments one should consider for describing a process is not known a priori, and $(iii)$ the boundary conditions associated with the Grad moment equations are unclear, which is apparently also the case for higher-order equations resulting from the CE expansion. Over the last couple of years, the two of the present authors have developed a computational methodology, which can compute the production terms associated with the Grad moment equations for $(i)$ a monatomic gas, $(ii)$ a mixture of monatomic gases, and $(iii)$ a dilute granular gas, all interacting with a general interaction potential (Gupta & Torrilhon, 2012, 2015a; Gupta, 2015). The order of magnitude method developed by Struchtrup (2004, 2005)—which regularizes the original G13 equations—identifies the required moments for describing a process in rarefied monatomic gases in a systematic way, although the method requires the Grad 26-moment (G26) equations which at present are not available for granular gases. The boundary conditions associated with the Grad moment equations for monatomic gases are typically obtained using the Maxwell accommodation model (Maxwell, 1879); however, those for granular gases are yet to be explored.

Within moment approximations the number of moments is crucial and has a strong influence on the results, especially when considering only few moments. It is known form studies in monatomic gases that, while moment approximations do converge to the solution of the Boltzmann equation for large number of moments, when using small systems oscillatory convergence pattern makes it difficult to judge the approximation quality, see Torrilhon (2015). Consequently, when using moment approximations, one should always compare a set of successive approximations with different number of moments. The result of a 14-moment theory alone is not conclusive and must be cross-checked with larger systems. A similar argument applies to using linearized models, where it is important to check the possible influence of non-linearity even if the outcome confirms its negligibility. Another important reason to consider moment equations beyond 14 moments is the long-term perspective of establishing a complete set of predictive moment equations for granular flows. For monatomic gases, the regularized 13-moment (R13) equations provide such a model. The R13 system is derived from the G26 model and only on that level relevant phenomena like Knudsen layers or non-gradient transport can be predicted.

In this paper, we derive the G26 equations for dilute granular gases through the Grad’s method of moments, although we neither consider the regularization of the Grad moment equations nor the required boundary conditions—these problems will be considered elsewhere in future. The fully nonlinear production terms associated with an even larger Grad moment system—which includes the fully contracted sixth moment as a field variable along with the $26$ moments and is referred to as the system of Grad 27-moment (G27) equations here—are computed with computer algebra software Mathematica® and presented in appendix A for hard-sphere interaction potential. The production terms for the G26 equations can readily be found from those for the G27 equations by simply discarding the terms containing the fully contracted sixth moment. Here, it is noteworthy to point out that the present work does not have any restriction on the coefficient of restitution other than being a constant and is, consequently, expected to be applicable not only for nearly elastic granular gases but also for those having large inelasticity. With the semi-linearized G26 and G27 equations, we investigate the HCS—which has also been studied a lot theoretically, numerically as well as experimentally—of a freely cooling granular gas to show that the decay of granular temperature closely follows Haff’s law (Haff, 1983) while the other higher-order moments relax on a faster time scale than the granular temperature. We further study the effects of nonlinear terms of scalar fourth moment and that of linear terms of scalar sixth moment on Haff’s law. Following the approach of Garzó (2013), we determine the constitutive relations for the stress and heat flux for dilute granular gases through the G26 equations and compare them with those obtained via CE expansion and via computer simulations. We further investigate the linear stability of HCS by scrutinizing the eigenmodes of longitudinal and transverse problems associated with the G26 system and other sub-systems. Similar problems related to the stability of eigenmodes have been investigated earlier by Brey et al. (1996) for five moment theory and by Kremer & Marques Jr. (2011) for thirteen moment theory. We find that one eigenmode from each Grad moment theory considered in this paper is unstable for inelastic gases while others are stable, while all the eigenmodes are stable for monatomic (elastic) gases. Moreover, the thirteen moment theory of Kremer & Marques Jr. (2011) shows that the unstable heat mode of the longitudinal system remains unstable below a certain value of the coefficient of restitution for all wavenumbers while the present work reveals that the unstable heat mode of the longitudinal system associated with all moment theories expect for the G13 theory becomes stable above a critical wavenumber for all values of the coefficient of restitution. By investigating the critical wavenumbers of the longitudinal and transverse systems from the G26 theory, we also study the critical system size, beyond which the system becomes unstable, and compare the results with the existing theories and simulations. The findings of the paper will be useful in better understanding of granular gases, in developing new mathematical models—such as regularized moment equations—and boundary conditions, and in capturing some intriguing features of granular gases theoretically.

The rest of the paper is structured as follows. Kinetic theory for granular gases is briefly reviewed in § 2. Grad’s method of moments is outlined and applied to derive the Grad moment systems of various orders—in particular, the system of G26 equations—for granular gases in § 3. The HCS of a freely cooling granular gas is investigated through the G26 and G27 equations in § 4 in order to discern the effects of nonlinear terms of scalar fourth moment and of linear terms of scalar sixth moment on HCS. The constitutive relations for the stress and heat flux are computed through the G26 equations in § 5. The linear stability of HCS is analyzed in § 6 through the G26 equations by decomposing them into longitudinal and transverse systems. The conclusions of the paper are given in § 7.

2 Short review of kinetic theory

We consider a dilute granular gas composed of smooth-inelastic-identical hard-spheres of mass $m$ and diameter $d$. The binary collision between two such spheres having pre-collisional velocities $\bm{c}$ and $\bm{c}_{1}$ leads to the following velocity transformation after the collision (see e.g., Brilliantov & Pöschel, 2004; Rao & Nott, 2008):

[TABLE]

where $\bm{c}^{\prime}$ and $\bm{c}_{1}^{\prime}$ are the post-collisional velocities of the respective spheres, $\bm{g}=\bm{c}-\bm{c}_{1}$ is the relative velocity, $\hat{\bm{k}}$ is the unit vector directed from the center of one sphere to that of other at the time of collision, and $e$ is the coefficient of normal restitution (also, referred to as the coefficient of restitution). The coefficient of restitution $e$, in principle, is not a constant and usually depends on impact velocity (Bizon et al., 1999; Brilliantov & Pöschel, 2004). Nevertheless, for simplicity, we assume that $e$ is constant with $0\leq e\leq 1$ in this work. The two limiting cases of $e=0$ and $e=1$ correspond to sticky and perfectly elastic collisions, respectively.

In kinetic theory, the state of a dilute granular gas can be described with single-particle velocity distribution function $f\equiv f(t,\bm{x},\bm{c})$ which obeys the (inelastic) Boltzmann equation (Goldshtein & Shapiro, 1995; Sela & Goldhirsch, 1998; Brilliantov & Pöschel, 2004)

[TABLE]

where $\Theta$ is the Heaviside step function, $\bm{F}$ is the external force per unit mass and usually do not depend on $\bm{c}$, $f_{1}\equiv f(t,\bm{x},\bm{c}_{1})$, $f^{\prime\prime}\equiv f(t,\bm{x},\bm{c}^{\prime\prime})$, $f_{1}^{\prime\prime}\equiv f(t,\bm{x},\bm{c}_{1}^{\prime\prime})$, the integration limits of $\hat{\bm{k}}$ extend over the unit sphere $S^{2}$ and Einstein summation convention is assumed over repeated indices. The velocities $\bm{c}^{\prime\prime}$ and $\bm{c}_{1}^{\prime\prime}$ are the pre-collisional velocities in an inverse collision and are related to the $\bm{c}$ and $\bm{c}_{1}$ as follows (see e.g., Brilliantov & Pöschel, 2004; Rao & Nott, 2008).

[TABLE]

Owing to brevity, hereafter, we shall omit the limits of integrations. Nevertheless, the integration over any velocity space will stand for the volume integral over $\mathbb{R}^{3}$ and that over $\hat{\bm{k}}$ will stand for the volume integral over the unit sphere $S^{2}$. The right-hand side (RHS) of (2) is referred to as the (inelastic) Boltzmann collision operator.

The hydrodynamic variables—number density $n\equiv n(t,\bm{x})$, macroscopic velocity $\bm{v}\equiv\bm{v}(t,\bm{x})$, granular temperature $T\equiv T(t,\bm{x})$—are directly related to the velocity distribution function as follows.

[TABLE]

where $\bm{C}(t,\bm{x},\bm{c})=\bm{c}-\bm{v}(t,\bm{x})$ is the peculiar velocity. Here, the definition of the granular temperature—which is a measure of fluctuating kinetic energy—is adopted following the references Brilliantov & Pöschel (2004); Garzó (2013), although some authors—for instance, Lun et al. (1984); Campbell (1990); Saha & Alam (2014)—also refer $\theta=T/m$ as granular temperature (defined as fluctuating kinetic energy per unit mass). It should be noted that in the case of monatomic gases, $T$ in (6) is replaced with $k_{B}T_{\textrm{th}}$, where $k_{B}$ denotes the Boltzmann constant and $T_{\textrm{th}}$ is the thermodynamic temperature.

The governing equations for the hydrodynamic variables can be derived from the Boltzmann equation (2) as follows. Let us first introduce the moments. For a particle property, $\psi\equiv\psi(t,\bm{x},\bm{c})$, its average value $\langle\psi\rangle\equiv\langle\psi\rangle(t,\bm{x})$ is defined in terms of the distribution function $f$ as

[TABLE]

The integral $\int\!\psi f\,\mathrm{d}\bm{c}$ on the RHS of (7) is referred to as the moment of the velocity distribution function with respect to $\psi$. Clearly,

[TABLE]

The governing equation for the moment $\int\!\psi f\,\mathrm{d}\bm{c}$—often, referred to as the moment equation for $n\langle\psi\rangle$ or the transfer equation for property $\psi$—is obtained by multiplying the Boltzmann equation (2) with $\psi$ and integrating over the velocity space $\bm{c}$. The transfer equation for the quantity $\psi$ reads

[TABLE]

where $\frac{\mathrm{D}}{\mathrm{D}t}\equiv\frac{\partial}{\partial t}+\bm{v}\cdot\bm{\nabla}$ is the material derivative, the underline denotes the flux term, and

[TABLE]

with $\psi^{\prime}\equiv\psi(t,\bm{x},\bm{c}^{\prime})$ etc., is the rate of change of $\langle\psi\rangle$ per unit volume due to collisions and referred to as the production term or the (inelastic) Boltzmann collision integral corresponding to moment $\int\!\psi f\,\mathrm{d}\bm{c}$. While writing (2), the symmetry properties of the Boltzmann collision operator (Brilliantov & Pöschel, 2004; Rao & Nott, 2008) have been employed. It should be noted from the transfer equation (9) that irrespective of the value of $\psi$ chosen for defining a moment by (7), the transfer equation (9) will always contain an additional moment of one more order in the flux term. On substituting $\psi$ in (9) with $1$, $c_{i}$ and $\frac{1}{3}m\,C^{2}$ successively, we get the mass, momentum and energy balance equations, respectively, which read

[TABLE]

[TABLE]

[TABLE]

In (12) and (13), $\sigma_{ij}\equiv\sigma_{ij}(t,\bm{x})$ and $q_{i}\equiv q_{i}(t,\bm{x})$ are the component of stress tensor and heat flux, respectively, and—following the notations of Struchtrup (2005)—these are given by

[TABLE]

where the angle brackets around the indices denote the symmetric and traceless part of the tensor (Struchtrup, 2005), and $\zeta$ in (13) is the cooling rate due to inelastic collisions:

[TABLE]

Note that the RHSs of (11) and (12) vanish due to conservation of mass and momentum; however, the energy is not conserved due to dissipative collisions resulting into non vanishing RHS in (13). It is worth pointing out that the stress and heat flux appearing in (12) and (13)—in general—have non-vanishing collisional contributions in addition to the usual kinetic contribution given by (14) (Brey et al., 1997; Garzó & Dufty, 1999; Rao & Nott, 2008; Garzó, 2013). However, for dilute granular gases, the kinetic contributions to the stress and heat flux dominate over their respective collisional transfer contributions, and therefore the collisional transfer contributions to the stress and heat flux can be neglected for dilute granular gases (Brilliantov & Pöschel, 2004). In any case, since (12) and (13) are obtained here by taking the velocity moments of the distribution function—which is essentially the idea of Grad’s method of moments too, only kinetic contributions to the stress and heat flux will emerge in the equations. The collisional contributions to stress and heat flux, in the case of dense granular gases, must be computed separately (Garzó, 2013).

Obviously, the system of equations (11)–(13) in hydrodynamic variable $n$, $v_{i}$ and $T$ is not closed, since it contains the unknowns $\sigma_{ij}$, $q_{i}$ and $\zeta$, and in order to deal with this system further one must close it. Here, we employ Grad’s method of moments (Grad, 1949b) in order to obtain a closed system of equations.

3 Grad’s method of moments

The function $\psi$ in (7) can be chosen in infinitely many ways. Therefore, the transfer equation (9) leads to an infinite hierarchy of moment equations. However, in practice, only a finite number of moment equations—obtained by truncating the infinite hierarchy of moment equations at a certain level—are used. Nonetheless, this truncated system of equation is not closed due to the (underlined) flux term in (9). In order to obtain a closed (finite) system of moment equations, Grad (1949b) expanded the velocity distribution function $f$ in a finite linear combination of the $N$-dimensional Hermite polynomials (Grad, 1949a) in peculiar velocity and computed the unknown coefficients in the expansion in terms of the considered moments by satisfying their definitions with the approximated distribution function. This method of obtaining a closed set of moment equations is referred to as Grad’s method of moments and its details can be found in Grad (1949b) and in many standard textbooks, see e.g., Struchtrup (2005); Kremer (2010).

On including the governing equations for stress ($\sigma_{ij}$) and heat flux ($q_{i}$) into the system of mass, momentum and energy balance equations (11)–(13), one obtains the well-known 13-moment equations (in three dimensions). Here, we want to derive and explore the G26 equations for dilute granular flows. To this end, let us first introduce the general form of a moment. The typical form of $\psi$ is $m\,C^{2a}C_{\langle i_{1}}C_{i_{2}}\dots C_{i_{n}\rangle}$, where $a,n\in\mathbb{N}_{0}$ and the angle brackets around the indices again denote the symmetric and traceless part of the corresponding quantity (Struchtrup, 2005). Thus, a general moment is given by

[TABLE]

Clearly,

[TABLE]

where $\rho=m\,n$ is the mass density and $\theta=T/m$.

3.1 26-moment equations

In addition to the well-known 13 moments, the system of 26-moment equations include full third rank tensor and one- and full-traces of the fourth moment (i.e., the 26-moment equations include the moments $n$, $v_{i}$, $T$, $\sigma_{ij}$, $q_{i}$, $u_{ijk}^{0}$, $u_{ij}^{1}$, $u^{2}$). The system of G26 equations is obtained by substituting $\psi$ in (9) with $1$, $c_{i}$, $\frac{1}{3}m\,C^{2}$, $m\,C_{\langle i}C_{j\rangle}$, $\frac{1}{2}m\,C^{2}C_{i}$, $m\,C_{\langle i}C_{j}C_{k\rangle}$, $m\,C^{2}C_{\langle i}C_{j\rangle}$ and $m\,C^{4}$ successively. The system of 26-moment equations consists of the mass, momentum and energy balance equations (11)–(13) and other higher-order moment equations, which on using the abbreviations

[TABLE]

read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

are the production terms. Notice that $\mathcal{P}^{0}=\mathcal{P}_{i}^{0}=0$ due to mass and momentum conservation, and $\mathcal{P}^{1}=-3\,n\,T\,\zeta$.

In (18), $w$ has been defined as the difference of the full trace of fourth moment to its value computed with the Maxwellian distribution function

[TABLE]

and $R_{ij}$ is defined as the difference of the one trace of fourth moment to its value computed with the Grad 13-moment (G13) distribution function so that the quantities $m_{ijk}$, $R_{ij}$ and $w$ vanish for the G13 theory. It may be noted that the value of the full trace of fourth moment $u^{2}$ when computed either with Maxwellian distribution function or with the G13 distribution function is $15\rho\theta^{2}$.

3.2 Grad 26-moment closure

Clearly, the system of 26-moment equations (eqs. (11)–(13) and (19)–(3.1)) is not closed, since it contains the unknown higher-order moments $u_{ijkl}^{0}$, $u_{i}^{2}$ and $u_{ijk}^{1}$, and on top of that the production terms are also not known. The system is closed with the Grad distribution function based on the 26 moments considered, which is referred to as the G26 distribution function and reads

[TABLE]

Insertion of the G26 distribution function (3.2) into the definitions of higher-order moments $u_{ijkl}^{0}$, $u_{i}^{2}$ and $u_{ijk}^{1}$, and into the production terms (3.1) expresses them in terms of the considered 26 moments. The unknown higher-order moments $u_{ijkl}^{0}$, $u_{i}^{2}$ and $u_{ijk}^{1}$ turn into

[TABLE]

where the subscript “$|\textrm{G26}$” just denotes that these moments are evaluated with the G26 distribution function (3.2). Unfortunately, the production terms (3.1) are still not easy to evaluate by hand. The two authors of the present paper implemented the strategy for computing the production terms into computer algebra software Mathematica® and obtained the fully nonlinear production terms associated with the G26 equations for dilute granular gases of smooth hard-spheres. The details of the computation can be found in Gupta & Torrilhon (2012) and the source code for the computation is provided as supplementary material with the present paper. Interested readers are also referred to Gupta & Torrilhon (2015b); Gupta (2015) for deriving the higher-order moment equations for monatomic gas mixtures and to Gupta & Torrilhon (2015a); Gupta (2015) for learning the computation of their associated production terms, which might be useful in developing higher-order moment theories for granular gas mixtures. For the sake of completeness, we provide the production terms associated with the G26 equations (eqs. (11)–(13) and (19)–(3.1)) in appendix A; the production terms associated with the G26 equations are (A)–(A) on taking $\Xi=0$. Notice from (A) that the cooling rate computed via the G26 equations is given by

[TABLE]

where

[TABLE]

is the collision frequency. From (3.2), one can see that the cooling rate computed via the G26 equations is also proportional to $(1-e^{2})$ as obtained in previous studies, e.g., Sela & Goldhirsch (1998); van Noije & Ernst (1998); Brilliantov & Pöschel (2004); Kremer & Marques Jr. (2011), and it vanishes identically for monatomic gases ($e=1$) resulting into the conservation of energy. Furthermore, on dropping the underlined nonlinear terms in (15), the expression for the dissipation matches with that in Kremer & Marques Jr. (2011).

3.3 Various Grad moment systems

Equations (11)–(13) and (19)–(3.1) along with (A)–(A) form the system of G26 equations. Various Grad moment systems may be obtained from the system of G26 equations.

The system of well-known G13 equations contains the balance equations for variables $n,v_{i},T,\sigma_{ij},q_{i}$, i.e., the system of G13 equations consists of equations (11)–(13), (19) and (3.1) along with (A)–(A) and $m_{ijk}=R_{ij}=w=0$. 2. 2.

The system of G14 equations contains the balance equations for variables $n,v_{i},T,\sigma_{ij},q_{i},w$, i.e., the system of G14 equations includes equations (11)–(13), (19), (3.1) and (3.1) along with (A)–(A), (A) and $m_{ijk}=R_{ij}=0$. 3. 3.

The system of G20 equations contains the balance equations for variables $n,v_{i},T,\sigma_{ij},q_{i},m_{ijk}$, i.e., the system of G20 equations includes equations (11)–(13) and (19), (3.1), (3.1) along with (A)–(A) and $R_{ij}=w=0$. 4. 4.

The system of G21 equations contains the balance equations for variables $n,v_{i},T,\sigma_{ij},q_{i},m_{ijk},w$, i.e., the system of G21 equations includes equations (11)–(13), (19)–(3.1) and (3.1) along with (A)–(A), (A) and $R_{ij}=0$.

It is commonly accepted that the dependence of heat flux on the density gradient in addition to the temperature gradient is akin to inclusion of the full trace of the fourth moment into the moment system (Kremer & Marques Jr., 2011; Garzó, 2013). In that sense, the G13 and G20 theories for granular flows may not lead to meaningful results. Nevertheless, we shall also include them in this study for comparison purposes.

4 Homogeneous cooling state of a freely cooling granular gas

The state of a granular gas when in the absence of any external forces its granular temperature decays continuously but its spatial homogeneity is maintained is termed as the homogeneous cooling state (Brilliantov & Pöschel, 2004). In this section, we study the HCS of a granular gas with Grad moment equations. We assume a spatially homogeneous state (i.e., $\frac{\partial(\cdot)}{\partial x_{i}}=0$) without any external force acting on the particles of the granular gas (i.e., $\bm{F}=\bm{0}$). Furthermore, following Kremer & Marques Jr. (2011), in the production terms (eqs. (A)–(A)), we shall retain only those nonlinear terms which are the product of the scalar fourth moment ($w$) and a vectorial or tensorial moment, and all other nonlinear terms are simply ignored. This means that we are focusing our attention on the early evolution stage of homogeneously cooling granular gas.

Often, It is more convenient to use a dimensionless variable $\Delta:=\frac{w}{\rho\,\theta^{2}}$ instead of the field variable $w$ in the Grad moment systems, where the governing equations for the former can be obtained by (11), (13) and (3.1). Therefore, in the following, we shall use this new variable $\Delta$ while writing the Grad moment system. The G26 equations in the zero external force case and with aforementioned simplification of the production terms—in terms of $\Delta$—read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

Notice that (4) contains $\Delta^{2}$ term, even though we have dropped all the nonlinear terms including $\Delta^{2}$ terms in the RHSs of (32) and (3.1) owing to aforementioned simplification.

We now introduce the following scaling

[TABLE]

With the scaling (42) and noting that $\Delta=w_{\ast}/(n_{\ast}T_{\ast}^{2})$, the G26 equations in the HCS (i.e., $\frac{\partial(\cdot)}{\partial x_{i}}=0$, $\bm{F}=\bm{0}$) reduce to

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

4.1 Haff’s law

By following heuristic approach, Haff (1983) discovered that in a freely cooling granular gas (with constant coefficient of restitution), the decay rate of granular temperature is given by

[TABLE]

where $\bar{n}$ is the average number density. The solution of (51) is given by

[TABLE]

where $\tau_{\circ}^{-1}\propto\bar{n}\,d^{2}(1-e^{2})\sqrt{T_{\circ}}$ is an inverse time scale, see e.g., Haff (1983); Brilliantov & Pöschel (2004). Equation (52) is termed as Haff’s law for evolution of the granular temperature of a freely cooling granular gas (with a constant coefficient of restitution).

On comparing (45) and (51), one readily perceives that Haff’s law can be obtained from (45), provided $\Delta$ is constant or, in other words,

[TABLE]

With a constant value of $\Delta$ ($=\alpha$, let us say)—which is obtained from condition (53)—equations (43) and (45) along with the initial conditions $n_{\ast}(0)=T_{\ast}(0)=1$ yield Haff’s law for evolution of the dimensionless granular temperature $T_{\ast}$:

[TABLE]

where

[TABLE]

is the inverse of a (dimensionless) time scale $\tau_{0}$.

It is worthwhile to note that, in the light of condition (53), the constant $\alpha$ in (55) is an equilibrium point (Strogatz, 1994) of differential equation (50). Here, we shall discard the underlined term (proportional to $\Delta^{2}$) in (50), although we shall analyze the effect of this term in § 4.3. On discarding the underlined term in (50), one readily obtains the equilibrium point—which is the constant $\alpha$ in (55)—as

[TABLE]

is same as the coefficient of second Sonine polynomial $S_{2}(v^{2})$ while performing the CE expansion on the inelastic Boltzmann equation, see e.g., van Noije & Ernst (1998); Brilliantov & Pöschel (2004).

4.2 Relaxation of moments in the homogeneous cooling state

It is clear from (43) and (44) that the number density $n_{\ast}$ remains constant while the macroscopic velocity $v_{i}^{\ast}$ remains zero (since there is no macroscopic velocity in the initial state) during the homogeneous cooling. Furthermore, equations (46)–(50) are coupled and, therefore need to be solved numerically for further analysis. It may also be noticed from the structure of (47) and (48) that vanishing initial conditions on $q_{i}^{\ast}$ and $m_{ijk}^{\ast}$ will result into vanishing solution for these quantities. Similarly, vanishing initial conditions on both $\sigma_{ij}^{\ast}$ and $R_{ij}^{\ast}$ will also result into zero solution for them; however, owing to the coupling on the RHSs of (46) and (4), a non-vanishing initial condition on any one of them will trigger the non-vanishing solution for both of them. On the contrary, vanishing initial condition on $\Delta$ results into non-vanishing solution for it for all values of the coefficient of restitution except for $e=1$ and $e=1/\sqrt{2}$. In the following analysis, we take the initial conditions as $n_{\ast}(0)=T_{\ast}(0)=\sigma_{ij}^{\ast}(0)=q_{i}^{\ast}(0)=m_{ijk}^{\ast}(0)=R_{ij}^{\ast}(0)=1$, $v_{i}^{\ast}(0)=0$ and $\Delta(0)=15$.

Figure 1 illustrates the numerical solution for the (dimensionless) granular temperature and other higher-order moments obtained by solving the system (43)–(50) (without the underlined term in (50)) along with the aforementioned initial conditions. The solutions for the number density and velocity are not shown since they are just constants. Figure 1(a) exhibits the granular temperature for two values of the coefficient of restitution $e=0.75$ (in blue color) and $e=0.95$ (in red color). The solid lines depict the numerical solution for the granular temperature obtained by solving the system (43)–(50) (without the underlined term in (50)) along with the aforementioned initial conditions while the symbols denote the corresponding granular temperature computed via Haff’s law (54). Clearly, the numerical results (denoted by solid lines) are in good agreement with those obtained via Haff’s law (denoted by symbols). Moreover, it is also clear that the granular temperature relaxes faster with decreasing the coefficient of restitution. This is due to the fact that more inelastic particles dissipate more energy during the collision in comparison to the less inelastic ones resulting into the faster decay of the granular temperature for the former in comparison to the latter. The relaxation of other moments—$\sigma_{ij}^{\ast}$, $q_{i}^{\ast}$, $m_{ijk}^{\ast}$, $R_{ij}^{\ast}$ and $(\Delta/15)$—with time $t_{\ast}$ for $e=0.75$ is displayed in figure 1(b). It turns out that all these moments decay with time much faster than the granular temperature. Among themselves, $m_{ijk}^{\ast}$ decays with time $t_{\ast}$ faster than any other moment, followed by $R_{ij}^{\ast}$, $\sigma_{ij}^{\ast}$, $q_{i}^{\ast}$ and $\Delta$ (in the order of fast to slow). That the stress (dashed red line) decays faster than the heat flux (dotted blue line) which decays faster than $\Delta/15$ (magenta line with symbols) is in agreement with the findings of Kremer & Marques Jr. (2011).

4.3 Effect of nonlinear terms of dimensionless scalar moment ($\Delta$) on Haff’s law

In § 4.1, we obtained Haff’s law by dropping all the nonlinear terms of non-equilibrium moments ($\sigma_{ij}^{\ast}$, $q_{i}^{\ast}$, $m_{ijk}^{\ast}$, $R_{ij}^{\ast}$, $\Delta$) in (45) and (50). In this subsection, we shall investigate the effect of nonlinear terms of (dimensionless) scalar moment ($\Delta$) on Haff’s law in detail.

4.3.1 Case 1: Effect of $\Delta^{2}$ term present in (50)

The constant $\alpha$ in characteristic time scale $\tau_{0}$ for the temperature decay via Haff’s law (54) was obtained by discarding $\Delta^{2}$ term in (50). We now consider (50) without dropping any term. Equation (50) has two equilibrium points for $e\neq 1$. Nevertheless, a simple stability analysis of these equilibrium points (see Chapter 2 of Strogatz, 1994) shows that only one equilibrium point is stable whereas the other one is unstable. In the limiting case of $e\to 1$, the stable equilibrium point tends to zero while the unstable equilibrium point tends to infinity. This makes sense as (50) in the case of $e=1$ has only one equilibrium point (which is zero). Furthermore, in the limiting case of $e\to 1$, the stable equilibrium point leads to infinite relaxation time $\tau_{0}$ which is meaningful for monatomic gases (i.e., for $e=1$); however the unstable equilibrium point in this limit leads to a finite relaxation time which is meaningless. Owing to these reasons, we shall neglect the unstable equilibrium point.

Figure 2(a) portrays the stable (black filled circle) and unstable (white filled circle) equilibrium points for $e=0.95$ on a phase portrait. We consider only the stable equilibrium point (black filled circle) and neglect the unstable equilibrium point (white filled circle). The main panel of figure 2(b) illustrates the temperature decay via Haff’s law (54) with the inverse time scale (55) for three values of the coefficient of restitution $e=0.5,0.75,0.95$, shown by green, blue and red colors, respectively. The solid lines depict Haff’s law when the constant $\alpha$ in (55) is taken as the stable equilibrium point of (50) while the symbols show that when $\alpha$ is taken as the single equilibrium point (56)1 (i.e., in the linear case). In the inset of figure 2(b), we also plot the corresponding absolute difference between the temperature values obtained in both the cases. It turns out that the difference is many order of magnitude smaller than the original temperature values. Moreover, the difference is even smaller for large values of the coefficient of restitution. Thus we conclude that the $\Delta^{2}$ term present in (50) do not play any significant role on Haff’s law.

4.3.2 Case 2: Effect of $\Delta^{2}$ terms present on the RHSs of (45) and (3.1)

To investigate the effect of nonlinear terms of scalar moment (i.e., $\Delta^{2}$ terms) on Haff’s law one needs to include $\Delta^{2}$ terms in the production terms of the energy balance equation (13) and $w$ balance equation (3.1). On including $\Delta^{2}$ terms along with linear terms on the RHSs of (13) and (3.1), and introducing the scaling (42), the new energy and $\Delta$ balance equations in the HCS read

[TABLE]

[TABLE]

Notice again that although we included the nonlinear terms proportional to $\Delta^{2}$ on the RHSs of (13) and (3.1), the RHS of (4.3.2) contains $\Delta^{3}$ terms as well. Similar to above, for obtaining Haff’s law form (57), $\Delta$ must be a constant or, in other words, condition (53) must be fulfilled. With a constant value of $\Delta$ ($=\beta$, let us say), equations (43) and (57) along with the initial conditions $n_{\ast}(0)=T_{\ast}(0)=1$ again yield Haff’s law (54). However, the inverse of dimensionless characteristic time scale in this case reads

[TABLE]

Again, in the light of condition (53), the constant $\beta$ in (59) is an equilibrium point of differential equation (4.3.2).

Equation (4.3.2) has three equilibrium points for $e\neq 1$. Again, a simple stability analysis of these equilibrium points (see Chapter 2 of Strogatz, 1994) shows that only one equilibrium point is stable whereas the other two are unstable. Figure 3(a) displays the stable (black filled circle) and unstable equilibrium points (white filled circles) for $e=0.95$ on a phase portrait.

In the limiting case of $e\to 1$, the stable equilibrium point tends to zero, the left unstable point tends to $-480$ and the right unstable equilibrium point tends to infinity. That one equilibrium point of (4.3.2) in the limiting case of $e\to 1$ tends to infinity makes sense as (4.3.2) in the case of $e=1$ has only two equilibrium points. Furthermore, in the limiting case of $e\to 1$, the stable equilibrium point and the left unstable equilibrium point lead to infinite relaxation time $\tau_{0}$ which is meaningful for monatomic gases (i.e., for $e=1$); however the right unstable equilibrium point in this limit leads to a vanishing relaxation time which is meaningless, therefore we neglect the right unstable equilibrium point. Although, the left unstable equilibrium point lead to meaningful infinite relaxation time $\tau_{0}$ in the limit $e\to 1$, its value in the limit $e\to 1$ itself is not meaningful because the equilibrium points are the steady state solutions of (4.3.2) and the distribution function for monatomic gases (i.e., for $e=1$) in the steady state is Maxwellian, and consequently, by definition, $w$ (and hence $\Delta$) must vanish in the limit $e\to 1$ and in the steady state.

Therefore, we again consider only the stable equilibrium point (marked by black filled circle) and neglect the unstable equilibrium points (marked by white filled circles). The main panel of figure 3(b) illustrates the temperature decay via Haff’s law (54) for three values of the coefficient of restitution $e=0.5,0.75,0.95$, shown by green, blue and red colors, respectively. The solid lines depict the results obtained with the inverse time scale (59) and the constant $\beta$ in (59) as the stable equilibrium point of (4.3.2) (marked by black filled circle in the phase portrait) whereas the symbols delineate the corresponding results obtained with the inverse time scale (55) and the constant $\alpha$ as given in (56)1 (i.e., in the linear case). In the inset of figure 3(b), we also plot the corresponding absolute difference between the temperature values obtained in both the cases. There is practically no difference in both the results as the difference is again many order of magnitude smaller than the original temperature values.

Thus, we conclude that the nonlinear terms of (dimensionless) scalar moment $\Delta$ do not play any significant role on Haff’s law.

4.4 Effect of higher-order scalar moment on Haff’s law

In this subsection, we shall analyze effect of the sixth order scalar moment $u^{3}$, which is the full trace of the sixth order moment, on Haff’s law. Nevertheless, for convenience, we introduce a new dimensionless variable $\Xi$ which is related to the sixth order scalar moment $u^{3}$ via

[TABLE]

so that $\Xi$ vanishes if it is computed with $f$ as either Maxwellian or G13 or G26 distribution function.

We now include the governing equation for the moment $\Xi$ to our existing 26-moment system (eqs. (11)–(13) and (19)–(3.1)) and close this system with the Grad distribution function (which we shall name as the G27 distribution function)

[TABLE]

We shall refer to the system of these 27-moment equations closed with the G27 distribution function as the system of G27 equations.

For our purposes, it is not necessary to derive the full G27 equations but we do need the production terms for the energy, $\Delta$ and $\Xi$ balance equations. All the production terms associated with the G27 equations are given in appendix A. Note that the production terms for the $\Delta$ and $\Xi$ balance equations are given by

[TABLE]

respectively. In the HCS, the dimensionless energy, $\Delta$ and $\Xi$ balance equations from the system of G27 equations—on employing (42) and and discarding all the nonlinear terms of non-equilibrium scalar, vector or tensor moments—read

[TABLE]

[TABLE]

[TABLE]

Clearly, for obtaining Haff’s law form (62), $\Delta$ and $\Xi$ must be constants or, in other words, conditions

[TABLE]

must be fulfilled. With constant values of $\Delta$ ($=\varkappa_{1}$, let us say) and $\Xi$ ($=\varkappa_{2}$, let us say), eqs. (43) and (62) along with the initial conditions $n_{\ast}(0)=T_{\ast}(0)=1$ again yield Haff’s law (54). However, the inverse of dimensionless characteristic time scale now reads

[TABLE]

In the light of conditions (65), the constants ($\varkappa_{1},\varkappa_{2}$) in (66) are the equilibrium points ($\Delta_{\mathrm{eq}},\Xi_{\mathrm{eq}}$) of the system of two first-order ordinary differential equations (4.4) and (4.4) (Strogatz, 1994). The system of equations (4.4) and (4.4) has a unique equilibrium point. Figure 4 illustrates the temperature decay via Haff’s law (54) for three values of the coefficient of restitution $e=0.5,0.75,0.95$ shown by green, blue and red colors, respectively. The solid lines depict the results obtained with the inverse time scale (66) and the constants $(\varkappa_{1},\varkappa_{2})$ in (66) as the unique equilibrium point of (4.4) and (4.4) whereas the symbols delineate the corresponding results obtained with the inverse time scale (55) and the constant $\alpha$ as given in (56)1. In the inset of figure 4(b), we also plot the corresponding absolute difference between the temperature values obtained in both the cases. Again, there is practically no difference in the results as the differences are many order of magnitude smaller than the original temperature values.

Furthermore, in view of § 4.3, it is expected that the nonlinear terms of scalar moments in (62)–(4.4) would also have only negligible effect on Haff’s law. Thus, it is concluded empirically that the nonlinear terms of non-equilibrium scalar moments as well as higher-order scalar moments do not have any significant effect on Haff’s law. Therefore, it suffices to consider only those moment systems which consists of the fourth order scalar moment ($\Delta$)—i.e., to consider 14- or 26-moment system. Moreover, in the production terms of these moment systems, it is sufficient to retain only those nonlinear terms which are products of $\Delta$ and a vector or a tensor.

Owing to this reason, equations (30)–(4) without the underlined term in (4) will be referred to as the G26 equations henceforth.

5 Constitutive relations for the stress and heat flux: Navier–Stokes, and Fourier laws

The mass, momentum and energy balance equations (11)–(13) are not closed since they contain the stress ($\sigma_{ij}$), heat flux ($q_{i}$) and cooling rate ($\zeta$) as additional unknowns. One of the main goal in the kinetic theory for granular gases is to derive the constitutive relations for these unknowns in order to close the system of mass, momentum and energy balance equations. To the linear approximation in the spatial gradients, these constitutive relations read (see e.g., Jenkins & Richman, 1985a, b; Garzó & Dufty, 1999; Garzó, 2013)

[TABLE]

Equations (67)1 and (67)2 are the Navier–Stokes law and Fourier law, respectively. The coefficients $\zeta_{0}$ and $\zeta_{1}$ in (67)3 are the zeroth- and first-order contributions, respectively, to the cooling rate. The transport coefficients $\eta$, $\kappa$ are referred to as the shear viscosity and thermal conductivity, respectively; and the transport coefficient $\lambda$ is a special coefficient for granular gases which vanishes identically for elastic gases. Typically, these transport coefficients are obtained by a formal CE expansion on the Boltzmann equation (see e.g., Brey et al., 1998a; Sela & Goldhirsch, 1998; Gupta, 2011), nevertheless some authors have also computed them via computer-aided kinetic theory (Noskowicz et al., 2007) and via Grad’s method of moment, (see e.g., Jenkins & Richman, 1985a; Kremer & Marques Jr., 2011; Garzó, 2013). Although Grad’s method of moments does not restrict the coefficient of restitution to be close to $1$, G13 theory of Jenkins & Richman (1985a) and the Maxwellian iteration procedure on the G14 equations of Kremer & Marques Jr. (2011) neglect the effect of cooling rate on the transport coefficients. On the other hand, the work of Garzó (2013) incorporates the effect of cooling rate on the transport coefficients and, consequently, the transport coefficients obtained in this work are in good agreement with those obtained via the classical CE expansion. In this section, we shall determine these transport coefficients through the G26 equations (eqs. (30)–(4) without the underlined term in (4)). To this end, we shall follow the approach of Garzó (2013) so that the effect of cooling rate on the transport coefficients is incorporated.

We are interested in the linear approximation in the spatial gradients. To this end, let us first analyze the zeroth-order contributions in the spatial gradients.

5.1 Zeroth-order contributions in the spatial gradients

To zeroth-order in the spatial gradients, $\Delta$ balance equation (4) yields

[TABLE]

Consequently, to zeroth-order in the spatial gradients, the mass, momentum and energy balance equations (30)–(32) reduce to

[TABLE]

with

[TABLE]

while the balance equations for the higher moments (eqs. (4)–(4)) simplify to

[TABLE]

respectively, with

[TABLE]

Equations (71)1,4 imply that

[TABLE]

Thus, to zeroth-order in the spatial gradients, $\sigma_{ij}$, $q_{i}$, $m_{ijk}$ and $R_{ij}$ are zero while $\Delta=15\,a_{2}$. This means that all the non-equilibrium vectorial and tensorial moments are at least of linear order in the spatial gradients.

5.2 First-order contributions in the spatial gradients

We now investigate the moments up to first-order in the spatial gradients, which will lead to the desired Navier–Stokes and Fourier laws.

To first-order in the spatial gradients, the mass, momentum and energy balance equations (30)–(32) yield

[TABLE]

where $\zeta$ needs to be considered up to linear order in the spatial gradients. Furthermore, the stress and heat flux balance equations (eqs. (4) and (4)), to first-order in the spatial gradients, reduce to

[TABLE]

with

[TABLE]

Notice that we do not need the balance equations for other moments except (4), the one for $R_{ij}$, for first-order approximation. Equation (4) is only needed because the contribution of $R_{ij}$ is required in (77). To first-order in the spatial gradients, the $R_{ij}$ balance equation (4) reduces to

[TABLE]

Here, we ignore the time derivative of $R_{ij}$ by assuming that it decreases with time faster than stress and heat flux. With this assumption, equation (80) yields

[TABLE]

With (81), equation (77) turns to

[TABLE]

The time derivatives of the stress and heat flux in (82) and (78) are now calculated as follows. From (67)1,2, we have

[TABLE]

For first-order approximation in the spatial gradients, the time derivatives on the RHSs of (83) are required only up to zeroth-order in spatial gradients. Hence, by virtue of (69), the underlined terms in (83) vanish while the underbraced term in (83)2 becomes

[TABLE]

By dimensional analysis, it turns out that $\eta\propto\sqrt{T}$, $\kappa\propto\sqrt{T}$ and $\lambda\propto T^{3/2}$. Therefore, to zeroth-order in spatial gradients—on using (69)3—the time derivatives of the transport coefficients read

[TABLE]

The time derivatives of the stress and heat flux (83) on using (84) and (85) become

[TABLE]

Inserting the time derivatives of the stress and heat flux from (86) into (82) and (78), and comparing the coefficients of the spatial gradient of each hydrodynamic variable, one obtains the transport coefficients in the Navier–Stokes and Fourier laws (67)1,2 as

[TABLE]

where

[TABLE]

are the elastic values—in first Sonine approximation—of the shear viscosity and thermal conductivity, respectively. They are also referred to as the Enskog viscosity and Enskog thermal conductivity, respectively (Brilliantov & Pöschel, 2004). Moreover, $\eta^{\ast}$, $\kappa^{\ast}$ and $\lambda^{\ast}$ in (87) are the reduced shear viscosity, reduced thermal conductivity and reduced coefficient corresponding to $\lambda$. These reduced transport coefficients are given by

[TABLE]

Note that if we take $\nu_{\sigma R}^{\ast}=0$, i.e., if we ignore the contribution of $R_{ij}$, the reduced transport coefficients (89) obtained here agree with those obtained via Grad’s method of moments in Garzó (2013) for dilute granular gases and also with those obtained by first-Sonine approximation via the CE expansion in Brey et al. (1998a). If we further neglect the cooling effect (i.e., $\zeta_{0}^{\ast}=0$) as well, the reduced transport coefficients obtained here concur with those derived in Kremer & Marques Jr. (2011).

In the elastic case ($e=1$), the reported (exact) value of the reduced shear viscosity $\eta^{\ast}$ is $1.016034$ which was obtained by Pekeris & Alterman (1957) by reducing the Boltzmann equation to an ordinary differential equation of order $4$ and subsequently integrating it numerically. Reinecke & Kremer (1990) obtained almost the same value $\eta^{\ast}\approx 1.016028$ for the reduced shear viscosity by employing Maxwell iteration procedure on Grad moment equations—at fifth order of approximation. The reduced shear viscosity for $e=1$ obtained from (89)1 is $\eta^{\ast}=205/202\,(\approx 1.014851)$, which is exactly same as that obtained at second order approximation in Reinecke & Kremer (1990), thanks to the coupling between $\sigma_{ij}$ and $R_{ij}$ in their balance equations. Note that the first Sonine approximation (Brey et al., 1998a), modified Sonine approximation (Garzó et al., 2007), the G14 theory of Kremer & Marques Jr. (2011) and Grad moment theory of Garzó (2013), they all in the elastic case lead to $\eta^{\ast}=1$, which is the value obtained at first order approximation in Reinecke & Kremer (1990). Furthermore, the reported value of the reduced thermal conductivity $\kappa^{\ast}$ in the elastic case is $1.025218$ (Pekeris & Alterman, 1957). Reinecke & Kremer (1990) again obtained almost the same value $\kappa^{\ast}\approx 1.025197$ at fifth order approximation. However, the reduced thermal conductivity obtained here as well as that obtained by all the four aforementioned approaches is $\kappa^{\ast}=1$, which is the value obtained at first order approximation in Reinecke & Kremer (1990). The reason for not obtaining an improved value of the reduced thermal conductivity here is that, unlike the RHS of the stress balance equation which is coupled with the one trace of the fourth moment ($R_{ij}$), the RHS of the heat flux balance equation does not have coupling with other vectorial or tensorial moments in its considered form.

5.3 Comparison with existing theories and computer simulations

In this subsection, we shall compare the reduced transport coefficients derived here with those obtained by CE expansion as well as with those obtained through computer simulations.

Figure 5 illustrates the reduced shear viscosity $\eta^{\ast}$ while figure 6 illustrates the reduced thermal conductivity $\kappa^{\ast}$ and reduced coefficient $\lambda^{\ast}$ for different values of the coefficient of restitution. The solid (black) lines in both the figures represent the results obtained with expressions (89) in the present work. The dashed (green) lines are the plots for the reduced transport coefficients derived at first Sonine approximation through CE expansion on the Boltzmann equation in Brey et al. (1998a). As the reduced transport coefficients obtained via Grad’s method of moments in Garzó (2013) coincide with those in Brey et al. (1998a), the dashed (green) lines also display the reduced transport coefficients obtained in Garzó (2013). The dash-dotted (magenta) lines display the results obtained with the theoretical expressions for the reduced transport coefficients deduced through the modified version of the first Sonine approximation in Garzó et al. (2007). The squares are the results obtained with the theoretical expressions derived via the so-called computer-aided method devised by Noskowicz et al. (2007) while the triangles denote the numerical solution of the Boltzmann equation obtained through Green–Kubo (GK) relations by means of the direct-simulation Monte Carlo (DSMC) method (Bird, 1994) in Brey et al. (2005). The circles in figure 5 also denote the DSMC simulation results from Montanero et al. (2005) obtained with another method—by the implementation of an external force which compensates for the collisional cooling. The DSMC simulation data from this method for $e=\{0.6,0.7,0.8,0.9,1\}$ were obtained by Montanero et al. (2005) while those for $e=\{0.2,0.3,0.4,0.5\}$ were obtained by Garzó et al. (2007).

In general, the reduced shear viscosity from the first Sonine approximation (dashed green line) agrees with that from the DSMC simulations using GK relations (triangles). Nevertheless, the results from the DSMC simulations with the implementation of external force (shown by circles) differ from those obtained with the GK relations (shown by triangles), especially for $e\lesssim 0.7$. The difference between the results from the two simulations could be due to velocity correlations in the correlation function present in the GK relation for the shear viscosity (Garzó et al., 2007). In comparison to the results from the modified version of the first Sonine approximation (dash-dotted magenta line), the present results for the reduced shear viscosity (solid black line) are in better agreement with the DSMC simulations of Garzó et al. (2007) even for small coefficient of restitution. Furthermore, the present results for the reduced shear viscosity (solid black line) coincide with those obtained from the computer-aided method devised by Noskowicz et al. (2007) (shown by squares) for almost all values of the coefficient of restitution. Nevertheless, in the elastic ($e=1$) case, the theoretical expressions of Noskowicz et al. (2007) lead to $\eta^{\ast}=1.01205$ while the DSMC simulations of Garzó et al. (2007) give $\eta^{\ast}\approx 1.016$. Thus, in the elastic ($e=1$) case, $\eta^{\ast}\approx 1.014851$ from the expression (89)1 obtained here is much closer to its true value in comparison to that from the theoretical expressions of Noskowicz et al. (2007).

The expressions for the reduced transport coefficients corresponding to the heat flux, (i.e., for $\kappa^{\ast}$ and $\lambda^{\ast}$) derived here are the same as those obtained via the first Sonine approximation (Brey et al., 1998a) and via Grad moment method of Garzó (2013) because in the considered linearized form, the RHS of heat flux balance equation (4) does not couple with any other vectorial or tensorial moments. Therefore, the curves of the reduced thermal conductivity $\kappa^{\ast}$ and the reduced coefficient $\lambda^{\ast}$ in figure 6 from (89)2,3 (solid black lines) coincide with those from the first Sonine approximation (Brey et al., 1998a) or from Grad moment method of Garzó (2013) (dashed green lines). It is clear from figure 6 that the results for $\kappa^{\ast}$ and $\lambda^{\ast}$ from the present work, from first Sonine approximation and from Garzó (2013) agree with the DSMC simulations as well as with those from the computer-aided method of Noskowicz et al. (2007) only for large coefficient of restitution ($e\gtrsim 0.7$). For $e\lesssim 0.7$, similar to the first Sonine approximation and Grad moment method of Garzó (2013), the present work also overestimates both the transport coefficients.

Motivated from the satisfactory result for the reduced shear viscosity through the G26 equations, it is expected that a suitable coupling on the RHS of the heat flux balance equation (4) would improve the reduced transport coefficients $\kappa^{\ast}$ and $\lambda^{\ast}$ significantly. The first such coupling of RHS of the heat flux balance equation in the semi-linearized setting as considered here will be introduced on including the full trace of fifth order moment into the moment system, i.e., one would need to consider the system of 29 moment equations. This is beyond the scope of the present paper and will be considered elsewhere in future.

Although, the modified Sonine approach (Garzó et al., 2007) is an excellent theoretical approach whose predictions for all the transport coefficients ($\eta^{\ast}$, $\kappa^{\ast}$ and $\lambda^{\ast}$) are in very good agreement with the simulation results, the G26 theory slightly improves the results for the reduced shear viscosity coefficient $\eta^{\ast}$. However, the latter could not yield the other two transport coefficients related to the heat flux ($\kappa^{\ast}$ and $\lambda^{\ast}$) correctly for $e\lesssim 0.7$, exactly in the same way as the first Sonine approach (Brey et al., 1998a). Here, it is worthwhile to note that while it is an important validity test for moment equations to compute the precise form of the transport coefficients, it is well-known that these coefficients have only a little relevance in non-equilibrium situations when a linear relation between fluxes and gradients does not hold anymore. Instead, it has been shown (for elastic gases) that moment equations allow for cross effects like thermal stresses and non-gradient heat fluxes (Torrilhon, 2016) that are expected to also have a systematic influence in granular flows. The actual goal of moment equations is not only to match the transport coefficients but also to go beyond these classical theories and provide an enhanced fluid dynamic theory for granular gases. Moreover, our long-term perspective is to establish a complete set of predictive moment equations—somewhat similar to the R13 equations for monatomic gases—for granular flows, which certainly requires the G26 equations or beyond. Thus this work is the first step towards this long-term goal.

6 Linear stability analysis

In this section, we investigate the stability of the HCS due to small perturbations through various moment theories—particularly, with the G26 theory (equations (30)–(4) without the underlined term in (4))—developed in the earlier sections. We assume that the amplitudes of these perturbations are sufficiently small so that the linear analysis remains valid.

For the linear stability analysis, we decompose all the field variables into their reference values—i.e., their respective solutions in the HCS—and into perturbations from their respective solutions in the HCS, i.e., we define

[TABLE]

where $n_{0}$ is the constant number density and $T_{H}(t)$ is the granular temperature in the HCS; $\Delta_{H}$ is the constant solution for $\Delta$ in the HCS; $v_{H}(t)=\sqrt{T_{H}(t)/m}$ is a reference speed in the HCS and is proportional to the adiabatic sound speed in the HCS; the reference values (HCS solutions) for all other field variables are denoted by the superscript ‘$(H)$’; the quantities with tilde denote the dimensionless perturbations in the field variables from their respective solutions in the HCS. Note that as the underlined term in (4) is discarded, the HCS solution for $\Delta$ is $\Delta_{H}=15\,a_{2}$.

Inserting the field variables from (90) into the G26 equations, and neglecting all the nonlinear terms of the perturbed field variables (denoted with tilde in (90)), one obtains the system of linear partial differential equations in (dimensionless) perturbed field variables with time-dependent coefficients, which is given in appendix B. These equations are further simplified by exploiting the fact—concluded in § 4.2—that in the HCS, the non-equilibrium vectorial and tensorial moments—$\sigma_{ij}^{(H)}(t)$, $q_{i}^{(H)}(t)$, $m_{ijk}^{(H)}(t)$, $R_{ij}^{(H)}(t)$—decay faster than the granular temperature $T_{H}(t)$; therefore, we drop the terms containing $\sigma_{ij}^{(H)}(t)$, $q_{i}^{(H)}(t)$, $m_{ijk}^{(H)}(t)$, $R_{ij}^{(H)}(t)$ in (127)–(B). Now, it is possible to convert this system of partial differential equations to a new system of partial differential equations having constant coefficients as follows. We introduce a length scale

[TABLE]

to make the space variables dimensionless (i.e., $\tilde{x}_{i}=x_{i}/\ell$, where tilde again denotes the dimensionless space variable), and a dimensionless time $\tilde{t}$ (see McNamara, 1993) such that

[TABLE]

This leads to

[TABLE]

where

[TABLE]

With these definitions of dimensionless space and time, system (127)–(B) now simplifies to

[TABLE]

where the coefficients

[TABLE]

with

[TABLE]

depend only on the parameter $e$, the coefficient of restitution.

Now, we assume a normal mode solution of the form

[TABLE]

for system (94)–(101). Here, quantities with hats denote the complex amplitudes of the perturbed field variables; $\mathbbm{i}$ is the imaginary unit; $\bm{k}$ and $\omega$ are the (dimensionless) wavevector and (dimensionless) frequency, respectively, of the disturbance. For the temporal stability analysis—to be analyzed here—the wavevector $\bm{k}$ is assumed to be real and the frequency $\omega$ is assumed to be complex. The real part of the complex frequency, $\mathrm{Re}(\omega)$, determines the phase velocity $\bm{v}_{\mathrm{ph}}=\mathrm{Re}(\omega)/\bm{k}$ of the corresponding wave whereas the imaginary part of the complex frequency, $\mathrm{Im}(\omega)$, determines whether the amplitude of the disturbance grows or decays in time. The imaginary part of the complex frequency, $\mathrm{Im}(\omega)$, is referred to as the growth rate. From the normal mode solution (104), it is clear that the solution will decay (or grow) in time if the growth rate is negative (or positive). Consequently, stability of the system requires the growth rate to be non-positive, i.e., $\mathrm{Im}(\omega)\leq 0$.

If we assume that the wavevector of the disturbance is parallel to the $x$-axis, i.e., $\bm{k}=k\,\hat{\bm{x}}$ where the wavenumber $k$ is the magnitude of the wavevector $\bm{k}$ and $\hat{\bm{x}}$ is the unit vector in $x$-direction, we get two independent eigenvalue problems—namely, the longitudinal problem and the transverse problem—for the amplitude of the disturbance in two dimensions. These problems read

[TABLE]

respectively, where

[TABLE]

and

[TABLE]

For nontrivial solution of each eigenvalue problem in (105), the determinant of each of the two matrices $\mathcal{A}$ and $\mathcal{B}$ must vanish, i.e., $\mathrm{det}(\mathcal{A})=0$ and $\mathrm{det}(\mathcal{B})=0$. These two conditions on the determinant of matrices $\mathcal{A}$ and $\mathcal{B}$ lead to the following dispersion relations

[TABLE]

for the longitudinal and transverse problems (105), respectively. The coefficients $\mathbbm{a}_{1},\mathbbm{a}_{2},\dots,\mathbbm{a}_{8}$ and $\mathbbm{b}_{1},\mathbbm{b}_{2},\dots,\mathbbm{b}_{5}$ in (108) and (109) are functions of the wavenumber $k$ and the coefficient of restitution $e$; although the explicit values of these coefficients are not given here for brevity.

6.1 Eigenmodes

Condition (108) results into eight eigenmodes for the longitudinal system (105)1. These eigenmodes for the coefficient of restitution $e=0.75$ and $e=1$ are shown in figures 7 and 8, respectively. The left panels in each of figures 7 and 8 depict the real part of frequency, $\mathrm{Re}(\omega)$, while the right panels display the growth rate, $\mathrm{Im}(\omega)$, all plotted over the wavenumber $k$. For very small wavenumbers ($k\lesssim 0.02$), all the eigenmodes of the longitudinal system (105)1 are stationary in the case of $e=0.75$ (see the left panel of figure 7); however, as the wavenumber starts increasing, a pair of traveling modes commences at $k\approx 0.02$, another pair of slow traveling modes commences at $k\approx 0.043$ and a third pair of even slower traveling modes starts at $k\approx 0.057$. These six eigenmodes (which are traveling for large $k$) are referred to as the sound modes (Brilliantov & Pöschel, 2004; Garzó, 2005). Each pair of sound modes propagates in the opposite directions (for large $k$) since the corresponding eigenvalues are in complex conjugate pairs. The remaining two eigenmodes continue to be stationary for all wavenumbers since the frequencies associated with them are purely imaginary. These eigenmodes (which remain stationary for all wavenumbers) are referred to as the heat modes (Brilliantov & Pöschel, 2004; Garzó, 2005). The left panel of figure 8 illustrates that there are two (heat) eigenmodes which remain stationary for all wavenumbers and three pairs of traveling (sound) eigenmodes in the elastic case (i.e., for $e=1$) as well. Nevertheless, two pairs of sound modes commence propagating in the opposite directions at $k=0$ itself in this case. These two pairs of sound modes travel almost with the same speeds for $0\leq k\lesssim 0.18$ whereas for $k\gtrsim 0.18$ one pair of sound modes travels faster than the other. It may also be interesting to note that out of these two pairs of sound modes, the speed of the one pair traveling with slower speed coincides with that of a third pair of sound modes—which starts propagating at $k\approx 0.095$ in the opposite directions with even slower speed—for $0.28\lesssim k\lesssim 0.45$. It can also be noticed by comparing the left panels of figures 7 and 8 that the speeds of sound modes in the case of $e=0.75$ (figure 7) are almost same as those of corresponding sound modes in the elastic case (figure 8) for wavenumbers $k\gtrsim 0.7$.

Since the complex eigenvalues of a matrix occur in pairs, the growth rates, $\mathrm{Im}(\omega)$, of a pair of sound modes coincide beyond the wavenumber at which they start propagating. This is clearly reflected in right panels of figures 7 and 8: the imaginary parts of frequencies, $\mathrm{Im}(\omega)$, of the three pairs of traveling waves—in the left panels $(i)$ starting at $k\approx 0.02$, $k\approx 0.043$ and $k\approx 0.057$ in the case of $e=0.75$ (figure 7) and $(ii)$ two pairs starting at $k=0$ and one starting at $k\approx 0.095$ in the elastic case (figure 8)—merge together beyond these wavenumber values (compare the corresponding left and right panels of figures 7 and 8). As we have discussed above that non-positive growth rate ($\mathrm{Im}(\omega)\leq 0$) implies stability and vice versa, it is clear from the right panel of figure 7 that for $e=0.75$, one heat mode and all the sound modes are stable for all wavenumbers whereas the other heat mode (which has $\mathrm{Im}(\omega)>0$ for some wavenumbers) is unstable for small wavenumbers. On the other hand, it is evident from the right panel of figure 8 that in the elastic case (i.e., for $e=1$), all the heat and sound modes of the longitudinal system (105)1 are stable for all wavenumbers. The unstable heat mode in the case of $e=0.75$ turns stable beyond $k\approx 0.2005$ (see the right panel of figure 7) since beyond this value of the wavenumber, $\mathrm{Im}(\omega)$ becomes negative. A value of wavenumber $k$ at which the growth rate, $\mathrm{Im}(\omega)$, switches its sign is referred to as the critical wavenumber. In other words, the corresponding eigenmode is unstable (or stable) for the wavenumbers below the critical wavenumber since $\mathrm{Im}(\omega)>0$ (or $\mathrm{Im}(\omega)\leq 0$) for them while it is stable (or unstable) for those above the critical wavenumber since $\mathrm{Im}(\omega)\leq 0$ (or $\mathrm{Im}(\omega)>0$) for them. Thus, $k\approx 0.2005$ is the critical wavenumber for the longitudinal system (105)1 in the case of $e=0.75$.

Condition (109) leads to five eigenmodes for the transverse system (105)2. The eigenmodes of the transverse system (105)2 are referred to as the shear modes (Brilliantov & Pöschel, 2004; Garzó, 2005). These shear modes for the coefficient of restitution $e=0.75$ and $e=1$ are illustrated in figures 9 and 10, respectively. The left panels in each of figures 9 and 10 again delineate the real part of frequency, $\mathrm{Re}(\omega)$, while the right panels portray the growth rate, $\mathrm{Im}(\omega)$, all plotted over the wavenumber $k$. For very small wavenumbers ($k\lesssim 0.065$), all the shear modes are stationary in the case of $e=0.75$ (see the left panel of figure 9). A pair of traveling shear modes commences at $k\approx 0.065$ and another pair of slow traveling shear modes commences at $k\approx 0.562$. Each pair of traveling shear modes propagates in the opposite directions, and one remaining shear mode continues to be stationary for all wavenumbers since the frequency associated with it is purely imaginary. For the same reason as discussed above, the imaginary parts of frequencies of each pair of traveling shear modes coincide beyond the wavenumbers at which they start propagating (one pair coincides for $k\gtrsim 0.066$ and another for $k\gtrsim 0.563$ in the right panel of figure 9). Furthermore, it is noted from the right panel of figure 9 that one shear mode in the case of $e=0.75$ is also unstable for wavenumber values below the critical wavenumber which is $k\approx 0.2827$ while it is stable for wavenumber values above the critical wavenumber; the remaining four shear modes are always stable for all values of wavenumber $k$.

In the elastic case, for $e=1$ (figure 10), all the shear modes are stationary for small wavenumbers ($k\lesssim 0.23$). A pair of traveling shear modes commences at $k\approx 0.23$ but turns back to become stationary on slight increase in the wavenumber (at $k\approx 0.241$), refer to the zoomed region shown in the insets of figure 10. On further increase in the wavenumber, a pair of traveling shear modes starts at $k\approx 0.25$ and another pair of slow traveling shear modes commences at $k\approx 0.574$. Each pair of traveling shear modes propagates in the opposite directions, and one shear mode remains stationary for all wavenumbers since the frequency associated with it is purely imaginary. For the same reason as discussed above, the imaginary parts of frequencies of each pair of traveling shear modes coincide for the wavenumbers for which they are propagating—one pair coincides for $0.231\lesssim k\lesssim 0.24$ (see the inset in the right panel of figure 10), one pair for $k\gtrsim 0.251$ and another for $k\gtrsim 0.575$ in the right panel of figure 10). Furthermore, it is noted from the right panel of figure 10 that all five shear modes are always stable in the elastic case. It can again be noticed from the left panels of figures 9 and 10 that the speeds of traveling shear modes in the case of $e=0.75$ (figure 9) are almost same as those of corresponding traveling shear modes in elastic case (figure 10) for wavenumbers $k\gtrsim 0.6$.

6.2 Eigenmodes in small wavenumber limit

There is another interesting classification of eigenmodes in the small wavenumber limit (or long wavelength limit), i.e., in the limit $k\to 0$: an eigenmode is referred to as a hydrodynamic mode if the frequency $\omega(k)$ of this eigenmode vanish in the limit $k\to 0$ while it is referred to as a kinetic mode if its frequency $\omega(k)$ attains a nonzero constant value in the limit $k\to 0$ (Kremer & Marques Jr., 2011). In order to explore the behavior of the eigenmodes in the small wavenumber limit (i.e., in the limit $k\to 0$), the frequency $\omega$ is expressed in powers of $k$ (see Kremer & Marques Jr., 2011):

[TABLE]

The unknown coefficients $\gamma_{0},\gamma_{1},\gamma_{2},\dots$ in the above expansion are determined by inserting expansion (110) for $\omega$ into (108) and (109), and solving the algebraic equations resulting from the comparison of coefficients of each power of $k$ on both sides of each equation. With this technique, it turns out that in the limit $k\to 0$, the frequencies of eight eigenmodes of the longitudinal system (105)1 in the inelastic (i.e., $e\neq 1$) case are related to the wavenumber $k$ via

[TABLE]

and those in the elastic (i.e., $e=1$) case are related to the wavenumber $k$ via

[TABLE]

The unknown constants $\xi_{-}$, $\xi_{+}$, $\vartheta$ and $\vartheta_{0}$ in (111) are given by

[TABLE]

while the constants $\vartheta_{1}$ and $\vartheta_{2}$ in (111) are too cumbersome to write here.

It is clear from (111) that, in the limit $k\to 0$, out of eight eigenmodes of the longitudinal system (105)1, one mode is hydrodynamic while the rest seven modes are kinetic in the inelastic case (i.e., $e\neq 1$); however, from (112), it is evident that three modes out of eight eigenmodes of the longitudinal system (105)1 are hydrodynamic in the elastic case while only the remaining five modes are kinetic. Figure 11 illustrates the real (in the left panel) and imaginary (in the right panel) parts of frequencies $\omega$—associated with the eigenmodes of the longitudinal system (105)1—plotted over the coefficient of restitution $e$ for $k=0$. The solid lines in figure 11 denote the frequencies obtained directly from the condition $\mathrm{det}(\mathcal{A})=0$ on substituting $k=0$ in the matrix $\mathcal{A}$ while the circles denote the frequencies obtained from the analytic expressions (111) and (112) in the limit $k\to 0$; moreover, the numbers in figure 11 represent the numbering of eigenmodes as given in (111) and (112). Clearly, the frequencies associated with the eigenmodes of the longitudinal system (105)1 in the limit $k\to 0$ from both the methods match perfectly well. From the right panel of figure 11, one may notice that the magnitude of the imaginary part of frequency, $\mathrm{Im}(\omega)$, for eigenmodes $2$ and $3$ decrease as the coefficient of restitution increases and they vanish for $e=1$; thus, it is concluded that the longitudinal system (105)1 has three hydrodynamic modes exclusively for the elastic case ($e=1$); in the inelastic case ($e\neq 1$), the longitudinal system (105)1 always has only a single hydrodynamic mode. The left panel of figure 11 shows that for the wavenumber $k=0$ only two eigenmodes ($7$ and $8$) have nonzero real part of frequency for the coefficient of restitution $0\leq e\lesssim 0.73$ while the rest six eigenmodes ($1,2,\dots,6$) are purely diffusive since the frequencies associated with them are purely imaginary. This can also be seen from (111): apart from frequency $\omega^{(1)}$ which is zero in the limit $k\to 0$ for all $e$, the frequencies $\omega^{(2)},\omega^{(3)},\dots,\omega^{(6)}$ are purely imaginary in the limit $k\to 0$; the frequencies $\omega^{(7)}$ and $\omega^{(8)}$ are, in general, complex with nonzero real part—even in the limit $k\to 0$—since $\xi_{-}$ and $\xi_{+}$ may be complex depending on the value of the coefficient of restitution; nonetheless, in the limit $k\to 0$, $\omega^{(7)}$ and $\omega^{(8)}$ also become purely imaginary for $0.73\lesssim e\leq 1$—it can be seen form (112) for $e=1$ as well as in the left panel of figure 11. One can further perceive from the right panel of figure 11 that the imaginary part of frequency of one eigenmode (denoted by number ‘$3$’) is always positive for $e\neq 1$; consequently, one eigenmode of longitudinal system (105)1 in the inelastic case ($e\neq 1$) is unstable in the limit $k\to 0$ whereas all the other eigenmodes are stable in this limit. Nevertheless, all the eigenmodes of longitudinal system (105)1 in the elastic case ($e=1$) are always stable for any wavenumber (see also the right panel of figure 8).

By the similar frequency expansion technique employed above, the frequencies of five eigenmodes of the transverse system (105)2 in the limit $k\to 0$ are related to the wavenumber $k$ in the inelastic ($e\neq 1$) case via

[TABLE]

and, in the elastic ($e=1$) case, these frequencies simplify to

[TABLE]

The constants $\vartheta_{3}$ and $\vartheta_{4}$ in (113) are also too cumbersome to write here.

Expressions (113) and (114) suggest that, in the limiting case of limit $k\to 0$, all five eigenmodes of the transverse system (105)2 are kinetic for all values of the coefficient of restitution except in the elastic case (i.e., for $e=1$) for which $\xi_{0}=0$ and consequently, one of the five eigenmodes is hydrodynamic in the elastic limit. Figure 12 illustrates the real (in the left panel) and imaginary (in the right panel) parts of frequencies $\omega$—associated with the eigenmodes of the transverse system (105)2—plotted over the coefficient of restitution $e$ for $k=0$. The solid lines in figure 12 denote the frequencies obtained directly from the condition $\mathrm{det}(\mathcal{B})=0$ on substituting $k=0$ in the matrix $\mathcal{B}$ while the circles denote the frequencies obtained from the analytic expressions (113) and (114) in the limit $k\to 0$; moreover, the numbers in figure 12 represent the numbering of eigenmodes as given in (113). Clearly, the frequencies associated with the eigenmodes of the transverse system (105)2 in the limit $k\to 0$ from both the methods match perfectly well. From the right panel of figure 12, one may notice that the imaginary part of frequency, $\mathrm{Im}(\omega)$, for one eigenmode (denoted by ‘$1$’) is non-negative and decreases as the coefficient of restitution increases and, finally, vanishes for $e=1$; therefore, it is concluded that one eigenmode (denoted by ‘$1$’) of the transverse system (105)2 in the elastic case (i.e., for $e=1$) is hydrodynamic since its real and imaginary parts are zero (cf. (114)1) and that this mode of the transverse system (105)2 in the inelastic case ($e\neq 1$) is unstable in the limit $k\to 0$ whereas all the other eigenmodes are stable in this limit. Nonetheless, all the eigenmodes of transverse system (105)2 in the elastic case ($e=1$) are always stable for any wavenumber (see also the right panel of figure 10). The left panel of figure 12 shows that for the wavenumber $k=0$, only two eigenmodes ($4$ and $5$) have nonzero real part of frequency for the coefficient of restitution $0\leq e\lesssim 0.73$ while the remaining three eigenmodes ($1,2$ and $3$) are purely diffusive since the frequencies associated with them are purely imaginary. This can also be seen from (113): the frequencies $\omega^{(1)},\omega^{(2)}$ and $\omega^{(3)}$ are purely imaginary in the limit $k\to 0$; the frequencies $\omega^{(4)}$ and $\omega^{(5)}$ are, in general, complex with nonzero real part—even in the limit $k\to 0$—since $\xi_{-}$ and $\xi_{+}$ may be complex depending on the value of the coefficient of restitution; nonetheless, in the limit $k\to 0$, $\omega^{(4)}$ and $\omega^{(5)}$ also become purely imaginary for $0.73\lesssim e\leq 1$—it can be seen form (114) for $e=1$ as well as in the left panel of figure 12.

6.3 Comparison among various Grad moment theories

As discussed in § 3.3, one can obtain a lower-level system of Grad moment equations by dropping the appropriate field variables in the system of G26 equations. In the same way, one can obtain the longitudinal and transverse systems associated with the G13, G14, G20 and G21 equations by dropping the appropriate variables and corresponding rows and columns of the matrices $\mathcal{A}$ and $\mathcal{B}$ in (105). Note that the transverse systems for the G13 and G14 equations are same since the field variable $\hat{\Delta}$ does not appear in the list of unknowns of the transverse system (see (105)2); for the same reason, the transverse systems for the G20 and G21 equations are also same. It is also important to note that the results for the G13 theory presented below are equivalent to those of Kremer & Marques Jr. (2011) since they assumed a constant value for the field variable $\Delta$ (given in (68)) and studied the eigenmodes of the G13 theory essentially.

We have demonstrated in § 6.1 that in the inelastic ($e\neq 1$) case, one eigenmode of both the longitudinal and transverse systems (105)—a heat mode of the longitudinal system (105)1 and a shear mode of the transverse system (105)2—is unstable below some critical wavenumber values (see the right panels of figures 7 and 9). To analyze the stability of the longitudinal and transverse systems (105), we investigate the complex frequency of the least stable eigenmode in each system. We have found (but not shown here) that the real part of the complex frequency of the least stable eigenmode in each system is always zero; therefore, the instabilities of both the longitudinal and transverse systems are stationary. This means that the complex frequency $\omega$ vanishes (i.e., $\omega=0$) for the least stable eigenmode in each of the longitudinal and transverse systems. Consequently, the principle of exchange of instabilities (Drazin & Reid, 1981) is valid for both the systems.

Figures 13 and 14 illustrate the zero contours of (a) the least stable heat mode of the longitudinal system and (b) the least stable shear mode of the transverse system associated with the G13 and G26 equations, respectively in $(e,k)$-plane. The black solid line in both the figures represents the critical wavenumber below which the system is unstable (region depicted with gray color) and above which the system is stable (region shown with white color). The left panels of figures 13 and 14 unveil that the longitudinal system associated with the G13 equations is unstable for $0\leq e\lesssim 0.21$ and for $0\leq k\leq 1$ (in fact, for any wavenumber apparently) whereas that associated with the G26 equations is stable above some critical wavenumber for all values of $e$. Thus, the G13 theory of Kremer & Marques Jr. (2011) may not be suitable to the granular flows having coefficient of restitution $e\lesssim 0.35$ (the point where the critical wavenumber from the G13 theory attains a sudden jump). The critical wavenumber in the transverse system associated with the G13 equations (right panel of figure 13)—for small values of the coefficient of restitution ($0\leq e\lesssim 0.4$)—increases with increasing $e$, however, that in the transverse system associated with the G26 equations (right panel of figure 14) decreases with increasing $e$ in the same region. Moreover, the instability region is more in the transverse system associated with the G26 equations than that in the transverse system associated with the G13 equations (cf. the right panels of figures 13 and 14).

In order to compare the critical wavenumbers for the longitudinal and transverse systems associated with various moment theories (G13, G14, G20, G21 and G26), we again plot the zero contours of the least stable eigenmode of the (a) longitudinal and (b) transverse systems associated with these moment theories in $(e,k)$-plane in figure 15. Recall that the transverse systems for the G13 and G14 equations are same and those for the G20 and G21 equations are also same. Therefore, the curves representing the G14 (in red color) and G21 (in black color) theories in right panel of figure 15 also represent the G13 and G20 theories, respectively. Apparently, all moment theories predict same critical wavenumber for the longitudinal system when $0.94\lesssim e\leq 1$ (see the left panel of figure 15) and same critical wavenumber for the transverse system when $0.75\lesssim e\leq 1$ (see the right panel of figure 15). From the right panel of figure 15, we see that for the transverse system, the critical wavenumber profiles predicted by all Grad moment theories are qualitatively similar except for that in the region where $0\leq e\lesssim 0.4$; in this region, the critical wavenumber predicted by the G14 equations increases with increasing $e$ while that predicted by the G21 equations remains more or less constant and that predicted by the G26 equations decreases with increasing $e$. Also, the stability region for the transverse system is decreasing as the number of moments in the system are increasing. On the other hand, the critical wavenumber profiles for the longitudinal system predicted by various moment theories are quite different from each other, especially for small values of the coefficient of restitution ($0\leq e\lesssim 0.4$). From the left panel of figure 15, it seems that based on the critical wavenumber profiles, one can classify the Grad moment theories into two groups: one containing the G13 and G20 theories and the other containing the G14, G21 and G26 theories. Notice that the first group of theories does not contain the perturbed part of the scalar moment $\tilde{\Delta}$ in contrast to the other group of theories. One can see that the critical wavenumber profiles predicted by the G13 and G20 theories, which do not contain $\tilde{\Delta}$, are very different from those predicted by the G14, G21 and G26 theories, which contain $\tilde{\Delta}$. The critical wavenumber from the G20 theory closely follows that from the G13 theory for $0.26\lesssim e\leq 1$ and they both have a sudden jump at around $e\approx 0.35$; nevertheless, for $0\leq e\lesssim 0.21$, the G13 theory does not give any critical wavenumber whereas the G20 theory does give critical wavenumbers, which means that the longitudinal system associated with G13 equations is always unstable for $0\leq e\lesssim 0.21$ while that associated with G20 equations is stable above a critical wavenumber for all value of $e$. This also means that the addition of more moments into the system is stabilizing the longitudinal system associated with the G13 equations. In the other group of theories (G14, G21 and G26), the critical wavenumber profiles from the longitudinal system associated with the G14 and G21 are qualitatively similar—including kinks at $e\approx 0.24$ and $e\approx 0.23$ in critical wavenumber profiles for G14 and G21 theories, respectively—since they both contain $\tilde{\Delta}$ and do not contain $\hat{R}_{ij}$. The instability region for the longitudinal system associated with the G26 equations is also more than that associated with G14 and G21 equations. It may also be stated from the left panel of figure 15 that the number of moments ought to be increased as the inelasticity increases; while stating this, we have ignored the G20 theory (magenta line in the figure) which is, any way, not meaningful for granular gases since it does not contain the scalar fourth order moment.

6.4 Critical length

It is well-known (and have also been discussed above) that the HCS of a freely cooling granular gas is unstable in general, and leads to the formation of velocity vortices (due to the instability of shear mode) and, subsequently, to the formation of density clusters (due to the instability of heat mode). Nevertheless, it is also known that these instabilities are confined to long wavelengths (or to small wavenumbers), and thus may not be observed in systems having small enough system size, (see e.g., Brey et al., 1998b; Brilliantov & Pöschel, 2004; Garzó, 2005). Here, we want to find the critical system size (length), above which the system will become unstable, through the G26 theory presented in this work.

Recall from § 6.3 that the complex frequency $\omega$ is zero for the least stable modes of both the longitudinal and transverse systems. Therefore the critical wavenumber for the least stable shear mode of the longitudinal system and that for the least stable heat mode of the transverse system can be obtained directly by inserting $\omega=0$ in (108) and (109), respectively. Inserting $\omega=0$ in (108) and (109), they yield $\mathbbm{a}_{8}=0$ and $\mathbbm{b}_{5}=0$, respectively. Note that $\mathbbm{a}_{8}$ is a six degree polynomial in $k$ while $\mathbbm{b}_{5}$ is a four degree polynomial in $k$. Nevertheless, the solution of $\mathbbm{a}_{8}=0$ leads to only one meaningful value of $k$, which is the critical wavenumber $k_{h}$ for the least stable heat mode of the longitudinal system (105)1; similarly, the solution of $\mathbbm{b}_{5}=0$ also leads to only one meaningful value of $k$, which is the critical wavenumber $k_{s}$ for the least stable shear mode of the transverse system (105)2. The explicit expressions for $k_{h}$ and $k_{s}$, in a compact form, can be written as

[TABLE]

where

[TABLE]

The plots of the critical wavenumbers $k_{h}$ and $k_{s}$ over the coefficient of restitution $e$ are exactly same as the solid line in figure 14 (a) and (b), respectively. Thus, a heat mode of the longitudinal system (105)1 in the regime $k>k_{h}$ will always decay while that in the regime $k<k_{h}$ will grow exponentially; similarly a shear mode of the transverse system (105)2 in the regime $k>k_{s}$ will always decay while that in the regime $k<k_{s}$ will grow exponentially. Notice from figure 14 that $k_{s}>k_{h}$ for all values of the coefficient of restitution.

In the normal mode analysis considered above, the wavevector is assumed to be in the $x$-direction; thus the $x$-direction is periodic with period $2\pi/k$. Consequently, the smallest admissible wavenumber in a system imposed with periodic boundary conditions is $2\pi/L$, where $L$ is the largest system size (length). Hence, corresponding to the highest critical wavenumber, one can determine a critical system size (length) $L_{c}$ such that the system is stable for $L<L_{c}$ while unstable for $L>L_{c}$. In other words, the critical length $L_{c}$ for the system considered in this work is determined by

[TABLE]

is the dimensionless critical length and $\ell$ is the length scale defined in (91) (recall that the wavenumber—and hence $k_{h}$ and $k_{s}$—are dimensionless and the length scale for making them dimensionless was $\ell$). Since we have concluded above that $k_{s}>k_{h}$, the critical length $L_{c}$ is given by

[TABLE]

where $\ell_{0}=1/(\sqrt{2}\,\pi n_{0}d^{2})$ is the mean free path of a hard-sphere dilute gas.

Figure 16 illustrates the critical length of the system in units of mean free path ($L_{c}/\ell_{0}$) as a function of the coefficient of restitution $e$. The solid (black) line denotes the critical length (in units of mean free path) computed from (118) with expression (116) for $k_{s}$, which has been obtained through the G26 theory presented above. The dashed (green) line denotes the critical length (in units of mean free path) computed from theoretical expression obtained by the linear stability analysis of Navier–Stokes and Fourier equations in Garzó (2005). It is evident from the figure that the results on critical length from the G26 (black solid line) and Navier–Stokes–Fourier (green dashed line) theories agree for $0.55\lesssim e\leq 1$ while differ significantly for more inelastic systems as the latter is not reliable for the systems with large inelasticity. The squares delineate the critical length (in units of mean free path) computed from the theoretical expression given in Brey et al. (1998b), which is obtained by the linear stability analysis of a kinetic model due to Brey et al. (1997). Clearly, the results from the G26 theory and those from the kinetic model by Brey et al. (1997) agree perfectly for all values of $e$.

The simulation results (depicted by circles and triangles) in figure 16 are included only to assess the results obtained from the G26 theory qualitatively. The circles depict the results of two-dimensional DSMC simulations carried out by Brey et al. (1998b). One can see in Brey et al. (1998b) that these simulation results are in excellent agreement with the theoretical results from the linear stability of the kinetic model by Brey et al. (1997), and the latter (shown by squares) are in perfect agreement with those from the G26 theory (black solid line). Therefore it can be expected that the results from full three-dimensional DSMC simulation of a dilute granular gas on critical length would agree with those predicted by the G26 theory (black solid line). The triangles represent the results of molecular dynamic simulations carried out by Mitrano et al. (2011) at solid fraction $\phi=0.1$. Note that the results from the G26 theory presented above are in the dilute limit ($\phi\to 0$). One may notice from Mitrano et al. (2011) that their results are also in excellent agreement with those obtained from the theoretical expression of Garzó (2005) at $\phi=0.1$ for $0.5\leq e\leq 1$, and the latter (green dashed line)—for $\phi=0$—again agree perfectly with those obtained from the G26 theory for $0.5\leq e\leq 1$. Therefore, it is stated that the G26 theory also provides the correct critical length for all values of the coefficient of restitution.

7 Conclusion

Higher-order Grad moment equations—up to first 26-moments—for granular gases have been derived by employing the Grad’s method of moments to the inelastic Boltzmann equation, and the production terms associated with the derived moment equations have been presented. The production terms have, in fact, been computed and presented for the G27 system (which consists of the first 26-moments and the fully contracted sixth moment). The production terms associated with the G26 equations can easily be computed by dropping the terms containing the fully contracted sixth moment in the given production terms associated with the G27 equations. The HCS of a freely cooling granular gas has been explored with the G26 equations in a semi-linearized setting and it has been found that the temperature decay in the HCS closely follows Haff’s law while the other higher-order moments decay on a faster time scale. Further, the nonlinear terms of fully contracted fourth moment have also been included and, by exploiting the stability analysis of fixed points in a dynamical system, it has been shown that some of the fixed points of the system are unstable and with the only stable fixed point it has been concluded that the nonlinear terms, indeed, have only negligible effect on Haff’s law. The G27 equations has also been scrutinized and the stability analysis of fixed points in a dynamical system has again been exploited to deduce that even the inclusion of scalar sixth order moment into the G26 system has negligible effect on Haff’s law. By following the approach of Garzó (2013), the transport coefficients in the Navier–Stokes and Fourier laws for dilute granular gases have been determined through the G26 equations and compared with those obtained from various theoretical methods and computer simulations. The reduced shear viscosity obtained in the present work has been found to be in good agreement with that from the computer-aided method of Noskowicz et al. (2007) and with the DSMC simulations of Garzó et al. (2007) for almost all values of the coefficient of restitution. In fact, the present work yields slightly better results for the reduced shear viscosity in comparison to the modified Sonine approach (Garzó et al., 2007). However, in contrast to the modified Sonine approach (Garzó et al., 2007), which yields good results for all the transport coefficients, the other two transport coefficients from the present work agree only with the first Sonine approximation, hence overestimate the corresponding coefficients for $e\lesssim 0.7$. This suggests to include the full trace of fifth moment into the G26 system in order to get coupling in the semi-linearized RHS of the heat flux balance equation. It is expected that the inclusion of the full trace of fifth moment—i.e., Grad 29-moment system—would improve the reduced transport coefficients $\kappa^{\ast}$ and $\lambda^{\ast}$, although this problem will be considered somewhere else in future.

The linear stability of the HCS has been analyzed through the G26 system and various sub-systems by decomposing them into the longitudinal and transverse systems. It has been found that a heat mode of the longitudinal system and a shear mode of the transverse system in the case of an inelastic gas are unstable whereas all the eigenmodes for both the systems are stable in the case of an elastic gas. By comparing the least stable eigenmodes from various theories, it has been established that the unstable heat mode of the longitudinal system obtained with the 13-field eigenmode theory of Kremer & Marques Jr. (2011) remains unstable for all wavenumbers below a certain coefficient of restitution whereas that obtained with any other higher-order moment theory becomes stable above some critical wavenumber for all values of the coefficient of restitution. This is apparently an artefact of assuming the dimensionless fourth moment as a constant in Kremer & Marques Jr. (2011) rather than considering it as a field variable while analyzing the eigenmodes. Out of all the theories considered, the G26 theory have produced the smoothest critical wavenumber profile. The critical wavenumber profiles from various Grad moment theories also suggested that the number of moments ought to be increased with increasing inelasticity. Further investigation of the critical wavenumbers from the G26 theory has unveiled that the value of the critical length—beyond which the system becomes unstable—is essentially driven by the unstable shear mode of the transverse system. The critical length profile obtained from the G26 theory is in excellent agreement with those obtained by the existing theories and also in qualitatively good agreement with those from the simulations.

Acknowledgments

The authors acknowledge anonymous reviewers for their valuable suggestions, which significantly improved the paper. VKG gratefully acknowledges financial support from the SRM University in the form of seed grant through the project “Grain”. PS acknowledges financial support from IIT Madras in the form of New Faculty Initiation Grant (MAT/15-16/833/NFIG/PRIY) and New Faculty Seed Grant (MAT/16-17/671/NFSC/PRIY).

Appendix A Production terms

The fully nonlinear production terms associated with the G27 equations using definition $\Delta=w/(\rho\,\theta^{2})$ are given below. The corresponding production terms associated with the G26 equations (eqs. (11)–(13) and (19)–(3.1)) are obtained from (A)–(A) on taking $\Xi=0$ and ignoring (A).

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

is the collision frequency.

Appendix B Governing equations for perturbed field variables

Inserting the field variables from (90) into the G26 equations (30)–(4) (without the underlined term) and neglecting all the nonlinear terms in perturbations, one obtains the following system of linear partial differential equations in the perturbed field variables.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]