Transport coefficients for granular suspensions at moderate densities

TL;DR

This paper derives explicit expressions for the transport coefficients of moderately dense granular suspensions using kinetic theory, accounting for gas-particle interactions, and finds the homogeneous steady state to be linearly stable unlike in dry granular materials.

Contribution

The study extends the Enskog kinetic theory to include gas effects in dense granular suspensions and provides explicit Navier-Stokes transport coefficients considering inelasticity and density.

Findings

Transport coefficients depend on inelasticity and density differently than in dry granular media.

Homogeneous steady state in suspensions is linearly stable, unlike in dry granular systems.

Explicit expressions for transport coefficients are obtained under steady-state conditions.

Abstract

The Enskog kinetic theory for moderately dense granular suspensions is considered as a model to determine the Navier-Stokes transport coefficients. The influence of the interstitial gas on solid particles is modeled by a viscous drag force term plus a stochastic Langevin-like term. The suspension model is solved by means of the Chapman--Enskog method conveniently adapted to dissipative dynamics. The momentum and heat fluxes as well as the cooling rate are obtained to first order in the deviations of the hydrodynamic field gradients from their values in the homogeneous steady state. Since the cooling terms (arising from collisional dissipation and viscous friction) cannot be compensated for by the energy gained by grains due to collisions with the interstitial gas, the reference distribution (zeroth-order approximation of the Chapman--Enskog solution) depends on time through its dependence on temperature. On the other hand, to simplify the analysis and given that we are interested in computing transport properties in the first order of deviations from the reference state, the steady-state conditions are considered. This simplification allows us to get explicit expressions for the Navier--Stokes transport coefficients. As expected, the results show that the dependence of the transport coefficients on both inelasticity and density is clearly different from that found in its granular counterpart (no gas phase). Finally, a linear stability analysis of the hydrodynamic equations with respect to the homogeneous steady state is performed. In contrast to the granular case (no gas-phase), no instabilities are found and hence, the homogeneous steady state is (linearly) stable.

Equations

$$\frac{\partial f}{\partial t}+\mathbf{v}\cdot \nabla f+\mathcal{F}f=J_{\text{E}}[\mathbf{r},\mathbf{v}|f,f],$$

$$J_{\text{E}}\left[\mathbf{r},{\mathbf{v}}_1|f,f\right]={\sigma }^{d-1}\int d{\mathbf{v}}_2\int d\hat{\mathbf{\sigma }}\Theta (\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})(\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})\left[{\alpha }^{-2}{f}_2(\mathbf{r},\mathbf{r}-\mathbf{\sigma },{\mathbf{v}}_1^{\prime \prime },{\mathbf{v}}_2^{\prime \prime },t)-f_2(\mathbf{r},\mathbf{r}+\mathbf{\sigma },{\mathbf{v}}_1,{\mathbf{v}}_2,t)\right]$$

$$f_2({\mathbf{r}}_1,{\mathbf{r}}_2,{\mathbf{v}}_1,{\mathbf{v}}_2,t)=\chi ({\mathbf{r}}_1,{\mathbf{r}}_2)f({\mathbf{r}}_1,{\mathbf{v}}_1,t)f({\mathbf{r}}_2,{\mathbf{v}}_2,t),$$

$${\mathbf{v}}_1^{\prime \prime }={\mathbf{v}}_1-\frac{1}{2}\left(1+{\alpha }^{-1}\right)(\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})\hat{\mathbf{\sigma }},{\mathbf{v}}_2^{\prime \prime }={\mathbf{v}}_2+\frac{1}{2}\left(1+{\alpha }^{-1}\right)(\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})\hat{\mathbf{\sigma }}.$$

$$n(\mathbf{r},t)=\int d\mathbf{v}f(\mathbf{r},\mathbf{v},t).$$

$${\mathcal{F}}^{\text{drag}}f\to -\gamma \frac{\partial }{\partial \mathbf{v}}\cdot \left(\mathbf{v}-{\mathbf{U}}_g\right)f,$$

$$⟨{\mathbf{F}}_i^{\text{st}}(t)⟩=0,⟨{\mathbf{F}}_i^{\text{st}}(t){\mathbf{F}}_j^{\text{st}}(t^{\prime })⟩=2m\gamma T_{\text{ex}}\mathsf{I}{\delta }_{ij}\delta (t-t^{\prime }),$$

$${\mathcal{F}}^{\text{st}}f\to -\frac{\gamma T_{\text{ex}}}{m}\frac{{\partial }^{2}f}{\partial v^2}.$$

$$\varphi =\frac{{\pi }^{d/2}}{2^{d-1}d\Gamma \left(\frac{d}{2}\right)}n{\sigma }^d.$$

$$\gamma ={\gamma }_{0}R(\varphi ),$$

$$R(\varphi )=1+3\sqrt{\frac{\varphi }{2}}.$$

$$R(\varphi )=1+\frac{3}{\sqrt{2}}{\varphi }^{1/2}+\frac{135}{64}\varphi \ln \varphi +11.26\varphi (1-5.1\varphi +16.57{\varphi }^2-21.77{\varphi }^3)-\varphi \chi (\varphi )\ln {ϵ}_m.$$

$$\frac{\partial f}{\partial t}+\mathbf{v}\cdot \nabla f-\gamma \Delta \mathbf{U}\cdot \frac{\partial f}{\partial \mathbf{v}}-\gamma \frac{\partial }{\partial \mathbf{v}}\cdot \mathbf{V}f-\gamma \frac{T_{\text{ex}}}{m}\frac{{\partial }^{2}f}{\partial v^2}=J_{\text{E}}[\mathbf{r},\mathbf{V}|f,f].$$

$$\mathbf{U}(\mathbf{r},t)=\frac{1}{n(\mathbf{r},t)}\int d\mathbf{v}\mathbf{v}f(\mathbf{r},\mathbf{v},t)$$

$$T(\mathbf{r},t)=\frac{m}{dn(\mathbf{r},t)}\int d\mathbf{v}{V}^{2}f(\mathbf{r},\mathbf{v},t).$$

$$D_{t}n+n\nabla \cdot \mathbf{U}=0,$$

$$D_t\mathbf{U}=-{\rho }^{-1}\nabla \cdot \mathsf{P}-\gamma \Delta \mathbf{U},$$

$$D_{t}T+\frac{2}{dn}\left(\nabla \cdot \mathbf{q}+\mathsf{P}:\nabla \mathbf{U}\right)=2\gamma \left(T_{\text{ex}}-T\right)-\zeta T.$$

$${\mathsf{P}}^{\text{k}}=\int d\mathbf{v}m\mathbf{V}\mathbf{V}f(\mathbf{r},\mathbf{v},t),{\mathbf{q}}^{\text{k}}=\int d\mathbf{v}\frac{m}{2}{V}^2\mathbf{V}f(\mathbf{r},\mathbf{v},t),$$

$${\mathsf{P}}^{\text{c}}=\frac{1+\alpha }{4}m{\sigma }^d\int d{\mathbf{v}}_1\int d{\mathbf{v}}_2\int d\hat{\mathbf{\sigma }}\Theta (\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})(\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12}{)}^2\hat{\mathbf{\sigma }}\hat{\mathbf{\sigma }}{\int }_0^{1}dx{f}_2\left[\mathbf{r}-x\mathbf{\sigma },\mathbf{r}+(1-x)\mathbf{\sigma },{\mathbf{v}}_1,{\mathbf{v}}_2,t\right],$$

$${\mathbf{q}}^{\text{c}}=\frac{1+\alpha }{4}m{\sigma }^d\int d{\mathbf{v}}_1\int d{\mathbf{v}}_2\int d\hat{\mathbf{\sigma }}\Theta (\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})(\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12}{)}^2({\mathbf{G}}_{12}\cdot \hat{\mathbf{\sigma }})\hat{\mathbf{\sigma }}{\int }_0^{1}dx{f}_2\left[\mathbf{r}-x\mathbf{\sigma },\mathbf{r}+(1-x)\mathbf{\sigma },{\mathbf{v}}_1,{\mathbf{v}}_2,t\right],$$

$$\zeta =\frac{\left(1-{\alpha }^2\right)}{4dnT}m{\sigma }^{d-1}\int d{\mathbf{v}}_1\int d{\mathbf{v}}_2\int d\hat{\mathbf{\sigma }}\Theta (\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})(\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12}{)}^3{f}_2(\mathbf{r},\mathbf{r}+\mathbf{\sigma },{\mathbf{v}}_1,{\mathbf{v}}_2,t).$$

$$\frac{\partial f}{\partial t}-\gamma \frac{\partial }{\partial \mathbf{v}}\cdot \mathbf{v}f-\gamma \frac{T_{\text{ex}}}{m}\frac{{\partial }^{2}f}{\partial v^2}=J_{\text{E}}[\mathbf{V}|f,f],$$

$$J_{\text{E}}\left[f,f\right]=\chi {\sigma }^{d-1}\int d{\mathbf{v}}_2\int d\hat{\mathbf{\sigma }}\Theta (\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})(\hat{\mathbf{\sigma }}\cdot {\mathbf{g}}_{12})\left[{\alpha }^{-2}f(v_1^{\prime \prime })f(v_2^{\prime \prime })-f(v_1)f(v_2)\right].$$

$${\partial }_{t}T=2\gamma \left(T_{\text{ex}}-T\right)-\zeta T.$$

$$\bigg{[}2\gamma\left(\theta^{-1}-1\right)-\zeta\bigg{]}T\frac{\partial f}{\partial T}-\gamma\frac{\partial}{\partial\mathbf{v}}\cdot\mathbf{v}f-\frac{\gamma T_{\text{ex}}}{m}\frac{\partial^{2}f}{\partial v^{2}}=J_{\text{E}}[f,f],$$

$$\zeta (t)=\frac{{\pi }^{(d-1)/2}}{4d\Gamma \left(\frac{d+3}{2}\right)}(1-{\alpha }^2)\frac{m{\sigma }^{d-1}}{nT}\chi \int d{\mathbf{v}}_1\int d{\mathbf{v}}_2{g}_{12}^{3}f({\mathbf{v}}_1,t)f({\mathbf{v}}_2,t).$$

$$f_0(v,t)=n{\left(\frac{m}{2\pi T(t)}\right)}^{d/2}\exp \left(-\frac{m{v}^2}{2T(t)}\right)$$

$${\partial }_{t}T=2\gamma \left(T_{\text{ex}}-T\right).$$

