Moment theories for a d-dimensional dilute granular gas of Maxwell molecules

TL;DR

This paper derives and analyzes moment equations for a d-dimensional dilute granular gas of Maxwell molecules, providing insights into the Navier-Stokes relations and stability of the homogeneous cooling state.

Contribution

It introduces a comprehensive set of moment equations for granular gases and demonstrates their sufficiency for Navier-Stokes relations and stability analysis.

Findings

Moment equations derived up to (d+3)(d^2+6d+2)/6 moments.

Navier-Stokes constitutive relations depend only on certain moments.

Critical system size for instability matches simulation results.

Abstract

Various systems of moment equations---consisting of up to (d+3)(d^2+6d+2)/6 moments---in a general dimension d for a dilute granular gas composed of Maxwell molecules are derived from the inelastic Boltzmann equation by employing the Grad moment method. The Navier--Stokes-level constitutive relations for the stress and heat flux appearing in the system of mass, momentum and energy balance equations are determined from the derived moment equations. It has been shown that the moment equations only for the hydrodynamic field variables (density, velocity and granular temperature), stress and heat flux---along with the time-independent value of the fourth cumulant---are sufficient for determining the Navier--Stokes-level constitutive relations in the case of inelastic Maxwell molecules, and that the other higher-order moment equations do not play any role in this case. The homogeneous cooling state of a freely cooling granular gas is investigated with the system of the Grad (d+3)(d^2+6d+2)/6-moment equations and its various subsystems. By performing a linear stability analysis in the vicinity of the homogeneous cooling state, the critical system size for the onset of instability is estimated through the considered Grad moment systems. The results on critical system size from the presented moment theories are found to be in reasonably good agreement with those from simulations.

Tables

Table 1. Table 1: Number of unknown field variables in Grad moment systems in d 𝑑 d dimensions

Field variables	Unknowns in $d$ dimensions	Unknowns in $3$ dimensions	Unknowns in $2$ dimensions
$\rho$	$1\text{rdelim}\}5-5mm[$d2+5d+22!$\left]\text{rdelim}\right\}8-5mm[$(d+1)(d2+8d+6)3!$]$	1 \rdelim}50pt[  13]  \rdelim}80mm[26]	1 \rdelim}50pt[  8]  \rdelim}80mm[13]
$v_i$	$d$	3	2
$\theta$	$1$	1	1
${\sigma }_{ij}$	$\frac{\left(d+2\right)\left(d-1\right)}{2!}$	5	2
$q_i$	$d$	3	2
$u_{ijk}^0$	$\frac{\left(d+4\right)d\left(d-1\right)}{3!}$	7	2
$u^2$	$1$	1	1
$u_{ij}^1$	$\frac{\left(d+2\right)\left(d-1\right)}{2!}$	5	2
$u_i^2$	$d$	3	2
Total	$\frac{\left(d+3\right)\left(d^2+6d+2\right)}{3!}$	29	15

Table 2. Table 2: Threshold values of the coefficient of restitution e th subscript 𝑒 th e_{\mathrm{th}} below which the longitudinal and transverse systems for IMM remain always unstable.

		NSF	G13	G14	G26	G29
Longitudinal	$d=2$	0.62798	0.60211	0.52174	0.37473	0.56356
	$d=3$	0.46551	0.46033	0.38608	0.06256	0.40157
Transverse	$d=2$	–	–	–	0.41360	0.32349
	$d=3$	–	–	–	0.23030	0.16867

Equations

$$\frac{\partial f}{\partial t}+c_i\frac{\partial f}{\partial x_i}+F_i\frac{\partial f}{\partial c_i}=J[\mathbf{c}|f,f],$$

$$J[\mathbf{c}|f,f]=\frac{\mathring{\nu }}{n{\Omega }_d}{\int }_{{\mathbb{R}}^d}{\int }_{S^{d-1}}\left[\frac{1}{e}f(t,\mathbf{x},{\mathbf{c}}^{\prime \prime })f(t,\mathbf{x},{\mathbf{c}}_1^{\prime \prime })-f(t,\mathbf{x},\mathbf{c})f(t,\mathbf{x},{\mathbf{c}}_1)\right]\mathrm{d}\hat{\mathbf{k}}\mathrm{d}{\mathbf{c}}_1.$$

$$\mathring{\nu }=\frac{d+2}{2}\nu ,\text{where}\nu =\frac{4{\Omega }_d}{\sqrt{\pi }(d+2)}n\text{mathbbm}{d}^{d-1}\sqrt{\frac{T}{m}}$$

$${\mathbf{c}}^{\prime \prime }=\mathbf{c}-\frac{1+e}{2e}(\hat{\mathbf{k}}\cdot \mathbf{g})\hat{\mathbf{k}}\text{and}{\mathbf{c}}_1^{\prime \prime }={\mathbf{c}}_1+\frac{1+e}{2e}(\hat{\mathbf{k}}\cdot \mathbf{g})\hat{\mathbf{k}},$$

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

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

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

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

$$\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$$

$$\frac{\mathrm{D}T}{\mathrm{D}t}+\frac{2}{d}\frac{1}{n}\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 =-\frac{m}{dnT}\int C^{2}J[\mathbf{c}|f,f]\mathrm{d}\mathbf{c}.$$

$${\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},$$

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

$${\mathcal{P}}_{i_1{i}_2\dots i_r}^a:=m\int C^{2a}{C}_{⟨i_1}{C}_{i_2}\dots C_{i_r⟩}J[\mathbf{c}|f,f]\mathrm{d}\mathbf{c},$$

$$\left(\genfrac{}{}{0ex}{}{d+r-1}{r}\right)-\left(\genfrac{}{}{0ex}{}{d+r-3}{r-2}\right)=\frac{d+2r-2}{d+r-2}\left(\genfrac{}{}{0ex}{}{d+r-2}{r}\right).$$

$$\begin{array}{rlrl}{m}_{ijk}& :=u_{ijk}^0,& \Delta & :=\frac{u^2}{d(d+2)\rho {\theta }^2}-1,\\ R_{ij}& :=u_{ij}^1-(d+4)\theta {\sigma }_{ij},& {\phi }_i& :=u_i^2-4(d+4)\theta q_i\end{array}\}$$

$$\frac{\mathrm{D}{\sigma }_{ij}}{\mathrm{D}t}+\frac{\partial m_{ijk}}{\partial x_k}+\frac{4}{d+2}\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}$$

$$+\theta \frac{\partial {\sigma }_{ij}}{\partial x_j}-\frac{{\sigma }_{ij}}{\rho }\left(\frac{\partial {\sigma }_{jk}}{\partial x_k}-\theta \frac{\partial \rho }{\partial x_j}\right)+m_{ijk}\frac{\partial v_j}{\partial x_k}$$

$$+\frac{d+4}{d+2}{q}_i\frac{\partial v_j}{\partial x_j}+\frac{d+4}{d+2}{q}_j\frac{\partial v_i}{\partial x_j}+\frac{2}{d+2}{q}_j\frac{\partial v_j}{\partial x_i}=\frac{1}{2}{\mathcal{P}}_i^1,$$

$$\frac{\mathrm{D}{m}_{ijk}}{\mathrm{D}t}$$

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

$$\frac{\mathrm{D}\Delta }{\mathrm{D}t}$$

$$+\frac{1}{d(d+2)}\frac{1}{\rho {\theta }^2}\left[\frac{\partial {\phi }_i}{\partial x_i}+4(d+2)q_i\frac{\partial \theta }{\partial x_i}-8\frac{q_i}{\rho }\left(\frac{\partial {\sigma }_{ij}}{\partial x_j}+\theta \frac{\partial \rho }{\partial x_i}\right)+4R_{ij}\frac{\partial v_i}{\partial x_j}\right]$$

$$\displaystyle=\frac{1}{d(d+2)}\frac{1}{\rho\theta^{2}}\Big{[}\mathcal{P}^{2}-2(d+2)(1+\Delta)\theta\,\mathcal{P}^{1}\Big{]},$$

$$\frac{\mathrm{D}{R}_{ij}}{\mathrm{D}t}$$

$$+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}{d}\theta {\sigma }_{ij}\frac{\partial v_k}{\partial x_k}-\frac{2(d+4)}{d}\frac{{\sigma }_{ij}}{\rho }\left(\frac{\partial q_k}{\partial x_k}+{\sigma }_{kl}\frac{\partial v_k}{\partial x_l}\right)$$

$$+\frac{\partial u_{ijk}^1}{\partial x_k}-(d+4)\theta \frac{\partial m_{ijk}}{\partial x_k}-2\frac{m_{ijk}}{\rho }\left(\frac{\partial {\sigma }_{kl}}{\partial x_l}+\rho \frac{\partial \theta }{\partial x_k}+\theta \frac{\partial \rho }{\partial x_k}\right)+2u_{ijkl}^0\frac{\partial v_k}{\partial x_l}$$

$$+\frac{d+6}{d+4}\left(R_{ij}\frac{\partial v_k}{\partial x_k}+2R_{k⟨i}\frac{\partial v_{j⟩}}{\partial x_k}\right)+\frac{4}{d+4}{R}_{k⟨i}\frac{\partial v_k}{\partial x_{j⟩}}+2(d+4)\Delta \rho {\theta }^2\frac{\partial v_{⟨i}}{\partial x_{j⟩}}$$

$$={\mathcal{P}}_{ij}^1-(d+4)\theta {\mathcal{P}}_{ij}^0-\frac{d+4}{d}\frac{{\sigma }_{ij}}{\rho }{\mathcal{P}}^1,$$

$$\frac{\mathrm{D}{\phi }_i}{\mathrm{D}t}$$

$$\displaystyle-4R_{ij}\frac{\partial\theta}{\partial x_{j}}-(d+4)\big{[}(d+6)+(d+2)\Delta\big{]}\theta^{2}\frac{\partial\sigma_{ij}}{\partial x_{j}}-2(d+4)^{2}\theta\sigma_{ij}\frac{\partial\theta}{\partial x_{j}}$$

$$-(d+2)(d+4)\left[2\rho {\theta }^3\frac{\partial \Delta }{\partial x_i}+(1+3\Delta ){\theta }^3\frac{\partial \rho }{\partial x_i}+(3+5\Delta )\rho {\theta }^2\frac{\partial \theta }{\partial x_i}\right]$$

$$-4\frac{R_{ij}}{\rho }\left(\frac{\partial {\sigma }_{jk}}{\partial x_k}+\theta \frac{\partial \rho }{\partial x_j}\right)+4u_{ijk}^1\frac{\partial v_j}{\partial x_k}-4(d+4)\theta m_{ijk}\frac{\partial v_j}{\partial x_k}$$

$$+\frac{8(d+4)}{d+2}\theta \left(q_i\frac{\partial v_j}{\partial x_j}+q_j\frac{\partial v_i}{\partial x_j}+q_j\frac{\partial v_j}{\partial x_i}\right)+\frac{d+6}{d+2}{\phi }_i\frac{\partial v_j}{\partial x_j}+\frac{d+6}{d+2}{\phi }_j\frac{\partial v_i}{\partial x_j}$$

$$+\frac{4}{d+2}{\phi }_j\frac{\partial v_j}{\partial x_i}={\mathcal{P}}_i^2-2(d+4)\theta {\mathcal{P}}_i^1-\frac{4(d+4)}{d}\frac{q_i}{\rho }{\mathcal{P}}^1.$$

$$f_{|G13}$$

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

$${\mathcal{P}}^1=$$

Full text

1 Introduction

Under strong excitation, granular materials resemble ordinary (molecular) gases and are referred to as rapid granular flows or granular gases (Campbell, 1990; Goldhirsch, 2003). The prototype model of a granular gas is a dilute system comprised of smooth (frictionless) identical hard spheres—with no interstitial fluid—colliding pairwise and inelastically with a constant coefficient of (normal) restitution $0\leq e\leq 1$, with $e=0$ referring to perfectly sticky collisions and $e=1$ to perfectly elastic collisions (Campbell, 1990; Goldhirsch, 2003; Brilliantov & Pöschel, 2004; Garzó, 2019). In the dilute limit, this system can be described by a single-particle velocity distribution function, which is the fundamental quantity in kinetic theory and obeys the Boltzmann equation suitably modified to incorporate energy dissipation due to inelastic collisions. The resemblance of granular gases to ordinary gases has motivated the development of several kinetic theory based tools for granular gases by suitably modifying these tools to account for energy dissipation due to inelastic collisions in the last three decades, and it is still an active area of research; see, e.g., Jenkins & Richman (1985a, b); Sela & Goldhirsch (1998); Brey et al. (1998a); Garzó & Santos (2003); Santos (2003); Brilliantov & Pöschel (2004); Bisi et al. (2004); Garzó & Santos (2011); Kremer & Marques Jr. (2011); Garzó (2013); Khalil et al. (2014); Kremer et al. (2014); Gupta & Shukla (2017); Gupta et al. (2018); Garzó (2019). Nevertheless, the non-conservation of energy in granular systems makes kinetic theory based tools much more involved and has profound consequences on their behaviour, leading to a raft of intriguing phenomena pertaining to granular matter.

Kinetic theory of classical (monatomic) gases offers systematic ways of deriving the transport equations for the field variables. The two notable approaches in kinetic theory, around which various solution techniques and some other models have been developed, are the Chapman–Enskog (CE) expansion (Chapman & Cowling, 1970) and the Grad moment method (Grad, 1949). While these approaches have been instrumental in understanding several problems from a theoretical point of view, both have their own shortcomings. The former, which is adequate for flows close to equilibrium, considers the transport equations only for the hydrodynamic field variables (density, velocity and temperature) and provides the constitutive relations for additional unknowns, namely the stress and heat flux, in these equations. Despite being successful in deriving the Euler equations (at zeroth order of expansion) and the classical Navier–Stokes and Fourier (NSF) equations (at first order of expansion), the usefulness of models (the Burnett equations and beyond) resulting from the higher-order CE expansion remains scarce mainly due to inherent instabilities (Bobylev, 1982). On the other hand, the Grad moment method (Grad, 1949) furnishes the governing equations—referred to as the moment equations—for more field variables than the hydrodynamic ones and employs a Hermite polynomial expansion to close the system of moment equations. A set of moment equations emanating from the Grad moment method is always linearly stable but suffers from the loss of hyperbolicity (Müller & Ruggeri, 1998; Cai et al., 2014b)—an essential property for the well-posedness of a system of partial differential equations (PDEs). The loss of hyperbolicity renders a Grad moment system to show some unphysical behaviour, e.g. unphysical sub-shocks within the shock profile above a critical Mach number. Yet, the Grad moment method has a clear advantage of linearly stable equations and, hence, is preferred over the CE expansion for describing nonequilibrium flows of monatomic gases.

