Kinetic theory of granular particles immersed in a molecular gas

Abstract

The transport coefficients of a dilute gas of inelastic hard spheres immersed in a gas of elastic hard spheres (molecular gas) are determined. We assume that the number density of the granular gas is much smaller than that of the surrounding molecular gas, so that the latter is not affected by the presence of the granular particles. In this situation, the molecular gas may be treated as a thermostat (or bath) of elastic hard spheres at a fixed temperature. The Boltzmann kinetic equation is the starting point of the present work. The first step is to characterise the reference state in the perturbation scheme, namely the homogeneous state. Theoretical results for the granular temperature and kurtosis obtained in the homogeneous steady state are compared against Monte Carlo simulations showing a good agreement. Then, the Chapman–Enskog method is employed to solve the Boltzmann equation to first order in spatial gradients. In dimensionless form, the Navier–Stokes–Fourier transport coefficients of the granular gas are given in terms of the mass ratio $m/m_g$ ($m$ and $m_g$ being the masses of a granular and a gas particle, respectively), the (reduced) bath temperature and the coefficient of restitution. Interestingly, previous results derived from a suspension model based on an effective fluid–solid interaction force are recovered in the Brownian limit ($m/m_g \to \infty$). Finally, as an application of the theory, a linear stability analysis of the homogeneous steady state is performed showing that this state is always linearly stable.

1. Introduction

A challenging problem in statistical physics is the understanding of multiphase flows, namely the flow of solid particles in two or more thermodynamic phases. Needless to say, these types of flows occur in many industrial settings (such as circulating fluidised beds) and can also affect our daily lives due to the fact that the comprehension of them may ensure vital needs of humans such as clean air and water (Subramaniam 2020). Among the different types of multiphase flows, a particularly interesting set corresponds to the so-called particle-laden suspensions in which small, immiscible and typically dilute particles are immersed in a carrier fluid. The dynamics of gas–solid flows is rich and extraordinarily complex (Gidaspow 1994; Jackson 2000; Koch & Hill 2001; Fox 2012; Tenneti & Subramaniam 2014; Fullmer & Hrenya 2017; Lattanzi et al. 2020) so their understanding poses a great challenge. Even the study of granular flows in which the effect of interstitial fluid is neglected (Campbell 1990; Goldhirsch 2003; Brilliantov & Pöschel 2004; Rao & Nott 2008; Garzó 2019) entails enormous difficulties.

In the case that the particle-laden suspensions are dominated by collisions (Subramaniam 2020), the extension of the classical kinetic theory of gases (Chapman & Cowling 1970; Ferziger & Kaper 1972; Résibois & de Leener 1977) to granular suspensions can be considered as an appropriate tool to model these systems. In this context and assuming nearly instantaneous collisions, the influence of gas-phase effects on the dynamics of solid particles is usually incorporated in the starting kinetic equation in an effective way via a fluid–solid interaction force (Koch 1990; Gidaspow 1994; Jackson 2000). Some models for granular suspensions (Louge, Mastorakos & Jenkins 1991; Tsao & Koch 1995; Sangani et al. 1996; Wylie et al. 2009; Parmentier & Simonin 2012; Heussinger 2013; Wang et al. 2014; Saha & Alam 2017; Alam, Saha & Gupta 2019; Saha & Alam 2020) only consider the Stokes linear drag law for gas–solid interactions. Other models (Garzó et al. 2012) include also an additional Langevin-type stochastic term.

For small Knudsen numbers, the Langevin-like suspension model mentioned above (Garzó et al. 2012) has been solved by means of the Chapman–Enskog method (Chapman & Cowling 1970) adapted to dissipative dynamics. Explicit expressions for the Navier–Stokes–Fourier transport coefficients have been obtained in terms of the coefficient of restitution and the parameters of the suspension model (Garzó et al. 2012; Gómez González & Garzó 2019). Knowledge of the forms of the transport coefficients has allowed an assessment of not only the impact of inelasticity on them (which was already analysed in the case of dry granular fluids Brey et al. 1998; Garzó & Dufty 1999) but also the influence of the interstitial gas on the momentum and heat transport. Beyond the Navier–Stokes domain, this type of suspension model has also been considered to compute the rheological properties in sheared gas–solid suspensions (see e.g. Tsao & Koch 1995; Sangani et al. 1996; Parmentier & Simonin 2012; Heussinger 2013; Seto et al. 2013; Kawasaki, Ikeda & Berthier 2014; Chamorro, Vega Reyes & Garzó 2015; Hayakawa, Takada & Garzó 2017; Saha & Alam 2017; Alam et al. 2019; Hayakawa & Takada 2019; Gómez González & Garzó 2020; Saha & Alam 2020; Takada et al. 2020).

The quantitative and qualitative accuracies of the (approximate) analytical results derived from the kinetic-theory two-fluid model (Garzó et al. 2012) have been confronted against computer simulations in several problems. In particular, the critical length for the onset of velocity vortices in the homogeneous cooling state of gas–solid flows obtained from a linear stability analysis presents an acceptable agreement with molecular dynamics (MD) simulations carried out for strong inelasticity (Garzó et al. 2016). Simulations using a computational fluid dynamics solver (Capecelatro & Desjardins 2013; Capecelatro, Desjardins & Fox 2015) of Radl & Sundaresan (2014) have shown a good agreement in the mean slip velocity with the kinetic-theory predictions (Fullmer & Hrenya 2016). On the other hand, kinetic theory has also been assessed for describing clustering instabilities in sedimenting fluid–solid systems; good agreement is found at high solid-to-fluid density ratios although the agreement is weaker for intermediate and low density ratios (Fullmer et al. 2017). In the case of non-Newtonian flows, the theoretical results (Saha & Alam 2017; Alam et al. 2019; Saha & Alam 2020) derived from the Stokes drag model for the ignited–quenched transition and the rheology of a sheared gas–solid suspension have been shown to compare very well with computer simulations. Regarding the Langevin-like model (Garzó et al. 2012), the rheological properties of a moderately dense inertial suspension computed by a simpler version of this model exhibit a quantitatively good agreement with MD simulations in the high-density region (Takada et al. 2020). In addition, the extension to binary mixtures of this suspension model has been tested against Monte Carlo data and MD simulations for both time-dependent and steady homogeneous states with an excellent agreement (Khalil & Garzó 2014; Gómez González, Khalil & Garzó 2020; Gómez González & Garzó 2021).

In spite of the reliability of the generalised Langevin and Stokes drag models for capturing in an effective way the impact of gas phase on grains, it would be desirable to propose a suspension model that considers the real collisions between solid and gas particles. In the context of kinetic theory and as already mentioned in previous works (Gómez González et al. 2020), a possibility would be to describe gas–solid flows in terms of a set of two coupled kinetic equations for the one-particle velocity distribution functions of the solid and gas phases. Nevertheless, the determination of the transport coefficients of the solid particles starting from the above suspension model is a very intricate problem. A possible way of overcoming the difficulties inherent to the description of gas–solid flows when one attempts to involve the different types of collisions is to assume that the properties of the gas phase are unaffected by the presence of solid particles. In fact, although sometimes not explicitly stated, this is one of the overarching assumptions in most of the suspension models reported in the granular literature. This assumption can be clearly justified in the case of particle-laden suspensions where the granular particles (or ‘granular gas’) are sufficiently rarefied (dilute particles), and hence the properties of the interstitial gas can be supposed to be constant. This means that the background gas can be treated as a thermostat at a constant temperature $T_g$.

Under these conditions and inspired by the work of Biben, Martin & Piasecki (2002), we propose here the following suspension model. We consider a set of granular particles immersed in a bath of elastic particles (molecular gas) at equilibrium at a certain temperature $T_g$. While the collisions between granular particles are inelastic (and characterised by a constant coefficient of normal restitution $\alpha$), the collisions between the granular and gas particles are considered to be elastic. In the homogeneous steady state (HSS), the energy lost by the solid particles due to their collisions among themselves is exactly compensated for by the energy gained by the grains due to their elastic collisions with particles of the molecular gas. In other words, the gas of inelastic hard spheres (granular gas) is thermostatted by a bath of elastic hard spheres. The dynamic properties of this system in HSSs were studied years ago independently by Biben et al. (2002) and Santos (2003). Our goal here is to go beyond the homogeneous state and determine the transport coefficients of a granular gas immersed in a molecular gas when the magnitude of the spatial gradients is small (Navier–Stokes domain).

It is quite apparent that this suspension model (granular particles plus molecular gas) can be seen as a binary mixture in which the concentration of one of the species (tracer species or granular particles) is much smaller than the other one (excess species or molecular gas). In these conditions, it is reasonable to assume that the state of the background gas (excess species) is not perturbed by the presence of the tracer species (granular particles). In addition, although the density of grains is very small, we will take into account not only the collisions between solid and gas particles, but also the grain–grain collisions in the kinetic equation of the one-particle distribution function $f(\boldsymbol {r},\boldsymbol {v};t)$ of the granular gas. In spite of the simplicity of the model, it can be considered sufficiently robust since it retains most of the basic features of granular suspensions such as the competition between the different spatial and time scales. As we show in this paper, in contrast to previous suspension models reported in the granular literature (Koch 1990; Gidaspow 1994; Jackson 2000), the present model incorporates a new parameter: the ratio between the mass $m$ of a granular particle and the mass $m_g$ of a particle of the molecular gas (modelled as an elastic gas of hard spheres).

The objective of the present paper is twofold. On the one hand, we want to determine the conditions under which the expressions for the transport coefficients derived from a Langevin-like suspension model (Garzó et al. 2012; Gómez González & Garzó 2019) are consistent with those achieved here from a collisional model. A careful analysis shows that the present results reduce to those previously found (Garzó et al. 2012; Gómez González & Garzó 2019) when $m \gg m_g$ (Brownian limit). Apart from assessing the consistency, the analysis allows one to express the drift or friction coefficient $\gamma$ (which is a free parameter in the Langevin-like model) in terms of the mass ratio $m/m_g$ and the bath temperature $T_g$. On the other hand, beyond the Brownian limit, we extend the expressions of transport properties to arbitrary values of the mass ratio. This allows us to offer a theory that can be employed not only in gas–solid systems for relatively massive particles (for which the analytical results obtained from the Langevin model are quite useful) but also in situations where the mass of granular particles is comparable to that of the elastic gas.

However, surprisingly, the results derived here for the transport coefficients are practically indistinguishable from those obtained from the Langevin model (Gómez González & Garzó 2019) for not relatively large values of the mass ratio (for a typical value of the reduced bath temperature $T_g^{*}=1000$, the Langevin results converge to those reported here for $m/m_g \simeq 50$). This means that the range of mass (or size) ratios where the collisional suspension model offers new results not covered by the Langevin-like model (Gómez González & Garzó 2019) is constrained to situations where the mass (or size) ratios of grains and gas particles are comparable. This is of course an important limitation of our model, especially if one is interested in real applications (fine aerosol particles in air) where the mass ratio $m/m_g$ is large.

Regarding the above point, some doubts are raised concerning the mechanisms that govern a gas–solid collision when the mass of the solid particles and that of the surrounding gas are comparable. Two different options are equally valid: either grains are no longer of a mesoscopic size, or the gas particles become granular ones. In the first case, we are dealing with a binary mixture of molecular gases, while the second case corresponds to a mixture of granular gases. In both cases, we assume that the concentration of one of the species is negligible (tracer limit), and so the state of the excess species is not affected by the other one. Although the limitation of comparable mass ratios reduces the applicability of the present model, it is well known that a vast number of parameters influence the dynamics of a collision. The sizes, surface properties and material that constitute the particles may perturb the processes of fracture, friction or internal vibrations that regulate the loss of energy in each collision to a greater or lesser extent. Through an appropriate selection of the particles’ features, the granular suspension can therefore still be modelled as a granular gas immersed in an ensemble of elastic hard spheres. Although this latter situation is not likely quite frequent in nature or in industrial set-ups (however, it seems feasible in protoplanetary disks Schneider et al. 2021), we think that the results reported in this paper for comparable masses or sizes may prove to be still useful for analysing computer simulation results where the granular gas may be thermostatted by a bath of elastic hard spheres (Biben et al. 2002).