$$2\gamma \left({\theta }_{\text{s}}^{-1}-1\right)-{\zeta }_{\text{s}}=0,$$

Full text

I Introduction

Although in nature granular matter is surrounded by an interstitial fluid (like the air, for instance), most of theoretical and computational studies have neglected the impact of the gas phase on the dynamics of solid particles. On the other hand, it is known that in many practical applications (like for instance species segregation in granular mixtures Möbius et al. (2001); Naylor et al. (2003); Sánchez et al. (2004); Wylie et al. (2008); Clement et al. (2010); Pastenes et al. (2014)) the effect of the surrounding fluid on grains cannot be ignored. Needless to say, at a kinetic theory level, the description of granular suspensions ( namely, a suspension of solid particles in a viscous gas) is a quite complex problem since a complete microscopic description of the gas-solid system involves the solution of a set of two coupled kinetic equations for each one of the velocity distribution functions of the different phases. Thus, due to the mathematical difficulties embodied in this approach and in order to gain some insight into this problem, an usual model for describing gas-solid flows Koch and Hill (2001) is to consider a kinetic equation for the solid particles where the influence of the surrounding fluid on them is modeled by means of an effective external force. As usual Garzó et al. (2012); Hayakawa et al. (2017), the external force modeling the effect of the gas phase is constituted by two terms: (i) a viscous drag force (via a term involving a drift or friction coefficient $\gamma$) accounting for the friction of grains on the interstitial fluid and (ii) a stochastic Langevin-like term (via a term involving the background or bath temperature $T_{\text{ex}}$) accounting for the energy gained by the grains due to their collisions with particles of the background fluid.

Recently the above suspension model has been employed to study the so-called discontinuous shear thickening in non-Newtonian gas-solid suspensions Hayakawa et al. (2017); Gómez González and Garzó (2019). The results show the transition from the discontinuous shear thickening (observed for very dilute gases) to the continuous shear thickening as the density of the system increases. These analytical results (approximately obtained by means of Grad’s moment method Hayakawa et al. (2017) and from an exact solution of the Boltzmann equation for inelastic Maxwell models Gómez González and Garzó (2019)) compare quite well with molecular dynamics simulations Hayakawa et al. (2017) for conditions of practical interest. This good agreement highlights again the good performance of kinetic theory tools in reproducing the transport properties of gas-solid flows.

On the other hand, to the best of our knowledge, most of the efforts in kinetic theory of granular suspensions has been mainly focused on non-Newtonian transport properties (which are directly related with the pressure tensor). In particular, much less is known about the energy transport in gas-solid flows. The knowledge of the transport coefficients associated with the heat flux is interesting by itself and also for possible practical applications in suspensions where temperature and density gradients are present in the system. In this context, it would be desirable to provide simulators with the appropriate expressions of the Navier–Stokes transport coefficients to work when studying gas-solid flows where collisions among particles are inelastic.

The aim of this paper is to determine the Navier–Stokes transport coefficients of granular suspensions in the framework of the Enskog kinetic equation. Since this equation applies for moderate densities (let’s say for instance, solid volume fraction $\phi\lesssim 0.25$ for hard spheres), the comparison between kinetic theory and molecular dynamics simulations becomes practical. Attempts on the evaluation of the Navier–Stokes transport coefficients for granular suspensions modeled by the Enskog equation have been previously published. Thus, in Ref. Garzó et al. (2012) the authors determined the transport coefficients of gas-solid flows starting from the suspension model constituted by the viscous drag force plus the stochastic Langevin term. Their results show that the effect of the gas phase on both the shear viscosity and the diffusive heat conductivity coefficients is non-negligible for industrially relevant portions of the parameter space. However, for the sake simplicity, the temperature dependence of the scaled friction coefficient $\gamma^{*}=\gamma/\nu(T)$ (where $\nu\propto T^{1/2}$ is an effective collision frequency for hard spheres and $T$ is the granular temperature) was implicitly neglected in the above calculations Garzó et al. (2012) to get analytic (explicit) expressions for the transport coefficients. The above temperature dependence of $\gamma^{*}$ was accounted for in a subsequent paper Garzó et al. (2016) but for a simplified model where only the drag force term was considered in the Enskog equation.

A more careful study was carried out later in Ref. Garzó et al. (2013a) where the transport coefficients were explicitly computed by considering both the temperature dependence of the reduced friction coefficient as well as the complete form of the suspension model. On the other hand, although computer simulations Koch and Sangani (1999) have clearly shown that the friction coefficient depends on the volume fraction, the calculations performed in Ref. Garzó et al. (2013a) were carried out by assuming that the driven parameters of the model are constant. Needless to say, the impact of the density dependence of $\gamma$ on transport properties is expected to be more relevant as the gas phase becomes denser. Apart from this simplification, although not explicitly stated, another limitation of the above theory Garzó et al. (2013a) is that it was obtained by neglecting contributions to the transport coefficients coming from an additional density dependence of the zeroth-order distribution $f^{(0)}$ (in fact, although this simplification was noted in a subsequent erratum Garzó et al. (2013b), it has not been implemented so far in the calculations). This extra density dependence of $f^{(0)}$ is expected to be involved in the evaluation of the heat flux transport coefficients.

The question arises then as to whether, and if so to what extent, the conclusions drawn from Ref. Garzó et al. (2013a) may be altered when the above two new ingredients (density dependence of both the distribution $f^{(0)}$ and the friction coefficient $\gamma$) are accounted for in the theory. In this paper we address this question by extending the results derived in Ref. Garzó et al. (2013a) to situations not covered by previous studies on granular suspensions. The present theory subsumes all previous analyses Garzó et al. (2012, 2016); Garzó et al. (2013a), which are recovered in the appropriate limits. In particular, a comparison between the results obtained here for the transport coefficients with those derived in Garzó et al. (2013a) shows that while the expression of the shear viscosity coefficient is formally equivalent to the one obtained before, the heat flux transport coefficients and the first-order contribution to the cooling rate differ from those reported in Ref. Garzó et al. (2013a).

As in previous works Garzó and Dufty (1999); Lutsko (2005); Garzó et al. (2012), the transport coefficients are obtained by solving the Enskog equation by means of the application of the Chapman–Enskog method Chapman and Cowling (1970). Since a reference equilibrium state is missing in granular gases, an important point in the Chapman–Enskog expansion is the choice of the zeroth-order solution $f^{(0)}$ (reference base state of the perturbation scheme). While in the dry granular case (no gas phase) the distribution $f^{(0)}$ is chosen to be the local version of the homogeneous cooling state, there is more flexibility in the choice of $f^{(0)}$ in driven granular gases ( or, equivalently in gas-solid flows). In the case of gas-solid flows Garzó et al. (2013a), for simplicity one possibility is to take a steady distribution $f^{(0)}$ at any point of the system Garzó and Montanero (2002); Garzó (2011). However, the presence of the interstitial fluid introduces the possibility of a local energy unbalance and hence, the zeroth-order distribution is not in general a stationary distribution. This fact introduces new contributions to the transport coefficients, which were not considered when a local steady state was assumed at zeroth-order Garzó and Montanero (2002); Garzó (2011). Thus, for general small deviations from the homogeneous steady state the energy gained by grains due to collisions with the background fluid cannot be compensated locally with the cooling terms (viscous friction plus inelastic collisions). Thus, although we are interested in determining the transport coefficients under steady state conditions, we have to start from an unsteady zeroth-order solution in order to achieve the integral equation verifying the first-order solution $f^{(1)}$. The solution to this equation under steady state conditions provides the explicit forms of the transport coefficients.

The plan of the paper is as follows. In section II, the Enskog kinetic equation for granular suspensions is introduced and the corresponding balance equations for the densities of mass, momentum, and energy are derived. Then, section III studies the homogeneous steady state where some theoretical predictions are compared against available computer simulation results. The comparison shows an excellent agreement for conditions of practical interest. Section IV addresses the Chapman–Enskog expansion around the unsteady reference distribution $f^{(0)}(\mathbf{r},\mathbf{v},t)$ up to first-order in spatial gradients. The explicit expressions of the Navier–Stokes transport coefficients and the cooling rate are displayed in section V for steady state conditions. In dimensionless form, these coefficients are given in terms of the coefficient of restitution $\alpha$, the volume fraction $\phi$, and the (reduced) background temperature $T_{\text{ex}}^{*}$. The dependence of the transport coefficients and the cooling rate on the parameter space is illustrated for several systems showing that the influence of the gas phase on them is in general quite significant. As an application of the results found here, a linear stability analysis of the Navier–Stokes hydrodynamic equations around the homogeneous steady state is carried out in section VI; the analysis shows that the homogeneous state is linearly stable. This finding agrees with the previous stability analysis performed in Ref. Garzó et al. (2013a). We close the paper in section VII with a brief discussion of the results reported here.

II Enskog kinetic equation for granular suspensions

We consider a set of solid particles of diameter $\sigma$ and mass $m$ immersed in a viscous gas. Collisions between grains are inelastic and are characterized by a (positive) constant coefficient of normal restitution $\alpha\leq 1$, where $\alpha=1$ corresponds to elastic collisions (ordinary gases). At moderate densities, the one-particle velocity distribution function of solid particles $f(\mathbf{r},\mathbf{v};t)$ obeys the Enskog kinetic equation

[TABLE]

where

[TABLE]

is the Enskog collision operator. Here,

[TABLE]

$d$ is the dimensionality of the system ($d=2$ for disks and $d=3$ for spheres), $\bm{\sigma}=\sigma\widehat{\bm{\sigma}}$, $\widehat{\bm{\sigma}}$ being a unit vector, $\Theta$ is the Heaviside step function, and ${\bf g}_{12}={\bf v}_{1}-{\bf v}_{2}$. The double primes on the velocities in Eq. (2) denote the initial values $\{\mathbf{v}_{1}^{\prime\prime},\mathbf{v}_{2}^{\prime\prime}\}$ that lead to $\{\mathbf{v}_{1},\mathbf{v}_{2}\}$ following a binary collision:

[TABLE]