To circumvent the problems associated with the Grad moment method, a number of moment methods have been proposed in the literature. Levermore (1996) propounded the maximum-entropy approach for closing a moment system. Although the maximum-entropy approach of Levermore (1996) produces hyperbolic systems of moment equations by construction, it is extremely difficult to obtain Levermore’s moment equations in an explicit form beyond the 10-moment case (which does not include the heat flux) because the fluxes associated with higher moments cannot be expressed in a closed form. As Levermore’s 10-moment system is not capable of describing heat conduction, it is not very useful for describing gaseous processes. In addition, larger moment systems resulting from the maximum-entropy approach are prone to serious mathematical issues (Junk, 1998; Junk & Unterreiter, 2002). To alleviate problems associated with the maximum-entropy approach, McDonald & Torrilhon (2013) proposed affordable numerical approximations to the maximum-entropy closures for problems involving heat transfer and presented a robust and affordable version of Levermore’s 14-moment system (that includes the heat flux). Although the 14-moment system proposed by McDonald & Torrilhon (2013) is capable of predicting accurate and smooth shock structures even for relatively large Mach numbers, it is not globally hyperbolic. In a completely different approach that focused on producing hyperbolic moment equations, Torrilhon (2010) introduced a novel closure computed with multi-variate Pearson-IV-distributions for the 13-moment system; however the approach seems unlikely to work for higher-order moment systems. In another approach, Struchtrup & Torrilhon (2003) introduced a model, termed as the regularised 13-moment (R13) equations, which regularises the original Grad 13-moment (G13) equations by employing a CE-like expansion around a pseudo-equilibrium state. Subsequently, the R13 equations were also derived in an elegant and clean way by Struchtrup (2005) via the order of magnitude approach. The R13 equations are linearly stable, predict smooth shock structure for all Mach numbers and can capture several rarefaction effects, such as Knudsen layers, with good accuracy for sufficiently small Knudsen numbers (Struchtrup & Torrilhon, 2003; Torrilhon & Struchtrup, 2004). To cover more transition-flow regime, Gu & Emerson (2009) employed the regularisation approach of Struchtrup & Torrilhon (2003) to derive the regularised 26-moment (R26) equations. It may be noted, however, that the R13 and R26 equations are also not hyperbolic. The Grad and regularised moment equations consisting of an arbitrary number of moments have also been implemented and solved numerically by Cai & Li (2010). In the last few years, various other regularisation techniques that yield globally hyperbolic moment equations have been introduced. Cai et al. (2013) proposed a regularisation of the Grad moment equations in one space dimension based on investigating the properties of the Jacobian matrix of fluxes in the system and derived globally hyperbolic moment equations (HME) in one space dimension. Further, Cai et al. (2014a) generalised the method to derive HME in multi-dimension. Koellermeier et al. (2014) employed quadrature-based projection methods, which alter the structure of a moment system in a desired way, to obtain hyperbolic systems of the so-called quadrature-based moment equations (QBME). Fan et al. (2016) proposed a generalised framework, which is capable of deriving various existing as well as some new systems of regularised hyperbolic moment equations, based on the so-called operator projection method. A remarkable drawback of HME and QBME is that they cannot be written in a conservative form (Koellermeier & Torrilhon, 2017). Consequently, the standard finite volume schemes cannot be applied to solve systems of HME and QBME numerically. Recently, some non-conservative numerical schemes have been proposed by Koellermeier & Torrilhon (2017, 2018) for the numerical solution of QBME in one and two dimensions. While numerical methods for solving general three-dimensional unsteady flow problems with moment equations are still intractable, the method of fundamental solution (MFS) enables us to develop efficient meshfree numerical methods for solving three-dimensional steady flow problems with the linearised moment equations. Recently, the MFS has been developed for the G13 and R13 equations in Lockerby & Collyer (2016) and Claydon et al. (2017), respectively.

The system of the R13 equations (also the system of the R26 equations), despite being non-hyperbolic, may be regarded as the most promising continuum model for describing rarefied monatomic gas flows since it is accompanied with appropriate boundary conditions (Gu & Emerson, 2007; Torrilhon & Struchtrup, 2008), and has already been successful in describing a number of canonical flows (see Torrilhon, 2016, and references therein). Motivated from the accomplishments of the moment method (in particular, the R13 equations) in the case of monatomic gases, the Grad and regularised moment equations have also been developed for monatomic gas mixtures (Gupta & Torrilhon, 2015; Gupta, 2015; Gupta et al., 2016). It is important to note here that the derivation of the regularised moment equations requires higher-order Grad moment equations, for instance, the derivations of the R13 and R26 equations require the Grad 26-moment (G26) and Grad 45-moment equations, respectively, and that most of the aforementioned works on the moment method employ either some simplified kinetic models to replace the Boltzmann collision operator or the Maxwell potential for molecular interactions. The latter, introduced by Maxwell (1867), is inversely proportional to the fourth power of the distance between the colliding molecules and makes the collision rate independent of the relative velocity between the colliding molecules, which greatly simplifies the original Boltzmann equation. Remarkably, for Maxwell molecules (i.e. for molecules interacting with the Maxwell interaction potential), the collisional production terms—the terms emanating from the Boltzmann collision operator in the moment equations—can be computed without the knowledge of explicit form of the distribution function and, moreover, they turn out to be bilinear combinations of moments of the same or lower order, resulting into a one-way coupling on the right-hand sides of a moment system. This makes the moment equations for Maxwell molecules tractable. For more details on the moment method for Maxwell molecules, the reader is referred to a review paper by Santos (2009).

The development of kinetic theory of granular gases started out with two seminal works by Jenkins & Savage (1983) and Lun et al. (1984), which introduced kinetic theory for smooth inelastic hard spheres (IHS), followed by the pioneering work of Jenkins & Richman (1985b) on kinetic theory for rough inelastic hard disks (IHD). The aforementioned methods, namely the CE expansion and the Grad moment method, in kinetic theory of classical gases have also been extended to granular gases, with the main goal of determining the NSF-level transport coefficients appearing in the expressions for the stress and heat flux, since the hydrodynamic equations closed with the NSF-level constitutive relations are sufficient to describe flows involving small spatial gradients. The CE expansion to zeroth order was first employed by Goldshtein & Shapiro (1995) to obtain the Euler-like hydrodynamic equations for rough granular flows. Subsequently, Brey et al. (1998a) and Garzó & Dufty (1999) determined the NSF-level transport coefficients for dilute and dense granular gases of IHS, respectively, by means of the first-order CE expansion in powers of a uniformity parameter that estimates the strength of spatial gradients of the hydrodynamic field variables. The derivation of Burnett equations (i.e. the second-order CE expansion) even for the prototype model of a granular gas is an arduous task. Yet, by performing a generalised CE expansion in powers of two small parameters, namely the Knudsen number and the degree of inelasticity, Sela & Goldhirsch (1998) determined the constitutive relations for the stress and heat flux up to Burnett order for a smooth granular gas of IHS. The requirement of the degree of inelasticity being small for performing asymptotic expansion limits the validity of Burnett equations derived by Sela & Goldhirsch (1998) to nearly elastic granular gases. Lutsko (2005) further extended the CE expansion to dense granular fluids with arbitrary energy loss models and determined the NSF-level constitutive relations. Not only did his work consider arbitrary inelasticity but also a velocity-dependent coefficient of restitution, providing the NSF-level constitutive relations for more realistic granular fluids.

Granular flows of interest often fall beyond the regime covered by Newtonian hydrodynamics since the strength of spatial gradients in flows of practical interest is not small due to the inherent coupling between the spatial gradients and inelasticity (Goldhirsch, 2003). Consequently, for such flows, the granular NSF equations obtained from the first-order CE expansion are not adequate and the Burnett equations for IHS are not meaningful due to their validity being restricted to nearly elastic granular gases besides Bobylev’s instability. Such granular flows can alternatively be modelled by the Grad moment method. The method was extended to granular fluids first by Jenkins & Richman who derived the G13 equations for a dense and smooth granular gas of IHS (Jenkins & Richman, 1985a) and the Grad 16-moment equations for a dense and rough granular gas of IHD (Jenkins & Richman, 1985b). It is well-established that the fourth cumulant (scalar fourth moment of the velocity distribution function) ought to be included in the list of the field variables for appropriate description of processes in granular gases; for instance, a theoretical description of the recently observed Mpemba effect in granular fluids requires the fourth cumulant as a field variable (Lasanta et al., 2017). Keeping that in mind, Risso & Cordero (2002) included the fourth cumulant in the list of the field variables to derive the Grad 9-moment equations for a bidimensional granular gas, and utilised them to investigate the homogeneous cooling state (HCS) and the steadily heated state of a bidimensional granular gas. Bisi et al. (2004) attempted to extend the Grad moment method to one-dimensional dilute granular flows of viscoelastic hard spheres. It may be noted that all the aforementioned works on moment method for granular fluids are also restricted to nearly elastic particles. The Grad 14-moment (G14) equations for a dilute granular gas of IHS were introduced by Kremer & Marques Jr. (2011) wherein the authors exploited the G14 equations to obtain the NSF-level constitutive relations for the stress and heat flux via the Maxwell iteration procedure and to investigate the linear stability of the HCS. Although the G14 equations introduced by Kremer & Marques Jr. (2011) were not restricted to nearly elastic particles, their procedure to obtain the constitutive relations did not incorporate the effect of the collisional dissipation. Consequently, the constitutive relations determined by them are valid only for nearly elastic granular gases. The issue was resolved by Garzó (2013) who proposed a procedure to determine the NSF-level constitutive relations incorporating the contributions through the collisional dissipation as well. Although the work of Garzó (2013) yielded the accurate NSF-level constitutive relations for moderately dense granular gases in a general dimension, it only computed the collisional contributions to stress and heat flux exploiting the G14 distribution function but did not provide the G14 equations explicitly. Very recently, Gupta et al. (2018) derived the fully nonlinear G26 equations for dilute granular gases of IHS. Following the approach of Garzó (2013), they determined the NSF-level constitutive relations for the stress and heat flux through the G26 equations. The coefficient of the shear viscosity found by them through the G26 equations turned out to be the best among those obtained via any other theory so far. Notwithstanding, the other transport coefficients related to the heat flux obtained through the G26 equations in Gupta et al. (2018) were exactly the same as those obtained via the CE expansion at the first Sonine approximation (Brey et al., 1998a) or via the G14 distribution function (Garzó, 2013), and the authors adduced that the Grad 29-moment (G29) theory, which includes the flux of the fourth cumulant as field variable, would be able to improve the transport coefficients related to the heat flux.

Despite these ever-improving developments, the fact is that the Boltzmann equation for IHS, and hence the models stemming from the Boltzmann equation for IHS, remains difficult to deal with. To circumvent the difficulties pertaining to models for IHS, a model of inelastic Maxwell molecules (IMM) was proposed at the beginning of this century (Ben-Naim & Krapivsky, 2000; Carrillo et al., 2000; Ernst & Brito, 2002). Similarly to the model of Maxwell molecules for monatomic gases, the IMM model also makes the collision rate of the inelastic Boltzmann equation independent of the relative velocity of the colliding molecules and thereby simplifies the inelastic Boltzmann equation greatly. In the past few years, the IMM model has received tremendous attention as the simple structure of the Boltzmann collision operator for IMM enables us to describe many properties of granular gases analytically, such as the high-velocity tails (Ben-Naim & Krapivsky, 2002; Ernst & Brito, 2002) and the fourth cumulant (Ernst & Brito, 2002; Santos, 2003) in the HCS, and the NSF-level transport coefficients (Santos, 2003). Moreover, the experimental results on the velocity distribution in driven granular gases composed of magnetic grains are well-described by the IMM model (Kohlstedt et al., 2005). The paper by Garzó & Santos (2011) presents a comprehensive review of the IMM model. The two relevant works here are by Santos (2003) and Khalil et al. (2014), which respectively derive the NSF- and Burnett-level transport coefficients for a $d$-dimensional dilute granular gas of IMM by means of the CE expansion. It is worthwhile to note that the work of Khalil et al. (2014), in contrast to that of Sela & Goldhirsch (1998), contains only one smallness parameter, proportional to the spatial gradient of a hydrodynamic field, for performing the CE expansion but is not restricted to nearly elastic granular gases. Nevertheless, as also pointed out in Khalil et al. (2014) as a cautionary note, a regularisation of Burnett equations for IMM is apparently necessary to extricate Bobylev’s instability. Furthermore, as mentioned above, for a proper description of many processes in granular fluids, it is imperative to include the scalar fourth moment as a field variable. Therefore, a moment-based modelling of granular gases seems to be necessary for proper description of processes involving large spatial gradients.

Aiming to the long-term perspective of establishing a complete set of predictive moment equations—for which appropriate boundary conditions, the MFS for steady flow problems and a general numerical framework for unsteady flow problems can be developed—for granular gases, the main objective of this paper is to derive the Grad moment equations—comprising of up to $(d+3)(d^{2}+6d+2)/6$ moments—for a $d$-dimensional dilute (unforced) granular gas of IMM. Here, $d=2$ refers to planar disk flows and $d=3$ to three-dimensional sphere flows. Following the procedure due to Garzó (2013), the NSF-level transport coefficients for a dilute granular gas of IMM are determined from the derived Grad moment equations for IMM. The Grad $(d+3)(d^{2}+6d+2)/6$-moment equations are then utilised to study the HCS of a freely cooling granular gas of IMM. As it is well-known that the HCS of a granular gas is unstable but the instabilities are confined to large systems (see, e.g., Brilliantov & Pöschel, 2004, and references therein), the linear stability of the HCS is investigated with the considered Grad moment systems and the results are employed to estimate the critical system size for the onset of instability.

It is worthwhile to note that a Grad moment system for a dilute granular gas differs from that for a rarefied monatomic gas only on the right-hand sides by virtue of different Boltzmann collision operators, therefore it is expected that a Grad moment system for dilute granular gases, similarly to a Grad moment system for monatomic gases, will also suffer from the loss of hyperbolicity. A detailed investigation of the hyperbolicity of the Grad moment systems derived in this paper and their regularisations will, however, be considered elsewhere in the future. From an application point of view, the Grad moment equations derived in the present work have limited applications at present due to unavailability of the associated boundary conditions, which are beyond the scope of the present paper and will also be considered elsewhere in the future. Nonetheless, the Grad moment equations developed in this paper can be utilised to investigate problems that do not require boundary conditions, e.g. the shock-tube problem, by employing numerical techniques specialised to moment equations developed, for example, in Torrilhon (2006); Cai & Li (2010); Koellermeier & Torrilhon (2017).

The layout of the paper is as follows. The Boltzmann equation for IMM and the basic transport equations (i.e. mass, momentum and energy balance equations) for granular gases of IMM are presented in § 2. The considered Grad moment systems are presented in § 3. The NSF transport coefficients for a dilute granular gas of IMM are determined from the Grad moment equations in § 4. The HCS of a freely cooling granular gas is explored through the Grad moment equations in § 5. The linear stability analysis of the HCS is performed in § 6. The paper ends with a short summary and conclusion in § 7.