The plan of the paper is as follows. The Boltzmann kinetic equation for a granular gas thermostatted by a bath of elastic hard spheres is presented in § 2 along with the corresponding balance equations for the densities of mass, momentum and energy. The Brownian limit ($m/m_g \to \infty$) is also considered; in this limit the Boltzmann–Lorentz operator (accounting for the rate of change of $f$ due to the elastic collisions between grains and gas particles) reduces to the Fokker–Planck operator (Résibois & de Leener 1977; McLennan 1989), which is the basis of the Langevin-like suspension model (Garzó et al. 2012). Section 3 is devoted to the study of the HSS. Although the HSS was already analysed by Santos (2003) for a three-dimensional system ($d=3$), we revisit here this study by extending the analysis to an arbitrary number of dimensions $d$. Section 4 addresses the application of the Chapman–Enskog-like expansion (Chapman & Cowling 1970) to the Boltzmann kinetic equation. Since the system is slightly disturbed from the HSS, the expansion is around the local version of the homogeneous state which is in general a time-dependent distribution. Explicit expressions for the Navier–Stokes–Fourier transport coefficients are obtained in § 5 by considering the leading terms in a Sonine polynomial expansion. As an application of the results reported in § 5, a linear stability analysis of the HSS is carried out in § 6. As expected, the analysis shows that the HSS is linearly stable regardless of the value of the mass ratio $m/m_g$. Finally, in § 7 we summarise our main conclusions.

2. Boltzmann kinetic equation for a granular gas surrounded by a molecular gas

We consider a gas of inelastic hard disks ($d=2$) or spheres ($d=3$) of mass $m$ and diameter $\sigma$. The spheres are assumed to be perfectly smooth, so that collisions between any two particles of the granular gas are characterised by a (positive) constant coefficient of normal restitution $\alpha \leq 1$. When $\alpha =1$ ($\alpha <1$), the collisions are elastic (inelastic). The granular gas is immersed in a gas of elastic hard disks or spheres of mass $m_g$ and diameter $\sigma _g$ (‘molecular gas’). A collision between a granular particle and a particle of the molecular gas is considered to be elastic. As discussed in § 1, we are interested here in describing a situation where the granular gas is sufficiently rarefied (the number density of granular particles is much smaller than that of the molecular gas) so that the state of the molecular gas is not affected by the presence of solid (grains) particles. In this sense, the background (molecular) gas may be treated as a thermostat, which is at equilibrium at a temperature $T_g$. Thus, the velocity distribution function $f_g$ of the molecular gas is the Maxwell–Boltzmann distribution:

(2.1)\begin{equation} f_g(\boldsymbol{V}_g)=n_g \left(\frac{m_g}{2{\rm \pi} T_g}\right)^{d/2} \exp \bigg(-\frac{m_g V_g^{2}}{2T_g}\bigg), \end{equation}

where $n_g$ is the number density of the molecular gas and $\boldsymbol {V}_g=\boldsymbol {v}-\boldsymbol {U}_g$, in which $\boldsymbol {U}_g$ is the mean flow velocity of the molecular gas. Note that here, for the sake of generality, we have assumed that $\boldsymbol {U}_g \neq \boldsymbol {U}$ ($\boldsymbol {U}$ being the mean flow velocity of the granular gas; see its definition in (2.7)). In addition, for the sake of simplicity, the Boltzmann constant $k_{B}=1$ throughout the paper.

In the low-density regime, the time evolution of the one-particle velocity distribution function $f(\boldsymbol {r}, \boldsymbol {v}, t)$ of the granular gas is given by the Boltzmann kinetic equation. Since the granular particles collide among themselves and with the particles of the molecular gas, in the absence of external forces the velocity distribution $f(\boldsymbol {r}, \boldsymbol {v}, t)$ verifies the kinetic equation:

(2.2)\begin{equation} \frac{\partial f}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot} \boldsymbol{\nabla} f=J[\boldsymbol{v}\,|\,f,f]+J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

Here, the Boltzmann collision operator $J[\boldsymbol {v}\,|\,f,f]$ gives the rate of change of the distribution $f(\boldsymbol {r}, \boldsymbol {v}, t)$ due to binary inelastic collisions between granular particles. On the other hand, the Boltzmann–Lorentz operator $J_g[\boldsymbol {v}\,|\,f,f_g]$ accounts for the rate of change of the distribution $f(\boldsymbol {r}, \boldsymbol {v}, t)$ due to elastic collisions between granular and molecular gas particles.

The explicit form of the nonlinear Boltzmann collision operator $J[\boldsymbol {v}\,|\,f,f]$ is (Garzó 2019)

(2.3)\begin{align} J[\boldsymbol{v}_1\,|\,f,f]=\sigma^{d-1}\int {\rm d}\boldsymbol{v}_{2} \int {\rm d}\widehat{\boldsymbol {\sigma }}\varTheta (\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12}) (\widehat{\boldsymbol {\sigma }} \boldsymbol{\cdot} \boldsymbol{g}_{12}) [\alpha^{{-}2}f(\boldsymbol{v}_{1}'')f(\boldsymbol{v}_{2}'')-f(\boldsymbol{v}_{1})f(\boldsymbol{v}_{2})], \end{align}

where $\boldsymbol {g}_{12}=\boldsymbol {v}_1-\boldsymbol {v}_2$ is the relative velocity, $\widehat {\boldsymbol {\sigma }}$ is a unit vector along the line of centres of the two spheres at contact and $\varTheta$ is the Heaviside step function. In (2.3), the double primes denote pre-collisional velocities. The relationship between pre-collisional $(\boldsymbol {v}_1'', \boldsymbol {v}_2'')$ and post-collisional $(\boldsymbol {v}_1, \boldsymbol {v}_2)$ velocities is

(2.4a,b)\begin{equation} \boldsymbol{v}_1''=\boldsymbol{v}_1-\frac{1+\alpha}{2\alpha}(\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12})\widehat{\boldsymbol {\sigma }},\quad \boldsymbol{v}_2''=\boldsymbol{v}_2+\frac{1+\alpha}{2\alpha}(\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12})\widehat{\boldsymbol {\sigma }}. \end{equation}

The form of the linear Boltzmann–Lorentz collision operator $J_g[\boldsymbol {v}\,|\,f,f_g]$ is (Résibois & de Leener 1977; Garzó 2019)

(2.5)\begin{align} J_g[\boldsymbol{v}_1\,|\,f,f_g]=\overline{\sigma}^{d-1}\int {\rm d}\boldsymbol{v}_{2} \int {\rm d}\widehat{\boldsymbol {\sigma }}\varTheta (\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12}) (\widehat{\boldsymbol {\sigma }}\boldsymbol{\cdot} \boldsymbol{g}_{12})[f(\boldsymbol{v}_{1}'')f_g(\boldsymbol{v}_{2}'')-f(\boldsymbol{v}_{1})f_g(\boldsymbol{v}_{2})], \end{align}

where $\overline{\sigma }=(\sigma +\sigma _g)/2$. In (2.5), the relationship between $(\boldsymbol {v}_1'',\boldsymbol {v}_2'')$ and $(\boldsymbol {v}_1,\boldsymbol {v}_2)$ is

(2.6a,b)\begin{equation} \boldsymbol{v}_1''=\boldsymbol{v}_1-2 \mu_g(\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12})\widehat{\boldsymbol {\sigma }},\quad \boldsymbol{v}_2''=\boldsymbol{v}_2+2 \mu (\widehat{{\boldsymbol {\sigma }}} \boldsymbol{\cdot} \boldsymbol{g}_{12})\widehat{\boldsymbol {\sigma }}, \end{equation}

where $\mu _g=m_g/(m+m_g)$ and $\mu =m/(m+m_g)$.

The relevant hydrodynamic fields of the granular gas are the number density $n(\boldsymbol {r}; t)$, the mean flow velocity $\boldsymbol {U}(\boldsymbol {r}; t)$ and the granular temperature $T(\boldsymbol {r}; t)$. They are defined, respectively, as

(2.7)\begin{equation} \left\{n, n \boldsymbol{U}, dn T\right\}=\int {\rm d}\boldsymbol{v} \{1, \boldsymbol{v}, m V^{2}\} f(\boldsymbol{v}), \end{equation}

where $\boldsymbol {V}=\boldsymbol {v}-\boldsymbol {U}$ is the peculiar velocity. As said before, in general the mean flow velocity $\boldsymbol {U}$ of solid particles is different from the mean flow velocity $\boldsymbol {U}_g$ of molecular gas particles. As we show later, the difference $\boldsymbol {U}-\boldsymbol {U}_g$ induces a non-vanishing contribution to the heat flux.

The macroscopic balance equations for the granular gas are obtained by multiplying (2.2) by $\{1, \boldsymbol {v}, m V^{2}\}$ and integrating over velocity. The result is

(2.8)$$\begin{gather} D_t n+n\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}=0, \end{gather}$$

(2.9)$$\begin{gather}\rho D_t\boldsymbol{U}={-}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{\mathsf{P}}+\boldsymbol{\mathcal{F}}[f], \end{gather}$$

(2.10)$$\begin{gather}D_tT+\frac{2}{dn}\left(\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{q}+\boldsymbol{\mathsf{P}}:\boldsymbol{\nabla}\boldsymbol{U}\right)={-} T \zeta-T \zeta_g. \end{gather}$$

Here, $D_t=\partial _t+\boldsymbol {U}\boldsymbol {\cdot }\boldsymbol {\nabla }$ is the material derivative, $\rho =m n$ is the mass density of solid particles and the pressure tensor $\boldsymbol{\mathsf{P}}$ and the heat flux vector $\boldsymbol {q}$ are given, respectively, as

(2.11a,b)\begin{equation} \boldsymbol{\mathsf{P}}=\int {\rm d}\boldsymbol{v} m \boldsymbol{V}\boldsymbol{V} f(\boldsymbol{v}),\quad \boldsymbol{q}=\int {\rm d}\boldsymbol{v} \frac{m}{2} V^{2} \boldsymbol{V} f(\boldsymbol{v}). \end{equation}

Since the Boltzmann–Lorentz collision term $J_g[\boldsymbol {v}\,|\,f,f_g]$ does not conserve momentum, then the production of momentum $\boldsymbol {\mathcal {F}}[f]$ is in general different from zero. It is defined as

(2.12)\begin{equation} \boldsymbol{\mathcal{F}}[f]=\int {\rm d}\boldsymbol{v} m \boldsymbol{V} J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

In addition, the partial production rates $\zeta$ and $\zeta _g$ are given, respectively, as

(2.13a,b)\begin{equation} \zeta={-}\frac{m}{d n T}\int {\rm d}\boldsymbol{v} V^{2} J[\boldsymbol{v}\,|\,f,f], \quad \zeta_g={-}\frac{m}{d n T}\int {\rm d}\boldsymbol{v} V^{2} J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

The cooling rate $\zeta$ gives the rate of kinetic energy loss due to inelastic collisions between particles of the granular gas. It vanishes for elastic collisions. The term $\zeta _g$ gives the transfer of kinetic energy between the particles of the granular and molecular gases. It vanishes when the granular and molecular gases are at the same temperature ($T_g=T$).