In addition, $\chi[{\bf r},{\bf r}\pm\bm{\sigma}|\{n(t)]$ is the equilibrium pair correlation function at contact as a functional of the nonequilibrium density field $n({\bf r},t)$ defined by

[TABLE]

In Eq. (1), the operator $\mathcal{F}$ represents the fluid-solid interaction force that models the effect of the viscous gas on solid particles. In order to fully account for the influence of the interstitial molecular fluid on the dynamics of grains, a instantaneous fluid force model is employed Garzó et al. (2012, 2016); Hayakawa et al. (2017). For low Reynolds numbers, it is assumed that the external force $\mathbf{F}$ acting on solid particles is composed by two independent terms. One term corresponds to a viscous drag force $\mathbf{F}^{\text{drag}}$ proportional to the (instantaneous) velocity of particle $\mathbf{v}$. This term takes into account the friction of grains on the viscous gas. Since the model attempts to mimic gas-solid flows, the drag force is defined in terms of the relative velocity $\mathbf{v}-\mathbf{U}_{g}$ where $\mathbf{U}_{g}$ is the (known) mean flow velocity of the surrounding molecular gas. Thus, the drag force $\mathbf{F}^{\text{drag}}=-m\gamma\left(\mathbf{v}-\mathbf{U}_{g}\right)$ is represented in the Enskog equation (1) by the term

[TABLE]

where $\gamma$ is the drag or friction coefficient. The second term in the total force corresponds to a stochastic force that tries to simulate the kinetic energy gain due to eventual collisions with the (more rapid) molecules of the background fluid. It does this by adding a random velocity to each particle between successive collisions Williams and MacKintosh (1996). This stochastic force $\mathbf{F}^{\text{st}}$ has the form of a Gaussian white noise with the properties van Kampen (1981)

[TABLE]

where $\mathsf{I}$ is the unit tensor and $i$ and $j$ refer to two different particles. Here, $T_{\text{ex}}$ can be interpreted as the temperature of the background (or bath) fluid. In the context of the Enskog kinetic equation, the stochastic external force is represented by a Fokker–Planck operator of the form van Kampen (1981); van Noije and Ernst (1998)

[TABLE]

Note that the strength of correlation in Eq. (8) has been chosen to be consistent with the fluctuation-dissipation theorem for elastic collisions van Kampen (1981).

Although the drift coefficient $\gamma$ is in general a tensor, here for simplicity we assume that this coefficient is a scalar proportional to the square root of $T_{\text{ex}}$ because the drag coefficient is proportional to the viscosity of the solvent Koch and Hill (2001). In addition, as usual in granular suspension models Koch (1990); Koch and Sangani (1999), $\gamma$ is a function of the solid volume fraction

[TABLE]

Thus, the drift coefficient $\gamma$ can be written as

[TABLE]

where $\gamma_{0}\propto\eta_{g}\propto\sqrt{T_{\text{ex}}}$, $\eta_{g}$ being the viscosity of the solvent or gas phase. In the case of hard spheres ($d=3$), for Stokes flow we can use the existing analytical closure derived by Koch Koch (1990) for the function $R(\phi)$ in the case of very dilute suspensions ($\phi\leq 0.1$):

[TABLE]

For $\phi>0.1$, Koch and Sangani Koch and Sangani (1999) used simulations based on multipole expansions to propose the $\phi$-dependence of $R$. It is given by

[TABLE]

Here, $\epsilon_{m}\sigma$ can be regarded as a length scale characterizing the impact of non-continuum effects on the lubrication forces between two smooth particles at contact. Typical values of $\epsilon_{m}$ are in the range 0.01–0.05. Since this term contributes to $R(\phi)$ through a weak logarithmic factor, the influence of its explicit value is not important in the final results. Here, we take $\epsilon_{m}=0.01$ as a typical value.

The suspension model defined by Eqs. (1), (6), and (8) is a simplified version of the model employed in Ref. Garzó et al. (2013a) to get the Navier–Stokes transport coefficients. In this latter model Garzó et al. (2012), the friction coefficient of the drag force ($\gamma_{b}$ in the notation of Ref. Garzó et al. (2013a)) and the strength of the correlation ($\xi_{b}^{2}$ in the notation of Ref. Garzó et al. (2013a)) are considered to be in general different. Here, as mentioned before, both coefficients are related as $\xi_{\text{b}}^{2}=2\gamma_{\text{b}}T_{\text{ex}}/m^{2}$ to be consistent with the fluctuation-dissipation theorem. Thus, some of the results derived in this paper (mainly those regarding homogeneous states) can be directly obtained from those reported in Ref. Garzó et al. (2013a) by making the changes $\gamma_{b}\to m\gamma$ and $\xi_{b}^{2}\to 2\gamma T_{\text{ex}}/m$ with $R(\phi)=1$. We have preferred in this paper to adopt the notation introduced in Eqs. (6) and (8) because this is the notation used in previous studies of sheared granular suspensions Hayakawa et al. (2017); Gómez González and Garzó (2019).

According to Eqs. (6) and (8), the Enskog equation (1) reads

[TABLE]

Here, $\Delta\mathbf{U}=\mathbf{U}-\mathbf{U}_{g}$, $\mathbf{V}=\mathbf{v}-\mathbf{U}$ is the peculiar velocity, and

[TABLE]

is the mean particle velocity. Another relevant hydrodynamic field is the granular temperature $T(\mathbf{r},t)$ defined as

[TABLE]

Note that in the model defined in Garzó et al. (2013a)) the mean flow velocity of the interstitial gas is assumed to be equal to the mean flow velocity of solid particles ($\mathbf{U}_{g}=\mathbf{U}$) for the sake of simplicity.

The macroscopic balance equations for the granular suspension are obtained when one multiplies the Enskog equation (13) by $\{1,m{\bf v},mv^{2}\}$ and integrates over velocity. After some algebra, one gets the balance equations Garzó and Dufty (1999); Garzó et al. (2012); Garzó et al. (2013a)

[TABLE]

[TABLE]

[TABLE]

In the above equations, $D_{t}=\partial_{t}+{\bf U}\cdot\nabla$ is the material derivative and $\rho=mn$ is the mass density. The cooling rate $\zeta$ is proportional to $1-\alpha^{2}$ and is due to dissipative collisions. The pressure tensor ${\sf P}({\bf r},t)$ and the heat flux ${\bf q}({\bf r},t)$ have both kinetic and collisional transfer contributions, i.e., ${\sf P}={\sf P}^{\text{k}}+{\sf P}^{\text{c}}$ and ${\bf q}={\bf q}^{\text{k}}+{\bf q}^{\text{c}}$. Their kinetic contributions are defined by

[TABLE]

and the collisional transfer contributions are Garzó and Dufty (1999)

[TABLE]

[TABLE]

where ${\bf G}_{12}=\frac{1}{2}({\bf V}_{1}+{\bf V}_{2})$ is the velocity of center of mass. Finally, the cooling rate $\zeta$ is given by

[TABLE]

Before closing this section, it is important to recall the range of validity of the suspension model (13). As already discussed before Garzó et al. (2012), the assumptions made in the model are relevant to the range of dimensionless physical parameters encountered in a circulating fluidized bed (low Reynolds numbers and moderate densities). A crucial aspect of the model is that the form of the Enskog collision operator $J_{\text{E}}[\mathbf{r},\mathbf{v}|f,f]$ is assumed to be the same as for a dry granular gas (i.e., when the influence of the interstitial gas is neglected). This means that the collision dynamics does not contain any parameter of the environmental gas. As it has been noted in several papers Koch (1990); Tsao and Koch (1995); Sangani et al. (1996); Koch and Hill (2001); Wylie et al. (2003), the above assumption requires that the mean-free time between collisions is assumed to be much less than the time needed by the fluid forces to significantly affect the dynamics of solid particles. Thus, we expect that the suspension model (3) may be reliable in situations where the gas phase has a weak influence on the motion of grains (solid particles immersed in air, for instance). Of course, this assumption fails for instance in the case of liquid flows (high density) where the presence of fluid must be taken into account in the collision process.

III Homogeneous steady state

Before computing the transport coefficients, it is instructive to analyze the homogeneous steady state. This state was widely analyzed in Refs. Garzó et al. (2013a); Chamorro et al. (2013). For homogeneous situations, the density $n$ and the temperature $T$ are spatially uniform, and with an appropriate selection of the frame of reference, the mean flow velocities vanish ($\mathbf{U}=\mathbf{U}_{g}=\mathbf{0}$). Consequently, Eq. (13) becomes

[TABLE]

where

[TABLE]

Here, $\chi$ is the pair correlation function evaluated at the (homogeneous) density $n$. The collision operator (23) can be recognized as the Boltzmann operator for inelastic collisions multiplied by the factor $\chi$. For homogeneous states, the only nontrivial balance equation is that of the temperature (18):

[TABLE]

As usual, for times longer than the mean free time, one expects that the system achieves a hydrodynamic regime where the distribution $f$ qualifies as a normal distribution Chapman and Cowling (1970) in the sense that $f$ depends on time only through its dependence on the temperature $T$. In this regime, $\partial_{t}f=(\partial_{T}f)(\partial_{t}T)$ and Eq. (22) reads

[TABLE]

where $\theta\equiv T/T_{\text{ex}}$ and use has been made of Eq. (24). In addition, for homogeneous states, Eq. (21) gives the following form for the cooling rate $\zeta$:

[TABLE]

For elastic collisions ($\alpha=1$ and so, $\zeta=0$), as expected Eq. (25) admits the solution

[TABLE]

where the temperature obeys the time-dependent equation

[TABLE]

The system therefore is in a time-dependent “equilibrium state” before reaching the asymptotic steady state where $T=T_{\text{ex}}$. For inelastic collisions, $\zeta\neq 0$ and to date the solution to Eq. (25) has not been found.

On the other hand, after a transient stage, the system achieves a steady state characterized by the steady temperature $T_{\text{s}}$. According to Eq. (28), $T_{\text{s}}$ is given by the condition

[TABLE]

where the subscript s means that the quantities are evaluated at $T=T_{\text{s}}$. At a given value of the environmental temperature $T_{\text{ex}}$ (which acts as a bath temperature in the sense that it is considered as a thermal energy reservoir), Eq. (29) implies that in the steady state the energy gained by grains due to their collisions with the interstitial fluid ($\gamma T_{\text{ex}}$) is exactly compensated by the cooling terms arising from collisional dissipation ($\zeta T$) and viscous friction ($\gamma T$). Moreover, as usual in the granular literature, the effects of the energy balance on the internal degrees of freedom of grains are not considered in the description.

As shown in previous works Garzó et al. (2013a); Chamorro et al. (2013); García de Soria et al. (2012, 2013), dimensionless analysis requires that $f_{\text{s}}$ has the scaled form

[TABLE]