2 The Boltzmann equation and the hydrodynamic equations for IMM

We consider a dilute granular gas composed of smooth-identical-inelastic $d$-dimensional spherical particles of mass $m$ and diameter $\mathbbm{d}$. The state of such a gas can be fully described by a single-particle velocity distribution function $f\equiv f(t,\bm{x},\bm{c})$—where $t$, $\bm{x}$ and $\bm{c}$ denote the time, position and instantaneous velocity of a particle, respectively—that obeys the inelastic Boltzmann equation (Brilliantov & Pöschel, 2004)

[TABLE]

where $\bm{F}$ is the external force per unit mass that does not usually depend on $\bm{c}$, $J[\bm{c}|f,f]$ is the (inelastic) Boltzmann collision operator and the Einstein summation applies over repeated indices throughout the paper (unless mentioned otherwise). For $d$-dimensional IMM, the Boltzmann collision operator has a simplified form given by (Ben-Naim & Krapivsky, 2002; Ernst & Brito, 2002; Garzó & Santos, 2007; Garzó, 2019)

[TABLE]

In the Boltzmann collision operator for IMM (2), $e$ is the (constant) coefficient of restitution and $\mathring{\nu}\equiv\mathring{\nu}(e)$—a free parameter in the model—is an effective collision frequency that is typically chosen in such a way that the results from the Boltzmann equations for IHS and IMM agree in an optimal way (Garzó & Santos, 2007). In particular, the agreement of cooling rates from the Boltzmann equations for IHS and IMM leads to (Garzó & Santos, 2007; Khalil et al., 2014)

[TABLE]

is the collision frequency associated with the Navier–Stokes shear viscosity of an elastic (monatomic) gas with $\Omega_{d}=2\pi^{d/2}/\Gamma(d/2)$ being the total solid angle in $d$ dimensions, $n\equiv n(t,\bm{x})$ the number density and $T\equiv T(t,\bm{x})$ the granular temperature, which is a measure of the fluctuating kinetic energy. The velocities $\bm{c}^{\prime\prime}$ and $\bm{c}_{1}^{\prime\prime}$ in (2) are the pre-collisional velocities of the colliding molecules that transform to the post-collisional velocities $\bm{c}$ and $\bm{c}_{1}$ in an inverse collision following the relations (Sela & Goldhirsch, 1998; Brilliantov & Pöschel, 2004):

[TABLE]

where $\bm{g}=\bm{c}-\bm{c}_{1}$ is the relative velocity of the colliding molecules, $\hat{\bm{k}}$ is the unit vector joining the centres of the colliding molecules at the time of collision. The integration limits of $\hat{\bm{k}}$ in (2) extend over the $d$-dimensional unit sphere $S^{d-1}$. Although the limits of integration will be dropped henceforth for the sake of succinctness, an integration over any velocity space will stand for the volume integral over $\mathbb{R}^{d}$ and that over $\hat{\bm{k}}$ will stand for the volume integral over the $d$-dimensional unit sphere $S^{d-1}$.

The hydrodynamic variables—number density $n\equiv n(t,\bm{x})$, macroscopic velocity $\bm{v}\equiv\bm{v}(t,\bm{x})$ and granular temperature $T\equiv T(t,\bm{x})$—relate to the velocity distribution function via

[TABLE]

where $\bm{C}(t,\bm{x},\bm{c})=\bm{c}-\bm{v}(t,\bm{x})$ is the peculiar velocity. The governing equations for the hydrodynamic variables—namely, the mass, momentum and energy balance equations—can be derived from the Boltzmann equation (1) by multiplying it with $1$, $c_{i}$ and $\frac{1}{d}m\,C^{2}$, and integrating each of the resulting equations over the velocity space successively. The mass, momentum and energy balance equations, respectively, read

[TABLE]

where $\mathrm{D}/\mathrm{D}t\equiv\partial/\partial t+v_{k}\,\partial/\partial x_{k}$ is the material derivative. The right-hand sides of the mass and momentum balance equations (8) and (9) vanish due to the conservation of mass and momentum. However, owing to dissipative collisions among grains, the energy is not conserved, yielding a nonzero right-hand side in the energy balance equation (10) with the nonzero cooling rate $\zeta$ being given by

[TABLE]

Furthermore, $\sigma_{ij}\equiv\sigma_{ij}(t,\bm{x})$ and $q_{i}\equiv q_{i}(t,\bm{x})$ in (9) and (10) are the stress tensor and heat flux, respectively, and are given by

[TABLE]

where the angle brackets around the indices denote the symmetric and traceless part of the corresponding tensor; see appendix A for its definition.

Needless to say, the system of mass, momentum and energy balance equations (8)–(10) for the hydrodynamic variables $n$, $v_{i}$ and $T$ is not closed since it encompasses the additional unknowns $\sigma_{ij}$, $q_{i}$ and $\zeta$, and in order to deal with this system any further, it is indispensable to close it. Typically, the closure for system (8)–(10) is obtained by means of the CE expansion, which yields the constitutive relations for $\sigma_{ij}$, $q_{i}$ and $\zeta$ to various orders of approximation (see, e.g., Brey et al., 1998a; Sela & Goldhirsch, 1998; Garzó & Dufty, 1999; Gupta, 2011; Khalil et al., 2014). However, as also stated in § 1, system (8)–(10) closed with the constitutive relations obtained at the zeroth and first orders of the CE expansion is not adequate for describing processes involving large spatial gradients while system (8)–(10) closed with the constitutive relations obtained at the second and higher orders of the CE expansion suffers from Bobylev’s instability. On the other hand, the Grad moment method is capable of yielding more accurate models that do not suffer from Bobylev’s instability and are expected to be valid for processes involving large spatial gradients. Therefore, in what follows, the Grad moment method will be employed for deriving a few closed sets of macroscopic transport equations for a $d$-dimensional dilute granular gas of IMM.

3 Grad moment method

The central goal of the moment method is to have reduced complexity while allowing for more accurate models for rarefied gases. It is well-known that the direct solutions of the Boltzmann equation are computationally expensive since the Boltzmann equation is solved for the velocity distribution function, which depends on total $2d+1$ variables (1 for time, $d$ for space and $d$ for velocity). The idea of moment method is to consider a finite number of equations for moments, instead of the Boltzmann equation, that depend only on $d+1$ variables (1 for time and $d$ for space); and the hope is that a sufficient number of moment equations would recover the solution from the Boltzmann equation (to a certain extent). The details of the Grad moment method are skipped here for the sake of brevity but they—for monatomic gases—can be found in Grad (1949) and in standard textbooks, e.g. Struchtrup (2005); Kremer (2010), and—for granular gases of three-dimensional hard spheres—in Gupta et al. (2018).

Inclusion of the governing equations for the stress ($\sigma_{ij}$) and heat flux ($q_{i}$) along with the system of mass, momentum and energy balance equations (8)–(10) leads to the well-known system of the 13-moment equations in three dimensions. In this paper, some Grad moment systems consisting of higher-order moments will also be derived and investigated. To this end, it is convenient to introduce a general symmetric-traceless moment

[TABLE]

and its associated collisional production term (or collisional moment)

[TABLE]

where the angle brackets around the indices again denote the symmetric and traceless part of the corresponding quantity; see appendix A for its definition. From definitions (13) and (14), it is straightforward to verify that $u^{0}=m\,n=\rho$, $u_{i}^{0}=0$, $u^{1}=d\,n\,T=d\,\rho\,\theta$, $u_{ij}^{0}=\sigma_{ij}$, $u_{i}^{1}=2\,q_{i}$, $\mathcal{P}^{0}=\mathcal{P}_{i}^{0}=0$ and $\mathcal{P}^{1}=-d\,n\,T\,\zeta$. Here $\rho=m\,n$ is the mass density and $\theta=T/m$.

3.1 Counting moments in $d$ dimensions

Before deriving the various moment systems, it is worthwhile to know how many moments a Grad moment system contains in a general dimension $d$. As it is more convenient to work with symmetric-traceless moments, the number of moments in a Grad moment system in a general dimension can be determined by knowing the number of independent components in a symmetric $r$-rank tensor and the number of traces in this tensor. Indeed, the number of independent components of a fully symmetric $r$-rank tensor in $d$ dimensions is $\binom{d+r-1}{r}=\frac{(d+r-1)!}{r!\,(d-1)!}$, and the number of traces in this tensor is [math] for $r\in\{0,1\}$ while $\binom{d+r-3}{r-2}=\frac{(d+r-3)!}{(r-2)!\,(d-1)!}$ for $r\in\mathbb{N}\setminus\{1\}$. Consequently, the number of independent components of a fully symmetric-traceless $r$-rank tensor ($r\in\mathbb{N}\setminus\{1\}$) in $d$ dimensions is

[TABLE]

Notably, any symmetric-traceless $r$-rank tensor ($r\in\mathbb{N}$) in two dimensions has only two independent components while any symmetric-traceless $r$-rank tensor ($r\in\mathbb{N}$) in three dimensions has $2r+1$ independent components. The counting of number of moments in some of the Grad moment systems considered in this paper is illustrated in table 1. Notwithstanding, any Grad moment system considered in this paper henceforth will be referred by its number of moments in three dimensions since Grad moment systems with the number of moments in three dimensions are more familiar to us (see, e.g., Jenkins & Richman, 1985a; Levermore, 1996; Struchtrup, 2005; Kremer & Marques Jr., 2011). For instance, the Grad $\frac{d^{2}+5d+2}{2!}$-, $\frac{(d+1)(d^{2}+8d+6)}{3!}$- or $\frac{(d+3)(d^{2}+6d+2)}{3!}$-moment systems will simply be referred to as the Grad $13$-, $26$- or $29$-moment systems, respectively, which we are more acquainted with.

3.2 The system of the 29-moment equations

The system of the 29-moment equations includes the governing equations for the third rank tensor, for one- and full-traces of the fourth rank tensor and for full-trace of the fifth rank tensor along with the governing equations for the well-known 13 moments. In other words, the system of the 29-moment equations consists of the governing equations for the moments $n$, $v_{i}$, $T$, $\sigma_{ij}$, $q_{i}$, $u_{ijk}^{0}$, $u^{2}$, $u_{ij}^{1}$, $u_{i}^{2}$, and is obtained by multiplying the Boltzmann equation (1) with $1$, $c_{i}$, $\frac{1}{d}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^{4}$, $m\,C^{2}C_{\langle i}C_{j\rangle}$ and $m\,C^{4}C_{i}$, and integrating each of the resulting equations over the velocity space successively. The detailed derivation of the 29-moment equations is provided as supplementary material. Here, they are presented directly. The system of the 29-moment equations consists of the mass, momentum and energy balance equations (8)–(10) and other higher-order moment equations, which on using the abbreviations

[TABLE]

read

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The abbreviations (15) are introduced in such a way that $m_{ijk}$, $\Delta$, $R_{ij}$ and $\varphi_{i}$ vanish if computed with the well-known G13 distribution function

[TABLE]

where

[TABLE]

is the Maxwellian distribution function (Garzó, 2013). In general, the computation of the collisional production terms $\mathcal{P}_{i_{1}i_{2}\dots i_{r}}^{a}$ requires the knowledge of the distribution function and is not easy for particles interacting with a general interaction potential. Nevertheless, for IMM (considered in this work), the collisional production terms can be evaluated easily—indeed, without the knowledge of the explicit form of the distribution function. A strategy for computing them for IMM in an automated way using the computer algebra software Mathematica® is demonstrated in appendix B. Using this strategy, the production terms associated with the G29 equations for $d$-dimensional IMM have been computed. They turn out to be

[TABLE]

where the coefficients $\zeta_{0}^{\ast}$, $\nu_{\sigma}^{\ast}$, $\nu_{q}^{\ast}$, $\nu_{m}^{\ast}$, $\nu_{R}^{\ast}$, $\nu_{\varphi}^{\ast}$, $\alpha_{0}$, $\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$, $\varsigma_{0}$, $\varsigma_{1}$, $\varsigma_{2}$ and $\varsigma_{3}$ depend only on the dimension $d$ and coefficient of restitution $e$, and are relegated to appendix C for better readability. Collisional production terms (24)–(29) for IMM agree with those obtained in Garzó & Santos (2007), wherein they have been computed till fourth order. Moreover, the coefficients $\zeta_{0}^{\ast}$, $\nu_{\sigma}^{\ast}$, $\nu_{q}^{\ast}$, $\nu_{R}^{\ast}$, $\nu_{\varphi}^{\ast}$, and $\alpha_{1}$ relate to the collisional rate $\nu_{2r|s}$—associated with the Ikenberry polynomial $Y_{2r|i_{1}i_{2}\dots i_{s}}(\bm{C})$—given in Santos & Garzó (2012) for $s\in\{0,1,2\}$ via $\zeta_{0}^{\ast}\,\nu=\nu_{2|0}$, $\nu_{\sigma}^{\ast}\,\nu=\nu_{0|2}$, $\nu_{q}^{\ast}\,\nu=\nu_{2|1}$, $\nu_{R}^{\ast}\,\nu=\nu_{2|2}$, $\nu_{\varphi}^{\ast}\,\nu=\nu_{4|1}$ and $\alpha_{1}\,\nu=d(d+2)\nu_{4|0}$. I could not find the full expression for the collisional production term (30) for granular gases in the existing literature. Nonetheless, for monatomic gases (i.e. for $d=3$ and $e=1$), it can be found, for instance, in Gu & Emerson (2009)—although not explicitly. The source code for computing the above collisional production terms is provided as supplementary material with the present paper. The collisional production terms associated with the G26 equations for three-dimensional IHS can be found in Gupta & Torrilhon (2012); Gupta et al. (2018).

The relation $\mathcal{P}^{1}=-d\,n\,T\,\zeta$ on exploiting (24) gives the cooling rate for IMM:

[TABLE]

where $\zeta_{0}^{\ast}=(d+2)(1-e^{2})/(4d)$ (see (133)). The cooling rate (31) is the same as that obtained in Santos (2003); Garzó & Santos (2011); Khalil et al. (2014), and vanishes identically for monatomic gases (i.e. for $e=1$), guaranteeing the conservation of energy for them. It is important to note from (31) that the cooling rate for IMM neither depends on the gradients of any field nor on any higher-order moment (in contrast to the cooling rate for IHS that also depends on the scalar fourth moment $\Delta$; see Gupta et al. (2018)).

3.3 Grad 29-moment closure

The system of the 29-moment equations for IMM (eqs. (8)–(10) and (16)–(3.2) along with collisional production terms (24)–(30)) is still not closed as it possesses the additional unknown moments $u_{ijkl}^{0}$, $u_{ijk}^{1}$, $u_{ij}^{2}$, $u^{3}$. The system is closed with the Grad distribution function based on the considered 29 moments, which is referred to as the G29 distribution function. The ($d$-dimensional) G29 distribution function $f_{|\mathrm{G29}}$ reads