The macroscopic size of grains entails that their gravitational potential energy is much larger than the usual thermal energy scale $k_{B}T$. For example, a grain of common sand at room temperature ($T=300$ K) would require a energy of the order of $10^{5}k_{B}T$ to rise a distance equal to its diameter when subjected to the action of gravity (Heinrich, Nagel & Behringer 1996). Therefore, thermal fluctuations have a negligible effect on the dynamics of grains, and so they are considered athermal systems. On the other hand, when particles are subjected to a strong excitation (e.g. vibrating walls or air-fluidised beds), the external energy supplied to the system can compensate for the energy dissipated by collisions and the effects of gravity. In this situation (rapid-flow conditions), particles’ velocities acquire some kind of random motion that looks much like the motion of the atoms or molecules in an ordinary or molecular gas. The rapid-flow regime opens up the possibility of establishing a relation between particles’ response to the external supply of energy and some kind of temperature. In this context, the granular temperature $T$ can be interpreted as a statistical quantity measuring the deviations (or fluctuations) of the velocities of grains with respect to its mean value $\boldsymbol {U}$. Just as in the ordinary case, the velocity fluctuations approach is the basic assumption in the construction of out-of-equilibrium theories such as the kinetic theory. Since granular gases are athermal, the granular temperature $T$ has no relation to the conventional thermodynamic temperature associated with the second thermodynamic law. On the other hand, although the classical ensemble averages provide a thermodynamic interpretation for the temperature $T_g$ of the bath (or molecular gas) through the definition of entropy, statistical mechanics shows that the thermodynamic temperature is the same as that obtained using the velocity fluctuations approach (up to a factor including the mass of the particles and the Boltzmann constant) (see for instance the review paper of Goldhirsch (2008) for a more detailed discussion of this issue). Moreover, coming back to the suspensions framework, the fluctuation–dissipation theorem supports the athermal statistical interpretation of $T_g$ since its value has been demonstrated to coincide with an effective temperature derived from the Einstein relation in both experiments and simulations (Puglisi, Baldassarri & Loreto 2002; Garzó 2004; Chen & Hou 2014). Thus, $T_g$ is treated also here in the same way as $T$. Namely, $T_g$ (defined through the distribution $f_g(\boldsymbol {v})$) is used to assess the deviations of the velocities of the molecular gas with respect to the mean value $\boldsymbol {U}_g$. In addition, we also assume that the values of $m$, $n$, $T$ and $T_g$ are such that a description based on classical mechanics is appropriate.

2.1. Brownian limit ($m/m_g \to \infty$)

The Boltzmann equation (2.2) applies in principle for arbitrary values of the mass ratio $m/m_g$. On the other hand, a physically interesting situation arises in the so-called Brownian limit, namely when the granular particles are much heavier than the particles of the surrounding molecular gas ($m/m_g \to \infty$). In this case, a Kramers–Moyal expansion (Résibois & de Leener 1977; Rodríguez, Salinas-Rodríguez & Dufty 1983; McLennan 1989) in the velocity jumps $\delta \boldsymbol {v}=(2/(1+m/m_g)) (\widehat {\boldsymbol {\sigma }}\boldsymbol {\cdot } \boldsymbol {g}_{12}) \boldsymbol {g}_{12}$ allows us to approximate the Boltzmann–Lorentz operator $J_g[\boldsymbol {v}\,|\,f,f_g]$ by the Fokker–Planck operator $J_g^{{FP}}[\boldsymbol {v}\,|\,f,f_g]$ (Résibois & de Leener 1977; Rodríguez et al. 1983; McLennan 1989; Brey, Dufty & Santos 1999a; Sarracino et al. 2010):

(2.14)\begin{equation} J_g[\boldsymbol{v}\,|\,f,f_g]\to J_g^{{FP}}[\boldsymbol{v}\,|\,f,f_g]=\gamma \frac{\partial}{\partial \boldsymbol{v}}\boldsymbol{\cdot} \left(\boldsymbol{v}+\frac{T_g}{m} \frac{\partial}{\partial \boldsymbol{v}}\right)f(\boldsymbol{v}), \end{equation}

where the drift or friction coefficient $\gamma$ is defined as

(2.15)\begin{equation} \gamma=\frac{4{\rm \pi}^{(d-1)/2}}{d\varGamma\left(\dfrac{d}{2}\right)} \biggl(\frac{m_g}{m}\biggr)^{1/2}\left(\frac{2T_g}{m}\right)^{1/2}n_g \overline{\sigma}^{d-1}. \end{equation}

While obtaining (2.14) and (2.15), it has been assumed that $\boldsymbol {U}_g=\boldsymbol {0}$ and that the distribution $f(\boldsymbol {v})$ of the granular gas is a Maxwellian distribution.

Most of the suspension models employed in the granular literature to fully account for the influence of an interstitial molecular fluid on the dynamics of grains are based on the replacement of $J_g[\boldsymbol{v}\,|\,f,f_g]$ by the Fokker–Planck operator (2.14) (Koch & Hill 2001). More specifically, for general inhomogeneous states, the impact of the background molecular gas on solid particles is through an effective force composed of three different terms: (i) a term proportional to the difference ${\rm \Delta} \boldsymbol {U}=\boldsymbol {U}-\boldsymbol {U}_g$, (ii) a drag force term mimicking the friction of grains on the viscous interstitial gas and (iii) a stochastic Langevin-like term accounting for the energy gained by grains due to their interactions with particles of the molecular gas (neighbouring particles effect) (Garzó et al. 2012). This yields the following kinetic equation for gas–solid suspensions:

(2.16)\begin{equation} \frac{\partial f}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} f-\gamma{\rm \Delta} \boldsymbol{U}\boldsymbol{\cdot} \frac{\partial f}{\partial \boldsymbol{v}}-\gamma \frac{\partial}{\partial \boldsymbol{v}}\boldsymbol{\cdot} \boldsymbol{V} f-\gamma \frac{T_{g}}{m}\frac{\partial^{2} f}{\partial v^{2}}=J[\boldsymbol{v}\,|\,f,f]. \end{equation}

Note that there are three different scalars $(\beta,\gamma,\xi )$ in the suspension model proposed by Garzó et al. (2012); each one of the coefficients is associated with the different terms of the fluid–solid force. For the sake of simplicity, the results derived by Gómez González & Garzó (2019) were obtained by assuming that $\beta =\gamma =\xi$.

Since the model attempts to mimic gas–solid flows where $\boldsymbol {U}\neq \boldsymbol {U}_g$, note that one has to make the replacement $\boldsymbol {v}\to \boldsymbol {v}-\boldsymbol {U}_g$ in (2.14) to obtain (2.16). The Boltzmann equation (2.16) has been solved by means of the Chapman–Enskog method (Chapman & Cowling 1970) to first order in spatial gradients. Explicit forms for the Navier–Stokes–Fourier transport coefficients have been obtained in steady-state conditions, namely when the cooling terms are compensated for by the energy gained by the solid particles due to their collisions with the bath particles (Garzó, Chamorro & Vega Reyes 2013; Gómez González & Garzó 2019). Thus, the results derived in the present paper must be consistent with those previously obtained by Gómez González & Garzó (2019) when the limit $m/m_g \to \infty$ is considered in our general results.

3. Homogeneous steady state

As a first step and before studying inhomogeneous states, we consider the HSS. The HSS is the reference base state (zeroth-order approximation) used in the Chapman–Enskog perturbation method (Chapman & Cowling 1970). Therefore, its investigation is of great importance. The HSS was widely analysed by Santos (2003) for a three-dimensional granular gas. Here, we extend these calculations to a general dimension $d$.

In the HSS, the density $n$ and temperature $T$ are spatially uniform, and with an appropriate selection of the frame reference, the mean flow velocities vanish ($\boldsymbol {U}=\boldsymbol {U}_g=\boldsymbol {0}$). Consequently, the Boltzmann equation (2.2) reads

(3.1)\begin{equation} \frac{\partial f(\boldsymbol{v};t)}{\partial t}=J[\boldsymbol{v}\,|\,f,f]+J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

Moreover, the velocity distribution $f(\boldsymbol {v};t)$ of the granular gas is isotropic in $\boldsymbol {v}$ so that the production of momentum $\boldsymbol {\mathcal {F}}[f]=\boldsymbol {0}$, according to (2.12). Thus, the only non-trivial balance equation is that of the temperature (2.10), namely $\partial _t \ln T=-(\zeta +\zeta _g)$. As mentioned in § 2, since collisions among granular particles are inelastic, the cooling rate $\zeta >0$. A collision between a particle of the granular gas and a particle of the molecular gas is elastic, and so the total kinetic energy of two colliding particles in such a collision is conserved. On the other hand, since in the steady state the background gas acts as a thermostat, the mean kinetic energy of granular particles is smaller than that of the molecular gas, and so $T< T_g$. This necessarily implies that $\zeta _g<0$. Therefore, in the steady state, the terms $\zeta$ and $|\zeta _g|$ exactly compensate each other and one gets the steady-state condition $\zeta +\zeta _g=0$. This condition allows one to get the steady granular temperature $T$. However, according to the definitions (2.13a,b), the determination of $\zeta$ and $\zeta _g$ requires knowledge of the velocity distribution $f(\boldsymbol {v})$. For inelastic collisions ($\alpha \neq 1$), to date the solution of the Boltzmann equation (3.1) has not been found. On the other hand, a good estimate of $\zeta$ and $\zeta _g$ can be obtained when the first Sonine approximation to $f$ is considered (Brilliantov & Pöschel 2004). In this approximation, $f(\boldsymbol {v})$ is given by

(3.2)\begin{equation} f(\boldsymbol{v})\simeq f_{{MB}}(\boldsymbol{v})\left\{1+\frac{a_2}{2} \left[\left(\frac{m v^{2}}{2T}\right)^{2}-(d+2)\frac{mv^{2}}{2T}+\frac{d(d+2)}{4}\right]\right\}, \end{equation}

where

(3.3)\begin{equation} f_{{MB}}(\boldsymbol{v})=n \biggl(\frac{m}{2{\rm \pi} T}\biggr)^{d/2} \exp \left(-\frac{m v^{2}}{2T}\right) \end{equation}

is the Maxwell–Boltzmann distribution and

(3.4)\begin{equation} a_2=\frac{1}{d(d+2)}\frac{m^{2}}{n T^{2}}\int {\rm d}\boldsymbol{v}v^{4} f(\boldsymbol{v})-1 \end{equation}

is the kurtosis or fourth cumulant. This quantity measures the departure of the distribution $f(\boldsymbol {v})$ from its Maxwellian form $f_{{MB}}(\boldsymbol {v})$. From experience with the dry granular case (van Noije & Ernst 1998; Garzó & Dufty 1999; Montanero & Santos 2000; Santos & Montanero 2009), the magnitude of the cumulant $a_2$ is expected to be very small, and so the Sonine approximation (3.2) to the distribution $f$ turns out to be reliable. In the case that $|a_2|$ does not remain small for high inelasticity, one should include cumulants of higher order in the Sonine polynomial expansion of $f$. However, the possible lack of convergence of the Sonine polynomial expansion for very small values of the coefficient of restitution (Brilliantov & Pöschel 2006a,b) puts in doubt the reliability of the Sonine expansion in the high-inelasticity region. Here, we restrict ourselves to values of $\alpha$ where $|a_2|$ remains relatively small.

The expressions of $\zeta$ and $\zeta _g$ can now be obtained by replacing in (2.13a,b) $f$ by its Sonine approximation (3.2). Retaining only linear terms in $a_2$, the forms of the dimensionless production rates $\zeta ^{*}=(\ell \zeta /v_{th})$ and $\zeta _g^{*}=(\ell \zeta /v_{th})$ can be written as (van Noije & Ernst 1998; Brilliantov & Pöschel 2006a)

(3.5a,b)\begin{equation} \zeta^{*}=\tilde{\zeta}^{(0)}+\tilde{\zeta}^{(1)} a_2, \quad \zeta_g^{*}=\tilde{\zeta}_g^{(0)}+\tilde{\zeta}_g^{(1)} a_2, \end{equation}

where