where $v_{0}=\sqrt{2T_{\text{s}}/m}$ is the thermal speed and the unknown scaled distribution $\varphi_{\text{s}}$ is a function of the dimensionless parameters $\mathbf{c}\equiv\mathbf{v}/v_{0}$ and $\gamma_{\text{s}}^{*}$ where

[TABLE]

Here, $T_{\text{ex}}^{*}\equiv T_{\text{ex}}/(m\sigma^{2}\gamma_{0}^{2})$ is the (reduced) background gas temperature. In the second relation of Eq. (31), $\ell=1/(n\sigma^{d-1})$ is proportional to the mean free path of hard spheres. The scaling given by Eq. (30) is equivalent to the one proposed in Refs. Garzó et al. (2013a); Chamorro et al. (2013) when one makes the mapping $\xi_{\text{s}}^{*}\to 2\lambda\theta_{\text{s}}^{-3/2}$ with $R(\phi)=1$. Here, $\xi_{\text{s}}^{*}$ is defined by Eq. (24) of Garzó et al. (2013a). This means that the results for homogeneous states can be directly obtained from those derived in Refs. Garzó et al. (2013a); Chamorro et al. (2013) by making the above change. On the other hand, we have preferred here to revisit the homogeneous state in order to check the previous results.

In terms of $\varphi_{\text{s}}$, in the steady state, Eq. (22) for $f_{\text{s}}$ can be rewritten as

[TABLE]

where we have introduced the dimensionless collision operator $J_{\text{E}}^{*}=\ell v_{0}^{d-1}J_{\text{E}}/n$. Although the exact form of $\varphi_{\text{s}}$ is not known, an indirect information on it can be obtained from the kurtosis or fourth cumulant

[TABLE]

The cumulant $a_{2,\text{s}}$ measures the deviation of $\varphi_{s}$ from its Maxwellian form $\pi^{-d/2}e^{-c^{2}}$. This coefficient can be obtained by multiplying Eq. (32) by $c^{4}$ and integrating over velocity. The result is

[TABLE]

where $\zeta_{\text{s}}^{*}\equiv\ell\zeta_{\text{s}}/v_{0}$ and

[TABLE]

Upon deriving Eq. (34) use has been made again of the steady state condition (29).

As expected, Eq. (34) cannot be solved unless one knows the collisional moments $\zeta_{\text{s}}^{*}$ and $\beta_{4}$. As in previous works van Noije and Ernst (1998); Garzó et al. (2013a); Chamorro et al. (2013), a good estimate of $\zeta_{\text{s}}^{*}$ and $\beta_{4}$ can be obtained by replacing $\varphi_{\text{s}}$ by its leading Sonine approximation van Noije and Ernst (1998):

[TABLE]

In this case, retaining only linear terms in $a_{2,\text{s}}$, one has

[TABLE]

where van Noije and Ernst (1998)

[TABLE]

[TABLE]

and

[TABLE]

With these results, Eq. (34) can be easily solved with the result

[TABLE]

Notice that in Eq. (41), $\gamma_{\text{s}}^{*}$ is consistently obtained from the steady state condition (29) by replacing $\zeta_{\text{s}}^{*}\to\zeta_{0}^{(0)}$. The expression (41) agrees with the one obtained in Ref. Garzó et al. (2013a) when one takes the steady state condition $\xi_{\text{s}}^{*}=2\gamma_{\text{s}}^{*}+\zeta_{0}^{(0)}$ in Eq. (31) of Garzó et al. (2013a).

Once $a_{2,\text{s}}$ is known, the dependence of the cooling rate on both the coefficient of restitution $\alpha$ and the (reduced) external temperature $T_{\text{ex}}^{*}$ can be obtained from the first relation of Eq. (37). Finally, the (reduced) steady temperature $\theta_{s}$ is determined by solving the cubic equation

[TABLE]

As expected, Eq. (42) is consistent with Eq. (33) of Ref. Garzó et al. (2013a) for the steady temperature when one takes $R(\phi)=1$ and makes the replacement $\xi_{\text{s}}^{*}\to 2\lambda\theta_{\text{s}}^{-3/2}$.

Figure 1 shows the $\alpha$-dependence of the fourth cumulant $a_{2,\text{s}}$ for hard disks ($d=2$) with the solid volume fraction $\phi=0.25$. In the case of hard disks, we have chosen the following form for $\chi(\phi)$ Torquato (1995):

[TABLE]

The theoretical results given by Eq. (41) are compared against the results obtained in Ref. Garzó et al. (2013a) by numerically solving the Enskog equation from the direct simulation Monte Carlo (DSMC) method Bird (1994). The parameters of the simulation are $m=1$, $\sigma=0.01$, $\gamma_{0}=1$, and $T_{\text{ex}}=1$. In addition, the function $R(\phi)=1$ in the simulations. Although this figure was already presented in Ref. Garzó et al. (2013a), we plot it again here to remark the excellent agreement between theory and simulations observed in the complete range of values of $\alpha$. Since the values of $a_{2,\text{s}}$ are very small (in fact their magnitude is smaller than the one found in the dry granular case van Noije and Ernst (1998); Montanero and Santos (2000)) then, the Sonine approximation (36) can be considered as a good representation of the scaled distribution $\varphi_{s}(c)$. As a complement of Fig. 1, Fig. 2 shows $a_{2,\text{s}}$ versus $\phi$ for two values of $\alpha$. It is quite apparent that the qualitative dependence of the fourth cumulant on the density depends strongly on the inelasticity since while $a_{2,\text{s}}$ decreases monotonically with $\phi$ at $\alpha=0.8$, the opposite happens at $\alpha=0.6$. We do not actually have an intuitive explanation for the change of behaviour of $a_{2,\text{s}}$ when the coefficient of restitution varies from 0.8 to 0.6. Next, the (reduced) temperature $\theta_{s}$ is considered. Figure 3 shows $\theta_{\text{s}}$ versus $\alpha$ for $d=2$, $\phi=0.25$, and the same parameters as the one considered in Figs. 1 and 2. First, as expected, $\theta_{\text{s}}=1$ for elastic collisions. Moreover, the steady granular temperature decreases with inelasticity. It is illustrated in Fig. 4 ( which was also plotted in Ref. Garzó et al. (2013a)) where $\theta_{\text{s}}$ is plotted against the density $\phi$ for two different values of $\alpha$. Figures 3 and 4 highlight again the excellent agreement between theory and simulations, even for extreme values of both inelasticity and/or density.

IV Transport around the homogeneous steady state. Chapman–Enskog expansion

As in previous studies Garzó and Dufty (1999); Garzó et al. (2013a); Khalil and Garzó (2013), we assume that we perturb the homogeneous steady state by small spatial gradients. These perturbations give rise to nonzero contributions to the pressure tensor and the heat flux, which are characterized by transport coefficients. The evaluation of the transport coefficients is the main objective of the present contribution. In order to get them, we will solve the Enskog equation (13) by means of the Chapman–Enskog method Chapman and Cowling (1970) conveniently adapted to granular fluids. As usual, the Chapman–Enskog method assumes the existence of a normal solution such that all space and time dependence of the velocity distribution function occurs through the hydrodynamic fields, namely,

[TABLE]

The notation on the right hand side indicates a functional dependence on the density, temperature and flow velocity. For small spatial variations (i.e., low Knudsen numbers), this functional dependence can be made local in space through an expansion in the gradients of the hydrodynamic fields. To generate it, $f$ is written as a series expansion in a formal parameter $\epsilon$ measuring the non-uniformity of the system,

[TABLE]

where each factor of $\epsilon$ means an implicit gradient of a hydrodynamic field. In contrast to the case of dry granular gases Garzó and Dufty (1999), in ordering the different level of approximations in the kinetic equation, one has to characterize the magnitude of the drift term $\gamma$ relative to the gradients as well as the term $\Delta\mathbf{U}$. With respect to the first term, since $\gamma$ does not induce any flux in the system, it is considered to be of zeroth-order in gradients. Regarding the term $\Delta\mathbf{U}$, since in the absence of gradients $\mathbf{U}$ tends to $\mathbf{U}_{g}$ after a transient period, then $\Delta\mathbf{U}$ is expected to be at least to first order in the spatial gradients.

According to the expansion (45), the Enskog operator $J_{\text{E}}$ and the time derivative $\partial_{t}$ are also given in the representations

[TABLE]

The coefficients in the time derivative expansion are identified by a representation of the fluxes and the cooling rate in the macroscopic balance equations as a similar series through their definitions as functionals of $f$. This is the usual Chapman–Enskog method Chapman and Cowling (1970); Garzó and Santos (2003) for solving kinetic equations. The expansions (46) yield similar expansions for the heat and momentum fluxes and the cooling rate when substituted into Eqs. (19)–(21):

[TABLE]

Here, we shall restrict our calculations to the first order in the uniformity parameter $\epsilon$.

IV.1 Zeroth-order approximation

To zeroth order in the expansion, the distribution $f^{(0)}$ obeys the kinetic equation

[TABLE]

where $J_{\text{E}}^{(0)}[f^{(0)},f^{(0)}]$ is given by Eq. (23) with the replacement $f_{\text{s}}\to f^{(0)}({\bf r},{\bf v},t)$. The conservation laws at this order are given by $\partial_{t}^{(0)}n=0$, $\partial_{t}^{(0)}\mathbf{U}=\mathbf{0}$, and

[TABLE]

where $\zeta^{(0)}$ is determined from Eq. (21) to zeroth order. In particular, as said in section III, a good approximation to $\zeta^{(0)}$ is given by the first relation of Eq. (37), namely,

[TABLE]

The kinetic equation (48) can be rewritten in terms of the derivative $\partial_{T}f^{(0)}$ when one takes into account the zeroth-order balance equations:

[TABLE]

Equation (51) has the same form as the corresponding Enskog equation (25) for a strictly homogeneous state. However, in Eq. (51) $f^{(0)}(\mathbf{r},\mathbf{v},t)$ is a local distribution. Therefore, as in the homogeneous state, the solution to Eq. (51) can be written in the form (30) (with the replacement $T_{\text{s}}\to T$) where the scaled distribution $\varphi(\mathbf{c},\lambda,\theta)$ obeys the unsteady equation

[TABLE]

where $\zeta_{0}^{*}\equiv\ell\zeta^{(0)}/v_{0}(T)$ and $\gamma^{*}=\lambda\theta^{-1/2}$. Upon deriving Eq. (52) use has been made of the property

[TABLE]