[TABLE]

The details of computing the G29 distribution function (3.3) can be found in appendix D. Insertion of the G29 distribution function (3.3) into the definitions of unknown moments $u_{ijkl}^{0}$, $u_{ijk}^{1}$, $u_{ij}^{2}$ and $u^{3}$ expresses them in terms of the considered 29 moments:

[TABLE]

where the subscript “$|\mathrm{G29}$” denotes that these moments are computed with the G29 distribution function (3.3).

3.4 The G29 system for IMM

Equations (8)–(10) and (16)–(3.2) closed with (33) and (24)–(30) form the system of the G29 equations for $d$-dimensional IMM. Combining all of them, the system of the G29 equations for $d$-dimensional IMM reads

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

and the other coefficients $\zeta_{0}^{\ast}$, $\nu_{\sigma}^{\ast}$, $\nu_{q}^{\ast}$, $\nu_{m}^{\ast}$, $\nu_{R}^{\ast}$, $\nu_{\varphi}^{\ast}$, $\varsigma_{0}$, $\varsigma_{1}$, $\varsigma_{2}$ and $\varsigma_{3}$ appearing on the right-hand sides of the G29 equations (34)–(3.4) depend only on the dimension $d$ and coefficient of restitution $e$ (see appendix C for their expressions). For $d=3$ and $e=1$, these coefficients become $\nu_{\sigma}^{\ast}=1$, $\nu_{q}^{\ast}=2/3$, $\nu_{m}^{\ast}=3/2$, $\nu_{\!\Delta}^{\ast}=2/3$, $\nu_{R}^{\ast}=7/6$, $\nu_{\varphi}^{\ast}=1$, $\nu_{R\sigma}^{\ast}=0$, $\nu_{\varphi q}^{\ast}=0$, $\varsigma_{0}=2/45$, $\varsigma_{1}=2/3$, $\varsigma_{2}=28/15$ and $\varsigma_{3}=2/3$, which are the same as the respective coefficients for monatomic gases of Maxwell molecules; see, e.g., Struchtrup (2005) and Gu & Emerson (2009). In particular, the vanishing coefficients $\nu_{R\sigma}^{\ast}=0$ and $\nu_{\varphi q}^{\ast}=0$ make the right-hand sides of the linearised G29 equations for monatomic gases of Maxwell molecules completely decoupled, which is not the case for granular gases. Furthermore, the underlined nonlinear terms in (3.4)–(3.4) will be discarded for simplicity while investigating the HCS of a granular gas in § 5.

3.5 Various Grad moment systems

The abbreviations (15) have been introduced in such a way that the smaller systems of the Grad moment equations can be obtained directly from the G29 system (34)–(3.4). The other Grad moment systems considered in this paper are as follows.

The G13 system: The system of the 13-moment equations contains the governing equations for variables $n$, $v_{i}$, $T$, $\sigma_{ij}$ and $q_{i}$, i.e. it consists of equations (34)–(3.4). However, equations (34)–(3.4) contain additional unknowns $m_{ijk}$, $\Delta$ and $R_{ij}$ that vanish on being computed with the G13 distribution function (22). Thus, the G13 system for $d$-dimensional IMM consists of equations (34)–(3.4) with $m_{ijk}=\Delta=R_{ij}=0$. 2. 2.

The G14 system: The system of the 14-moment equations contains the governing equations for variables $n$, $v_{i}$, $T$, $\sigma_{ij}$, $q_{i}$ and $\Delta$, i.e. it consists of equations (34)–(3.4) and (3.4). However, equations (34)–(3.4) and (3.4) contain additional unknowns $m_{ijk}$, $R_{ij}$ and $\varphi_{i}$ that also vanish on being computed with the G14 distribution function

[TABLE]

which can be obtained easily by following a similar procedure presented in appendix D. Thus, the G14 system for $d$-dimensional IMM consists of equations (34)–(3.4) and (3.4) with $m_{ijk}=R_{ij}=\varphi_{i}=0$. 3. 3.

The G26 system: The system of the 26-moment equations contains the governing equations for variables $n$, $v_{i}$, $T$, $\sigma_{ij}$, $q_{i}$, $m_{ijk}$, $\Delta$ and $R_{ij}$, i.e. it consists of equations (34)–(3.4) and (3.2)–(3.2) with the right-hand sides computed using the collisional production terms (24)–(30). However, equations (3.2)–(3.2) contain additional unknowns $u_{ijkl}^{0}$, $\varphi_{i}$ and $u_{ijk}^{1}$ that are computed with the G26 distribution function

[TABLE]

which can also be obtained easily by following a similar procedure presented in appendix D. With the G26 distribution function (3), $u_{ijkl}^{0}$ and $\varphi_{i}$ vanish and $u_{ijk}^{1}$ turns out to be $u_{ijk|\mathrm{G26}}^{1}=(d+6)\,\theta\,m_{ijk}$, which is exactly the same as the value of $u_{ijk}^{1}$ obtained with the G29 distribution function (3.3). Therefore, inserting the G26 closure (i.e. $u_{ijkl}^{0}=\varphi_{i}=0$ and $u_{ijk}^{1}=(d+6)\,\theta\,m_{ijk}$), equations (3.2)–(3.2) turn to (3.4)–(3.4) in which $\varphi_{i}=0$. Thus, the G26 system for $d$-dimensional IMM consists of equations (34)–(3.4) with $\varphi_{i}=0$.

It is worthwhile to note that the G13 and G26 theories belong to the category of ordered moment theories, which always include the neglected fluxes of a moment theory at the previous level (Torrilhon, 2015). Also, there are other moment theories, which consider complete (traces and traceless) moments of a given order; such moment theories are referred to as full moment theories (Torrilhon, 2015). The first few examples of full moment theories are the Grad 10-, 20- and 35-moment theories (in three dimensions). In this sense, the G14 and G29 theories considered in the present work neither belong to the category of ordered moment theories nor to that of full moment theories.

4 Transport coefficients in the NSF laws

Recall that system (8)–(10) of the mass, momentum and energy balance equations was not closed due to the presence of additional unknowns: the stress $\sigma_{ij}$, heat flux $q_{i}$ and cooling rate $\zeta$. One of the major goals of kinetic theory is to furnish a closure for the system of the mass, momentum and energy balance equations in the form of constitutive relations. Traditionally, these constitutive relations are derived by performing the CE expansion on the Boltzmann equation. An alternative, but relatively much easier, way to determine the constitutive relations is by means of a CE-like expansion—in powers of a small parameter (usually, the Knudsen number)—performed on the Grad moment system. For monatomic gases of Maxwell molecules, it can be shown via the order of magnitude approach that a CE-like expansion on the G13 equations yields the Euler, NSF and Burnett constitutive relations at the zeroth, first and second orders of expansion, respectively (Struchtrup, 2005). Thus, for monatomic gases of Maxwell molecules, the G13 equations already contain the Burnett equations. Such a CE-like expansion procedure of Struchtrup (2005) on the Grad moment equations for IMM is much more involved due to non-conservation of energy, and—at its present understanding—does not yield the correct transport coefficients appearing in the constitutive relations. I still believe that a formal CE-like expansion procedure based on the order of magnitude of moments, which would yield the correct transport coefficients for granular gases (in particular, for IMM) can be devised; although it will be a topic for future research. Here, I follow the approach of Garzó (2013) to determine the transport coefficients in the NSF laws for a dilute granular gas of IMM through the Grad moment equations developed above.

The cooling rate in the energy balance equation (10) for IMM is given by (31) while the constitutive relations for the stress and heat flux for closing the system of the mass, momentum and energy balance equations (8)–(10)—to the linear approximation in spatial gradients—read (Jenkins & Richman, 1985a, b; Garzó & Dufty, 1999; Garzó, 2013)

[TABLE]

where $\eta$, $\kappa$ and $\lambda$ are the transport coefficients. The coefficient $\eta$ is referred to as the shear viscosity and $\kappa$ as the thermal conductivity; the coefficient $\lambda$ is a Dufour-like coefficient (Alam et al., 2009; Kremer et al., 2014; Shukla et al., 2019) that vanishes identically for monatomic gases. Equations (46a) and (46b) are the Navier–Stokes’ law and the Fourier’s law for granular gases, respectively. Equations (46a) and (46b) together are referred to as the NSF laws for granular gases.

4.1 Zeroth-order contributions in spatial gradients

To zeroth order in the spatial gradients, the mass, momentum and energy balance equations (34)–(36) reduce to

[TABLE]

while the balance equations for the other higher moments (eqs. (37)–(3.4)) reduce to

[TABLE]

Equations in (50) readily imply that

[TABLE]

where

[TABLE]

is the same as the value of the fourth cumulant for IMM reported in previous studies (Santos, 2003; Khalil et al., 2014). Thus, to zeroth order in spatial gradients, $\sigma_{ij}$, $q_{i}$, $m_{ijk}$, $R_{ij}$ and $\varphi_{i}$ are zero while $\Delta=a_{2}$.

4.2 First-order contributions in spatial gradients

Now, the terms having first-order spatial derivatives are also retained in the moment equations. To first order in spatial gradients, moment equations (34)–(3.4), read

[TABLE]

Notice that, unlike IHS (see Gupta et al. (2018)), here none of the balance equations (3.4)–(3.4) is required for determining the transport coefficients for IMM up to first-order accuracy in spatial gradients, since the stress and heat flux balance equations (eqs. (56) and (57)) have no coupling with the higher moments. The balance equations (3.4)–(3.4) will only be needed for computing the transport coefficients beyond the first-order accuracy in spatial gradients, which is not the focus of the present work.

The time derivatives of the stress and heat flux in (56) and (57) are computed using dimensional analysis and using the zeroth-order accurate mass, momentum and energy balance equations (47). They turn out to be (Garzó, 2013; Gupta et al., 2018)

[TABLE]

Now, in the first-order accurate stress and heat flux balance equations ((56) and (57)), one replaces $\sigma_{ij}$ and $q_{i}$ using (46) and their time derivatives using (58). Subsequent comparison of the coefficients of each hydrodynamic gradient in both the resulting equations leads to the transport coefficients in the NSF laws (46):

[TABLE]

where

[TABLE]

are the shear viscosity and thermal conductivity, respectively, in the elastic limit; and $\eta^{\ast}$, $\kappa^{\ast}$ and $\lambda^{\ast}$ are the reduced shear viscosity, reduced thermal conductivity and reduced Dufour-like coefficient, respectively. These reduced transport coefficients are given by

[TABLE]

Expressions (61) for the reduced transport coefficients agree with those obtained at first order of the CE expansion for IMM, e.g. in Santos (2003); Khalil et al. (2014); Garzó & Santos (2011). Indeed, the structure of these transport coefficients is very similar to those for IHS (Brey et al., 1998a; Garzó, 2013) except for the fact that $a_{2}$, $\zeta_{0}^{\ast}$, $\nu_{\sigma}^{\ast}$ and $\nu_{q}^{\ast}$ for IMM and IHS are different. Despite the structural similarity, the transport coefficients $\kappa^{\ast}$ and $\lambda^{\ast}$ for IMM ((61b) and (61c)) diverge at a certain value of the coefficient of restitution and do not yield meaningful values below it—in contrast to the transport coefficients for IHS which are meaningful for all values of the coefficient of restitution and are in reasonably good agreement with the simulations. This issue pertaining to IMM can readily be appreciated by inspecting the explicit dependence of the reduced transport coefficients on the coefficient of restitution and dimension as follows. Inserting $a_{2}$, $\zeta_{0}^{\ast}$, $\nu_{\sigma}^{\ast}$ and $\nu_{q}^{\ast}$ from (52), (133), (134) and (135) in the reduced transport coefficients (61), they are expressed as a function of the coefficient of restitution and dimension (Garzó & Santos, 2011):

[TABLE]

Clearly, both $\kappa^{\ast}$ and $\lambda^{\ast}$ have singularities at $e=(4-d)/(3d)$, for which the denominators of both of them vanish. In particular, the denominators of both $\kappa^{\ast}$ and $\lambda^{\ast}$ vanish at $e=1/3$ for $d=2$ and at $e=1/9$ for $d=3$. Moreover, below the singularities, i.e. for $e<(4-d)/(3d)$, both $\kappa^{\ast}$ and $\lambda^{\ast}$ are negative, which is unphysical. The existence of these singularities is apparently attributed to the breakdown of hydrodynamics in granular gases of IMM for $e\leq(4-d)/(3d)$ due to the lack of time scale separation between the kinetic and hydrodynamic parts of the distribution function (Brey et al., 2010). Owing to these singularities, it is customary to write the heat flux as a linear combination of the gradients of $T$ and $n\sqrt{T}$ instead of its usual representation (46b) (see, e.g., Garzó et al., 2007; Garzó & Santos, 2011). The heat flux in the new form reads

[TABLE]

is referred to as the modified thermal conductivity (Garzó et al., 2007). The reduced modified thermal conductivity $\kappa^{\prime\ast}=\kappa^{\prime}/\kappa_{0}$—using (61)–(63)—is given by

[TABLE]

The reduced modified thermal conductivity $\kappa^{\prime\ast}$ does not possess the above singularity and hence is finite for all $0\leq e\leq 1$ and for $d=2,3$.

4.3 Comparison with existing theories and computer simulations

I have not found any simulation data on the transport coefficients for IMM. Therefore, in this subsection, I compare the reduced transport coefficients for IMM obtained above with those for IHD ($d=2$) and IHS ($d=3$) obtained through various theoretical and simulation methods.

The reduced transport coefficients $\eta^{\ast}$, $\kappa^{\ast}$, $\lambda^{\ast}$ and $\kappa^{\prime\ast}$ for a dilute granular gas are plotted over the coefficient of restitution $e$ in figures 1–4, respectively. The left and right panels in each figure exhibit the results for $d=2$ and $d=3$, respectively. The thick solid (red) line in each figure denotes the result for IMM obtained from (61) or (62) and (64), which have been obtained in this paper through the Grad moment equations. Recall that the reduced transport coefficients for IMM obtained through the moment method above are exactly the same as those obtained at first order of the CE expansion (see Santos, 2003; Garzó & Santos, 2011). Therefore, the thick solid (red) line in each figure also represents the results for IMM from the first-order CE expansion. The remaining lines and symbols in figures 1–4 depict the results for IHD (in case of $d=2$) or for IHS (in case of $d=3$). The thin solid (green) and dash-dotted (magenta) lines are the plots for the reduced transport coefficients from Garzó et al. (2007) obtained at the first Sonine and modified first Sonine approximations, respectively, in the CE expansion. The dashed (black) lines depict the reduced transport coefficients obtained through the G14 distribution function for $d$-dimensional IHS in Garzó (2013).