(3.6a,b)$$\begin{gather} \tilde{\zeta}^{(0)}=\frac{\sqrt{2}{\rm \pi}^{(d-1)/2}}{d\varGamma\left(\dfrac{d}{2}\right)}(1-\alpha^{2}),\quad \tilde{\zeta}^{(1)}=\frac{3}{16} \tilde{\zeta}^{(0)}, \end{gather}$$

(3.7a,b)$$\begin{gather}\tilde{\zeta}_g^{(0)}=2x(1-x^{2}) \left(\frac{\mu T}{T_g}\right)^{1/2}\gamma^{*},\quad \tilde{\zeta}_g^{(1)}=\frac{\mu_g}{8}x^{{-}3}[x^{2}\left(4-3\mu_g\right)-\mu_g] \left(\frac{\mu T}{T_g}\right)^{1/2}\gamma^{*}. \end{gather}$$

Here, $\ell =1/(n\sigma ^{d-1})$ is proportional to the mean free path of hard spheres, $v_{th}=\sqrt {2T/m}$ is the thermal velocity and we have introduced the auxiliary parameters

(3.8)\begin{equation} x=\left(\mu_g+\mu \frac{T_g}{T}\right)^{1/2} \end{equation}

and

(3.9a,b)\begin{equation} \gamma^{*}=\varepsilon \left(\frac{T_g}{T}\right)^{1/2},\quad \varepsilon=\frac{\ell \gamma}{\sqrt{2T_g/m}}= \frac{\sqrt{2}{\rm \pi}^{d/2}}{2^{d} d \varGamma\left(\dfrac{d}{2}\right)}\frac{1}{\phi \sqrt{T_g^{*}}}. \end{equation}

Here, $\phi =[{\rm \pi} ^{d/2}/2^{d-1}d\varGamma (d/2)]n\sigma ^{d}$ is the solid volume fraction and $T_g^{*}=T_g/(m\sigma ^{2} \gamma ^{2})$ is the (reduced) bath temperature. The dimensionless coefficient $\gamma ^{*}$ characterises the rate at which the collisions between grains and molecular particles occur. Equations (3.6a,b) and (3.7a,b) agree with those obtained by Santos (2003) for $d=3$.

To close the problem, we have to determine the kurtosis $a_2$. In this case, one has to compute the collisional moments

(3.10a,b)\begin{equation} \varLambda\equiv \int {\rm d}\boldsymbol{v}v^{4} J[\boldsymbol{v}\,|\,f,f],\quad \varLambda_g\equiv \int {\rm d}\boldsymbol{v} v^{4} J_g[\boldsymbol{v}\,|\,f,f_g]. \end{equation}

In the steady state, one has the additional condition $\varLambda +\varLambda _g=0$. The moments $\varLambda$ and $\varLambda _g$ have been obtained in previous works (van Noije & Ernst 1998; Brilliantov & Pöschel 2006a; Garzó, Vega Reyes & Montanero 2009; Garzó 2019) by replacing $f$ by its first Sonine form (3.2) and neglecting nonlinear terms in $a_2$. In terms of $\gamma ^{*}$, the expressions of $\{\varLambda ^{*},\varLambda _g^{*}\}=(\ell /(n v_{th}^{5}))\{\varLambda,\varLambda _g\}$ are given by

(3.11a,b)\begin{equation} \varLambda^{*}=\varLambda^{(0)}+\varLambda^{(1)} a_2,\quad \varLambda_g^{*}=\varLambda_g^{(0)}+\varLambda_g^{(1)} a_2, \end{equation}

where

(3.12)$$\begin{gather} \varLambda^{(0)}={-}\frac{{\rm \pi}^{(d-1)/2}}{\sqrt{2}\varGamma\left(\dfrac{d}{2}\right)} \left(d+\frac{3}{2}+\alpha^{2}\right)(1-\alpha^{2}), \end{gather}$$

(3.13)$$\begin{gather}\varLambda^{(1)}={-}\frac{{\rm \pi}^{(d-1)/2}}{\sqrt{2}\varGamma\left(\dfrac{d}{2}\right)} \left[\frac{3}{32}\left(10d+39+10\alpha^{2}\right)+\frac{d-1}{1-\alpha}\right](1-\alpha^{2}), \end{gather}$$

(3.14)$$\begin{gather}\varLambda_g^{(0)}={{d}x}^{{-}1}\left(x^{2}-1\right) [8 \mu_g x^{4}+x^{2}\left(d+2-8\mu_g\right)+\mu_g]\left(\frac{\mu T}{T_g}\right)^{1/2}\gamma^{*}, \end{gather}$$

(3.15)\begin{gather} \hspace{-2pc}\varLambda_g^{(1)} = \frac{d}{8} x^{{-}5}\{4x^{6}[30\mu_g^{3}-48\mu_g^{2}+3(d+8)\mu_g-2(d+2)]+ \mu_g x^{4} [{-}48\mu_g^{2} \nonumber\\ \hspace{4pc} +\,3(d+26)\mu_g-8(d+5)]+\mu_g^{2} x^{2}(d+14-9\mu_g)-3\mu_g^{3}\} \left(\frac{\mu T}{T_g}\right)^{1/2}\gamma^{*}. \end{gather}

For hard spheres ($d=3$), equations (3.12)–(3.15) are consistent with those previously obtained by Santos (2003).

Inserting (3.5a,b)–(3.7a,b) and (3.11a,b)–(3.15) into the steady-state conditions ($\zeta +\zeta _g=0$, $\varLambda +\varLambda _g=0$), one gets a set of coupled equations:

(3.16)$$\begin{gather} \tilde{\zeta}^{(0)}+\tilde{\zeta}_g^{(0)}+(\tilde{\zeta}^{(1)}+\tilde{\zeta}_g^{(1)})a_2=0, \end{gather}$$

(3.17)$$\begin{gather}\varLambda^{(0)}+\varLambda_g^{(0)}+(\varLambda^{(1)}+\varLambda_g^{(1)})a_2=0. \end{gather}$$

Eliminating $a_2$ in (3.16) and (3.17), one achieves the following closed equation for the temperature ratio $T/T_g$:

(3.18)\begin{equation} (\tilde{\zeta}^{(1)}+\tilde{\zeta}_g^{(1)})(\varLambda^{(0)}+\varLambda_g^{(0)})= (\tilde{\zeta}^{(0)}+\tilde{\zeta}_g^{(0)})(\varLambda^{(1)}+\varLambda_g^{(1)}). \end{equation}

For given values of $\alpha$, $\phi$ and $T_g^{*}$, the solution of (3.18) gives $T/T_g$. Once the temperature ratio is determined, the cumulant $a_2$ is simply given by

(3.19)\begin{equation} a_2={-}\frac{\tilde{\zeta}^{(0)}+\tilde{\zeta}_g^{(0)}}{\tilde{\zeta}^{(1)}+\tilde{\zeta}_g^{(1)}}={-}\frac{\varLambda^{(0)}+\varLambda_g^{(0)}}{\varLambda^{(1)}+\varLambda_g^{(1)}}. \end{equation}

The set of dimensionless control parameters of the problem considered here ($m/m_g$, $\alpha$, $\phi$, $T_g^{*}$) has been essentially chosen to perform a close and clean comparison with the previous results obtained by Gómez González & Garzó (2019) by means of the suspension model (2.16). On the other hand, another possible set of parameters are $m/m_g$, $\alpha$ and $\omega \equiv (n\sigma ^{d-1})/(n_g\overline {\sigma }^{d-1})$. The parameter $\omega$ represents the mean free path associated with the grain–gas collisions relative to that associated with grain–grain collisions. In fact, $\omega$ was considered independently by Biben et al. (2002) and Santos (2003) in their study of the HSS. The relationship between $\omega$ and $T_g^{*}$ is

(3.20)\begin{equation} \omega=\frac{2^{d+{3}/{2}}}{\sqrt{\rm \pi}}\left(\frac{m_g}{m}\right)^{1/2}\phi \sqrt{T_g^{*}}. \end{equation}

Note that $\omega$ encompasses the dependence on the volume fraction $\phi$ and the (reduced) bath temperature $T_g^{*}$ through the combination $\phi \sqrt {T_g^{*}}$. For this reason, according to (3.9a,b), the number of independent parameters in this set is $\alpha$, $m/m_g$ and $\omega$.

3.1. Brownian limit

Before illustrating the dependence of $T/T_g$ and $a_2$ on $\alpha$ for given values of $m/m_g$, $\phi$ and $T_g^{*}$, it is interesting to consider the Brownian limit $m/m_g\to \infty$. In this limiting case, $\mu _g\to 0$, $\mu \to 1$, $x\to \sqrt {T_g/T}$ and so

(3.21a,b)$$\begin{gather} \tilde{\zeta}_g^{(0)}\to 2\left(1-\frac{T_g}{T}\right)\gamma^{*}, \quad \tilde{\zeta}_g^{(1)}\to 0, \end{gather}$$

(3.22a,b)$$\begin{gather}\varLambda_g^{(0)}\to d(d+2)\left(\frac{T_g}{T}-1\right)\gamma^{*}, \quad \varLambda_g^{(1)}\to -d(d+2)\gamma^{*}. \end{gather}$$

Taking into account these results, the set of (3.16) and (3.17) can be written in the Brownian limit as

(3.23a,b)\begin{equation} 2\gamma^{*}\left(\frac{T_g}{T}-1\right)=\zeta^{*},\quad d(d+2)\left(\gamma^{*} a_2-\frac{1}{2}\zeta^{*}\right)=\varLambda^{*}. \end{equation}

These equations are the same as those derived by Gómez González & Garzó (2019) (see (29) and (34) of that paper) by using the suspension model (2.16). This shows the consistency of the present results in the HSS with those obtained in the Brownian limit.

3.2. Direct simulation Monte Carlo simulations

The previous analytical results have been obtained by using the first Sonine approximation (3.2) to $f$. Thus, it is worth solving the Boltzmann kinetic equation by means of an alternative method to test the reliability of the theoretical predictions for $T/T_g$ (3.18) and $a_2$ (3.19). The direct simulation Monte Carlo (DSMC) method developed by Bird (1994) is considered here to numerically solve the Boltzmann equation in the homogeneous state. Some technical details of the application of the DSMC method to the system studied in this paper are provided as supplementary material available at https://doi.org/10.1017/jfm.2022.410.

The dependence of the temperature ratio $\chi \equiv T/T_g$ on the coefficient of restitution $\alpha$ is plotted in figure 1 for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and several values of the mass ratio $m/m_g$. The value $T_g^{*}=1000$ has been chosen to guarantee that the grain–grain collisions play a relevant role in the dynamics of granular gas. Namely, that the value of the friction coefficient $\gamma$ is comparable to the value of the grain–grain collision frequency ($\nu =v_{th}/\ell$), so that the dimensionless coefficient $\gamma ^{*}\equiv \gamma /\nu$ is of the order of unity. This is fulfilled for the system studied here ($d=3$, $\phi =0.001$ and $T_g^{*}=1000$) since $\gamma ^{*}\simeq 10\sqrt {T_g/T}$, $T_g/T$ being of the order of unity (see figure 1). Thus, although the results displayed in this paper significantly differ from those found in the dry case (no gas phase), the effects of inelastic collisions still have importance for the dynamics of grains. The value $T_g^{*}=1000$ will therefore be maintained throughout this work.

Figure 1. Temperature ratio $\chi \equiv T/T_g$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=50, 10, 5$ and 1). The solid lines are the theoretical results obtained by numerically solving (3.18) and the symbols are the Monte Carlo simulation results. The dotted line is the result obtained by Gómez González & Garzó (2019) using the Langevin-like suspension model (2.16) while black circles refer to DSMC simulations implemented using the time-driven approach (see the supplementary material).