where the derivative $\partial\varphi/\partial\theta$ is taken at constant $\mathbf{c}$.

The velocity distribution function $f^{(0)}$ is isotropic in $\mathbf{V}$ so that, according to Eqs. (19)–(II), the heat flux to zeroth-order vanishes as expected ($\mathbf{q}^{(0)}=\mathbf{0}$) and the pressure tensor $P_{ij}^{(0)}=p\delta_{ij}$, where the hydrostatic pressure is

[TABLE]

As discussed in section III, although the explicit form of $\varphi$ is not known, a good approximation is given by the Sonine approximation (36). In particular, the equation for the unsteady fourth cumulant $a_{2}$ can be easily obtained from Eq. (52) as

[TABLE]

where $\Lambda^{(0)}\equiv 2\gamma^{*}\left(\theta^{-1}-1\right)-\zeta_{0}^{*}$ and $\beta_{4}$ is defined in Eq. (35). In the steady state, $\Lambda^{(0)}=0$ and the solution to Eq. (55) is given by Eq. (41) once one expands $\zeta_{0}^{*}$ and $\beta_{4}$ in powers of $a_{2}$. Beyond the steady state, Eq. (55) must be numerically solved to get the dependence of $a_{2}$ on the (reduced) temperature. On the other hand, as we will show in section V, in order to get the transport coefficients in the steady state we need to know the derivatives $\Delta_{\theta}\equiv(\partial a_{2}/\partial\theta)_{\text{s}}$, $\Delta_{\lambda}\equiv(\partial a_{2}/\partial\lambda)_{\text{s}}$, and $\Delta_{\chi}\equiv(\partial a_{2}/\partial\chi)_{\text{s}}$. These derivatives provide an indirect information (through the fourth cumulant) on the departure of the time-dependent solution $f^{(0)}$ from its stationary form $f_{\text{s}}$. According to Eq. (55), the former derivative is given by

[TABLE]

where here the expansions (37) have been considered and as usual nonlinear terms in $a_{2}$ have been neglected. In the steady state, the numerator and denominator of Eq. (56) vanish, hence the quantity $\Delta_{\theta}$ becomes indeterminate. As in Ref. Garzó et al. (2013a), this problem can be solved by applying l’Hôpital’s rule. The final result yields a quadratic equation for $\Delta_{\theta}$. However, given that the magnitude of $\Delta_{\theta}$ is quite small, one can neglect the term proportional to $\Delta_{\theta}^{2}$ in the above quadratic equation and obtain the simple expression

[TABLE]

Equation (57) is consistent with Eq. (A6) of Ref. Garzó et al. (2013a) when one neglects nonlinear terms in $(\partial a_{2}/\partial\xi^{*})_{\text{s}}$ and takes $\beta=\frac{1}{2}$. The derivatives $\Delta_{\lambda}$ and $\Delta_{\chi}$ can be easily derived from Eq. (55) with the results

[TABLE]

[TABLE]

Note that in Eqs. (57)–(59), $\theta_{\text{s}}$ is obtained from Eq. (42) by neglecting $a_{2,\text{s}}$. The dependence of the derivatives $\Delta_{\theta}$, $\Delta_{\lambda}$, and $\Delta_{\chi}$ on the coefficient of restitution $\alpha$ is plotted in Fig.  5 for $d=3$ and $\phi=0.25$. Here, $T_{\text{ex}}^{*}=0.9$; this is a typical value for the (reduced) background temperature used in previous simulations Hayakawa et al. (2017). It is seen that while the magnitude of $\Delta_{\lambda}$, and $\Delta_{\chi}$ is much smaller than that of the kurtosis $a_{2,\text{s}}$, $\Delta_{\theta}$ is of the same order of magnitude as $a_{2,\text{s}}$.

IV.2 First-order approximation

The mathematical steps involved in the derivation of the first-order distribution function $f^{(1)}$ are quite similar to those carried out in Ref. Garzó et al. (2013a). On the other hand, given that the calculations performed in this paper take into account some additional density dependencies not accounted for in the previous derivation Garzó et al. (2013a), we have preferred here to perform an independent calculation where most of the technical details are provided in the Appendix A for the sake of completeness. To first-order in spatial gradients, $f^{(1)}$ is given by

[TABLE]

where, in the steady state ($\Lambda^{0)}=0$), the quantities $\bm{\mathcal{A}}$, $\bm{\mathcal{B}}$, $\mathcal{C}_{ij}$, and $\mathcal{D}$ verify the following set of coupled linear integral equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In Eq. (64), $\zeta_{1}^{(1)}$ is a functional of $\mathcal{D}$ defined by Eq. (131). Moreover, in Eqs. (61)–(64), $\mathcal{L}$ is the linearized collision operator

[TABLE]

$R$ is defined by Eqs. (10)–(12) and the coefficients $\mathbf{A}$, $\mathbf{B}$, $C_{ij}$, and $D$ are functions of the peculiar velocity $\mathbf{V}$ and the hydrodynamic gradients. They are defined by Eqs. (111)–(114). Note that all the quantities appearing in Eqs. (61)–(64) are evaluated in the steady state (the subscript s has been omitted here for the sake of simplicity). Thus, the transport coefficients obtained by solving Eqs. (10)–(12) will be provided in terms of the steady temperature $T_{\text{s}}$. It is worthwhile to remark that since we are here interested in obtaining the momentum and heat fluxes in the first order of the deviations from the steady state, we only need to know the transport coefficients to zeroth order in the deviations. This means that the solution to the integral equations (61)–(64) will provide us the forms of the transport coefficients and the cooling rate in steady state conditions.

According to the Chapman–Enskog scheme Chapman and Cowling (1970), acceptable solutions to Eqs. (61)–(64) must obey

[TABLE]

These are necessary conditions for the solution to the integral equations to exist (the so-called Fredholm alternative MM56 ). Since $\bm{\mathcal{A}}(\mathbf{V})\propto\mathbf{A}(\mathbf{V})$, $\bm{\mathcal{B}}(\mathbf{V})\propto\mathbf{B}(\mathbf{V})$, $\mathcal{C}_{ij}(\mathbf{V})\propto C_{ij}(\mathbf{V})$, and $\mathcal{D}(\mathbf{V})\propto D(\mathbf{V})$, then the solubility conditions (66) can be proved when one takes into account the explicit forms of $\mathbf{A}$, $\mathbf{B}$, $C_{ij}$, and $D$.

V Navier–Stokes transport coefficients

To first order in the spatial gradients, the constitutive equations for the pressure tensor $P_{ij}^{(1)}$ and the heat flux $\mathbf{q}^{(1)}$ are

[TABLE]

[TABLE]

Here, $\eta$ is the shear viscosity, $\eta_{\text{b}}$ is the bulk viscosity, $\kappa$ is the thermal conductivity, and $\mu$ is the diffusive heat conductivity. This latter coefficient vanishes for ordinary gases ($\alpha=1$). While the coefficients $\eta$, $\kappa$, and $\mu$ have kinetic and collisional contributions, the bulk viscosity $\eta_{\text{b}}$ has only collisional contributions and hence, it vanishes for dilute gases. In addition, as already mentioned in Ref. Garzó et al. (2012), the forms of the collisional contributions to the transport coefficients are exactly the same as those obtained in the dry granular case (namely, in the absence of the gas phase) Garzó and Dufty (1999); Lutsko (2005), except that $a_{2,\text{s}}$ depends on $\gamma^{*}$. Thus, we will focus here our attention on the kinetic contributions to the transport coefficients and the cooling rate. Some technical details on this calculation are provided in the Appendix B.

V.1 Shear and bulk viscosities

The bulk viscosity $\eta_{\text{b}}$ is given by

[TABLE]

where

[TABLE]

is the low density value of the shear viscosity for an ordinary gas of hard spheres ($\alpha=1$). The shear viscosity $\eta$ can be written as

[TABLE]

where $K^{\prime}=(d+2)/8K$, $K$ is defined by Eq. (40) and the (reduced) collision frequency $\nu_{\eta}^{*}$ is Garzó et al. (2007)

[TABLE]

where $a_{2,\text{s}}$ is defined by Eq. (41). The expression (71) for the shear viscosity agrees with the one obtained in Ref. Garzó et al. (2013a) when $R(\phi)=1$. This is because the new contributions to the fluxes coming from the extra density dependencies not accounted for in Garzó et al. (2013a) do not affect the form of the pressure tensor.

V.2 Thermal conductivity and diffusive heat conductivity

The thermal conductivity is given by

[TABLE]

where

[TABLE]

is the low density value of the thermal conductivity for an ordinary gas of hard spheres ($\alpha=1$) and $\kappa_{\text{k}}$ denotes the kinetic contribution to the thermal conductivity. Its explicit expression is

[TABLE]

where $\zeta_{0}^{(1)}$ is defined by Eq. (38) and the derivative $\Delta_{\theta}$ is given by Eq. (57). In addition, the (reduced) collision frequency $\nu_{\kappa}^{*}$ is Garzó et al. (2007)

[TABLE]

To compare the expression (75) with the one derived in Ref. Garzó et al. (2013a) (see Eq. (65) of this reference), one has to make the mapping $\xi_{\text{s}}^{*}(\partial a_{2}/\partial\xi^{*})_{\text{s}}\to-(2/3)\theta_{\text{s}}\Delta_{\theta}$ and takes $R(\phi)=1$. In this case, we find that the form (75) of the thermal conductivity coefficient is consistent with the one obtained in Garzó et al. (2013a) except for the last term of the numerator (i.e., the term proportional to $\frac{1}{2}\left(1+\alpha\right)\theta_{\text{s}}\Delta_{\theta}$). This term comes from the collision integral (128). We have checked that Eq. (75) gives the correct result and hence it fixes the slight mistake of Eq. (65) of Ref. Garzó et al. (2013a).

The diffusive heat conductivity $\mu$ is

[TABLE]

where the kinetic contribution $\mu_{\text{k}}$ is given by

[TABLE]

Here, the derivatives $\Delta_{\lambda}$ and $\Delta_{\chi}$ are defined by Eqs. (58) and (59), respectively. The expression (78) agrees with Eq. (69) of Ref. Garzó et al. (2013a) when one neglects (i) the density dependence of the function $R$ (i.e., $\partial_{\phi}R=0$) and (ii) all the derivatives of $a_{2}$ with respect to $\theta$, $\lambda$, and $\chi$ in the steady state (i.e., $\Delta_{\theta}=\Delta_{\lambda}=\Delta_{\chi}=0$). In addition, as in the case of a dry granular gas Brey et al. (1998a); Garzó and Dufty (1999); Lutsko (2005), the coefficient $\mu$ vanishes for elastic collisions.