The right panel of figure 1 also displays a solid cyan line, which is not present in the other figures. This solid cyan line is the result for the reduced shear viscosity obtained with the G26 equations for IHS very recently by Gupta et al. (2018). Indeed, the other transport coefficients from the G14 or G26 equations remain the same; consequently, the solid cyan line in the right panels of each of figures 2–4 coincides with the dashed black line, and hence has not been shown separately. The squares are the results from the theoretical expressions derived via the so-called computer-aided method devised by Noskowicz et al. (2007) and the circles are the numerical solutions of the Boltzmann equation obtained via the direct simulation Monte Carlo (DSMC) method in Brey & Ruiz-Montero (2004); Montanero et al. (2005); Brey et al. (2005); Montanero et al. (2007). The paper by Garzó et al. (2007) summarises and presents the DSMC results in the aforementioned references that are computed with two approaches: (i) through Green–Kubo relations in Brey & Ruiz-Montero (2004) (for $d=2$) and Brey et al. (2005) (for $d=3$), and (ii) by implementing an external force in Montanero et al. (2005) and Montanero et al. (2007). The external force in the latter compensates for collisional cooling and yields somewhat better results. Therefore, the DSMC results in figures 1–4 are shown with the latter for whichever coefficient they are available else they are shown with the former—figure 1 and the right panel of figure 4 display the DSMC results with the latter while figures 2 and 3 show the DSMC results with the former. The DSMC results for the reduced shear viscosity from the latter approach were obtained by Montanero et al. (2005) for $e=0.6,0.7,0.8,0.9,1$ in the case of $d=3$ while that for $e=0.2,0.3,0.4,0.5$ in the case of $d=3$ and that for all $e$ in the case of $d=2$ were obtained by Garzó et al. (2007). The DSMC data for $\kappa^{\prime\ast}$ in two dimensions (left panel of figure 4) are apparently unavailable.

It is clear from figures 1–4 that the IMM model overpredicts all the transport coefficients significantly in comparison to the IHS model, despite the transport coefficients for IMM computed from the Grad moment method in the present work being exactly the same as those obtained through CE expansion for IMM, for example, in Santos (2003); Garzó & Santos (2011). The discrepancies between the results for IMM and IHS are apparently linked to the choice of the effective collision frequency $\mathring{\nu}$ (see (3$a$)) in the IMM model (Santos, 2003), which is chosen in such a way that the cooling rates from the Boltzmann equation for IHS and IMM remain exactly the same. Furthermore, the reduced transport coefficients $\kappa^{\ast}$ and $\lambda^{\ast}$ for IMM (shown by the thick solid red lines in figures 2 and 3) diverge at $e=1/3$ for $d=2$ (left panels of figures 2 and 3) and at $e=1/9$ for $d=3$ (right panels of figures 2 and 3), and remain unphysical below these values of the coefficient of restitution. On the other hand, the reduced modified thermal conductivity $\kappa^{\prime\ast}$ for IMM remains positive for all values of the coefficient of restitution in both dimensions (see figure 4). Nevertheless, $\kappa^{\prime\ast}$ for IMM is also much higher than that for IHS. In the case of $d=2$ (left panel of figure 4), $\kappa^{\prime\ast}$ for IHS from any theory first decreases then increases on increasing the coefficient of restitution (although the profiles of $\kappa^{\prime\ast}$ from the modified first Sonine approximation and first Sonine approximation/G14 theory differ significantly) whereas that for IMM decreases monotonically on increasing the coefficient of restitution. However, as the DSMC data are not available in this case, it is difficult to discern which theory for IHS yields better results in this case.

Among fully theoretical methods, the modified version of the first Sonine approximation (dash-dotted magenta lines) proposed by Garzó et al. (2007) for IHS seems to be the best model, which captures all the transport coefficient very well, although the G26 model of Gupta et al. (2018) was able to capture the coefficient of the reduced shear viscosity (but not the other transport coefficients) better than the modified first Sonine approximation.

5 The HCS of a freely cooling granular gas of IMM

The state of a force-free granular gas when its granular temperature decays constantly while its spatial homogeneity is maintained is termed as the HCS (Brilliantov & Pöschel, 2004). For studying the HCS, one considers a force-free (i.e. $\bm{F}=\bm{0}$) granular gas having an initial number density as $n(0,\bm{x})=n_{0}$ and initial granular temperature $T(0,\bm{x})=T_{0}$ at time $t=0$ when the gas is left to cool down freely due to inelastic collisions while maintaining the spatial homogeneity (i.e. $\partial(\cdot)/\partial x_{i}=0$).

In this section, I investigate the HCS of a $d$-dimensional granular gas of IMM with the Grad moment equations (34)–(3.4) presented above. The nonlinear (underlined) contributions on the right-hand sides of the Grad moment equations (34)–(3.4) are discarded in this section for simplicity. This means that our focus is on the early evolution stage of homogeneously cooling granular gas. Hence, the possibility of increase in the granular temperature in a cooling granular gas (reported recently for granular gases of aggregating particles by Brilliantov et al. (2018)) is disregarded, which possibly occurs at large times.

It is convenient to study the HCS with dimensionless variables obtained by introducing the following scaling:

[TABLE]

where $\theta_{0}=T_{0}/m$ and $\nu_{0}=\nu(t=0)=4\,\Omega_{d}\,n_{0}\,\mathbbm{d}^{d-1}\sqrt{T_{0}/m}/[\sqrt{\pi}(d+2)]$. With scaling (67), the G29 equations (34)–(3.4)—without the underlined terms—in the HCS (i.e. with $\partial(\cdot)/\partial x_{i}=0$, $\bm{F}=\bm{0}$) reduce to

[TABLE]

5.1 Haff’s law

Equations (68) and (69) with the initial conditions of the HCS imply $n_{\ast}(t_{\ast})=1$ and $v_{i}^{\ast}(t_{\ast})=0$. Therefore, equation (70) using the initial conditions of the HCS yields Haff’s law (Haff, 1983) for the evolution of the granular temperature:

[TABLE]

where

[TABLE]

Here, $\tau_{0}$ is the time scale in Haff’s law for IMM and $\tau_{\ast}$ is the corresponding dimensionless time scale. Haff’s law (77) with time scale (78) is exactly the same as that obtained in Garzó & Santos (2011) for IMM. It is worthwhile to note that, unlike the energy balance equation in the case of IHS that also contains the scalar fourth moment $\Delta$ (see, e.g., Kremer & Marques Jr., 2011; Gupta et al., 2018), equation (70) does not contain any other moment except $n_{\ast}$ and $T_{\ast}$. Consequently, Haff’s law for IMM does not depend on higher moments; or in other words, Haff’s law remains unchanged for IMM, no matter how large a moment system it is determined from.

Note that the dimensionless time scale $\tau_{\ast}$ in Haff’s law (77) for $d$-dimensional IHS obtained at first approximation of the Sonine expansion is given by (van Noije & Ernst, 1998)

[TABLE]

with an excellent estimate for the fourth cumulant (van Noije & Ernst, 1998)

[TABLE]

Figure 5 illustrates the decay of the dimensionless granular temperature in the HCS via Haff’s law (77$a$) for the coefficients of restitution $e=0.75$ (depicted by solid lines and squares) and $e=0.95$ (depicted by dotted lines and circles). The lines denote Haff’s law for IMM, where the associated (dimensionless) time scale $\tau_{\ast}$ for the decay is given by (78$a$), while the symbols denote that for IHS at first approximation of the Sonine expansion, where the associated (dimensionless) time scale $\tau_{\ast}$ is given by (79). Although the time scale $\tau_{\ast}$ for IMM does not contain the fourth cumulant $a_{2}$ while that for IHS does contain it, Haff’s law from IMM (lines) is still in very good agreement with that from IHS (symbols) in both two and three dimensions. The granular temperature relaxes faster with increasing inelasticity due to the fact that more inelastic particles dissipate more energy during a collision in comparison to the less inelastic ones.

5.2 Relaxation of other moments in the HCS

For monatomic gases, it can be shown through the order of magnitude method devised by Struchtrup (2004) that all the nonequilibrium moments ($\sigma_{ij}$, $q_{i}$, $m_{ijk}$, $\Delta$, $R_{ij}$ and beyond) are at least of first order in spatial gradients (see Struchtrup, 2004, 2005). In contrast, the order of magnitude method in the case of granular gases is not straightforward due to non-conservation of energy and, to the best of my knowledge, has never been attempted so far. Notwithstanding, I would expect that all the higher vectorial and tensorial moments ($\sigma_{ij}$, $q_{i}$, $m_{ijk}$, $R_{ij}$, $\varphi_{i}$ and beyond) for granular gases are also at least of first order in spatial gradients; this conjecture is well known for $\sigma_{ij}$ and $q_{i}$ in the case of granular gases as well. This means that all the vectorial and tensorial moments remain zero in the HCS because spatial gradients are zero in this state. Nonetheless, it is interesting to analyse how these higher-order moments relax in the HCS if started with non-vanishing initial conditions.

Equations (71)–(76) are coupled with (68) and (70), but can be solved analytically. In order to compare the decay rates of the moments, the initial conditions for all the variables in (71)–(76) are taken as unity. With these initial conditions, equations (71)–(76) yield the following solution for the other variables:

[TABLE]

where $\varkappa_{R}=\nu_{R\sigma}^{\ast}\tau_{\ast}/[(\nu_{R}^{\ast}-\nu_{\sigma}^{\ast})\tau_{\ast}-2]$ and $\varkappa_{\varphi}=\nu_{\varphi q}^{\ast}\tau_{\ast}/[(\nu_{\varphi}^{\ast}-\nu_{q}^{\ast})\tau_{\ast}-2]$. It is important to note that for dilute monatomic gases (i.e. in the case of $d=3$ and $e=1$) of Maxwell molecules, equations (68)–(76) with the same initial conditions yield the solution

[TABLE]

From solution (82), it is clear that, for dilute monatomic gases of Maxwell molecules, the third-order moment $m_{ijk}^{\ast}$ relaxes faster than all other higher-order moments considered in the present work; $R_{ij}^{\ast}$ relaxes slower than $m_{ijk}^{\ast}$ but faster than the remaining moments ($\sigma_{ij}^{\ast}$, $q_{i}^{\ast}$, $\Delta$ and $\varphi_{i}^{\ast}$); $\sigma_{ij}^{\ast}$ and $\varphi_{i}^{\ast}$ relax with equal relaxation rates but faster than $q_{i}^{\ast}$ and $\Delta$, which also relax with equal relaxation rates.

Figure 6 illustrates the relaxation of the other (dimensionless) moments—$\sigma_{ij}^{\ast}$, $q_{i}^{\ast}$, $m_{ijk}^{\ast}$, $\Delta$, $R_{ij}^{\ast}$ and $\varphi_{i}^{\ast}$ over the (dimensionless) time $t_{\ast}$ (via analytical solutions (81)) in two and three dimensions for granular gases (i.e. for $e\neq 1$): $e=0.95$ (main panels) and $e=0.25$ (insets). It turns out that all these moments relax with time much faster than the granular temperature. It is interesting to note that all the vectorial and tensorial moments relax to zero whereas the scalar moment $\Delta$ relaxes to a nonzero value $a_{2}$ for all $e\neq 1$. These can also be verified from analytical solutions (81) in the limit $t_{\ast}\to\infty$, since all the exponents in (81) are negative for all $e\neq 1$ (note that $\tau_{\ast}$, $\nu_{\sigma}^{\ast}$, $\nu_{q}^{\ast}$, $\nu_{m}^{\ast}$, $\nu_{\Delta}^{\ast}$, $\nu_{R}^{\ast}$, $\nu_{R\sigma}^{\ast}$, $\nu_{\varphi}^{\ast}$, $\nu_{\varphi q}^{\ast}$, $\alpha_{0}$, $\alpha_{1}$, $\alpha_{2}$ and $\alpha_{3}$ are positive for $e\neq 1$). For large values of the coefficient of restitution (main panels), all these moments decay monotonically and their decay rates are, apparently, proportional to their tensorial orders, i.e. the third-order moment ($m_{ijk}^{\ast}$) decays faster than the second-order moments ($\sigma_{ij}^{\ast}$ and $R_{ij}^{\ast}$), which decay faster than the vectorial moment ($q_{i}^{\ast}$ and $\varphi_{i}^{\ast}$), which decay faster than the scalar moment ($\Delta$). However, for sufficiently small values of the coefficient of restitution (insets), the higher-order moments ($R_{ij}^{\ast}$ and $\varphi_{i}^{\ast}$) do not decay monotonically. This is attributed to the fact that higher-order (sixth-order and beyond) moments for IMM are prone to diverge for sufficiently small values of the coefficient of restitution in the HCS (Santos & Garzó, 2012) and it is manifested already through the non-monotonic relaxation of $R_{ij}^{\ast}$ and $\varphi_{i}^{\ast}$ for small coefficients of restitution (insets), although they themselves do not diverge.

6 Linear stability analysis of the HCS

In this section, the temporal stability of the HCS of a freely cooling granular gas of IMM due to small perturbations will be analysed through the G29 and other Grad moment theories described in § 3.5 with $F_{i}=0$. The amplitudes of these perturbations are assumed to be sufficiently small in order to ensure the validity of the linear analysis.

For the linear stability analysis, all the field variables in (34)–(3.4) are decomposed into their reference state values (i.e. their solutions in the HCS) plus perturbations from their respective reference state values. In other words, the field variables in the G29 system (34)–(3.4) are written as

[TABLE]

where $T_{H}(t)$ is the granular temperature in the HCS and $v_{H}(t)=\sqrt{T_{H}(t)/m}$ is a reference speed proportional to the adiabatic sound speed in the HCS; and the quantities with tilde denote the dimensionless perturbations in the field variables from their respective solutions in the HCS.

Inserting expressions (83) for the field variables into (34)–(3.4) and discarding all the nonlinear terms of the perturbations, one obtains the system of linear PDEs in (dimensionless) perturbed field variables with time-dependent coefficients. This system of PDEs is transformed to a new system of PDEs having time-independent coefficients by introducing a length scale

[TABLE]

that is employed to make the space variables dimensionless (i.e. $\tilde{x}_{i}=x_{i}/\ell$, where tilde denotes the dimensionless space variable), and a dimensionless time

[TABLE]

that measures time as the number of effective collisions per particle (Brey et al., 1998b; Brilliantov & Pöschel, 2004; Garzó & Santos, 2007). The resulting system of PDEs having time-independent coefficients reads

[TABLE]

where

[TABLE]

System (86)–(94), admits a normal mode solution of the form

[TABLE]