Theoretical results are compared against DSMC simulations in figure 1, which ensures the reliability of the results derived in this section for two different reasons: (i) a good agreement between theory and simulation is found and (ii) the convergence towards the Brownian limit can be clearly observed. Surprisingly, this convergence is fully reached for relatively small values of the mass ratio ($m/m_g\approx 50$). We also find that the departure of $\chi$ from unity increases as the masses of the granular and gas particles are comparable. However, this unexpected result is only due to the way of scaling the variables. This is illustrated in figure 2 where we take $\omega$ instead of $T_g^{*}$ as input as in figure 2 of Santos (2003). In contrast to figure 1, as expected (Barrat & Trizac 2002; Dahl et al. 2002), it is quite apparent that the lack of energy equipartition is more noticeable as $m\gg m_g$. In addition, we also see that the impact of the mass ratio on the temperature ratio is apparently more significant when one fixes $\omega$ instead of $T_g^{*}$. The fact that the difference $1-\chi$ increases with decreasing mass ratio when $T_g^{*}$ is fixed (see figure 1) can be easily understood. According to (2.6a,b), the transmission of energy per individual collision from a molecular particle to a grain is greater when their masses are similar. Nonetheless, the constraint imposed by the way of scaling $\gamma$ leads to a dependence of $N/N_g$ on the mass ratio $m/m_g$ for fixed $\sigma _g$. Thus, $N_g/N\propto m/m_g$, and so the number density of the molecular gas increases with increasing mass ratio. In this way, the mean force exerted by the molecular particles on the grains is greater, and therefore the thermalisation caused by the presence of the interstitial fluid is much more effective. The steady temperature ratio $\chi$ is reached when the energy lost by collisions is compensated for by the energy provided by the bath. Hence, the non-equipartition of energy turns out then to be remarkable to small values of $m/m_g$ and $\alpha$. Finally, figure 3 shows the $\alpha$ dependence of the cumulant $a_2$ for the same parameters as in figure 1. As expected, we find that the magnitude of $a_2$ is in general small for not quite large inelasticity (for instance, $\alpha \gtrsim 0.5$); this result supports the assumption of a low-order truncation (first Sonine approximation) in the polynomial expansion of the distribution function. Figure 3 highlights the excellent agreement between theory and simulations, except for $m/m_g=1$ where small differences are present for very strong inelasticities. However, these discrepancies are of the same order as those found for dry granular gases (Montanero & Garzó 2002).

Figure 2. Temperature ratio $\chi \equiv T/T_g$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $\omega =0.1$ and four different values of the mass ratio $m/m_g$: $m/m_g=1$ (solid line), $m/m_g=10$ (dashed line), $m/m_g=100$ (dotted line) and $m/m_g=1000$ (dash-dotted line). The (reduced) bath temperature $T_g^{*}=61.36 (m/m_g)$.

Figure 3. Plot of the fourth cumulant $a_2$ as a function of the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=1, 5, 10$ and 50). The solid lines are the theoretical results obtained from (3.18) and the symbols are the Monte Carlo simulation results. The dotted line is the result obtained by Gómez González & Garzó (2019) using the Langevin-like suspension model (2.16) while black circles refer to DSMC simulations implemented using the time-driven approach (see the supplementary material).

4. Chapman–Enskog expansion. First-order approximation

We perturb now the homogeneous state by small spatial gradients. These perturbations will give non-zero contributions to the pressure tensor and the heat flux vector. The determination of these fluxes will allow us to identify the Navier–Stokes–Fourier transport coefficients of the granular gas. For times longer than the mean free time, we assume that the system evolves towards a hydrodynamic regime where the distribution function $f(\boldsymbol {r}, \boldsymbol {v};t)$ adopts the form of a normal or hydrodynamic solution. This means that all space and time dependence of $f$ only occurs through the hydrodynamic fields $n$, $\boldsymbol {U}$ and $T$:

(4.1)\begin{equation} f(\boldsymbol{r},\boldsymbol{v},t)=f\left[\boldsymbol{v}\,|\,n (t), \boldsymbol{U}(t), T(t) \right]. \end{equation}

The notation on the right-hand side indicates a functional dependence on the density, flow velocity and temperature. For low Knudsen numbers (i.e. small spatial variations), the functional dependence (4.1) can be made local in space by means of an expansion in powers of the gradients $\boldsymbol {\nabla } n$, $\boldsymbol {\nabla } \boldsymbol {U}$ and $\boldsymbol {\nabla } T$ (Chapman & Cowling 1970). In this case, $f$ can be expressed in the form $f=f^{(0)}+f^{(1)}+f^{(2)}+\cdots$, where the approximation $f^{(k)}$ is of order $k$ in spatial gradients. Here, since we are interested in the Navier–Stokes hydrodynamic equations, only terms up to first order in gradients are considered in the constitutive equations for the momentum and heat fluxes.

As has been noted in previous works on granular suspensions (Garzó et al. 2013; Gómez González & Garzó 2019), although one is interested in computing transport in steady conditions, the presence of the background molecular gas may induce a local energy unbalance between the energy lost due to inelastic collisions and the energy transfer via elastic collisions. Thus, we have to consider first the time-dependent distribution $f^{(0)}(\boldsymbol {r},\boldsymbol {v}; t)$ in order to arrive at the linear integral equations obeying the Navier–Stokes–Fourier transport coefficients. Then, to get explicit forms for the transport coefficients, the above integral equations are (approximately) solved under steady-state conditions.

To collect the different approximations in (2.2), one has to characterise the magnitude of the velocity difference ${\rm \Delta} \boldsymbol {U}$ relative to the gradients as well. In the absence of spatial gradients, the production of momentum term $\mathcal {F}[f]\propto {\rm \Delta} \boldsymbol {U}$ (see (6.5)) and hence, according to the momentum balance equation (2.9), the mean flow velocity $\boldsymbol {U}$ of the granular gas relaxes towards that of the molecular gas $\boldsymbol {U}_g$ after a transient period. Thus, the term ${\rm \Delta} \boldsymbol {U}$ must be considered to be at least of first order in the spatial gradients. In this case, the Maxwellian distribution $f_g(\boldsymbol {v})$ must also be expanded as

(4.2)\begin{equation} f_g(\boldsymbol{v})=f_g^{(0)}(\boldsymbol{V})+f_g^{(1)}(\boldsymbol{V})+\cdots, \end{equation}

where

(4.3)\begin{equation} f_g^{(0)}(\boldsymbol{V})=n_g \left(\frac{m_g}{2{\rm \pi} T_g}\right)^{d/2} \exp \bigg(-\frac{m_g V^{2}}{2T_g}\bigg) \end{equation}

and

(4.4)\begin{equation} f_g^{(1)}(\boldsymbol{V})={-}\frac{m_g}{T_g}\boldsymbol{V}\boldsymbol{\cdot} {\rm \Delta} \boldsymbol{U} f_g^{(0)}(\boldsymbol{V}). \end{equation}

According to the expansion (4.1), the pressure tensor $P_{ij}$, the heat flux $\boldsymbol {q}$ and the partial production rates $\zeta$ and $\zeta _g$ must also be expressed according to the perturbation scheme in the forms

(4.5a–d)\begin{equation} \left.\begin{gathered} P_{ij}=P_{ij}^{(0)}+P_{ij}^{(1)}+\cdots,\quad \boldsymbol{q}=\boldsymbol{q}^{(0)}+\boldsymbol{q}^{(1)}+\cdots,\\ \zeta=\zeta^{(0)}+\zeta^{(1)}+\cdots,\quad \zeta_g=\zeta_g^{(0)}+\zeta_g^{(1)}+\cdots . \end{gathered}\right\} \end{equation}

In addition, the time derivative $\partial _t$ is also given as $\partial _t=\partial _t^{(0)}+\partial _t^{(1)}+\cdots$. The action of the operators $\partial _t^{(k)}$ on the hydrodynamic fields can be identified when the expansions (4.5a–d) of the fluxes and the production rates are considered in the macroscopic balance equations (2.8)–(2.10). This is the conventional Chapman–Enskog method (Chapman & Cowling 1970; Garzó 2019) for solving the Boltzmann kinetic equation.

As usual in the Chapman–Enskog method (Chapman & Cowling 1970), the zeroth-order distribution function $f^{(0)}$ defines the hydrodynamic fields $n$, $\boldsymbol {U}$ and $T$:

(4.6)\begin{equation} \left\{n, n \boldsymbol{U},d n T\right\}=\int {\rm d}\boldsymbol{v}\{1,\boldsymbol{v}, m V^{2}\}f^{(0)}(\boldsymbol{V}). \end{equation}

The requirements (4.6) must be fulfilled at any order in the expansion, and so the distributions $f^{(k)}$ ($k\geq 1$) must thus obey the orthogonality conditions

(4.7)\begin{equation} \int {\rm d}\boldsymbol{v}\{1,\boldsymbol{v}, m V^{2}\}f^{(k)}(\boldsymbol{V})=\{0,\boldsymbol{0},0\}. \end{equation}

These are the usual solubility conditions of the Chapman–Enskog scheme.

The mathematical steps involved in the determination of the zeroth- and first-order distributions are quite similar to those made in previous works (Brey et al. 1998; Garzó & Dufty 1999; Garzó et al. 2012, 2013; Gómez González & Garzó 2019). Some of the technical details involved in this derivation are given in the supplementary material.

4.1. Navier–Stokes transport coefficients

To first order in spatial gradients and based on symmetry considerations, the pressure tensor $P_{ij}^{(1)}$ and the heat flux $\boldsymbol {q}^{(1)}$ are given, respectively, by

(4.8)$$\begin{gather} P_{ij}^{(1)}={-}\eta \left(\frac{\partial U_i}{\partial r_j}+ \frac{\partial U_j}{\partial r_i}-\frac{2}{d}\delta_{ij}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{U}\right), \end{gather}$$

(4.9)$$\begin{gather}\boldsymbol{q}^{(1)}={-}\kappa \boldsymbol{\nabla} T-\bar{\mu} \boldsymbol{\nabla} n-\kappa_U {\rm \Delta} \boldsymbol{U}. \end{gather}$$

Here, $\eta$ is the shear viscosity, $\kappa$ is the thermal conductivity, $\bar {\mu }$ is the diffusive heat conductivity and $\kappa _U$ is the velocity conductivity. While $\eta$, $\kappa$ and $\mu$ are the coefficients of proportionality between fluxes and hydrodynamic gradients, the coefficient $\kappa _U$ connects the heat flux with the velocity difference ${\rm \Delta} \boldsymbol {U}$ (‘convection’ current). Although this contribution is not present in dry granular gases (Garzó 2019), it also appears in the case of driven granular mixtures (Khalil & Garzó 2013, 2018). The coefficient $\kappa _U$ can be seen as a measure of the contribution to the heat flow due to ‘diffusion’ (in the sense that we have a ‘binary mixture’ of granular and gas particles where both species have different mean velocities). In this context, $\kappa _U$ can be regarded as an effect inverse to thermal diffusion (diffusion thermo-effect) (Chapman & Cowling 1970). It is important to recall that ${\rm \Delta} \boldsymbol {U}$ vanishes for HSSs, and hence one expects that ${\rm \Delta} \boldsymbol {U}$ can be expressed in terms of $\boldsymbol {\nabla } T$ and $\boldsymbol {\nabla } n$ in particular inhomogeneous situations.

The Navier–Stokes–Fourier transport coefficients $\eta$, $\kappa$ and $\mu$ are defined, respectively, as

(4.10)$$\begin{gather} \eta={-}\frac{1}{(d-1)(d+2)}\int {\rm d}\boldsymbol{v} R_{ij}(\boldsymbol{V})\mathcal{C}_{ij}(\boldsymbol{V}), \end{gather}$$