V.3 Cooling rate

The cooling rate is

[TABLE]

where $\zeta_{\text{s}}^{(0)}$ is given by Eq. (50) with the replacement $T\to T_{\text{s}}$. The coefficient $\zeta_{U}$ can be written as

[TABLE]

where

[TABLE]

[TABLE]

Here, we have introduced the quantities

[TABLE]

[TABLE]

It is quite apparent that $\zeta_{U}=0$ for elastic collisions ($\alpha=1$). As in the case of the diffusive heat conductivity, to compare Eq. (82) with the expression (73) for $\zeta_{1}^{(1)}$ obtained in Ref. Garzó et al. (2013a) one has to make the replacement $\theta\Delta_{\theta}\to-(3/2)\xi_{\text{s}}^{*}(\partial a_{2}/\partial\xi^{*})_{\text{s}}$, take $R(\phi)=1$, and neglect the derivatives of $a_{2}$ with respect to $\lambda$ and $\chi$ ($\Delta_{\lambda}=\Delta_{\chi}=0$). After these changes, we see that both results agree except for a misprint we have found in Eq. (73) of Ref. Garzó et al. (2013a). Note also that $\zeta_{U}\neq 0$ for dilute granular suspensions García de Soria et al. (2013).

V.4 Some illustrative examples

In summary, the Navier–Stokes transport coefficients $\eta_{\text{b}}$, $\eta$, $\kappa$, and $\mu$ are given by Eqs. (69), (71), (73), and (77), respectively, while the first-order contribution $\zeta_{U}$ to the cooling rate is given by Eqs. (80)–(82). As expected, all these coefficients present an intricate dependence on the coefficient of restitution $\alpha$, the density $\phi$, and the (reduced) background temperature $T_{\text{ex}}^{*}$. In addition, their dimensionless forms are defined in terms of the steady temperature $\theta_{\text{s}}$ and the derivatives $\Delta_{\theta}$, $\Delta_{\lambda}$, and $\Delta_{\chi}$. While these derivatives are explicitly given by Eqs. (57)–(59), the granular temperature is given in terms of the physical solution of the cubic equation (42).

As in previous works Garzó et al. (2012); Garzó et al. (2013a), it is quite apparent that one of the principal new features of the present paper lies on the dependence of the Navier–Stokes transport coefficients of granular suspensions on the coefficient of restitution $\alpha$. Therefore, to illustrate the differences between granular ($\alpha\neq 1$) and ordinary ($\alpha=1$) suspensions, the transport coefficients are scaled with respect to their values for elastic collisions. In addition, we consider a three-dimensional system ($d=3$) with $T_{\text{ex}}^{*}=0.9$ and three different values of the volume fraction $\phi$: $\phi=0.01$ (very dilute system), $\phi=0.1$, and $\phi=0.2$ (moderately dense system).

In Figs. 6–8, the above Navier–Stokes transport coefficients are plotted as functions of $\alpha$. While in the case of the shear viscosity and thermal conductivity coefficients we observe that their deviation from their forms for elastic collisions is in general significant, no happens the same in the case of the diffusive heat conductivity since the magnitude of the scaled coefficient $n\mu(\alpha)/T\kappa(1)$ is much smaller than that of the (scaled) coefficient $\kappa(\alpha)/\kappa(1)$. Since both $\kappa$ and $\mu$ characterize the heat flux, one could neglect the term proportional to the density gradient in the heat flux. Thus, for practical purposes and analogously to ordinary (elastic) suspensions, one could assume that the heat flux verifies Fourier’s law $\mathbf{q}^{(1)}=-\kappa\nabla T$. With respect to the $\alpha$-dependence of $\eta$ and $\kappa$, Figures 6 and 7 highlight that both transport coefficients are decreasing functions of the inelasticity regardless of the density of the system. In addition, the influence of collisional dissipation on momentum and heat transport increases with density, being very tiny in the limit of dilute suspensions. A comparison with the results obtained for dry granular fluids (see for instance, Fig. 1 of Ref. Garzó (2005)) shows significant differences between dry (no gas phase) and granular suspensions. In particular, both theory Garzó and Dufty (1999); Lutsko (2005); Garzó (2013) and simulations Montanero et al. (2005) show that for relatively dilute dry granular gases ($\phi\lesssim 0.1$) $\eta$ increases with inelasticity, while the opposite occurs for sufficiently dense dry granular fluids ($\phi\gtrsim 0.1$). The same qualitative behavior is observed for the thermal conductivity coefficient Garzó and Dufty (1999); Lutsko (2005); Garzó (2013). This non-monotonic behavior contrasts with the predictions found here for granular suspensions where $\eta$ and $\kappa$ always decreases with decreasing $\alpha$. Regarding the coefficient $\mu$, we see that the impact of density on it is significant since while $\mu$ is always positive for dilute suspensions, it can be negative for moderately dense suspensions. It is worthwhile to note that the behavior of the shear viscosity and thermal conductivity on both density and coefficient of restitution found here is qualitatively similar to that of a confined quasi-two-dimensional granular fluid Garzó et al. (2018).

Finally, the dependence of the magnitude of the first-order contribution $|\zeta_{U}|$ to the cooling rate is plotted in Fig. 9 for the same parameters employed in Figs. 6–8. As the coefficient $\mu$, $\zeta_{U}=0$ for elastic collisions. On the other hand, in contrast to the diffusive heat conductivity, we observe that the influence of inelasticity on $\zeta_{U}$ is important, specially at large densities. This means that the contribution of $\zeta_{U}$ to the cooling rate must be considered as the inelasticity increases.

VI Stability of the homogeneous steady state

The knowledge of the Navier–Stokes transport coefficients and the cooling rate opens up the possibility of solving the hydrodynamic equations for $n$, $\mathbf{U}$, and $T$ for situations close to the homogeneous steady state. The solution of the linearized hydrodynamic equations allows us to study the stability of the homogeneous steady state. This is likely one of the nicest applications of the Navier–Stokes equations. In order to obtain them, one has to substitute the equation of state (54), the Navier–Stokes constitutive equations (67) and (68) for the pressure tensor and heat flux, respectively, and Eq. (79) for the cooling rate into the exact balance equations (16)–(18). The Navier–Stokes hydrodynamic equations read

[TABLE]

[TABLE]

[TABLE]

As mentioned in several previous papers Garzó (2005); Garzó et al. (2006), the general form of the cooling rate $\zeta$ should include second-order gradient contributions of the form $\zeta_{n}\nabla^{2}n$ and $\zeta_{T}\nabla^{2}T$ in Eq. (87). Nevertheless, as shown for a dilute (dry) granular gas Brey et al. (1998a), given that the ratios $\zeta_{n}/\mu$ and $\zeta_{T}/\kappa$ were shown to be very small for not very inelastic particles, the terms $\zeta_{n}\nabla^{2}n$ and $\zeta_{T}\nabla^{2}T$ were neglected in the Navier–Stokes transport equations. We assume that the same happens for dense gases and hence, these second-order contributions can be neglected for practical purposes. Apart from this approximation, Eqs. (85)–(87) are exact to second order in the spatial gradients for a granular suspension at moderate densities.

The stability analysis of the homogeneous steady state was also carried out in Ref. Garzó et al. (2013a). On the other hand and as mentioned in section I, the present work generalizes the results derived before Garzó et al. (2013a) since it takes into account both an extra density dependence of the zeroth-order distribution $f^{(0)}$ and the dependence of the friction coefficient $\gamma$ on the volume fraction $\phi$ ($R(\phi)\neq 1$). Thus, it is worth to assess to what extent the previous theoretical results Garzó et al. (2013a) are indicative of what happens in the stability analysis of the homogeneous state when the above density dependencies for the transport coefficients and the cooling rate are considered. This is the main motivation of this Section.

To analyze the stability of the homogeneous solution, Eqs. (85)–(87) must be linearized around the homogeneous steady state. In this state, the hydrodynamic fields take the homogeneous steady values $n\equiv\text{const.}$, $T_{\text{s}}\equiv\text{const.}$, and $\mathbf{U}_{g}=\mathbf{U}\equiv\mathbf{0}$. For small spatial gradients, we assume that the deviations $\delta y_{\beta}(\mathbf{r},t)=y_{\beta}(\mathbf{r},t)-y_{\beta,\text{s}}$ are small, where $\delta y_{\beta}(\mathbf{r},t)$ denotes the deviations of the hydrodynamic fields $\left\{y_{\beta};\beta=1,\cdots,d+2\right\}=\left\{n,\mathbf{U},T\right\}$ from their values in the homogeneous steady state. Moreover, as usual we also suppose that the interstitial fluid is not perturbed and hence, $\mathbf{U}_{g}=\mathbf{U}=\mathbf{0}$.

It must be recalled that here, in contrast to the linear stability analysis for dry granular gases Garzó (2005); Brilliantov and Pöschel (2004); G19 , the reference state is stationary and so one does not have to eliminate the time dependence of the transport coefficients. On the other hand, in order to compare our results with those obtained for granular fluids Garzó (2005), the following space and time variables are introduced:

[TABLE]

The dimensionless time scale $\tau$ measures the average number of collisions per particle in the time interval between 0 and $t$. The unit length $\mathbf{r}^{\prime}$ is proportional to the mean free path of solid particles. As usual, a set of Fourier transformed dimensionless variables are then introduced by

[TABLE]

where $\delta y_{\mathbf{k}\beta}\equiv\left\{\rho_{\mathbf{k}}(\tau),\mathbf{w}_{\mathbf{k}}(\tau),\theta_{\mathbf{k}}(\tau)\right\}$ is defined as

[TABLE]

where here the wave vector $\mathbf{k}$ is dimensionless.

In terms of the above dimensionless variables, as expected, the $d-1$ transverse velocity components $\mathbf{w}_{\mathbf{k}\perp}=\mathbf{w}_{\mathbf{k}}-\left(\mathbf{w}_{\mathbf{k}}\cdot\widehat{\mathbf{k}}\right)\widehat{\mathbf{k}}$ (orthogonal to the wave vector $\mathbf{k}$) decouple from the other three modes. Their evolution equation is

[TABLE]