where $\tilde{\bm{\Psi}}=(\tilde{n},\tilde{v}_{i},\tilde{T},\tilde{\sigma}_{ij},\tilde{q}_{i},\tilde{m}_{ijk},\tilde{\Delta},\tilde{R}_{ij},\tilde{\varphi}_{i})^{\mathsf{T}}$ is the vector containing all the dimensionless perturbations and $\check{\bm{\Psi}}=(\check{n},\check{v}_{i},\check{T},\check{\sigma}_{ij},\check{q}_{i},\check{m}_{ijk},\check{\Delta},\check{R}_{ij},\check{\varphi}_{i})^{\mathsf{T}}$ the vector containing their corresponding complex amplitudes. Furthermore, in the normal mode solution (98), $\mathbbm{i}$ is the imaginary unit, $\bm{k}$ the dimensionless wavevector of the disturbance and $\omega$ the dimensionless frequency of the associated wave. For temporal stability analysis, the wavevector $\bm{k}$ is assumed to be real and the frequency $\omega$ is assumed to be complex. The real part of the frequency, $\mathrm{Re}(\omega)$, measures the phase velocity $\bm{v}_{\mathrm{ph}}$ of the wave via $\bm{v}_{\mathrm{ph}}=\mathrm{Re}(\omega)\bm{k}/k^{2}$, where $k=|\bm{k}|$ is the wavenumber, and the imaginary part of the frequency, $\mathrm{Im}(\omega)$, controls the growth/decay of the disturbance in time and is referred to as the growth rate. Form (98) of the normal mode solution deduces that the disturbance will grow (or decay) in time if the growth rate is positive (or negative). Consequently, for stability of the system, the growth rate must be non-positive, i.e. $\mathrm{Im}(\omega)\leq 0$.

Assuming that the wavevector of the disturbance is parallel to the $x$-axis, i.e. $\bm{k}=k\,\hat{\bm{x}}$, where $\hat{\bm{x}}$ is the unit vector in the $x$-direction, system (86)–(94), using relations in appendix E, can be decomposed into two independent eigenvalue problems, namely the longitudinal problem and the transverse problem, for the amplitude of the disturbance. It is worthwhile to note that in two dimensions, there is only one transverse direction along the $y$-axis while in three dimensions, there are two transverse directions along the $y$- and $z$-axes. Consequently, there is one transverse problem in two dimensions and two transverse problems in three dimensions. Nevertheless, the coefficient matrices associated with the both transverse problems in three dimensions are essentially the same; therefore it is sufficient to analyse only one transverse problem (let us say, that associated with the $y$-direction) in three dimensions. Thus, the longitudinal and transverse problems read

[TABLE]

respectively, where the matrices $\mathscr{L}\equiv\mathscr{L}(k,\omega,d,e)$ and $\mathscr{T}\equiv\mathscr{T}(k,\omega,d,e)$ are presented in appendix F.

For nontrivial solutions of the longitudinal and transverse problems (99), the determinants of both matrices $\mathscr{L}$ and $\mathscr{T}$ must vanish, i.e. $\det{(\mathscr{L})}=0$ and $\det{(\mathscr{T})}=0$. These conditions are the dispersion relations for the longitudinal and transverse systems and can, respectively, be written as

[TABLE]

where the coefficients $\mathbbm{a}_{r}$ ($r=1,2,\dots,9$) and $\mathbbm{b}_{s}$ ($s=1,2,\dots,6$) are functions of the wavenumber $k$, the dimension $d$ and the coefficient of restitution $e$; although the explicit values of these coefficients are not given here for conciseness. The solutions of the longitudinal and transverse problems (99) for each root $\omega\equiv\omega(k)$ of the corresponding dispersion relations (100) are referred to as the eigenmodes for the longitudinal and transverse problems (99), respectively.

6.1 Eigenmodes from the longitudinal and transverse systems

The nine roots of dispersion relation (100$a$) lead to nine eigenmodes for the longitudinal system (99$a$) associated with the G29 equations while the six roots of dispersion relation (100$b$) yield six eigenmodes for the transverse system (99$b$) associated with the G29 equations. The real part of the frequency ($\mathrm{Re}(\omega)$) associated with each eigenmode and its growth rate ($\mathrm{Im}(\omega)$) from both the longitudinal and transverse systems are illustrated in figures 7–10 for $d=2$ (figures 7 and 9), $d=3$ (figures 8 and 10)—with figures 7 and 8 being for the longitudinal system (99$a$) and figures 9 and 10 for the transverse system (99$b$). The top and bottom rows in each figure depict the eigenmodes for the inelastic ($e=0.75$) and elastic ($e=1$) cases, respectively while the left and right columns in each figure delineate the real part of the frequency and growth rate, respectively. It is apparent from the left columns of the figures 7 and 8 that four pairs out of the nine eigenmodes from the longitudinal system have nonzero $\mathrm{Re}(\omega)$, i.e. the four pairs of associated eigenmodes are travelling whereas one eigenmode is purely imaginary, i.e. it has $\mathrm{Re}(\omega)=0$ for all wavenumbers and hence always remains stationary. Similarly, it is clear from the left columns of figures 9 and 10 that two pairs out of the six eigenmodes from the transverse system have nonzero $\mathrm{Re}(\omega)$, i.e. the two pairs of associated eigenmodes are travelling, whereas two eigenmodes are purely imaginary and hence remain stationary for all wavenumbers. A travelling eigenmode is commonly referred to as a sound mode and a stationary eigenmode as a heat mode (Brilliantov & Pöschel, 2004; Garzó, 2005).

6.1.1 Longitudinal systems (figures 7 and 8)

For $e=0.75$ (top rows of figures 7 and 8), the first pair of sound modes originates at $k\approx 0.00343$ in the case of $d=2$ (at $k\approx 0.00328$ in the case of $d=3$) (see the insets on the left columns of figures 7 and 8), followed by a second pair of sound modes commencing at $k\approx 0.0587$ in the case of $d=2$ (at $k\approx 0.0453$ in the case of $d=3$) travelling slower than the first pair, followed by a third pair of sound modes starting at $k\approx 0.1992$ in the case of $d=2$ (at $k\approx 0.1114$ in the case of $d=3$) travelling even slower than the second pair (in general), followed by a fourth pair of sound modes commencing at $k\approx 0.3791$ in the case of $d=2$ (at $k\approx 0.461$ in the case of $d=3$) travelling even slower than the third pair. It should be noted, however, that below these wavenumbers, the respective eigenmodes are heat modes since their real parts are zero. Each pair of the sound modes starts propagating in opposite directions at the aforesaid wavenumbers as the eigenvalues corresponding to each pair of the sound modes are a pair of complex conjugates. Beyond the aforesaid wavenumbers, the imaginary parts of each (respective) pair of sound modes coincide due to the same reason. This can be clearly seen in the top rows and right columns of figures 7 and 8, in which the imaginary parts of the first, second, third and fourth pairs of sound modes coincide beyond $k\approx 0.00343,0.0587,0.1992,0.3791$, respectively, in the case of $d=2$ (beyond $k\approx 0.00328,0.0453,0.1114,0.461$ respectively, in the case of $d=3$). From the top rows and right columns of figures 7 and 8, it is evident that all the eigenmodes except one heat mode from the fourth pair (for which $\mathrm{Im}(\omega)$ coincide beyond $k\approx 0.3791$ in the case of $d=2$ and beyond $k\approx 0.461$ in the case of $d=3$) are stable as $\mathrm{Im}(\omega)\leq 0$ for them. The unstable heat mode remains unstable for wavenumbers $k<k_{c}$ but becomes stable for wavenumbers $k\geq k_{c}$, where $k_{c}$ is called the critical wavenumber, the wavenumber at which $\mathrm{Im}(\omega)$ flips its sign. For $e=0.75$, $k_{c}\approx 0.179$ in the case of $d=2$ and $k_{c}\approx 0.18$ in the case of $d=3$. The general behaviour of the eigenmodes of the longitudinal system for moderate to large values of $e$ is similar to those for $e=0.75$ (as shown in the top rows of figures 7 and 8), although $k_{c}\to 0$ as $e\to 1$. Thus, for moderately to nearly elastic granular gases, there exists a critical wavenumber $k_{c}$, below which the unstable heat mode renders the longitudinal system unstable. On the other hand, for sufficiently small values of $e$, the growth rates of some of the eigenmodes of the longitudinal system remain positive even for large wavenumbers and, hence, there does not exist a critical wavenumber in this case. This means that the longitudinal system remains always unstable for all coefficients of restitution below a certain value, which is not true (see, e.g., Gupta et al., 2018). Let us refer to this value of the coefficient of restitution—below which a system remains always unstable—as the threshold coefficient of restitution $e_{\mathrm{th}}$. For the longitudinal system associated with the G29 equations, $e_{\mathrm{th}}\approx 0.56356$ in the case of $d=2$ and $e_{\mathrm{th}}\approx 0.40157$ in the case of $d=3$.

For $e=1$ (bottom rows of figures 7 and 8), two pairs (one faster and the other slower) of sound modes start propagating in opposite directions already at $k=0$, followed by a third pair of sound modes appearing at $k\approx 0.2104$ in the case of $d=2$ (at $k\approx 0.1937$ in the case of $d=3$) travelling slower than both the first and second pairs, followed by an even slower fourth pair of sound modes commencing at $k\approx 0.3383$ in the case of $d=2$ (at $k\approx 0.4714$ in the case of $d=3$). Accordingly, the imaginary parts of the frequencies for the two pairs of sound modes commencing at $k=0$ coincide for all wavenumbers, that for the third pair coincide beyond $k\approx 0.2104$ in the case of $d=2$ (beyond $k\approx 0.1937$ in the case of $d=3$) and that for the fourth pair coincide beyond $k\approx 0.3383$ in the case of $d=2$ (beyond $k\approx 0.4714$ in the case of $d=3$); see the bottom rows and right columns of the figures. From the bottom rows and right columns of figures 7 and 8, it can be perceived that $\mathrm{Im}(\omega)\leq 0$ for all eigenmodes in this case. This means that the longitudinal system remains stable for all wavenumbers in the elastic case ($e=1$).

6.1.2 Transverse system (figures 9 and 10)

For $e=0.75$ (top rows of figures 9 and 10), the first pair of sound modes starts propagating in opposite directions at $k\approx 0.0458$ in the case of $d=2$ (at $k\approx 0.07216$ in the case of $d=3$), followed by a second pair of sound modes commencing at $k\approx 0.1977$ in the case of $d=2$ (at $k\approx 0.1882$ in the case of $d=3$) travelling slower than the first pair. Beyond these wavenumbers, the imaginary parts of each (respective) pair of travelling sound modes merge due to the aforementioned reason, which is clearly reflected in the top rows and right columns of figures 9 and 10: the first pair coincides beyond $k\approx 0.0458$ in the case of $d=2$ (beyond $k\approx 0.07216$ in the case of $d=3$) and the second beyond $k\approx 0.1977$ in the case of $d=2$ (beyond $k\approx 0.1882$ in the case of $d=3$). The remaining two out of six eigenmodes remain stationary for all wavenumbers, i.e. these modes have no oscillations since their frequencies are purely imaginary. The non-oscillatory eigenmodes of the transverse system (99$b$) are referred to as the shear modes (Brilliantov & Pöschel, 2004; Garzó, 2005). Clearly, there are two shear modes in the transverse system but one of them is unstable for wavenumbers $k<k_{c}$ as its growth rate is positive (see the top rows and right columns of figures 9 and 10), where $k_{c}$ is the critical wavenumber for the transverse system. For $e=0.75$, $k_{c}\approx 0.341$ in the case of $d=2$ and $k_{c}\approx 0.296$ in the case of $d=3$. The general behaviour of the eigenmodes of the transverse system for moderate to large values of $e$ is also similar to that for $e=0.75$ (as shown in the top rows of figures 9 and 10) with $k_{c}\to 0$ as $e\to 1$. Thus, for moderately to nearly elastic granular gases, the unstable shear mode renders the transverse system unstable for wavenumbers $k<k_{c}$; nevertheless, the system becomes stable for all wavenumbers $k\geq k_{c}$. However, analogously to the longitudinal system, the growth rates of some of the eigenmodes of the transverse system also remain positive for large wavenumbers for all $e$ below a threshold coefficient of restitution $e_{\mathrm{th}}$, implying that there does not exist a $k_{c}$ for all $e<e_{\mathrm{th}}$ and that the transverse system remains always unstable for all $e<e_{\mathrm{th}}$, which is also not true (see, e.g., Gupta et al., 2018). For the transverse system associated with the G29 equations, $e_{\mathrm{th}}\approx 0.32349$ in the case of $d=2$ and $e_{\mathrm{th}}\approx 0.16867$ in the case of $d=3$.

For $e=1$ (bottom rows of figures 9 and 10), the first pair of travelling sound modes starts propagating in opposite directions at $k\approx 0.1056$ in the case of $d=2$ (at $k\approx 0.08398$ in the case of $d=3$), followed by a second pair of sound modes commencing at $k\approx 0.2296$ in the case of $d=2$ (at $k\approx 0.22396$ in the case of $d=3$) travelling slower than the first pair. Accordingly, the imaginary parts of frequencies for the first pair of travelling sound modes merge beyond $k\approx 0.1056$ in the case of $d=2$ (beyond $k\approx 0.08398$ in the case of $d=3$) and that for the second pair merge beyond $k\approx 0.2296$ in the case of $d=2$ (beyond $k\approx 0.22396$ in the case of $d=3$); see the bottom rows and right columns of the figures. The remaining two eigenmodes are stable shear modes in the elastic case since the real parts of their associated frequencies are zero and imaginary parts (growth rates) are non-positive for all wavenumbers. Consequently, the transverse system also remains stable for all wavenumbers in the elastic case ($e=1$).

6.2 Comparison among various Grad moment theories

As discussed in § 3.5, a lower-level Grad moment system can be obtained from the G29 equations by discarding the appropriate field variables. Accordingly, the longitudinal and transverse problems associated with the G13, G14 and G26 systems can be obtained by eliminating the appropriate variables and corresponding rows and columns of the matrices $\mathscr{L}$ and $\mathscr{T}$ in (99). For comparison purpose, I shall also include the results obtained from the linear stability analysis of the system of the NSF equations for IMM along with these Grad moment systems. The linear-dimensionless NSF equations in the perturbed field variables are (86)–(88) with

[TABLE]

where the reduced transport coefficients $\eta^{\ast}$, $\kappa^{\ast}$ and $\lambda^{\ast}$ for IMM are given by (61). From the linear-dimensionless NSF equations, it is straightforward to obtain the longitudinal and transverse problems associated with the system of the NSF equations by following a similar procedure as above. The longitudinal and transverse problems associated with the system of the NSF equations read

[TABLE]