(4.11a,b)$$\begin{gather}\kappa={-}\frac{1}{dT}\int {\rm d}\boldsymbol{v}\boldsymbol{S}(\boldsymbol{V})\boldsymbol{\cdot} \boldsymbol{\mathcal{A}}(\boldsymbol{V}),\quad \bar{\mu}={-}\frac{1}{dn}\int {\rm d}\boldsymbol{v}\boldsymbol{S}(\boldsymbol{V})\boldsymbol{\cdot} \boldsymbol{\mathcal{B}}(\boldsymbol{V}), \end{gather}$$

while $\kappa _U$ is

(4.12)\begin{equation} \kappa_U={-}\frac{1}{d}\int {\rm d}\boldsymbol{v}\boldsymbol{S}(\boldsymbol{V})\boldsymbol{\cdot} \boldsymbol{\mathcal{E}}(\boldsymbol{V}). \end{equation}

In (4.10)–(4.12), we have introduced the quantities

(4.13a,b)\begin{equation} R_{ij}(\boldsymbol{V})=m\left(V_i V_j-\frac{1}{d}V^{2} \delta_{ij}\right),\quad \boldsymbol{S}(\boldsymbol{V})=\left(\frac{m}{2}V^{2}-\frac{d+2}{2}T\right)\boldsymbol{V}. \end{equation}

5. Sonine polynomial approximation to the transport coefficients in steady-state conditions

So far, all the results displayed in § 4 for the transport coefficients $\eta$, $\kappa$, $\bar {\mu }$ and $\kappa _U$ are exact. More specifically, their expressions are given by (4.10)–(4.12), where the unknowns $\boldsymbol {\mathcal {A}}$, $\boldsymbol {\mathcal {B}}$, $\mathcal {C}_{ij}$, $\mathcal {D}$ and $\boldsymbol {\mathcal {E}}$ are the solutions of a set of coupled linear integral equations displayed in the supplementary material. However, it is easy to see that the solution for general unsteady conditions requires one to solve numerically a set of coupled differential equations for $\eta$, $\kappa$, $\bar {\mu }$ and $\kappa _U$. Thus, in a desire of achieving analytical expressions of the transport coefficients, we consider steady-state conditions. In this case, the constraint $\zeta ^{(0)}+\zeta _g^{(0)}=0$ applies locally, and so the transport coefficients can be explicitly obtained. The procedure for deriving the expressions of the transport coefficients is described in the supplementary material and only their final forms are provided here.

5.1. Shear viscosity

The shear viscosity coefficient $\eta$ is given by

(5.1)\begin{equation} \eta=\frac{\eta_0}{\nu_\eta^{*}+K'\tilde{\nu}_\eta\gamma^{*}}, \end{equation}

where $\gamma ^{*}$ is defined in (3.9a,b), $\eta _0=[(d+2)\varGamma (d/2)/(8{\rm \pi} ^{(d-1)/2})]\sigma ^{1-d}\sqrt {m T}$ is the low-density value of the shear viscosity of an ordinary gas of hard spheres ($\alpha =1$) and $K'= \sqrt {2}(d+2)\varGamma (d/2)/(8{\rm \pi} ^{(d-1)/2})$. Moreover, we have introduced the (reduced) collision frequencies

(5.2)\begin{equation} \nu_\eta^{*}=\frac{3}{4d}\left(1-\alpha+\frac{2}{3}d\right)(1+\alpha), \end{equation}

(5.3)\begin{align} \tilde{\nu}_\eta &= \frac{1}{(d-1)(d+2)}\left(\frac{m}{m_g}\right)^{3} \mu_g \left(\frac{T_g}{T}\right)^{2}\theta^{{-}1/2} \nonumber\\ &\quad \times\{2(d+3)(d-1)\left(\mu-\mu_g \theta\right)\theta^{{-}2}(1+\theta)^{{-}1/2} \nonumber\\ &\quad +2d(d-1)\mu_g \theta^{{-}2}(1+\theta)^{1/2}+2(d+2)(d-1)\theta^{{-}1}(1+\theta)^{{-}1/2}\}, \end{align}

where $\theta =m T_g/(m_g T)$ is the ratio of the mean square velocities of granular and molecular gas particles. It is important to recall that all the quantities appearing in (5.1) are evaluated at the steady-state conditions.

5.2. Thermal conductivity, diffusive heat conductivity and velocity conductivity

We consider here the transport coefficients associated with the heat flux. The thermal conductivity coefficient $\kappa$ is

(5.4)\begin{equation} \kappa=\frac{d-1}{d}\frac{\kappa_0}{\nu_\kappa^{*}+K'\left(\tilde{\nu}_\kappa+\beta\right)\gamma^{*}}, \end{equation}

where $\kappa _0=[d(d+2)/2(d-1)](\eta _0/m)$ is the low-density value of the thermal conductivity for an ordinary gas of hard spheres and

(5.5)\begin{equation} \beta=\left(x^{{-}1}-3x\right)\mu^{3/2}\left(\frac{T_g}{T}\right)^{1/2}. \end{equation}

In (5.4), we have introduced the (reduced) collision frequencies

(5.6)$$\begin{gather} \nu_\kappa^{*}=\frac{1+\alpha}{d}\left[\frac{d-1}{2}+\frac{3}{16}(d+8)(1-\alpha)\right](1+\alpha), \end{gather}$$

(5.7)$$\begin{gather}\tilde{\nu}_\kappa=\frac{1}{2(d+2)}\mu \frac{\theta}{1+\theta} \left[G-(d+2)\frac{1+\theta}{\theta}F\right], \end{gather}$$

where

(5.8)\begin{align} F &= (d+2)(2 \delta +1)+4(d-1)\mu_g \delta \theta^{{-}1}(1+\theta)+3(d+3)\delta^{2}\theta^{{-}1} \nonumber\\ &\quad +(d+3)\mu_g^{2} \theta^{{-}1}(1+\theta)^{2} -(d+2)\theta^{{-}1}(1+\theta), \end{align}

(5.9)\begin{align} G &= (d+3)\mu_g^{2} \theta^{{-}2}(1+\theta)^{2}\left[d+5+(d+2)\theta\right] \nonumber\\ &\quad -\mu_g(1+\theta)\{4(1-d)\delta \theta^{{-}2}\left[d+5+(d+2)\theta\right]-8(d-1)\theta^{{-}1}\} \nonumber\\ &\quad +3(d+3)\delta^{2}\theta^{{-}2}\left[d+5+(d+2)\theta\right]+2\delta\theta^{{-}1} [24+11d+d^{2}+(d+2)^{2}\theta] \nonumber\\ &\quad +(d+2)\theta^{{-}1}\left[d+3+(d+8)\theta\right]-(d+2)\theta^{{-}2}(1+\theta)\left[d+3+(d+2)\theta\right] \end{align}

and $\delta \equiv \mu -\mu _g \theta$.

The diffusive heat conductivity $\bar {\mu }$ can be written as

(5.10)\begin{equation} \bar{\mu}=\frac{K' T}{n}\frac{\kappa \zeta^{*}}{\nu_\kappa^{*}+K'\tilde{\nu}_\kappa \gamma^{*}}. \end{equation}

Finally, the velocity conductivity $\kappa _U$ is given by

(5.11)\begin{equation} \kappa_U={-}\frac{n T}{2}\frac{K' \mu (1+\theta)^{{-}1/2}\theta^{{-}1/2}H}{\nu_\kappa^{*}+K'\tilde{\nu}_\kappa \gamma^{*}}\gamma^{*}, \end{equation}

where

(5.12)\begin{equation} H=(d+2)(1+2\delta)+4(1-d)\mu_g (1+\theta) \delta-3(d+3)\delta^{2}-(d+3)\mu_g^{2} (1+\theta)^{2}. \end{equation}

5.3. Brownian limit

Equations (5.1), (5.4), (5.10) and (5.11) provide the expressions of the transport coefficients $\eta$, $\kappa$, $\bar {\mu }$ and $\kappa _U$, respectively, for arbitrary values of the mass ratio $m/m_g$. As before in the homogeneous state, it is quite interesting to consider the limiting case $m/m_g \to \infty$ (Brownian limit). In this limit case, $\mu _g\to 0$, $\mu \to 1$, $\chi \equiv \textrm {finite}$, and so $\theta \to 0$, $x\to \chi ^{-1/2}$, $\delta \to 1-\chi ^{-1}$ and $\beta \to 1-3\chi ^{-1}$. This yields the results $\tilde {\nu }_\eta \to 2$ and $\tilde {\nu }_\kappa \to 3$, so that in the Brownian limit (5.1), (5.4), (5.10) and (5.11) reduce to (note that there is a typographical error in expression (78) of Gómez González & Garzó (2019) for the coefficient $\bar {\mu }$ since the denominator should be $\nu _\kappa ^{*}+3K'\gamma ^{*}$)

(5.13a,b)$$\begin{gather} \eta\to \frac{\eta_0}{\nu_\eta^{*}+2K'\gamma^{*}},\quad \kappa\to \frac{d-1}{d}\frac{\kappa_0}{\nu_\kappa^{*}+K'\left(\gamma^{*}- \frac{3}{2}\zeta^{*}\right)}, \end{gather}$$

(5.14a,b)$$\begin{gather}\bar{\mu}\to \frac{\kappa T}{n}\frac{K'\zeta^{*}}{\nu_\kappa^{*}+3K' \gamma^{*}},\quad \kappa_U \to 0. \end{gather}$$

Equations (5.13a,b) and (5.14a,b) agree with the results obtained by Gómez González & Garzó (2019) by using the suspension model (2.16). This confirms the self-consistency of the results obtained in this paper for general values of the mass ratio.

5.4. Some illustrative systems

In the steady state, the expressions of the Navier–Stokes–Fourier transport coefficients $\eta$, $\kappa$, $\bar {\mu }$ and $\kappa _U$ are provided by (5.1), (5.4), (5.10) and (5.11), respectively. As in previous works on transport in granular gases (Brey et al. 1998; Garzó & Dufty 1999; Garzó et al. 2012; Gómez González & Garzó 2019), to highlight the $\alpha$ dependence of the transport coefficients, they are scaled with respect to their values for elastic collisions. This scaling cannot be made in the case of the diffusive heat conductivity $\bar {\mu }$ since this coefficient vanishes for $\alpha =1$. In this case, we consider the scaled coefficient $n\bar {\mu }/T\kappa (1)$, where $\kappa (1)$ refers to the value of the thermal conductivity (5.4) for elastic collisions. All these scaled coefficients exhibit a complex dependence on the coefficient of restitution $\alpha$, the mass ratio $m/m_g$, the volume fraction $\phi$ (through the parameter $\varepsilon$ defined by (3.9a,b)) and the reduced temperature $T_g^{*}$ of the molecular gas.

Figures 4–7 show $\eta (\alpha )/\eta (1)$, $\kappa (\alpha )/\kappa (1)$, $n\bar {\mu }/T\kappa (1)$ and $\kappa _U(\alpha )/\kappa _U(1)$, respectively, as functions of the coefficient of restitution $\alpha$. Here, $\eta(1)$, $\kappa(1)$ and $\kappa_U(1)$ correspond to the values of $\eta$, $\kappa$ and $\kappa_U$ for elastic collisions. Moreover, in those plots we consider a three-dimensional system ($d=3$) with $\phi =0.001$ (very dilute granular gas), $T_g^{*}=1000$ and four different values of the mass ratio: $m/m_g=1, 5, 10$ and 50. We have also plotted the results obtained by Gómez González & Garzó (2019) using the suspension model (2.16). The results obtained from this model are expected to apply when $m \gg m_g$ (Brownian limit).

Figure 4. Plot of the (scaled) shear viscosity coefficient $\eta (\alpha )/\eta (1)$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=50, 10, 5$ and 1). The solid lines are the results derived in this paper while the dotted line is the result obtained by Gómez González & Garzó (2019) using the suspension model (2.16). Here, $\eta (1)$ refers to the shear viscosity coefficient when collisions between grains are elastic ($\alpha =1$).

Figure 5. Plot of the (scaled) thermal conductivity coefficient $\kappa (\alpha )/\kappa (1)$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=1, 5, 10$ and 50). The solid lines are the results derived in this paper while the dotted line is the result obtained by Gómez González & Garzó (2019) using the suspension model (2.16). Here, $\kappa (1)$ refers to the thermal conductivity coefficient when collisions between grains are elastic ($\alpha =1$).

Figure 6. Plot of the (scaled) diffusive heat conductivity coefficient $n\bar {\mu }(\alpha )/T\kappa (1)$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=1, 5, 10$ and 50). The solid lines are the results derived in this paper while the dotted line is the result obtained by Gómez González & Garzó (2019) using the suspension model (2.16). Here, $\kappa (1)$ refers to the thermal conductivity coefficient when collisions between grains are elastic ($\alpha =1$).