where $\eta^{*}=\eta/\sigma^{1-d}\sqrt{mT_{\text{s}}}$. The solution to Eq. (91) is

[TABLE]

Since both the (reduced) friction coefficient $\gamma_{\text{s}}^{*}$ and the (reduced) shear viscosity coefficient $\eta^{*}$ are positive, then the transversal shear modes of the granular suspension are linearly stable.

The remaining (longitudinal) modes correspond to $\rho_{\mathbf{k}}$, $\theta_{\mathbf{k}}$, and the longitudinal velocity component of the velocity field, $w_{\mathbf{k}\parallel}=\mathbf{w}_{\mathbf{k}}\cdot\widehat{\mathbf{k}}$ (parallel to $\mathbf{k}$). These modes are coupled and obey the equation

[TABLE]

where $\delta y_{\mathbf{k}\beta}(\tau)$ denotes now the set $\left\{\rho_{\mathbf{k}},w_{\mathbf{k}\parallel},\theta_{\mathbf{k}}\right\}$ and $\mathsf{M}$ is the square matrix

[TABLE]

Here, the (reduced) transport coefficient $\nu_{\ell}^{*}$, $\mu^{*}$, and $D_{\text{T}}^{*}$ are defined as

[TABLE]

while $p^{*}\equiv p_{\text{s}}/nT_{\text{s}}=1+2^{d-2}(1+\alpha)\chi\phi$, $\rho=mn$, and the quantities $C_{p}$, $C_{\chi}$, $C_{n}$, and $C_{\gamma}$ are given by

[TABLE]

[TABLE]

In the above equations, it is understood that the transport coefficients $\eta^{*}$, $\nu_{\ell}^{*}$, $D_{T}^{*}$, and $\mu^{*}$ are evaluated in the homogeneous steady state.

The longitudinal three modes have the form $\text{exp}\left[\Lambda_{\ell}(k)\tau\right]$ for $\ell=1,2,3,$ where $\Lambda_{\ell}(k)$ are the eigenvalues of the matrix $\mathsf{M}$, namely, they are the solutions of the cubic equation

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

In general, one of the longitudinal modes can be unstable for $k<k_{\text{h}}$, where $k_{\text{h}}$ is obtained from Eq. (98) when $\Lambda=0$, namely, $Z(k_{\text{h}})=0$. The result is

[TABLE]

At a fixed value of the background temperature $T_{\text{ex}}^{*}$, a careful analysis of the dependence of $k_{\text{h}}^{2}$ on both the coefficient of restitution $\alpha$ and the volume fraction $\phi$ shows that $k_{\text{h}}^{2}$ is always negative. This means that there are no physical values of the wave numbers for which the longitudinal modes become unstable. Therefore, as in the case of the transversal shear modes, we can conclude that all the eigenvalues of the dynamical matrix $\mathsf{M}$ have a positive real part and no instabilities are found in the homogeneous steady state of a granular suspension.

In summary, the stability analysis performed here by including the extra density dependencies of the transport coefficients shows no surprises relative to the earlier analysis Garzó et al. (2013a): the homogenous steady state of a moderately dense granular suspension is linearly stable. On the other hand, the dispersion relations derived here are different from those obtained in Ref. Garzó et al. (2013a) since for instance the functional form of the heat flux transport coefficients differs in both approaches.

VII Conclusions

In this paper we have undertaken a rather complete study on the transport properties of granular suspensions in the Navier–Stokes domain (first-order in the spatial gradients). The starting point of our study has been the Enskog kinetic equation where the effect of the gas phase on the solid particles is via the introduction of two additional terms: (i) a viscous drag force term proportional to the velocity of particle and (ii) a stochastic Langevin-like term. While the first term attempts to model the friction of solid particles on the viscous surrounding gas, the second term mimics the kinetic energy gained by grains due to eventual collisions with the more rapid molecules of the interstitial gas. Both terms are characterized by the friction coefficient $\gamma$ (which is a function of the volume fraction $\phi$) and the background temperature $T_{\text{ex}}$ (which is a known quantity of the model).

A previous attempt on the derivation of the Navier–Stokes transport coefficients of dense granular suspensions was worked out by Garzó et al. Garzó et al. (2013a) by starting from a similar suspension model. However, the above work has two deficiencies: (i) it neglects an additional density dependence of the zeroth-order distribution $f^{(0)}$ through the parameter $\lambda(\phi)$ (defined in Eq. (31)), and (ii) it assumes that the friction coefficient $\gamma$ is constant. While the former simplification may be relevant in the evaluation of the diffusive heat conductivity coefficient (the transport coefficient associated to the density gradient in the heat flux), the latter simplification may be not reliable as the suspension becomes denser. The present analysis incorporates both extra new ingredients (the density dependence of $\lambda$ in $f^{(0)}$ and $\gamma=\gamma_{0}R(\phi)$, $\gamma_{0}$ being constant) in the Chapman–Enskog solution. The results show that while these two new density dependencies do not formally affect the expression of the shear viscosity coefficient obtained in Ref. Garzó et al. (2013a), the forms of the heat flux transport coefficients and the cooling rate obtained here differ from those derived before. These findings are likely the most significant contributions of the present work. In this context, this paper complements and extends previous papers on transport properties in granular suspensions Garzó et al. (2012); Garzó et al. (2013a); Garzó et al. (2016).

Before considering inhomogeneous situations, the homogeneous steady state has been analyzed. As expected, after a transient period, the steady distribution function $f_{\text{s}}$ adopts the form (30) where the temperature dependence of the scaled distribution $\varphi_{\text{s}}$ is encoded through the dimensionless velocity $\mathbf{c}=\mathbf{v}/v_{0}$ ($v_{0}=\sqrt{2T_{\text{s}}/m}$ being the thermal speed) and the (scaled) friction coefficient $\gamma_{\text{s}}^{*}=\lambda(\phi)\theta_{\text{s}}^{-1/2}$ ($\theta_{\text{s}}=T_{\text{s}}/T_{\text{ex}}$ being the reduced steady temperature). As in previous works on granular fluids driven by thermostats García de Soria et al. (2012); Garzó et al. (2013a), the above scaling differs from the one assumed for undriven granular fluids van Noije and Ernst (1998); Brilliantov and Pöschel (2004); G19 where $\varphi_{\text{s}}$ depends on $T$ only through the scaled velocity $\mathbf{c}$. Although the exact form of $\varphi_{\text{s}}$ is not known, a good approximation of this distribution (at least in the thermal velocity region $c\sim 1$) is provided by the leading Sonine approximation (36). By using this distribution, we have explicitly obtained the fourth cumulant $a_{2,\text{s}}$; this coefficient provides an indirect information on the deviation of $\varphi_{\text{s}}$ from its Maxwellian form $\pi^{-d/2}e^{-c^{2}}$. Once $a_{2,\text{s}}$ is known, the steady temperature $\theta_{\text{s}}$ is obtained by solving the cubic equation (42). In spite of the above approximations, the theoretical predictions for $\theta_{\text{s}}$ and $a_{2,\text{s}}$ show an excellent agreement with Monte Carlo simulation results. As expected, the results obtained for homogeneous systems agree with those derived in Ref. Garzó et al. (2013a) when one makes the mapping $\xi_{\text{s}}^{*}\to 2\lambda\theta_{\text{s}}^{-3/2}$ with $R(\phi)=1$.

Once the steady reference state is well characterized, we have considered the transport processes occurring in granular suspensions with small spatial gradients of the hydrodynamic fields. In this situation, the Enskog kinetic equation has been solved by means of the Chapman–Enskog method Chapman and Cowling (1970) where only terms up to the first order in the spatial gradients have been retained (Navier–Stokes hydrodynamic order). As in previous papers on the application of the Chapman–Enskog method to granular systems Brey et al. (1998a); Garzó and Dufty (1999); Lutsko (2005); Garzó et al. (2013a), the spatial gradients have been assumed to be independent of the coefficient of restitution $\alpha$. Thus, although the constitutive equations for the irreversible fluxes are limited to first order in spatial gradients, the corresponding transport coefficients appearing in these equations apply a priori to arbitrary degree of collisional dissipation. This type of expansion differs from the ones considered by other authors Goldhirsch and Sela (1996); Sela et al. (1996); Sela and Goldhirsch (1998); Goldhirsch et al. (2005) where the Chapman–Enskog solution is given in powers of both the Knudsen number (or spatial gradients as in the conventional scheme) and the degree of collisional dissipation $\delta\equiv 1-\alpha^{2}$. The results reported here are consistent with the ones obtained in those papers Goldhirsch and Sela (1996); Sela et al. (1996); Sela and Goldhirsch (1998); Goldhirsch et al. (2005) in the limit $\delta\to 0$.

As in the Chapman–Enskog solution obtained in Ref. Garzó et al. (2013a), a subtle but important point is the choice of the zeroth-order approximation $f^{(0)}$ in the perturbation expansion. Although we are interested in obtaining the transport coefficients in steady state conditions, for general small perturbations around the homogeneous steady state, the density and temperature are specified separately in the local reference state $f^{(0)}$ and hence, it is not expected that the temperature is stationary at any point of the system. This means that $\partial_{t}^{(0)}T\neq 0$ in the reference base state and consequently, the complete determination of the Navier–Stokes transport coefficients requires to know for instance the temperature dependence of the fourth cumulant $a_{2}$ of the unsteady reference state. This of course involves the numerical integration of the differential equation (56). This is quite an intricate problem that goes beyond the objective of this paper. Since we are essentially motivated by a desire for analytic expressions, the steady state conditions have been considered. On the other hand, given that $\partial_{t}^{(0)}T\neq 0$ in the Chapman–Enskog scheme, the transport coefficients are defined not only in terms of the hydrodynamic fields in the steady state but also there are contributions to the transport coefficients [such as the derivatives $\Delta_{\theta}$, $\Delta_{\lambda}$, and $\Delta_{\chi}$ defined by Eqs. (57)–(59), respectively] accounting for the vicinity of the perturbed state to the steady state.

As usual, in order to obtain explicit expressions for the transport coefficients, the leading terms in a Sonine polynomial expansion have been considered. These forms have been displayed along the section V: the bulk $\eta_{\text{b}}$ and shear $\eta$ viscosities are given by Eqs. (69) and (71), respectively, the thermal conductivity $\kappa$ is given by Eqs. (73) and (75), the heat diffusive conductivity $\mu$ is given by Eqs. (77) and (78) and the first-order contribution $\zeta_{U}$ to the cooling rate is given by Eqs. (81) and (82). As said before, the expressions of $\eta_{\text{b}}$ and $\eta$ agree with those derived in Garzó et al. (2013a) (once one takes $R(\phi)=1$) while the expressions of $\kappa$, $\mu$, and $\zeta_{U}$ reduce to those obtained in Garzó et al. (2013a) when the contributions coming from the derivatives $\Delta_{\theta}$, $\Delta_{\lambda}$, and $\Delta_{\chi}$ are neglected.