respectively, where the matrix $\mathscr{L}_{\mathrm{NSF}}\equiv\mathscr{L}_{\mathrm{NSF}}(k,\omega,d,e)$ is also presented in appendix F. As far as the linear stability of the HCS is concerned, the NSF theory and all Grad moment theories—although not shown here explicitly for the NSF, G13, G14 and G26 equations—predict a similar behaviour for moderate to large values of the coefficient of restitution in the sense that one heat mode from the longitudinal system associated with each theory and one shear mode from the transverse system associated with each theory are unstable below some critical wavenumbers for granular gases while all the modes remain stable in the elastic case ($e=1$). The stability of a (longitudinal or transverse) system is regulated by its least stable eigenmode. Therefore, to analyse the stability of the longitudinal and transverse systems associated with different moment theories, the critical wavenumbers for the least stable mode (unstable shear mode for the longitudinal system and unstable heat mode for the transverse system) from each moment theory are plotted in the $(e,k_{c})$-plane in figure 11. The figure also includes the critical wavenumber profiles (shown by thin dashed black lines) for the least stable modes of the longitudinal and transverse systems associated with the NSF theory for IMM. The top and bottom rows of figure 11 display the critical wavenumbers for the longitudinal and transverse systems, respectively, while the left and right columns exhibit the results for $d=2$ and $d=3$, respectively. Since the critical wavenumber is that wavenumber where the growth rate ($\mathrm{Im}(\omega)$) changes its sign, the curves in figure 11 are essentially the zero contours—in the $(e,k_{c})$-plane—of the growth rate of the least stable mode in each system. Each curve for the critical wavenumber corresponds to $\mathrm{Im}(\omega)=0$ and hence divides the $(e,k_{c})$-plane into two parts demarcating the stable and unstable regions: the region below a curve is unstable as $\mathrm{Im}(\omega)>0$ in this region whereas that above this curve is stable as $\mathrm{Im}(\omega)<0$ in this region. In general, the NSF and all the moment theories considered here predict qualitatively similar critical wavenumber profiles. In particular, the critical wavenumber is zero in the elastic ($e=1$) case since both the longitudinal and transverse systems are stable in this case, and it increases with increasing inelasticity, in general. For nearly elastic granular gases ($0.9\leq e\leq 1$), the respective critical wavenumber profiles from the NSF and all moment theories coincide for both the longitudinal and transverse systems.

For the longitudinal system (top row in figure 11), the NSF and all the moment theories yield smooth critical wavenumber profiles for moderate to large values of the coefficient of restitution discerning the stability regions. However, each theory gives a threshold value of the coefficient of restitution below which the longitudinal system remains unstable (because for $e<e_{\mathrm{th}}$, the growth rates of some of the eigenmodes of the longitudinal system remain positive even for large wavenumbers), which is not true (see, e.g., Gupta et al., 2018). This simply means that the IMM model is not adequate for granular gases with moderate to large inelasticity. The threshold values of the coefficient of restitution $e_{\mathrm{th}}$ for the longitudinal systems associated with the NSF and different Grad moment theories are given in table 2. Furthermore, the G13, G14 and G26 theories lead to kinks at $e\approx 0.7,0.648,0.608$ in the case of $d=2$ (at $e\approx 0.562,0.522,0.446$ in the case of $d=3$), respectively, indicating that the region of applicability increases on increasing the number of moments.

For the transverse system (bottom row in figure 11), the critical wavenumbers for the G13 and G14 theories are exactly the same due to the fact that the scalar fourth moment $\check{\Delta}$ does not enter the transverse system associated with any moment theory (see (99$b$)). Hence the critical wavenumbers from both the theories are depicted by a single dot-dashed magenta line. Among the theories considered, the transverse systems associated with the NSF and G13 (or G14) theories give the critical wavenumbers for all coefficients of restitution whereas those associated with the G26 and G29 theories again lead to threshold values of the coefficient of restitution below which the transverse systems associated with them remain unstable for all wavenumbers (due to the same reason as above). This again restricts the employability of the IMM model to moderately to nearly elastic granular gases. The threshold values of the coefficient of restitution $e_{\mathrm{th}}$ for the transverse systems associated with the G26 and G29 theories are also given in table 2. The critical wavenumbers from the G26 (shown by dashed blue line) and G29 (shown by solid red line) theories closely follow each other for $0.7\lesssim e\leq 1$ in the case of $d=2$ and for $0.6\lesssim e\leq 1$ in the case of $d=3$. Similarly, the critical wavenumbers from the NSF (shown by thin dashed black line) and G13 or G14 (shown by dot-dashed magenta line) theories closely follow each other for $0.7\lesssim e\leq 1$ in the case of $d=2$ and for all values of $e$ in the case of $d=3$.

From figures 7–11, it is concluded that the instabilities of the longitudinal and transverse systems above—for moderately to nearly elastic granular gases—are confined to small wavenumbers (or long wavelengths), i.e. these instabilities are long-wave instabilities. Thus, there is a minimum system size, referred to as the critical system size, such that the instabilities will not appear in a system having size smaller than the critical system size.

6.3 Critical system size

It is well-established—through the linear stability analysis of hydrodynamic models, through the DSMC method as well as through molecular dynamics (MD) simulations—that the HCS of a freely cooling granular gas is unstable but a minimum critical system size is necessary for the onset of instabilities (see, e.g., Brey et al., 1998b; Brilliantov & Pöschel, 2004; Garzó, 2005; Gupta et al., 2018). Moreover, it is also known that during the instability phenomenon of the HCS, the instability of the unstable shear mode first engenders the formation of vortices in the system through linear effects (Brilliantov & Pöschel, 2004; Garzó, 2005) and subsequently clustering set in due to the instability of the shear mode through nonlinear effects (Brey et al., 1999; Goldhirsch, 2003). In order to verify through the moment theories that it is the unstable shear mode which is responsible for the onset of instabilities in a freely cooling granular gas, the critical wavenumbers for the longitudinal and transverse systems are plotted together in the $(e,k_{c})$-plane (but only for that range of $e$, which apparently leads to meaningful critical wavenumber profiles) in figure 12. Denoting the critical wavenumbers associated with the unstable heat and shear modes for any moment system by $k_{h}$ and $k_{s}$, respectively, it can be easily deduced that a heat (shear) mode of the longitudinal (transverse) system for wavenumbers $k>k_{h}$ ($k>k_{s}$) will always decay while that for wavenumbers $k<k_{h}$ ($k<k_{s}$) will grow exponentially. For the ranges of the coefficient of restitution shown in figure 12, $k_{s}\geq k_{h}$. Since the wavelength, and hence the critical system size, is inversely proportional to the wavenumber, it is the unstable shear mode from the transverse system which becomes unstable first. From the above discussion, the critical system size can be determined with the knowledge of the critical wavenumber for the transverse system itself. Indeed, it is possible to obtain the analytical expressions for the critical wavenumbers from the moment theories as follows.

It has been verified (although shown here only for the G29 system in the top rows of figures 7–10 in the case of $e=0.75$ for the sake of succinctness) that the real part of the frequency for the least stable eigenmode is either always zero (for longitudinal systems associated with the NSF, G14 and G26 theories and for transverse systems associated with all the theories) or is nonzero only for wavenumbers above the critical wavenumber (for longitudinal systems associated with the G13 and G29 theories). Therefore, for the least stable mode, $\omega=0$ at critical wavenumber. Consequently, the critical wavenumbers for the longitudinal and transverse systems associated with a moment theory can also be determined by substituting $\omega=0$ in their dispersion relations and solving the resulting equations for $k$. For instance, the critical wavenumbers for the longitudinal and transverse systems associated with the G29 theory can also be determined by inserting $\omega=0$ in (100) and solving the resulting equations, namely $\mathbbm{a}_{9}=0$ and $\mathbbm{b}_{6}=0$. Hence the roots of $\mathbbm{a}_{9}=0$ and $\mathbbm{b}_{6}=0$, respectively, yield the critical wavenumbers for the longitudinal and transverse systems associated with the G29 equations. It is worthwhile to note that the coefficients $\mathbbm{a}_{9}$ and $\mathbbm{b}_{6}$ are (even-degree) polynomials of degree eight and four in $k$, respectively. Similarly, the coefficients of $\omega^{0}$ in the dispersion relations for the longitudinal (transverse) systems associated with the NSF, G13, G14 and G26 are (even-degree) polynomials of degree four, four, four and six (two, two, two and four) in $k$, respectively. Consequently, each of them leads to more than one value of the critical wavenumber. Nevertheless, only one of them in each case is meaningful (positive) and that is the analytical expression for the corresponding critical wavenumber. After some algebra, the explicit expressions for $k_{h}$ and $k_{s}$ for each of the NSF and Grad moment systems, in a compact form, can be written as

[TABLE]

where the coefficients appearing in these expressions are relegated to appendix G for better readability; and “$|\mathrm{NSF}$” and “$|\mathrm{G\dots}$” in subscripts denote the NSF and Grad moment systems which the critical wavenumbers belong to. The plots of $k_{h}$ and $k_{s}$ from the analytical expressions agree completely with those in figure 12 (at least for the depicted ranges of the coefficient of restitution; for instance, $k_{h|\mathrm{G29}}$ becomes complex for smaller values of $e$).

The critical system size $L_{c}$ is estimated by $\tilde{L}_{c}=2\pi/\max{\{k_{h},k_{s}\}}$ (Garzó, 2005; Gupta et al., 2018), where $\tilde{L}_{c}:=L_{c}/\ell$ is the dimensionless critical system size and $\ell$ is the length scale defined in (84). From figure 12, $k_{s}>k_{h}$ for the NSF and all the moment theories considered in this work (for the shown ranges of $e$ in the figure); therefore the critical system size $L_{c}$ is given by

[TABLE]

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

The critical system size in units of the mean free path ($L_{c}/\ell_{0}$) is illustrated in figure 13 as a function of the coefficient of restitution $e$ for $d=2$ (left panel) and $d=3$ (right panel). The dot-dashed (magenta), dashed (blue) and solid (red) lines in the figure depict the critical system size from the G13 (or G14), G26 and G29 theories, respectively, computed with formula (112) using the analytical expressions for $k_{s}$ given in (107), (109) and (111), respectively. For comparison purpose, the figure also includes thin solid black lines depicting the critical system size from the NSF theory for IMM computed with formula (112) using the analytical expressions for $k_{s}$ given in (104). The dashed (green) lines with symbols delineate the critical system size computed from the theoretical expression given in Brey et al. (1998b), which was obtained via the linear stability analysis of a kinetic model for granular gases of IHS due to Brey et al. (1997). It should be noted that the $l_{0}$ used in Brey et al. (1998b) relates to the mean free path $\ell_{0}$ used in the present work via

[TABLE]

The dotted (black) line on the right panel of figure 13 denotes the critical system size determined from the theoretical expression obtained by the linear stability analysis of the granular NSF equations for IHS in Garzó (2005). The (red) circles on the left panel of figure 13 are the results from two-dimensional DSMC simulations carried out by Brey et al. (1998b). The triangles on the right panel of figure 13 delineate the results from MD simulations of IHS carried out by Mitrano et al. (2011) at solid fraction $\phi=0.1$ and are included only for qualitative comparison, since the present work deals with dilute granular gases for which $\phi\to 0$. Note that the critical system size in Mitrano et al. (2011) is scaled with the diameter of a particle while that in the present work is scaled with the mean free path. Therefore the MD simulations results of Mitrano et al. (2011) have been multiplied by a factor $6\sqrt{2}\,\phi\,\chi(\phi)$ while displaying them in figure 13. Here $\chi(\phi)=(2-\phi)/[2(1-\phi)^{3}]$ is the pair correlation function. The right panel of figure 13 also illustrates a solid cyan line, which depicts the critical system size for $\phi=0.1$ computed with the theoretical expressions derived in Garzó (2005), and is included just to show the good agreement between the theoretical results of Garzó (2005) and MD simulations results of Mitrano et al. (2011).

Figure 13 reveals that the critical system size from all the theories and simulations decreases with increasing inelasticity. For disk flows ($d=2$, left panel of figure 13), the critical system size from the NSF theory for IMM (thin solid black line) agrees well with that from the theoretical expression in Brey et al. (1998b). However, the critical system size from the G13 or G14 theory (dashed magenta line) seems to be slightly better than that from the NSF theory and agrees perfectly with that from the theoretical expression in Brey et al. (1998b) (dashed green line with symbols), and also agrees reasonably well with the DSMC results of Brey et al. (1998b) (red circles)—for $0.65\leq e\leq 1$. Nevertheless, the G26 and G29 theories somewhat underpredict the critical system size for all coefficients of restitution $e\gtrsim 0.65$, although their predictions are also close to the other theories for $e\gtrsim 0.9$.

For sphere flows ($d=3$, right panel of figure 13), I could not find any simulation data for the dilute limit, i.e. for solid fraction $\phi\to 0$. Therefore the data from MD simulations carried out by Mitrano et al. (2011) for solid fraction $\phi=0.1$ are included for comparison. It is important to note that the results from MD simulations by Mitrano et al. (2011) are in good agreement with those from theoretical expression of Garzó (2005) not only for $\phi=0.1$ but also for $\phi=0.4$ (see Mitrano et al., 2011, figure 9). Therefore, in the dilute limit ($\phi\to 0$), the results from the theoretical expression of Garzó (2005), shown by the dotted black line in figure 13, can be treated as a benchmark. Additionally, in this limit, the results from the theoretical expressions of Garzó (2005) and Brey et al. (1998b) are also in good agreement with each other. Clearly, the critical system size from the NSF theory for IMM is again in reasonably good agreement with that from Garzó (2005). Moreover, the critical system sizes from all the moment theories are also in good agreement with that from Garzó (2005) for $e\gtrsim 0.85$, but deviate slightly from the results of Garzó (2005) for $0.55\lesssim e\lesssim 0.85$, where the G26 and G29 theories again underpredict the critical system size while the G13 or G14 theory overpredicts it. Between the G26 and G29 theories, the latter seems to perform slightly better at moderate values of the coefficient of restitution for both $d=2$ and $d=3$.

From figure 13, it is apparent that the critical system size obtained from the NSF and Grad moment theories for IMM is in qualitatively good agreement with that obtained from the NSF-level theories and simulations for IHD/IHS, although it is also noticeable from the figure that some lower-order Grad moment theories (e.g. the G13 or G14 theory) perform better than some higher-order Grad moment theories (e.g. the G26 theory). A possible reason for this could be the choice of the effective collision frequency $\mathring{\nu}$ in the IMM model (see (3$a$)) that was chosen in such a way that the cooling rates for IHS and IMM remain exactly the same while the collisional production terms in the other moment equations for IMM follow accordingly based on this choice of the effective collision frequency. Consequently, with this choice of the effective collision frequency, even the NSF equations for IMM seem to perform better than some higher-order moment models. Therefore it would be interesting to explore other possible choices for the effective collision frequency in the future in such a way that the results from the Boltzmann equations for IHS and IMM agree in an optimal way so that a higher-order moment model would perform better than a lower-order moment model in the case of IMM, similarly to moment models for IHS.