Figure 7. Plot of the (scaled) velocity conductivity coefficient $\kappa _U(\alpha )/\kappa _U(1)$ versus the coefficient of normal restitution $\alpha$ for $d=3$, $\phi =0.001$, $T_g^{*}=1000$ and four different values of the mass ratio $m/m_g$ (from top to bottom: $m/m_g=1, 5, 10$ and 50). Here, $\kappa _U(1)$ refers to the velocity conductivity coefficient when collisions between grains are elastic ($\alpha =1$).

We observe that the deviations of the transport coefficients from their elastic forms are in general significant, especially when $m=m_g$. While the (scaled) shear viscosity and thermal conductivity coefficients exhibit a non-monotonic dependence on inelasticity, the (scaled) heat diffusive and velocity conductivity coefficients increase with increasing inelasticity, regardless of the value of the mass ratio considered. In addition, while $\eta (\alpha )<\eta (1)$, the opposite happens for the thermal conductivity since $\kappa (\alpha )>\kappa (1)$. With respect to the dependence on the mass ratio $m/m_g$, at a fixed value of the coefficient of restitution, it is quite apparent that while the (scaled) shear viscosity increases with increasing mass ratio, the (scaled) thermal conductivity decreases with increasing mass ratio. The same happens for the (scaled) coefficients $n\bar {\mu }/T\kappa (1)$ and $\kappa _U(\alpha )/\kappa _U(1)$ since both scaled coefficients decrease as the mass ratio increases. We also see that in the case $m/m_g=50$, the results derived here for $\eta (\alpha )/\eta (1)$, $\kappa (\alpha )/\kappa (1)$ and $n\bar {\mu }/T\kappa (1)$ practically coincide with those obtained in the Brownian limit by Gómez González & Garzó (2019). However, in the case $m/m_g=50$, the (scaled) velocity conductivity coefficient $\kappa _U(\alpha )/\kappa _U(1)$ ($\kappa _U=0$ for any value of $\alpha$ in the Brownian limit) is still clearly different from zero.

Although the results obtained here for the (scaled) transport coefficients depend on the values of the mass ratio and the (reduced) temperature of the molecular gas, it is worthwhile comparing the present results with those obtained for dry granular gases (i.e. in the absence of the molecular gas). In the case of the shear viscosity, a comparison between both systems (with and without the gas phase) shows significant discrepancies (see, for instance, figure 3.1 of Garzó 2019) even at a qualitative level since while $\eta$ increases with inelasticity for dry granular gases, the opposite happens here whatever the mass ratio considered. On the other hand, a more qualitative agreement is found for the the thermal conductivity (see, for instance, figure 3.2 of Garzó 2019) since $\kappa$ increases with decreasing $\alpha$ in both systems. In any case, important quantitative differences appear at strong dissipation since the influence of inelasticity on $\kappa$ is more relevant in the dry case than in the presence of the molecular gas. A similar conclusion is reached for the heat diffusive coefficient $\bar {\mu }$ (see, for instance, figure 3.3 of Garzó 2019) where the magnitude of this (scaled) coefficient for dry granular gases is much larger than that found here for granular suspensions.

6. Linear stability analysis of the HSS

Once the transport coefficients of the granular gas are known, the corresponding Navier–Stokes hydrodynamic equations can be explicitly displayed. To derive them, one has to take into account first that the the production of momentum $\boldsymbol {\mathcal {F}}[f]$ (defined in (2.12)) to first order in spatial gradients can be written as

(6.1)\begin{align} \boldsymbol{\mathcal{F}}^{(1)}[f^{(1)}] &= \int {\rm d}\boldsymbol{v} m\boldsymbol{V} \{J_g[f^{(0)},f_g^{(1)}]+J_g[f^{(1)},f_g^{(0)}]\} \nonumber\\ &={-}\xi {\rm \Delta} \boldsymbol{U}+\int {\rm d}\boldsymbol{v} m\boldsymbol{V}J_g[f^{(1)},f_g^{(0)}], \end{align}

where

(6.2)\begin{equation} \xi=\rho \mu \theta^{{-}1/2}(1+\theta)^{1/2}\gamma. \end{equation}

By symmetry reasons, the contributions of the first-order distribution $f^{(1)}(\boldsymbol {V})$ to the second integral in (6.1) come from the terms $\boldsymbol {\mathcal {A}}(\boldsymbol {V})$, $\boldsymbol {\mathcal {B}}(\boldsymbol {V})$ and $\boldsymbol {\mathcal {E}}(\boldsymbol {V})$. In the leading Sonine approximation (see the supplementary material), one gets

(6.3)\begin{align} \int {\rm d}\boldsymbol{v} m\boldsymbol{V}J_g[f^{(1)},f_g^{(0)}] &={-}\frac{2}{d(d+2)} \left(\frac{m}{nT^{2}}\kappa \boldsymbol{\nabla} \ln T+\frac{m}{T^{3}}\bar{\mu} \boldsymbol{\nabla} \ln n +\frac{m}{nT^{3}}\kappa_U {\rm \Delta} \boldsymbol{U}\right) \nonumber\\ &\quad \times \int {\rm d}\boldsymbol{v} m \boldsymbol{V} \boldsymbol{\cdot} J_g[\boldsymbol{S}f^{(0)},f_g^{(0)}]. \end{align}

The collision integral appearing in (6.3) can be computed with the result

(6.4)\begin{equation} \int {\rm d}\boldsymbol{v} m \boldsymbol{V}\boldsymbol{\cdot} J_g[\boldsymbol{S}f^{(0)},f_g^{(0)}]={-}\frac{d}{2}n T^{2} \mu \gamma \theta^{{-}1/2}\left(1+\theta\right)^{{-}1/2}. \end{equation}

Substitution of (6.4) into (6.3) leads to the final expression of $\boldsymbol {\mathcal {F}}^{(1)}[f^{(1)}]$:

(6.5)\begin{equation} \boldsymbol{\mathcal{F}}^{(1)}[f^{(1)}]={-}\xi {\rm \Delta} \boldsymbol{U}+\frac{\rho}{d+2}\mu\gamma\left(\frac{\kappa}{n}\boldsymbol{\nabla} \ln T+\frac{\bar{\mu}}{T}\boldsymbol{\nabla} \ln n+\frac{\kappa_U}{nT}{\rm \Delta} \boldsymbol{U}\right)X, \end{equation}

where

(6.6)\begin{equation} X(\theta)=\theta^{{-}1/2}\left(1+\theta\right)^{{-}1/2}. \end{equation}

Thus, when the constitutive equations (4.8)–(4.9) and (6.5) are substituted into the (exact) balance equations (2.8)–(2.10), one gets the Navier–Stokes hydrodynamic equations for a granular gas immersed in a molecular gas:

(6.7)\begin{gather} D_t n+n\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}=0, \end{gather}

(6.8)\begin{gather} \hspace{-3pc}\rho D_t U_i+\frac{\partial p}{\partial r_i} = \frac{\partial}{\partial r_j} \left[\eta\left(\frac{\partial U_j}{\partial r_i}+\frac{\partial U_i}{\partial r_j}-\frac{2}{d}\delta_{ij} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}\right)\right]-\xi {\rm \Delta} U_i \nonumber\\ \hspace{5pc}+\,\frac{\rho}{d+2}\mu\gamma X\left(\frac{\kappa}{n}\frac{\partial \ln T}{\partial r_i}+ \frac{\bar{\mu}}{T}\frac{\partial \ln n}{\partial r_i}+\frac{\kappa_U}{nT}{\rm \Delta} U_i\right), \end{gather}

(6.9)\begin{gather} \hspace{-8.9pc}D_tT+T(\zeta^{(0)}+\zeta_g^{(0)}) = \frac{2}{dn}\boldsymbol{\nabla}\boldsymbol{\cdot} \left(\kappa\boldsymbol{\nabla} T+\bar{\mu}\boldsymbol{\nabla} n+\kappa_U {\rm \Delta} \boldsymbol{U}\right) \nonumber\\ \hspace{10pc} +\,\frac{2}{dn}\left[\eta\left(\frac{\partial U_j}{\partial r_i}+ \frac{\partial U_i}{\partial r_j}-\frac{2}{d}\delta_{ij} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{U}\right) \frac{\partial U_i}{\partial r_j}-\frac{2}{d}T\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{U}\right]. \end{gather}

As said in § 5, we have not considered in (6.9) the first-order contributions to $\zeta$ and $\zeta _g$ since they vanish when non-Gaussian corrections to the distribution function $f^{(0)}$ are neglected. In addition, as already mentioned in several previous works (Garzó 2005; Garzó, Montanero & Dufty 2006), the above production rates should also include second-order contributions in spatial gradients. However, in the case of a dry dilute granular gas (Brey et al. 1998), it has been shown that these contributions are very small, and hence they can be neglected in the hydrodynamic equations. We expect here that the same happens for a granular suspension. Apart from the above approximations, the Navier–Stokes hydrodynamic equations (6.7)–(6.9) are exact to second order in the spatial gradients of $n$, $\boldsymbol {U}$ and $T$.

A simple solution of (6.7)–(6.9) corresponds to the HSS studied in § 2. A natural question is whether actually the HSS may be unstable with respect to long enough wavelength perturbations, as occurs for dry granular fluids in freely cooling flows (Goldhirsch & Zanetti 1993; McNamara 1993). This is one of the most characteristic features of granular gases; its origin is associated with the inelasticity of collisions. On the other hand, the stability of the HSS was also analysed by Gómez González & Garzó (2019) in the Brownian limit case showing that the HSS is always linearly stable. The objective of this section is to check whether the HSS is still linearly stable for arbitrary values of the mass ratio $m/m_g$.

Given that the present analysis is quite similar to that previously made by Gómez González & Garzó (2019), only the final results are provided here. Some mathematical steps are displayed in the supplementary material. In Fourier space, as expected (Brey et al. 1998; Garzó 2005), the $d-1$ transverse velocity components $\boldsymbol {w}_{\boldsymbol {k}\perp }=\boldsymbol {w}_{\boldsymbol {k}}-(\boldsymbol {w}_{\boldsymbol {k}} \boldsymbol {\cdot }\hat {\boldsymbol {k}})\hat {\boldsymbol {k}}$ (orthogonal to the wave vector $\boldsymbol {k}$) decouple from the other three modes. In terms of the dimensionless time $\tau$, their evolution equation is $\partial _\tau \boldsymbol {w}_{\boldsymbol {k}\perp }- \lambda _\perp (k)\boldsymbol {w}_{\boldsymbol {k}\perp }=0$. The explicit form of $\lambda _\perp (k)$ can be found in the supplementary material. A systematic analysis of the dependence of $\lambda _\perp (k)$ on the parameter space of the system shows that $\lambda _\perp (k)$ is always negative, and hence the transversal shear modes $\boldsymbol {w}_{\boldsymbol {k}\perp }(\boldsymbol {k},\tau )$ are linearly stable.