In reduced forms, it is quite apparent that the Navier–Stokes coefficients of the granular suspension exhibit a complex dependence on the (steady) temperature $\theta_{\text{s}}$, the coefficient of restitution $\alpha$, the solid volume fraction $\phi$, and the (reduced) background temperature $T_{\text{ex}}^{*}$. In addition, Figs. (6)–(8) highlight the significant impact of the gas phase on the Navier–Stokes transport coefficients $\eta$, $\kappa$, and $\mu$ since their $\alpha$-dependence is clearly different from the one previously found for dry granular gases Brey et al. (1998a); Garzó and Dufty (1999).

As an application of the previous results, the stability of the special homogeneous steady state solution has been analyzed. This has been achieved by solving the linearized Navier–Stokes hydrodynamic equations for small perturbations around the homogeneous steady state. The linear stability analysis performed here shows no new surprises relative to the earlier work Garzó et al. (2013a): the homogeneous steady state is linearly stable with respect to long enough wavelength excitations (namely, long enough small spatial gradients). On the other hand, it is worthwhile to recall that the conclusion reached here for the reference homogeneous steady state differs from the one found for freely cooling granular fluids where it was shown Brey et al. (1998a); Garzó (2005) that the resulting hydrodynamic equations exhibit a long wavelength instability for three of the hydrodynamic modes. This shows again the influence of the interstitial fluid on the dynamics of solid particles.

It is quite apparent that the theoretical results obtained in this paper under certain approximations should be tested against computer simulations. This would allow us to gauge the degree of accuracy of the theoretical predictions. As happens for dry granular gases Brey et al. (1998b, 1999); Brey and Ruiz–Montero (2004); Brey et al. (2005); Montanero et al. (2005, 2007); Mitrano et al. (2011, 2012); Brey and Ruiz–Montero (2013); Mitrano et al. (2014), we expect that the present results stimulate the performance of appropriate simulations where the kinetic theory calculations reported here can be assessed. We also plan to undertake such kind of simulations for the case of the shear viscosity. More specifically, we want to perform simulations of granular suspensions under uniform shear flow where the Navier–Stokes shear viscosity might be measured in the Newtonian regime (very small shear rates). Another possible project for the next future is the extension of the present results to the relevant subject of multicomponent granular suspensions. Work along these lines will be worked out in the near future.

Acknowledgements.

We want to thank Moisés García Chamorro for providing us the simulation data included in Figs. 1–4. The present work has been supported by the Spanish Government through Grant No. FIS2016-76359-P and by the Junta de Extremadura (Spain) Grant Nos. IB16013 (V.G.) and GR18079, partially financed by “Fondo Europeo de Desarrollo Regional” funds. The research of Rubén Gómez González has been supported by the predoctoral fellowship BES-2017-079725 from the Spanish Government.

Appendix A Some technical details on the first-order solution

Up to the first order in the expansion, the velocity distribution function $f^{(1)}$ verifies the Enskog kinetic equation

[TABLE]

where $D_{t}^{(1)}\equiv\partial_{t}^{(1)}+\mathbf{U}\cdot\nabla$ and $J_{\text{E}}^{(1)}[f,f]$ denotes the first-order contribution to the expansion of the Enskog collision operator in powers of the spatial gradients. In order to explicitly determine $J_{\text{E}}^{(1)}[f,f]$ we need the results

[TABLE]

[TABLE]

where $\chi$ is obtained from the functional $\chi(\mathbf{r},\mathbf{r}\pm\bm{\sigma}|n)$ by evaluating all density fields at $n(\mathbf{r},t)$. Taking into account Eqs. (104) and (105), $J_{\text{E}}^{(1)}$ reads Garzó et al. (2013a)

[TABLE]

where $\mathcal{L}$ is defined by Eq. (65) and the operator $\bm{\mathcal{K}}[X]$ is given by

[TABLE]

As already noted in Ref. Garzó et al. (2013a), upon obtaining Eq. (106) use has been made of the symmetry property $\mathcal{K}_{i}[\partial_{V_{j}}f^{(0)}]=\mathcal{K}_{j}[\partial_{V_{i}}f^{(0)}]$ that follows from the isotropy of the zeroth-order solution. Thus we are able to separate the contributions from the flow field gradients into independent traceless and diagonal components.

The macroscopic balance equations to first order in the gradients are

[TABLE]

where $\zeta^{(1)}$ is the first order contribution to the cooling rate. Since the cooling rate is a scalar, corrections to first-order in the gradients can arise only from $\nabla\cdot\mathbf{U}$ since $\nabla n$ and $\nabla T$ are vectors and the tensor $\partial_{j}U_{i}+\partial_{i}U_{j}-\frac{2}{d}\delta_{ij}\nabla\cdot\mathbf{U}$ is a traceless tensor. Thus, $\zeta^{(1)}$ can be written as

[TABLE]

The unknown quantity $\zeta_{U}$ is a functional of the first-order distribution $f^{(1)}$. A more explicit form for $\zeta_{U}$ is obtained by expanding Eq. (21) to first-order in gradients. This yields Eq. (80) where $\zeta_{1}^{(0)}$ and $\zeta_{1}^{(1)}$ are defined by Eqs. (81) and (131), respectively.

The use of the balance equations (108) allows us to evaluate the right-hand side of Eq. (103). The combination of these results with the forms (106) of the Enskog collision operator $J_{\text{E}}^{(1)}$ and (80) of $\zeta_{U}$ leads to the expression

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here, $p^{*}\equiv p/(nT)$. The structure of Eqs. (110)–(114) is formally equivalent to the ones derived for driven granular gases Garzó et al. (2013a). The only difference lies on the dependence of the zeroth-order solution $f^{(0)}$ on density and temperature.

As for dry granular gases Garzó and Dufty (1999), the solution to the kinetic equation (110) is given by Eq. (60) where the unknown functions $\bm{\mathcal{A}}$, $\bm{\mathcal{B}}$, $\mathcal{C}_{ij}$, and $\mathcal{D}$ are determined by solving Eq. (110). Since the gradients of the hydrodynamic fields are all independent, substitution of (60) into Eq. (110) yields a set of linear, inhomogeneous integral equations. In order to obtain them, one needs the result

[TABLE]

The integral equations (61)–(64) can be easily obtained after taking into account Eq. (115) and the steady state condition $\Lambda^{(0)}=0$.

Appendix B Kinetic contributions to the transport coefficients

In this Appendix we give some details on the determination of the kinetic contributions to the transport coefficients $\eta$, $\kappa$, and $\mu$ as well as the first-order contribution $\zeta_{U}$ to the cooling rate. Since all these quantities are obtained int he steady state, the subscript s appearing along the main text will be omitted here for the sake of brevity.

The kinetic part of the shear viscosity $\eta_{\text{k}}$ is defined as

[TABLE]

where

[TABLE]

As usual, to get $\eta_{\text{k}}$ one multiplies both sides of Eq. (61) by $D_{ij}$ and integrates over velocity. The result is

[TABLE]

where

[TABLE]

and Garzó and Dufty (1999); Lutsko (2005); Garzó (2013)

[TABLE]

The expression of $\eta_{\text{k}}$ can be easily obtained when one takes into account Eq. (120) and the explicit form (72) of $\nu_{\eta}$. This latter expression is obtained from Eq. (119) by considering the leading terms in a Sonine polynomial expansion of the unknown $\mathcal{C}_{ij}(\mathbf{V})$.

The kinetic parts $\kappa_{\text{k}}$ and $\mu_{\text{k}}$ are defined, respectively, as

[TABLE]

[TABLE]

where

[TABLE]

As in the case of $\eta_{\text{k}}$, $\kappa_{\text{k}}$ is obtained by multiplying both sides of Eq. (61) by $\mathbf{S}(\mathbf{V})$ and integrating over $\mathbf{v}$. The result is

[TABLE]

where use has been made of the steady state condition (29) and

[TABLE]

The right hand side of Eq. (124) can be computed when one takes into account Eq. (111) and the relationship (53). After some algebra, one gets

[TABLE]

where $f_{\text{M}}(\mathbf{c})=n\pi^{-d/2}v_{0}^{-d}\text{e}^{-c^{2}}$ and use has been made of the Sonine approximation (36) to $f^{(0)}$ and the property (53). The first collision integral involving the operator $\bm{\mathcal{K}}$ has been calculated in previous works Garzó and Dufty (1999); Lutsko (2005); Garzó (2013) and the result is

[TABLE]

The second collision integral in (126) has not been evaluated before. After some algebra, one gets

[TABLE]

With the above results, $\kappa_{\text{k}}$ can be finally written in the form (75). As in the case of $\nu_{\eta}$, the (reduced) collision frequency $\nu_{\kappa}$ can be well estimated by considering the leading Sonine approximation to $\bm{\mathcal{A}}$.

The evaluation of $\mu_{\text{k}}$ follows similar mathematical steps to those made for $\kappa_{\text{k}}$ since one has to multiply both sides of Eq. (62) by $\mathbf{S}(\mathbf{V})$ and integrate over $\mathbf{v}$. In order to get its explicit form (78), one needs the partial results

[TABLE]

[TABLE]

The expression (78) can be derived by using Eqs. (129) and (130).

Finally, the contribution $\zeta_{1}^{(1)}$ to the cooling rate $\zeta_{U}$ is defined as

[TABLE]

where the unknown function $\mathcal{D}(\mathbf{V})$ is the solution of the linear integral equation (64). As before, an approximate solution to (64) can be obtained by taking the Sonine approximation

[TABLE]

where

[TABLE]

The coefficient $e_{D}$ is given by

[TABLE]

Substitution of Eq. (132) into Eq. (131) gives

[TABLE]

where $\nu_{0}=nT/\eta_{0}$. The coefficient $e_{D}$ is obtained by substituting the Sonine solution (132) into the integral equation (64), multiplying it by the polynomial $F(\mathbf{V})$ and integrating over velocity. After some algebra one gets the expression (82) for $\zeta_{1}^{(1)}$.