7 Conclusion

Grad moment equations—consisting of up to 29 moments—for a $d$-dimensional dilute granular gas composed of IMM have been derived from the Boltzmann equation for IMM via the Grad moment method. A strategy for computing the collisional production terms associated with these moment equations in an automated way has been presented. Although the Maxwell interaction potential had been devised in such a way that the explicit form of the distribution function is not required to be known for determining the collisional production terms, and therefore the collisional production terms for IMM can, in principle, be evaluated using pen and paper, yet the complexity increases with an increase in the number of moments. Thus the presented strategy for computing the collisional production terms associated with the moment equations would really be useful when considering even more moments.

The transport coefficients in the NSF laws for dilute granular gases of IMM have been determined by following a procedure due to Garzó (2013) and it has been shown that the G13 equations for IMM are sufficient to derive the first-order (i.e. the NSF-level) transport coefficients and that the higher moments do not play any role in determining the NSF transport coefficients for IMM since the stress and heat flux balance equations at this order do not have any dependence on the higher moments. The higher-order moment equations will be required only for computing the transport coefficients beyond the NSF level. However, since the present work already provides some higher-order Grad moment equations, it would be interesting to compute the transport coefficients beyond the NSF level in the future by relating the Grad moment equations to the Burnett equations for IMM. Although the Grad moment theories for IMM presented in this work seem to overestimate all the transport coefficients, the NSF-level transport coefficients for IMM obtained with the Grad moment theories and with the CE expansion (Santos, 2003) are in complete agreement.

The HCS of a freely cooling granular gas has then been investigated and it has been found that the decay of the granular temperature in the HCS obeys Haff’s law but, in contrast to the case of IHS, does not depend on the higher moments since the fourth cumulant does not enter the energy balance equation for IMM. Yet, Haff’s laws for IMM and IHS have been found to be in good agreement with each other. Furthermore, the other higher moments have been found to relax much faster than the granular temperature.

As an application of the derived moment models, a linear stability analysis has been performed to scrutinise the stability of the HCS due to small perturbations. By decomposing each moment system into the longitudinal and transverse systems, it has been shown that a heat mode from the longitudinal system and a shear mode from the transverse system associated with each moment system are unstable for (moderately to nearly elastic) granular gases and that the unstable shear mode from the transverse system initiates instability in a homogeneously cooling granular gas. To assess the linear stability results, the critical system size for the onset of instability is investigated, and it has been found that the Grad moment theories for IMM yield a reasonably good estimate of the critical system size for granular gases with moderate to large coefficients of restitution.

It is important to note that the only assumption on the coefficient of restitution in this work is that it is a constant. Therefore the present work should, in principle, be applicable to granular gases with any degree of inelasticity, which is not the case unfortunately. Nevertheless, this should not be thought of as a problem with moment models presented here, rather it is a problem with the IMM model itself; for instance, the IMM model yields negative values for the coefficient of thermal conductivity below a certain value of the coefficient of restitution (Santos, 2003; Garzó & Santos, 2011).

It is anticipated that the Grad moment systems for dilute granular gases of IMM, similarly to those for monatomic gases, will also suffer from the loss of hyperbolicity. Therefore the hyperbolicity of these systems needs to be investigated in the future, which will also be useful in developing suitable numerical methods for solving them. Moreover, to overcome the undesirable consequences of the loss of hyperbolicity, a regularisation of Grad moment equations might also be necessary. The usefulness of the derived Grad moment systems is substantially limited by the unavailability of boundary conditions. Hence the development of boundary conditions complementing these Grad moment systems should be an immediate follow-up to the present work. Notwithstanding, the Grad moment systems presented in this work will be useful in describing granular processes involving large spatial gradients and are expected to pave the way to further developments including the developments of regularised moment models, required boundary conditions, efficient numerical frameworks including the MFS.

Acknowledgments

The author appreciates the anonymous referees for their informative suggestions, which have helped to improve the paper. The author gratefully acknowledges the financial supports through the “MATRICS” project MTR/2017/000693 funded by the SERB, India and through the Commonwealth Rutherford Fellowship availed at the University of Warwick, UK. The author is thankful to Prof. Vicente Garzó and Prof. Andrés Santos for some helpful suggestions on this work during the RGD31 symposium in Glasgow, to Dr. James Sprittles, Prof. Duncan Lockerby, Dr. Anirudh Singh Rana and Dr. Priyanka Shukla for some fruitful discussions, and to Prof. Manuel Torrilhon for implementing in Mathematica the Einstein summation, which has been used in this work.

Declaration of interests

The author reports no conflict of interest.

Appendix A Symmetric and traceless part of a tensor

Needless to say, the symmetric and traceless part of a tensor of rank zero or one is that tensor itself. The symmetric and traceless part of a rank two tensor $A_{ij}$ in $d$ dimensions is given by

[TABLE]

where the round brackets around the indices in (113) and in what follows always denote the symmetric part of the corresponding tensor. Here, $A_{(ij)}=(A_{ij}+A_{ji})/2$ is the symmetric part of $A_{ij}$.

The symmetric and traceless part of a rank three tensor $A_{ijk}$ in $d$ dimensions is given by

[TABLE]

where $A_{(ijk)}=(A_{ijk}+A_{ikj}+A_{jik}+A_{jki}+A_{kji}+A_{kij})/6$ is the symmetric part of $A_{ijk}$.

The symmetric and traceless part of a rank four tensor $A_{ijkl}$ in $d$ dimensions is given by

[TABLE]

where $A_{(ijkl)}=(A_{ijkl}+A_{ijlk}+A_{ikjl}+A_{iklj}+A_{iljk}+A_{ilkj}+A_{jikl}+A_{jilk}+A_{jkil}+A_{jkli}+A_{jlik}+A_{jlki}+A_{kijl}+A_{kilj}+A_{kjil}+A_{kjli}+A_{klij}+A_{klji}+A_{lijk}+A_{likj}+A_{ljik}+A_{ljki}+A_{lkij}+A_{lkji})/24$ is the symmetric part of $A_{ijkl}$.

The symmetric and traceless part of a rank $n$ tensor $A_{i_{1}i_{2}\dots i_{n}}$ in $d$ dimensions is given by

[TABLE]

where the coefficients $\beta_{n,k}$ are given by

[TABLE]

and the symmetric part of $A_{i_{1}i_{2}\dots i_{n}}$ is given by

[TABLE]

Appendix B Computation of the collisional production terms

For an arbitrary function $\psi(t,\bm{x},\bm{c})$, the collisional production term (or collisional moment) associated with it—on using the symmetry properties of the Boltzmann collision operator—reads (Garzó & Santos, 2007)

[TABLE]

where the velocity with single prime denotes the post-collisional velocity in a direct collision that transforms the pre-collision velocities $\bm{c}$ and $\bm{c}_{1}$ of the colliding molecules to the post-collisional velocities $\bm{c}^{\prime}$ and $\bm{c}_{1}^{\prime}$ via the relations (Brilliantov & Pöschel, 2004)

[TABLE]

Typically, $\psi$ is of the tensorial form: $\psi=m\,C^{2a}C_{\langle i_{1}}C_{i_{2}}\cdots C_{i_{n}\rangle}$ and hence the general form of the collisional production term is

[TABLE]

However, since the squared velocities in (B) can be easily expressed in index notation using the Einstein summation convention, for instance $C^{2}=C_{i_{0}}C_{i_{0}}$, and the indices in each term of (B) can be adjusted accordingly, it is convenient to first compute

[TABLE]

instead of computing (B) directly.

To compute the right-hand side of (122), the post-collisional peculiar velocities (marked with primes) in (122) are replaced with the pre-collisional peculiar velocities by exploiting the definition of the peculiar velocity and relation (120$a$). This changes the product of the post-collisional peculiar velocities in (122) to

[TABLE]

where $w_{0}=(1+e)/2$ and the round brackets around the indices again denote the symmetric part of the corresponding tensor (see appendix A for its definition). Substituting (123) into (122), one obtains

[TABLE]

where

[TABLE]

is termed as the scattering vector integral with $\hat{\bm{g}}=\bm{g}/g$. The structure of the integrand in (125) suggests that $I_{i_{1}i_{2}\dots i_{n}}$ will have the form

[TABLE]

where the unknown coefficients $a_{\beta}^{(n)}$ depend only on the dimension $d$, and are computed separately in § B.1. Insertion of (126) into (124) transpires the tensorial structure of $P_{i_{1}\ldots i_{n}}$:

[TABLE]

Expression (B), without the prefactor $m\,\mathring{\nu}/(n\,\Omega_{d})$, is entered as the starting point in Mathematica® script.

For general interaction potentials, specific forms of the distribution functions $f(\bm{c})$ and $f(\bm{c}_{1})$ must be provided in order to compute $P_{i_{1}\ldots i_{n}}$ further. Nevertheless, for IMM, specific forms of $f(\bm{c})$ and $f(\bm{c}_{1})$ are not required since the coefficients $a_{\beta}^{(j)}$ for IMM do not depend on the relative velocity $\bm{g}$. Indeed, for IMM, now all the components and the magnitude of the relative velocity $\bm{g}$ are replaced in terms of the peculiar velocities $\bm{C}$ and $\bm{C}_{1}$ by using the relation $\bm{g}=\bm{C}-\bm{C}_{1}$. It may be noted that the exponent of the magnitude of $\bm{g}$ is even in (B), which makes it easier to replace $g^{2\beta}$ in (B) using the relation $g^{2}=C^{2}+C_{1}^{2}-2C_{k}C_{1k}$. At this step, some vanishing integrals, such as

[TABLE]

are automatically taken care of in the Mathematica® script. Now, the double integrals in each term under the summation in (B) can be written as a product of two independent integrals, one over $\bm{c}$ and the other over $\bm{c}_{1}$, and each of these integrals can be expressed in terms of the considered moments. Note that the present work deals with the traceless moments and all the terms should be expressed as traceless moments. This is not very straightforward for tensors of rank more than three. Nevertheless, this step has also been incorporated in the Mathematica® script to express all the results in terms of traceless tensors. Finally, taking the traceless part of each term in the result, one obtains the required collisional production term. All the collisional production terms obtained in this work agree with those obtained in Garzó & Santos (2007) till fourth order, which validates the code. In principle, this Mathematica® script would be able to compute the collisional production terms for moments of any order. Nonetheless, as the code is not optimised, it takes a significantly long computation time in computing the collisional production terms associated with more than sixth-order moments.

B.1 Computation of the coefficients $a_{\beta}^{(n)}$

The unknown coefficients $a_{\beta}^{(n)}$ follow by appropriately contracting the two forms of $I_{i_{1}\ldots i_{n}}$ in (125) and (126) with combinations of $\hat{g}_{i}=g_{i}/g$ and with combinations of Kronecker deltas, successively. This results into linear systems of algebraic equations, which yield the coefficients $a_{\beta}^{(n)}$ as functions of the scalar integrals given by (van Noije & Ernst, 1998; Garzó & Santos, 2007)

[TABLE]

From (125) and (126), one has

[TABLE]

Contracting the above equation with $\hat{g}_{i}$ and using the fact that $\hat{k}_{i}\hat{g}_{i}=\hat{\bm{k}}\cdot\hat{\bm{g}}$, it readily follows that

[TABLE]

Again, from (125) and (126), one has

[TABLE]

Contracting the above equation with $\hat{g}_{i}\hat{g}_{j}$ and with $\delta_{ij}$ successively, one obtains

[TABLE]

and, thus

[TABLE]

The next integrals are treated analogously and one, eventually, finds

[TABLE]

[TABLE]

and so on; see Mathematica® file “kintegrals.nb” provided as supplementary material.

Appendix C Coefficients in the collisional production terms

The coefficients in the collisional production terms (24)–(30) associated with the G29 equations for IMM are as follows.

[TABLE]

Appendix D The G29 distribution function

The computation of the G29 distribution function is comparatively easier with the dimensionless variables. Let us introduce the dimensionless variables (denoted with bars) as follows:

[TABLE]

In the dimensionless variables, the definitions of the 29 moments can be recast as

[TABLE]

Let the G29 distribution function in the dimensionless form be given by

[TABLE]

where the angle brackets again denote the symmetric-traceless tensors and $\lambda$’s are the unknown coefficients that are determined by replacing $\bar{f}$ with $\bar{f}_{|\mathrm{G29}}$ in definitions (157), and solving the resulting system of algebraic equations for $\lambda$’s.

The integrals over velocity space are typically evaluated by transforming the integral from a $d$-dimensional Cartesian coordinate system to a $d$-dimensional spherical coordinate system. A useful identity, which employs this transformation, for evaluating the integral of an even function $h(C)$ in $C$ over the velocity space $\bm{C}$ is

[TABLE]

where the following identities have been used: for $n\geq 0$,

[TABLE]

Note that a similar integral for an odd function $h(C)$ in $C$ over the velocity space $\bm{C}$ vanishes. Furthermore, for an even function $h(C)$ in $C$, it can be shown that

[TABLE]

and, in general,

[TABLE]

for an even $n$ while the integral vanishes for an odd $n$. As a consequence of (160)–(162), it is straightforward to show that

[TABLE]

Now, replacing $\bar{f}$ with $\bar{f}_{|\mathrm{G29}}$ in definitions (157), and using identities (D)–(167), one obtains

[TABLE]

These equations yield

[TABLE]

Inserting these coefficients in (D), the dimensionless G29 distribution function reads

[TABLE]

which on introducing the dimensions using (151) yields the G29 distribution function (3.3).

Appendix E Explicit components of the traceless gradients

The explicit components of the traceless gradients in (86)–(94) are computed as follows. Using (113),

[TABLE]

where $\Phi\in\{v,q,\varphi\}$. From the above equation, it follows that

[TABLE]

The other components, if needed, can be computed analogously.

For a symmetric-traceless rank two tensor $\tilde{\Phi}_{ij}$, definition (114) gives

[TABLE]

where $\Phi\in\{\sigma,R\}$ in the present work. From the above equation, it follows that

[TABLE]

The other components, if needed, can be computed analogously.

Appendix F Coefficient matrices in (99) and (102)

The matrix $\mathscr{L}$ in the longitudinal problem (99$a$) can be written in the form of block matrices for better readability as

[TABLE]

with the block matrices being

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and the matrix $\mathscr{T}$ in the transverse problem (99$b$) reads

[TABLE]

where $I_{6}$ is the identity matrix of dimensions $6\times 6$.

The matrix $\mathscr{L}_{\mathrm{NSF}}$ in (102) reads

[TABLE]

where $I_{3}$ is the identity matrix of dimensions $3\times 3$.

Appendix G Coefficients in the analytical expressions of the critical wavenumbers

Using the abbreviations given in (97), the coefficients appearing in the analytical expressions (105)–(111) of the critical wavenumbers computed from various moment systems can be written as follows.

[TABLE]

with

[TABLE]