A careful analysis of the corresponding eigenvalues $\lambda _\ell (k)$ of the matrix $\boldsymbol{\mathsf{M}}$ obeying the remaining three longitudinal modes ($\rho _{\boldsymbol {k}}$, $\theta _{\boldsymbol {k}}$ and the longitudinal velocity component of the velocity field, $w_{\boldsymbol {k}\parallel }=\boldsymbol {w}_{\boldsymbol {k}}\boldsymbol {\cdot }\hat {\boldsymbol {k}}$) shows that one of the longitudinal modes could be unstable for values of the wavenumber $k< k_{\textrm {h}}$. The expression of $k_{\textrm {h}}$ is

(6.10)\begin{equation} k_{{h}}^{2}=\frac{\sqrt{2}(2\overline{\zeta}_g\gamma^{*}-\zeta^{*})+ \dfrac{2d}{d+2}\mu X \gamma^{*}[2\left(\zeta^{*}D_T^{*}-\overline{\zeta}_g \overline{\mu}^{*}\gamma^{*}\right)-\zeta^{*}\overline{\mu}^{*}]}{\overline{\mu}^{*}-D_T^{*}}, \end{equation}

where all the quantities appearing in (6.10) are defined in the supplementary material. As in the case of $\lambda _\perp (k)$, a study of the dependence of $k_{{h}}$ on the parameters of the system shows that $k_{{h}}^{2}$ is always negative. Consequently, there are no physical values of the wavenumber for which the longitudinal modes become unstable, and hence the longitudinal modes are also linearly stable.

In summary, the linear stability analysis of the HSS carried out here for a dilute granular gas surrounded by a molecular gas shows no surprises relative to the earlier study performed in the Brownian limit: the HSS is linearly stable for arbitrary values of the mass ratio $m/m_g$. However, the dispersion relations defining the dependence of the eigenvalues $\lambda _\perp (k)$ and $\lambda _{||}(k)$ on the parameter space are very different from those previously obtained when $m/m_g\to \infty$ (Gómez González & Garzó 2019). As an illustration, figure 8 shows the real parts of the eigenvalues $\lambda _{i,||}$ ($i=1,2,3$) and $\lambda _{\perp }$ as a function of the wavenumber $k$ for $\phi =0.001$, $T_g^{*}=1000$, $m/m_g=1$ and $\alpha =0.8$. It is quite apparent that all the eigenvalues are negative, as expected. In particular, although the longitudinal mode $\lambda _{1,||}$ is quite close to 0, the inset clearly shows that it is always negative. In addition, we also observe that in general the eigenvalues exhibit a very weak dependence on $k$; this contrasts with the results obtained for dry granular fluids (see, for instance, figure 4.7 of Garzó 2019).

Figure 8. Dispersion relations for a three-dimensional granular gas with $\phi =0.001$, $T_g^{*}=1000$, $m/m_g=1$ and $\alpha =0.8$. From top to bottom, the curves correspond to the longitudinal mode $\lambda _{1,||}$, the two degenerate shear (transversal) modes $\lambda _\perp$ (dotted line) and the two remaining longitudinal modes $\lambda _{3,||}$ and $\lambda _{2,||}$. The dependence of $\lambda _{1,||}$, $\lambda _{3,||}$ and $\lambda _\perp$ on $k$ is shown more clearly in the insets. Only the real parts of the eigenvalues are plotted.

7. Summary and concluding remarks

The main goal of this paper has been to determine the Navier–Stokes–Fourier transport coefficients of a granular gas (modelled as a gas of inelastic hard spheres) immersed in a bath of elastic hard spheres (molecular gas). We have been interested in a situation where the solid particles are sufficiently dilute, and hence one can assume that the state of the bath (molecular gas) is not affected by the presence of the granular gas. Under these conditions, the molecular gas can be considered as a thermostat kept at equilibrium at a temperature $T_g$. This system (granular gas thermostatted by a gas of elastic hard spheres) was originally proposed years ago by Biben et al. (2002) and it can be considered as an idealised collisional model for particle-laden suspensions. Thus, in contrast to the previous suspension models employed in the granular literature (Tsao & Koch 1995; Sangani et al. 1996; Garzó et al. 2012; Saha & Alam 2017) where the effect of the interstitial fluid on grains is accounted for via an effective fluid–solid force, the model considered here takes into account not only the inelastic collisions among grains themselves but also the elastic collisions between particles of the granular and molecular gases. Moreover, we have also assumed that the volume fraction occupied by the suspended solid particles is very small (low-density regime). In this case, the one-particle velocity distribution function $f(\boldsymbol {r}, \boldsymbol {v};t)$ of grains verifies the Boltzmann kinetic equation.

Before analysing inhomogeneous states, we have considered first homogeneous situations. The study of this state is important because it plays the role of the reference base state in the Chapman–Enskog solution to the Boltzmann equation. In this simple situation, the granular ‘temperature’ $T$ tends to thermalise to the bath temperature $T_g$. In the steady state, both competing effects (characterised by the cooling rates $\zeta$ and $\zeta _g$) cancel each other and a breakdown of energy equipartition appears ($T< T_g$). In the HSS, the temperature ratio $T/T_g$ and the kurtosis $a_2$ (measuring the deviation of the distribution function from its Maxwellian form) have been here estimated by considering the so-called first Sonine approximation (3.2) to the distribution function $f(\boldsymbol {v})$. In this approximation, the temperature ratio is obtained by solving the steady-state condition $\zeta +\zeta _g=0$ while the kurtosis is given by (3.19). These equations provide the dependence of $T/T_g$ and $a_2$ on the parameter space of the system. Our theoretical results extend to arbitrary dimensions the results obtained by Santos (2003) for hard spheres ($d=3$). To assess the accuracy of the (approximate) analytical results, a suite of Monte Carlo simulations have been also performed. Comparison between theory and simulations shows in general a very good agreement, especially in the case of the temperature ratio.

Once the homogeneous state is characterised, the next step has been to solve the Boltzmann equation by means of the Chapman–Enskog-like expansion to first order in spatial gradients (Chapman & Cowling 1970; Brilliantov & Pöschel 2004; Garzó 2019). The Navier–Stokes–Fourier transport coefficients have been explicitly determined by considering the leading terms in a Sonine polynomial expansion when steady-state conditions apply. Their forms are given by (5.1) for the shear viscosity $\eta$, (5.4) for the thermal conductivity $\kappa$, (5.10) for the diffusive heat conductivity $\bar {\mu }$ and (5.11) for the so-called velocity conductivity coefficient $\kappa _U$. This latter coefficient takes into account a contribution to the heat flux coming from the velocity difference ${\rm \Delta} \boldsymbol {U}$. It is quite apparent that the expressions for the transport coefficients show a complex dependence on $\alpha$, $m/m_g$, $T_g^{*}$ and $\phi$ (see figures 4–7). The dependence on the latter two parameters appears via the dimensionless drift coefficient $\gamma ^{*}$; this coefficient provides a characteristic rate for the elastic collisions between granular and bath particles. In general, we find that significant quantitative differences between dry granular gases (Brey et al. 1998; Garzó & Dufty 1999; Garzó 2013, 2019; Gupta 2020) and granular suspensions appear as the inelasticity in collisions increases. Thus, the impact of the gas phase on the transport coefficients of the granular gas cannot be in general neglected.

Interestingly, in the Brownian limit ($m/m_g \to \infty$), a careful analysis shows that the expressions of $\eta$, $\kappa$ and $\bar {\mu }$ reduce to those previously derived by Gómez González & Garzó (2019) using the Langevin-like model (2.16) for the instantaneous gas–solid force. In this limiting case, the present results show that the coefficient $\kappa _U$ vanishes, in agreement with Gómez González & Garzó (2019). Therefore, the results reported in this paper extend to arbitrary values of the mass ratio $m/m_g$ the results derived in previous works (Garzó et al. 2012; Gómez González & Garzó 2019). Nonetheless, the convergence to the Brownian limit is reached for relatively small values of the mass ratio ($m/m_g\simeq 50$ for $T_g^{*}=1000$). Thus, from a practical point of view, it seems that in most cases the Langevin-like model is able to assess the effect of the interstitial gas on transport properties of granular gas. As discussed in § 3.2, the convergence of the present results to the Brownian limit depends on the way of scaling the variables. Here, since we have been interested in recovering the results provided by Gómez González & Garzó (2019) in the Brownian limit, $T_g^{*}$ has been used as the control parameter regarding the force exerted by the molecular gas on grains. This fact implies that the influence of the mass ratio $m/m_g$ appearing in $\gamma$ (see (2.15)) on the thermalisation process is absorbed in the selection of $T_g^{*}$. On the other hand, following the scaling proposed by Biben et al. (2002) and Santos (2003), another possibility could have been to choose $\omega$ (see (3.20)) as the bath parameter. In such a way, an extra dependence on the mass ratio in the coefficient $\gamma ^{*}$ emerges and leads to a different dependence of both the temperature ratio (see figure 2) and the transport coefficients on the mass ratio.

Knowledge of the forms of the transport coefficients opens up the possibility of performing a linear stability analysis on the resulting continuum hydrodynamic equations. As in the Brownian limit case (Gómez González & Garzó 2019), the analysis shows that the HSS is always linearly stable whatever the mass ratio considered.

Although in real experiments the energy input for keeping rapid-flow conditions is usually done either by driving through the boundaries (Yang et al. 2002) or by bulk driving (as in air-fluidised beds; Schröter, Goldman & Swinney 2005; Abate & Durian 2006), these ways of supplying energy produce in many cases strong spatial gradients in the bulk domain, and so the Navier–Stokes hydrodynamics does not hold. To overcome this difficulty, it is quite common in computer simulations (see e.g. Puglisi et al. 1998; Paganobarraga et al. 2002; Prevost, Egolf & Urbach 2002; Fiege, Aspelmeier & Zippelius 2009; Shaebani, Sarabadani & Wolf 2013) to heat the system homogeneously by the presence of an external driving force or thermostat (Evans & Morriss 1990). A different (and likely more realistic) way of thermostatting a granular gas is by means of a sea of elastic hard spheres (Biben et al. 2002). In this context, given that in most of the computer simulation works the effects of the thermostat on the properties of the granular gas are ignored, the results derived in this paper may be useful for simulators when studying problems in granular fluids thermostatted by a bath of elastic hard spheres. Regardless of practical applications, needless to say a complete comprehension of gas–solid flows is still missing (Fullmer & Hrenya 2017; Morris 2020; Han et al. 2021). For this reason, the development of theoretical models for granular suspensions where collisions play a significant role becomes of great relevance to understand from a more fundamental view the results derived by means of effective models such as the Langevin-like model (Garzó et al. 2012).

One of the main limitations of the results derived in this paper is its restriction to the low-density regime. The extension of the present theory to a moderately dense granular suspension described by the Enskog kinetic equation is an interesting project for the future. These results could stimulate the performance of MD simulations to assess the reliability of the theory for finite densities. Furthermore, in an attempt to model gas–solid flows in fluidised beds where the fluid–solid interactions are governed not only by the drag force (Stokes drag force) but also by the Archimedes force (gas pressure gradient), the present model could be also extended to account for these new ‘additional forces’. Another challenging work could be the determination of the non-Newtonian rheological properties of a granular suspension under simple shear flow. This study would allow the extension of previous studies (Tsao & Koch 1995; Sangani et al. 1996; Chamorro et al. 2015; Saha & Alam 2017; Alam et al. 2019; Takada et al. 2020) to arbitrary values of the mass ratio $m/m_g$. The accuracy of the results for the shear viscosity can be tested by means of DSMC simulations. For instance, by performing local deviations of the velocity field from its value in the homogeneous state and analysing the relaxation process (Brey, Ruiz-Montero & Cubero 1999b; Brey & Cubero 2001). Another possible project could be to revisit the results obtained in this paper by considering the charge transport equation recently considered by Ceresiat, Kolehmainen & Ozel (2021). Work along these lines will be carried out in the future.