Shear-rate dependent transport coefficients in granular suspensions

TL;DR

This paper investigates how shear rate influences transport coefficients in granular suspensions using a kinetic model based on the inelastic Boltzmann equation, providing explicit forms and analyzing shear-rate dependence.

Contribution

It introduces a Chapman-Enskog-like expansion for granular suspensions near simple shear flow, deriving shear-rate dependent transport coefficients with explicit solutions.

Findings

Transport coefficients depend on shear rate and restitution coefficient.

Explicit forms of generalized transport coefficients are obtained under steady-state conditions.

Shear-rate dependence of transport coefficients is illustrated for various restitution coefficients.

Abstract

A recent model for monodisperse granular suspensions is used to analyze transport properties in spatially inhomogeneous states close to the simple (or uniform) shear flow. The kinetic equation is based on the inelastic Boltzmann (for low density gases) with the presence of a viscous drag force that models the influence of the interstitial gas phase on the dynamics of grains. A normal solution is obtained via a Chapman-Enskog-like expansion around a (local) shear flow distribution which retains all the hydrodynamic orders in the shear rate. To first-order in the expansion, the transport coefficients characterizing momentum and heat transport around shear flow are given in terms of the solutions of a set of coupled linear integral equations which are approximately solved by using a kinetic model of the Boltzmann equation. To simplify the analysis, the steady-state conditions when viscous heating is compensated by the cooling terms arising from viscous friction and collisional dissipation are considered to get the explicit forms of the set of generalized transport coefficients. The shear-rate dependence of some of the transport coefficients of the set is illustrated for several values of the coefficient of restitution.

Equations

$${\mathbf{F}}_{\text{fluid}}=-m\gamma (\mathbf{v}-{\mathbf{U}}_g),$$

$$\frac{\partial f}{\partial t}+\mathbf{v}\cdot \nabla f+\frac{\partial }{\partial \mathbf{v}}\cdot \left(\frac{{\mathbf{F}}_{\text{fluid}}}{m}f\right)=J[\mathbf{v}|f,f],$$

$$J\left[{\mathbf{v}}_1|f,f\right]$$

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

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

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

$$\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=J[\mathbf{v}|f,f],$$

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

$$D_t\mathbf{U}+(mn{)}^{-1}\nabla \cdot \mathsf{P}=-\gamma \Delta \mathbf{U},$$

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

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

$$P_{ij}(\mathbf{r},t)=m\int \text{d}\mathbf{v}{V}_i{V}_{j}f(\mathbf{r},\mathbf{v},t),$$

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

$$\zeta (\mathbf{r},t)=-\frac{m}{dn(\mathbf{r},t)T(\mathbf{r},t)}\int \text{d}\mathbf{v}{V}^{2}J[\mathbf{v}|f,f]$$

$$\frac{\partial f}{\partial t}-\gamma \frac{\partial }{\partial \mathbf{v}}\cdot \mathbf{v}f=J[\mathbf{v}|f,f],$$

$$\frac{\partial \ln T}{\partial t}=-2\gamma .$$

$$T(t)=T(0)e^{-2\gamma t},$$

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

$$U_i=a_{ij}{r}_j,a_{ij}=a{\delta }_{ix}{\delta }_{jy},$$

$$-a{V}_y\frac{\partial f}{\partial V_x}-\gamma \frac{\partial }{\partial \mathbf{V}}\cdot \mathbf{V}f=J[\mathbf{V}|f,f].$$

$$\frac{2}{dn}{P}_{xy}a=-2T\gamma -\zeta T.$$

$$P_{xy}=-\eta (\alpha ,a)a.$$

$$a_{ik}{P}_{kj}+a_{jk}{P}_{ki}+2\gamma P_{ij}={\Lambda }_{ij},$$

$${\Lambda }_{ij}\equiv \int \text{d}\mathbf{v}m{V}_i{V}_{j}J[\mathbf{V}|f,f].$$

$$f(\mathbf{V})\to f_{\text{M}}(\mathbf{V})\left(1+\frac{m}{2n{T}^2}{V}_i{V}_j{\Pi }_{ij}\right),$$

$$f_{\text{M}}(\mathbf{V})=n{\left(\frac{m}{2\pi T}\right)}^{d/2}{e}^{-m{V}^2/2T}$$

$${\Pi }_{ij}=P_{ij}-p{\delta }_{ij}$$

$${\Lambda }_{ij}=-\nu \left(\beta {\Pi }_{ij}+{\zeta }^{\ast }{P}_{ij}\right),$$

$$\nu =\frac{8}{d+2}\frac{{\pi }^{(d-1)/2}}{\Gamma \left(\frac{d}{2}\right)}n{\sigma }^{d-1}\sqrt{\frac{T}{m}}$$

$${\zeta }^{\ast }=\frac{\zeta }{\nu }=\frac{d+2}{4d}\left(1-{\alpha }^2\right)$$

Full text

I Introduction

Although in nature granular materials are usually immersed in a gas or liquid phase (like the air, for instance), the influence of the latter on the transport properties of solid particles is generally neglected in most theoretical and computational studies. However, high-velocity, gas-solid flows occur in a wide range of practical applications (like circulating fluidized beds, for instance) and hence, the impact of the gas phase on grains should be accounted for in many circumstances. An example corresponds to species segregation problems where several works 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) have shown that the presence of the interstitial fluid may significantly change the segregation phase-diagrams obtained in previous studies for (dry) granular flows (namely, when the role of the gas phase is neglected).

At a kinetic theory level, the description of such multiphase flows is quite intricate since the system involves two different phases (solid particles and interstitial fluid) and hence, one would need to solve a set of two coupled kinetic equations for each one of the velocity distribution functions of the different phases. On the other hand, in order to gain some insight into this complex problem, most of the models proposed in the literature for gas-solid flows have considered a single kinetic equation for the solid particles where the effect of the surrounding fluid on them is taken into account through an effective external force $\mathbf{F}_{\text{fluid}}$ models .

A simple and realistic way of modeling the fluid-solid interaction force $\mathbf{F}_{\text{fluid}}$ is by means of a viscous drag force given by

[TABLE]

where $m$ and $\mathbf{v}$ are the mass and the velocity of the particles, respectively, $\gamma$ is the friction coefficient (assumed to be proportional to the gas viscosity $\mu_{g}$), and $\mathbf{U}_{g}$ is the (known) mean velocity of the gas phase. The model defined in Eq. (1) has been recently considered in different papers to study the shear rheology of frictional hard-sphere suspensions H13 ; SMMD13 ; WGZS14 ; MSMD15 . In addition, model (1) can be seen as a simplified version of the more general particle acceleration model proposed in Ref. GTSH12 where the effect of the gas phase is not only accounted for by the drag force (1) but also by means of a Langevin-like term. This latter term takes into account the added effects coming from neighboring particles and can be neglected when the mean velocity of the solid particles follows the mean flow velocity of the gas ($\mathbf{U}\simeq\mathbf{U}_{g}$). Here, $\mathbf{U}$ [defined below in Eq. (5)] denotes the mean flow velocity of the solid particles. Thus, the results derived from this simple version of the model can be considered of practical interest to analyze linear transport in dilute gas-solid flows when the mean flow velocity of the solid and gas phases are practically the same (like, for instance, in the simple or uniform shear flow (USF) state TK95 ; SMTK96 ; ChVG15 ).

An interesting problem is to assess the impact of the interstitial fluid on the transport properties of solid particles under USF. As usual, solid particles are modeled as a gas of inelastic smooth hard spheres with a constant coefficient of restitution $0<\alpha\leq 1$. The USF state is likely the simplest flow problem since the only nonzero hydrodynamic gradient is $\partial U_{x}/\partial y\equiv a$, where $a$ is the constant shear rate. Due to its simplicity, this state has been widely studied in the past for dry elastic GS03 and inelastic C90 ; G03 gases as an ideal testing ground to shed light on the response of the system to large shear rates. Years ago, two independent papers L06 ; G06 analyzed momentum and heat transport around USF for a dry dilute granular gas in spatially inhomogeneous states close to the USF. The heat and momentum fluxes were determined to first order in the deviations of the hydrodynamic field gradients from their values in the reference USF state. Given that the granular gas is strongly sheared, the corresponding transport coefficients are nonlinear functions of both the shear rate and the coefficient of restitution $\alpha$. This is one of the main new added values of these constitutive equations. On the other hand, in order to get explicit results and due to the mathematical difficulties involved in the general non-stationary problem, a particular sort of perturbations were considered to obtain the generalized transport coefficients under steady state conditions. Given that the (scaled) shear rate $a^{*}\equiv a/\nu$ [$\nu$ is a collision frequency for hard spheres, see Eq. (29)] and $\alpha$ are coupled in the steady state, then the generalized transport coefficients are only functions of the coefficient of restitution $\alpha$.

The aim of this paper is to study transport around USF in dilute granular suspensions. As said before, the starting point is the inelastic Boltzmann equation Brilliantov and Pöschel (2004); P15 with the presence of the viscous drag force (1). As in Refs. L06 ; G06 , the Boltzmann equation is solved by means of a Chapman-Enskog-like expansion CC70 around the reference USF distribution $f^{(0)}$. Since the latter applies for arbitrary values of the shear rate $a$, the successive approximations $f^{(k)}$ in the perturbation expansion retain all the hydrodynamic orders in $a$. Consequently, the problem deals with two kinds of spatial gradients: small gradients due to perturbations of the USF and arbitrary shear rates due to the background shear flow. As in Refs. L06 ; G06 , the study here is restricted to first order in the spatial gradients in the density, temperature, and flow velocity. The question arises then as to whether, and if so to what extent, the conclusions drawn from Refs. L06 ; G06 may be altered when the new ingredient associated with the presence of the gas phase is accounted for in the theory.

In the first-order approximation, the momentum transport is characterized by the viscosity tensor $\eta_{ijk\ell}$ while the heat flux is characterized by the thermal conductivity tensor $\kappa_{ij}$ and the Dufour-like tensor $\mu_{ij}$. As in the case of dry granular gases, to get explicit analytical results, the steady state conditions are considered and hence, the (scaled) friction coefficient $\gamma^{*}$ (which characterizes the amplitude of the drag force) is given in terms of the (independent) relevant parameters $a^{*}$ and $\alpha$. This contrasts with the results offered in Refs. L06 ; G06 since the transport coefficients are now explicitly obtained as nonlinear functions of both the shear rate and the coefficient of restitution.

For ordinary fluids (elastic collisions), several previous works studied the shear-rate dependence of the thermal conductivity tensor under shear flow. Thus, Evans E91 derived years ago a Green-Kubo formula for the thermal conductivity in a strongly shearing fluid. In a similar way as in the equilibrium case, the thermal conductivity of a shearing steady state is expressed in terms of fluctuations in steady heat flux. This formula was subsequently employed to calculate the shear-rate dependence of the thermal conductivity of a Lennard-Jones fluid via nonequilibrium molecular dynamics simulations methods DE93 . In the context of kinetic theory, an explicit expression of the thermal conductivity tensor was derived G93 ; G95 by solving the Boltzmann equation by means of an expansion around the shear flow state. These analytical results were shown to compare qualitatively well with the computer simulations performed in Ref. DE93 . It must be noted that the calculations carried out in Refs. G93 ; G95 slightly differ from the ones carried out in this paper since the former require an additional external force to reach a steady state with constant pressure and linear shear field. Apart from these papers, a more recent paper SA14 for dry granular gases has determined the thermal conductivity tensor via an expansion around an anisotropic Gaussian distribution function. The authors derived a generalized Fourier law for the granular heat flux where the thermal conductivity is characterized by an anisotropic second rank tensor. A comparison between the results obtained here with those reported before for ordinary G93 ; G95 and granular SA14 sheared gases will be made in sec. VII.

The plan of the paper is as follows. In Sec. II the Boltzmann kinetic equation is introduced and its corresponding balance equations derived. Section III deals with the relevant results derived in the (unperturbed) USF problem by solving the Boltzmann equation by means of Grad’s moment method G49 . In Sec. IV the problem we are interested in is described and the set of coupled linear equations defining the generalized coefficients $\eta_{ijk\ell}$, $\kappa_{ij}$, and $\mu_{ij}$ are provided. Explicit expressions for these shear-rate dependent transport coefficients are then obtained in Sec. V by employing a kinetic model of the Boltzmann equation. The details of the calculations are displayed along several Appendices. The shear-rate dependence of some transport coefficients is illustrated in Sec. VI for different values of the coefficient of restitution. Finally, in Sec. VII the paper is closed with some concluding remarks.

II Boltzmann kinetic equation for monodisperse granular suspensions

We consider a granular suspension of solid particles of mass $m$ and diameter $\sigma$ immersed in a gas of viscosity $\mu_{g}$. Under rapid flow conditions, particles are modeled as a gas of smooth hard spheres or disks with inelastic collisions. The inelasticity of collisions is characterized by a constant (positive) coefficient of normal restitution $\alpha\leq 1$. As said in the Introduction, a simple and usual usual way of modeling the effect of the interstitial gas on the dynamic properties of the solid particles is through the presence of nonconservative external forces. These forces are incorporated into the corresponding Boltzmann kinetic equation of the solid particles. Thus, in the low-density regime, the one-particle velocity distribution function $f(\mathbf{r},\mathbf{v},t)$ of grains obeys the kinetic equation Brilliantov and Pöschel (2004)

[TABLE]

where the Boltzmann collision operator $J\left[{\bf v}|f,f\right]$ is given by

[TABLE]

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

[TABLE]

As mentioned in the Introduction, a simplest way of modeling the fluid-solid interaction force $\mathbf{F}_{\text{fluid}}$ is through the drag force (1) where

[TABLE]

is the mean flow velocity of the solid particles, and

[TABLE]

is the number density of particles. Thus, according to Eqs. (1) and (2), the Boltzmann equation becomes

[TABLE]

where $\Delta\mathbf{U}=\mathbf{U}-\mathbf{U}_{g}$, and ${\bf V}={\bf v}-{\bf U}$ is the peculiar velocity. Note that in the case of very dilute suspensions, $\gamma$ is assumed to be a constant Koch (1990); KS99 ; Wylie et al. (2003).

The macroscopic balance equations for the densities of mass, momentum and energy can be obtained by multiplying Eq. (7) by $1$, $m\mathbf{V}$, and $\frac{1}{2}mV^{2}$, respectively, and integrating over velocity. The result is

[TABLE]

[TABLE]

[TABLE]

Here, $D_{t}\equiv\partial_{t}+\mathbf{v}\cdot\nabla$ is the material derivative,

[TABLE]

is the granular temperature,

[TABLE]

is the pressure tensor,

[TABLE]

is the heat flux, and

[TABLE]

is the cooling rate characterizing the rate of energy dissipated due to collisions.

Notice that the interaction of solid particles with the gas phase is modeled solely by the friction term (1) since the term accounting for the momentum transferred from the gas (bath) to the granular particles (which is modeled by a stochastic force) has been neglected for the sake of simplicity. This stochastic force contributes to the Boltzmann equation (7) with a Langevin-like term of the form $-\frac{1}{2}\xi\partial^{2}f/\partial V^{2}$, where $\xi$ is the strength of the noise term. As said in sec. I, this stochastic term was considered in the complete suspension model proposed in Ref. GTSH12 . For elastic collisions and zero shear rate, the inclusion of the above stochastic term yields the balance equation $\partial_{t}T=-2T\gamma+m\xi$ and so, the Boltzmann equation (7) admits a stable steady equilibrium state. Indeed, it is precisely the condition of admitting an equilibrium state that gives rise to a fluctuation-dissipation theorem M89 fixing the strength of the noise term [i.e., $\xi=2\gamma T/m$, where $T$ is the steady equilibrium temperature]. The omission of this Langevin-like term could be justified if the bath temperature is very low compared to the granular temperature or if the mean flow velocities of solid and gas phases are quite similar GTSH12 .

On the other hand, in spite of the absence of the Langevin-like term, the Boltzmann equation (7) still admits a simple solution in the homogenous state (zero shear rate) for elastic collisions ($\alpha=1$). This solution is given by a time-dependent Maxwellian distribution. For homogeneous states, Eq. (7) becomes

[TABLE]

where an appropriate selection of the frame of reference where the mean flow velocity vanishes ($\mathbf{U}=\mathbf{U}_{g}=\mathbf{0}$) has been chosen. The only relevant balance equation is that of the temperature (10) which reads

[TABLE]

Since $\gamma\equiv\text{const.}$, then the solution to Eq. (16) is simply

[TABLE]

where $T(0)$ is the initial temperature. Under these conditions, it is easy to see that the Boltzmann equation (15) has the solution GSB90 ; PG14

[TABLE]

where $T(t)$ is given by (17). An H-theorem has been also proved GSB90 for the distribution $f_{0}$ in the sense that, starting from any initial condition and in the presence of the viscous drag force $\gamma\mathbf{v}$, the velocity distribution function $f(\mathbf{r},\mathbf{v},t)$ reaches in the long time limit the Maxwellian form (18) with a time-dependent temperature.

Before closing this Section, it is interesting to remark the situations in which the suspension model (7) is expected to provide reliable predictions. As has been previously discussed in several papers models ; TK95 ; SMTK96 ; Wylie et al. (2003), since the form of the Boltzmann collision operator (3) is the same as for a dry granular gas, one expects that the model (7) is appropriate for problems where the stresses applied by the gas phase on particles have only a weak influence on the dynamics of grains. This necessarily requires that the mean-free time between collisions is much shorter than the viscous relaxation time due to the viscous drag force. For other kind of systems (e.g., glass beads in liquid water), one should take into account the influence of the interstitial fluid on the Boltzmann collision operator.

III Simple shear flow problem in dilute granular suspensions

We assume now that the suspension is in steady USF. This state is macroscopically defined by a constant density $n$ and temperature $T$ and the mean velocity $\mathbf{U}$ is

[TABLE]

where $a$ is the constant shear rate. In addition, as usual in uniform sheared suspensions TK95 ; SMTK96 ; ChVG15 , the average velocity of particles follows the velocity of the gas phase and so, $\mathbf{U}=\mathbf{U}_{g}$. In this case, $\Delta{\bf U}=\textbf{0}$ and the Boltzmann equation (7) becomes

[TABLE]

Upon writing Eq. (20) use has been made of the fact that the USF state becomes spatially uniform when one expresses the Boltzmann equation in terms of the peculiar velocity $V_{i}=v_{i}-a_{ij}r_{j}$ DSBR86 . In the USF problem, the heat flux vanishes and the only relevant balance equation is that of the temperature (10). In the steady state and for the geometry of the USF, Eq. (10) reads

[TABLE]

Equation (21) implies that the viscous heating term ($-aP_{xy}>0$) is exactly canceled by the cooling terms arising from viscous friction ($\gamma T$) and collisional dissipation ($\zeta T$). Thus, in stationary conditions, for a given value of $\gamma$, the (steady) temperature is a function of the shear rate $a$ and the coefficient of restitution $\alpha$. Equivalently, one might chose $\gamma$ and $\alpha$ as independent parameters instead of $a$ and $\alpha$. This was the choice made in Refs. TK95 ; SMTK96 ; ChVG15 . Since we are mainly interested here in obtaining the shear-rate dependence of the transport coefficients, the former choice will be considered in this paper. A remarkable point is that a steady state is still possible for suspensions when the collisions are elastic ($\alpha=1$ and so, $\zeta=0$) provided $\gamma=-P_{xy}a/(dp)$, where $p=nT$ is the hydrostatic pressure.

The USF state is non-Newtonian. This can be characterized by generalized transport coefficients measuring the departure of transport coefficients from their Navier-Stokes forms. Thus, one can define a non-Newtonian shear viscosity coefficient $\eta(\alpha,a)$ by

[TABLE]

Moreover, while $P_{xx}=P_{yy}=P_{zz}$ in the Navier-Stokes domain, normal stress differences are present in the USF state.

The elements of the pressure tensor $P_{ij}$ can be obtained by multiplying both sides of Eq. (20) by $mV_{i}V_{j}$ and integrating over velocity. The result is

[TABLE]

where

[TABLE]

So far, the hierarchy (23) is still exact. However, the exact expression of the collision integral $\Lambda_{ij}$ is not known (even for elastic collisions). A good estimate of $\Lambda_{ij}$ can be obtained by using Grad’s approximation to $f$ G49 , namely,

[TABLE]

where

[TABLE]

is the local equilibrium distribution function, and

[TABLE]

is the traceless part of the pressure tensor. When Eq. (25) is substituted into the definition of $\Lambda_{ij}$ and nonlinear terms in $\Pi_{ij}$ are neglected, one gets the result G02

[TABLE]

where

[TABLE]

is an effective collision frequency,

[TABLE]

is the dimensionless cooling rate evaluated in the local equilibrium approximation and

[TABLE]

As will show below, the determination of the collisional moment $\Lambda_{ij}$ by considering only linear terms yields $P_{xx}\neq P_{yy}$ but $P_{yy}=P_{zz}$. This latter identity disagrees with computer simulation results TK95 ; SMTK96 ; ChVG15 . The evaluation of $\Lambda_{ij}$ by retaining all the quadratic terms in the pressure tensor $P_{ij}$ has been recently carried out in Ref. ChVG15 . As expected, the addition of these nonlinear terms allows to evaluate the normal stress differences in the plane normal to the laminar flow (e.g., $P_{yy}-P_{zz}$). However, given that this difference is quite small, the expression (28) can be considered as a reliable approximation. Apart from its simplicity, the linear Grad solution is also essentially motivated by the desire of analytic expressions that show in a clean way the shear-rate dependence of the rheological properties.

Once the collisional moment $\Lambda_{ij}$ is known, the set of coupled equations for $P_{ij}$ can be easily solved. In terms of the reduced shear rate $a^{*}=a/\nu$ and the coefficient of restitution $\alpha$, the expressions for the (scaled) elements $P_{ij}^{*}=P_{ij}/p$ are

[TABLE]

[TABLE]

where $\widetilde{a}=a^{*}/\beta$, and $\chi$ is the real root of the cubic equation

[TABLE]

namely,

[TABLE]

The (scaled) friction coefficient $\gamma^{*}=\gamma/\nu$ is defined as

[TABLE]

In the case of elastic collisions ($\alpha=1$), Eqs. (32)–(36) agree with those obtained GS03 for a thermostatted dilute gas under USF. Moreover, the analytical results given by Eqs. (32)–(36) compare quite well with Monte Carlo simulations of the Boltzmann equation ChVG15 , even for strong inelasticity.

Since $\gamma^{*}\geq 0$, then necessarily $2\beta\chi-\zeta^{*}\geq 0$, according to Eq. (36). This means that, at a given value of the coefficient of restitution, there is a threshold value of the (scaled) shear rate $a_{\text{th}}^{*}$ such that the steady state condition (21) admits a physical solution for $a^{*}\geq a_{\text{th}}^{*}$. This physical solution yields a positive granular temperature and is related to what Sangani et al. SMTK96 call ignited state. The value of $a_{\text{th}}^{*}$ is determined from the condition

[TABLE]

In particular, for elastic collisions, $\zeta^{*}=0$ and so, $a_{\text{th}}^{*}=0$. However, for inelastic collisions, $\zeta^{*}\neq 0$, and $a_{\text{th}}^{*}>0$. Thus, the rheological properties are only well-defined for shear rates beyond the nonvanishing $a_{\text{th}}^{*}$ in the case of granular suspensions ($\alpha\neq 1$). The $\alpha$-dependence of $a_{\text{th}}^{*}$ is plotted in Fig. 1 for $d=2$ and $d=3$. For strong inelasticity, the curves highlight that the granular suspension is in general beyond the Navier-Stokes domain (non-Newtonian regime) since the (reduced) threshold shear rate $a_{\text{th}}^{*}$ is not small in general. Thus, for instance $a_{\text{th}}^{*}\simeq 0.512$ at $\alpha=0.8$ in the physical three-dimensional case.

The fact that for granular suspensions ($\alpha\neq 1$) a steady state is only possible for sufficiently high shear rates can be easily understood from a physical point of view. For $\gamma^{*}=0$, the balance equation (37) establishes an intrinsic connection between the shear field [through the nonlinear function $\chi(\alpha,a^{*})$] and the collisional dissipation [through the cooling rate $\zeta^{*}(\alpha)$] in the system. Thus, the magnitude of the (scaled) shear rate $a^{*}$ is set by the coefficient of restitution $\alpha$. Since $\zeta^{*}\propto 1-\alpha^{2}$, then the cooling rate increases with inelasticity. Moreover, the rheological function $\chi$ increases with increasing $a^{*}$. Consequently, one needs to consider higher values of $a^{*}$ as $\alpha$ decreases to verify the condition (37) and achieve a steady state.

The (reduced) nonlinear shear viscosity $\eta^{*}=\eta/\eta_{0}$ can be easily identified from Eqs. (22) and (33). Here, $\eta_{0}=p/\nu$ is the Navier-Stokes shear viscosity of an ordinary (elastic) gas of hard spheres. The expression of $\eta^{*}$ is given by

[TABLE]

Since $\chi\sim a^{*2/3}$ for very large shear rates, then $\eta^{*}\sim a^{*-4/3}$ and goes to zero in the limit $a^{*}\to\infty$. To illustrate the shear-rate dependence of $\eta^{*}$, Fig. 2 shows the ratio $\eta^{*}(\alpha,a^{*})/\eta^{*}(\alpha,0)$ versus $a^{*}$ for $d=3$ and three different values of the coefficient of restitution $\alpha$. As mentioned before, except for elastic collisions and although $\eta^{*}$ is well defined for shear rates smaller than the threshold value $a_{\text{th}}^{*}$, the curves in Fig.  2 start from the point $a^{*}=a_{\text{th}}^{*}$ for $\alpha\neq 1$. It appears that shear thinning (viscosity decreases with increasing shear rate) is always present, regardless of the value of the coefficient of restitution. We also observe that, at a given value of $a^{*}$, inelasticity inhibits the momentum transport. However, the influence of inelasticity on the (scaled) shear viscosity is not quantitatively significant.

IV Transport coefficients for states close to USF

Let us assume that we perturb the USF by small spatial gradients. This will give rise to new contributions to the momentum and heat fluxes that can be characterized by generalized transport coefficients. Since the system is strongly sheared, the corresponding transport coefficients are highly nonlinear functions of the shear rate. The evaluation of these coefficients is the main objective of the present paper.

As in previous papers L06 ; G06 ; G07 , in order to analyze this problem one has to start from the Boltzmann equation (7) with a general time and space dependence. First, it is convenient to continue using the relative velocity ${\bf V}={\bf v}-{\bf U}_{0}$, where ${\bf U}_{0}={\sf a}\cdot{\bf r}$ is the flow velocity of the unperturbed USF state. As said before, the only nonzero element of the tensor $\mathsf{a}$ is $a_{ij}=a\delta_{ix}\delta_{jy}$. On the other hand, in the perturbed state the true velocity ${\bf U}$ is in general different from ${\bf U}_{0}$ since ${\bf U}={\bf U}_{0}+\delta{\bf U}$, $\delta{\bf U}$ being a small perturbation to ${\bf U}_{0}$. As a consequence, the true peculiar velocity is now ${\bf c}\equiv{\bf v}-{\bf U}={\bf V}-\delta{\bf U}$. In addition, for the sake of simplicity, we also assume that the interstitial gas is not perturbed and hence, $\mathbf{U}_{g}=\mathbf{U}_{0}$. Thus, in the Lagrangian frame moving with velocity $\mathbf{U}_{0}$, the convective operator $\mathbf{v}\cdot\nabla$ can be written as

[TABLE]

where the derivative $\nabla f$ is taken now at constant $\mathbf{V}$. In this case, the Boltzmann equation (7) reads

[TABLE]

The corresponding macroscopic balance equations associated with this disturbed USF state follows from the general equations (10)–(13) when one takes into account that ${\bf U}={\bf U}_{0}+\delta{\bf U}$. The result is

[TABLE]

[TABLE]

[TABLE]

where the pressure tensor ${\sf P}$, the heat flux ${\bf q}$ and the cooling rate $\zeta$ are defined by Eqs. (12)–(14), respectively, with the replacement ${\bf V}\rightarrow{\bf c}$.

Since we are interested here in states close to the USF state, it is assumed that the deviations from the USF state are small and hence, the spatial gradients of $n$, $\delta\mathbf{U}$, and $T$ are small. In this case, Eq. (40) can be solved by means of a generalization of the conventional Chapman-Enskog method CC70 , where the velocity distribution function is expanded around a local shear flow reference state in terms the small spatial gradients of the hydrodynamic fields relative to those of USF. This type of Chapman-Enskog-like expansion has been carried out for elastic gases to obtain the set of shear-rate dependent transport coefficients GS03 ; LD97 in a thermostatted shear flow problem and it has also been employed in the context of dry granular gases L06 ; G06 ; G07 .

The Chapman-Enskog method assumes the existence of a normal solution in which all space and time dependence of the distribution function occurs through a functional dependence of the hydrodynamic fields

[TABLE]

This solution expresses the fact that the space dependence of the shear flow is absorbed in $\mathbf{V}$ and the remaining space and time dependence is through a functional dependence on the fields $A(\mathbf{r},t)$. As in the conventional Chapman-Enskog method, this functional dependence can be made local by an expansion of $f$ in powers of spatial gradients:

[TABLE]

where the reference zeroth-order distribution function corresponds to the USF distribution function but taking into account the local dependence of the density and temperature and the change ${\bf V}\rightarrow{\bf V}-\delta{\bf U}({\bf r},t)$. The successive approximations $f^{(k)}$ are of order $k$ in the gradients of $n$, $T$, and $\delta{\bf U}$ but retain all the orders in the shear rate $a$. This is the main feature of this expansion. In addition, as in previous works GFHY16 , since the friction coefficient $\gamma$ does not induce any flux in the system, it is assumed then to be at least of zeroth order in the gradients. In this paper, only the first order approximation will be considered.

The expansion (45) yields the corresponding expansion for the fluxes and the cooling rate when one substitutes (45) into their definitions (12)–(14):

[TABLE]

[TABLE]

Finally, as in the usual Chapman-Enskog method, the time derivative is also expanded as

[TABLE]

where the action of each operator $\partial_{t}^{(k)}$ is obtained from the hydrodynamic equations (41)–(IV). These results provide the basis for generating the Chapman-Enskog solution to the inelastic Boltzmann equation (40).

IV.1 Zeroth-order approximation

Substituting the expansions (45)–(48) into Eq. (40), the kinetic equation for $f^{(0)}$ is given by

[TABLE]

To lowest order in the expansion the conservation laws are

[TABLE]

[TABLE]

As discussed in previous works L06 ; G06 ; G07 , for given values of $a$, $\gamma$ and $\alpha$, the steady state condition (21) establishes a mapping between the density and temperature so that every density corresponds to one and only one temperature. Since the density $n({\bf r},t)$ and temperature $T({\bf r},t)$ are specified separately in the local USF state, the viscous heating only partially compensates for the collisional cooling and friction viscous dissipation and so, $\partial_{t}^{(0)}T\neq 0$. Consequently, the zeroth-order distribution $f^{(0)}$ depends on time through its dependence on the temperature and the (dimensionless) parameters $a^{*}$, $\gamma^{*}$, and $\alpha$ must be considered as independent parameters for general infinitesimal perturbations around the USF state. The fact that the temperature must be considered as a time-dependent parameter has been already accounted for in previous perturbation solutions around driven non-steady states GMT13 ; GChV13 .

Since $f^{(0)}$ is a normal solution, then

[TABLE]

Upon deriving the last step in Eq. (52) use has been made of the fact that $f^{(0)}$ depends on $\delta{\bf U}$ only through the peculiar velocity ${\bf c}$. Substituting Eq. (52) into Eq. (49) yields the following kinetic equation for $f^{(0)}$:

[TABLE]

The zeroth-order solution leads to ${\bf q}^{(0)}={\bf 0}$ by symmetry. The closed set of equations defining the zeroth-order pressure tensor $\mathsf{P}^{(0)}$ can be obtained from Eq. (53) by taking into account Eq. (28). The result is

[TABLE]

where $\zeta_{0}^{*}\equiv\zeta^{(0)}/\nu$ is defined by Eq. (30).

The steady state solution of Eq. (IV.1) is given by Eqs. (32)– (34). However, for non-steady conditions, in general Eqs. (IV.1) must be solved numerically to get the dependence of the zeroth-order pressure tensor $P_{ij}^{(0)}(T)$ on temperature. In the hydrodynamic regime, it is expected that $P_{ij}^{(0)}$ adopts the form

[TABLE]

where the temperature dependence of the (dimensionless) pressure tensor $P_{ij}^{*}$ is through its dependence on $\gamma^{*}$ and $a^{*}$. Since $\gamma^{*}\propto T^{-1/2}$ and $a^{*}\propto T^{-1/2}$, then

[TABLE]

As it will be shown below, to determine the generalized transport coefficients in the steady state, one needs to know the derivatives $\partial_{\gamma^{*}}P_{ij}^{*}$ and $\partial_{a^{*}}P_{ij}^{*}$ in this state. These derivatives are evaluated in Appendix A. In what follows, $P_{ij}^{(0)}(T)$ will be considered as a known function of $T$.

The shear-rate dependence of the derivatives of the (reduced) pressure tensor with respect to $a^{*}$ and $\gamma^{*}$ in the steady state are illustrated in Fig. 3 for a three-dimensional suspension with elastic collisions ($\alpha=1$). Although we could not analytically prove the identity $\partial_{a^{*}}P_{xy}^{*}=\partial_{\gamma^{*}}P_{xy}^{*}$, numerical results systematically show this result. Since the magnitude of these derivatives is not in general quite small, it appears that their influence on transport cannot be in principle neglected.

IV.2 First-order approximation

The first order approximation is worked out in Appendix B. Only the final results are given here. The velocity distribution function $f^{(1)}$ is

[TABLE]

where the vectors ${\bf X}_{n}$ and ${\bf X}_{T}$ and the tensor ${\sf X}_{u}$ are the solutions of the following set of coupled linear integral equations:

[TABLE]

[TABLE]

[TABLE]

where ${\bf Y}_{n}({\bf c})$, ${\bf Y}_{T}({\bf c})$, and ${\sf Y}_{u}({\bf c})$ are defined by Eqs. (108)–(110), respectively, and $\zeta_{u,k\ell}$ is defined by Eq. (111). An approximate expression of $\zeta_{u,k\ell}$ is given by Eq. (112). In addition, ${\cal L}$ is the linearized Boltzmann collision operator around the USF state, namely,

[TABLE]

Note that, due to the presence of $P_{xy}^{(1)}$ in Eq. (104), the unknown coefficients $\eta_{xyk\ell}$ appear in the quantity $Y_{u,k\ell}$ of Eq. (IV.2). In the particular case of $\gamma^{*}=0$, Eqs. (IV.2)–(IV.2) are consistent with the results derived in Ref. G06 for dry granular gases.

With the distribution $f^{(1)}$ determined by Eq. (57), the first-order corrections to the fluxes are given by

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

Upon writing Eqs. (62)–(66) use has been made of the symmetry properties of $X_{n,i}$, $X_{T,i}$, and $X_{u,ij}$.

In the absence of gas phase ($\gamma^{*}=0$), for $a^{*}=0$ and $\alpha=1$, the conventional Navier-Stokes constitutive equations for ordinary gases are reobtained, namely,

[TABLE]

[TABLE]

Here, $\eta_{0}=p/\nu$ and $\kappa_{0}=d(d+2)\eta_{0}/2(d-1)m$ are the expressions of the shear viscosity and thermal conductivity coefficients, respectively, of an ordinary gas of disks ($d=2$) or hard ($d=3$) spheres CC70 . In the absence of shear rate, the expressions of the Navier-Stokes coefficients of a granular suspension have been recently derived in Ref. GFHY16 .

In general, the set of generalized transport coefficients $\eta_{ijk\ell}$, $\kappa_{ij}$, and $\mu_{ij}$ are nonlinear functions of the coefficient of restitution $\alpha$, the reduced shear rate $a^{*}$ and the reduced friction coefficient $\gamma^{*}$. The anisotropy induced in the system by the shear flow gives rise to new transport coefficients, reflecting broken symmetry. Since $P_{ij}^{(1)}$ is a symmetric and traceless tensor, then the viscosity tensor $\eta_{ijk\ell}$ is symmetric and traceless in $ij$, namely,

[TABLE]

The heat flux is expressed in terms of a thermal conductivity tensor $\kappa_{ij}$ and a Dufour-like tensor $\mu_{ij}$. While the diagonal elements of both tensors can be interpreted as generalizations of the Navier-Stokes transport coefficients, the off-diagonal elements $\kappa_{xy}$, $\kappa_{yx}$, $\mu_{xy}$, and $\mu_{yx}$ are generalizations of Burnett coefficients that, for small shear rates, are proportional to $a^{*}$. In addition, because of symmetry reasons, the off-diagonal elements $xz$, $zx$, $yz$, and $zy$ of the tensors $\kappa_{ij}$ and $\mu_{ij}$ are identically zero. This is consistent with Eqs. (IV.2) and (IV.2). The above behavior implies that if the thermal gradient is parallel to the $z$ axis ($\nabla T\parallel\widehat{\boldsymbol{z}}$), then $\mathbf{q}^{(1)}\parallel\widehat{\boldsymbol{z}}$, while if $\nabla T\perp\widehat{\boldsymbol{z}}$, then $\mathbf{q}^{(1)}\perp\widehat{\boldsymbol{z}}$. Similarly, many of the elements of the viscosity tensor $\eta_{ijk\ell}$ are zero. For instance, if the only nonzero velocity gradient is $\partial\delta U_{x}/\partial z$, then $P_{ij}^{(1)}=P_{xz}^{(1)}(\delta_{ix}\delta_{jz}+\delta_{jx}\delta_{iz})$.

IV.3 Steady state conditions

As in the case of dry granular gases ($\gamma^{*}=0$), the evaluation of the transport coefficients $\eta_{ijk\ell}$, $\kappa_{ij}$ and $\mu_{ij}$ for general unsteady conditions is quite intricate. This is due essentially to the fact that the temperature dependence of the velocity moments of the distribution $f^{(0)}$ must be numerically determined. Thus, since we want to get analytical expressions for those coefficients, the present study is limited to steady state conditions. This means that the relation (21) is considered at the end of the calculations. In this state, the (scaled) shear rate $a^{*}$ is coupled to the (reduced) friction coefficient $\gamma^{*}$ and the coefficient of restitution $\alpha$ so that, only two of the three parameters are independent. Here, as alluded to in Sec. II, $a^{*}$ and $\alpha$ are chosen as the independent (input) parameters of the problem. This allows us to independently assess the influence of shearing and inelasticity on momentum and heat transport. This contrasts with the analysis of dry granular gases G06 where both $a^{*}$ and $\alpha$ are considered as dependent parameters in the steady state.

Since the relation (21) holds in the steady state, the first term on the left hand side of the integral equations  (IV.2)–(IV.2) vanishes. In this case, these equations become

[TABLE]

[TABLE]

[TABLE]

In Eqs. (IV.3)–(IV.3) it is understood that all the quantities are evaluated in the steady state. Moreover, the dependence of $P_{ij}^{(0)}$ on the temperature $T$ is given by Eq. (56) while the dependence of $P_{ij}^{(0)}$ on the density can be written as

[TABLE]

V Results from a BGK-like kinetic model

Needless to say, the explicit form of the generalized transport coefficients $\eta_{ijk\ell}$, $\kappa_{ij}$, and $\mu_{ij}$ requires to solve the integral equations (IV.3)–(IV.3). Apart from the mathematical difficulties embodied in the Boltzmann collision operator ${\cal L}$, it is quite apparent that the fourth-degree velocity moments of the zeroth-order distribution $f^{(0)}$ are also needed to determine the heat flux transport coefficients $\mu_{ij}$ and $\kappa_{ij}$. Although these moments could in principle be determined from Grad’s moment method by including them in the trial distribution (25), their evaluation would be an intricate task.

A possible alternative could be the use of the so-called inelastic Maxwell models BCG00 ; BK00 ; EB02 , i.e., models for which the collision rate is independent of the relative velocity of the two colliding particles. The use of these models allows to obtain the velocity moments of the Boltzmann collision operator without the explicit knowledge of the velocity distribution function. This was the route followed in Ref. G07 to determine the shear-rate dependent transport coefficients in a dry granular sheared gas. However, apart from the difficulties associated with the evaluation of the fourth-degree moments and their derivatives, the results obtained for inelastic Maxwell models G07 show significant discrepancies from those obtained for inelastic hard spheres G06 .

Therefore, as in the previous study carried out for dry granular gases G06 , a model kinetic equation of the Boltzmann equation is considered to achieve explicit results. As for elastic collisions, the idea is to replace the true Boltzmann collision operator with a simpler, more tractable operator that retains the most relevant physical properties of the Boltzmann operator. Here, we consider a kinetic model BDS99 based on the well-known Bhatnagar-Gross-Krook (BGK) GS03 for ordinary gases where the operator $J[f,f]$ is note

[TABLE]

Here, $\nu$ is the effective collision frequency defined by Eq. (29), $f_{\text{M}}(\mathbf{c})$ is the Maxwellian distribution (26), $\beta$ is given by Eq. (31) and $\zeta$ is the cooling rate. It is easy to see that the BGK model yields the same expressions for the pressure tensor in the steady USF state than those derived from Grad’s method [Eqs. (32)–(34)]. Moreover, the fourth-degree velocity moments obtained from the BGK model compare quite well with Monte Carlo simulations ChVG15 ; AS05 of the Boltzmann equation. This confirms again the reliability of kinetic models to evaluate the velocity moments of the true Boltzmann equation GS03 .

In the perturbed USF problem, Eqs. (IV.3)–(IV.3) still apply with the replacements

[TABLE]

in the case of $X_{n,i}$ and $X_{T,i}$ and

[TABLE]

in the case of $X_{u,ij}$. In the above equations, $\zeta^{(0)}$ is the zeroth-order approximation to $\zeta$ which is given by Eq. (30). With the changes (75) and (76) all the generalized transport coefficients can be easily evaluated from Eqs. (IV.3)–(IV.3). Details of these calculations are given in Appendix C.

VI Shear-rate dependence of the generalized transport coefficients

The general results derived in the previous sections clearly show that the dependence of the generalized transport coefficients on both $a^{*}$ and $\alpha$ is quite complex. Since the main goal of the present paper is to assess the shear-rate dependence of $\eta_{ijk\ell}$, $\kappa_{ij}$, and $\mu_{ij}$ for given values of $\alpha$, we illustrate here this dependence for some relevant elements of the above tensors by two different values of $\alpha$: $\alpha=1$ (ordinary suspensions) and $\alpha=0.8$ (granular suspensions). Moreover, a three-dimensional system ($d=3$) is considered in all the plots and hence, $a_{\text{th}}^{*}\simeq 0.512$ for $\alpha=0.8$.

To analyze the shear-rate dependence of the transport coefficients, it is convenient first to introduce the dimensionless coefficients $\eta_{ijk\ell}^{*}\equiv\eta_{ijk\ell}/\eta_{0}$, $\kappa_{ij}^{*}\equiv\kappa_{ij}/\kappa_{0}$, and $\mu_{ij}^{*}\equiv n\mu_{ij}/T\kappa_{0}$. Here, $\eta_{0}=p/\nu$ and $\kappa_{0}=((d+2)/2)nT/(m\nu)$ are the elastic values of the shear viscosity and thermal conductivity coefficients, respectively, for a dilute gas given by the BGK kinetic model.

VI.1 Viscosity tensor

The (reduced) elements of the viscosity tensor $\eta_{ijk\ell}^{*}$ are determined by solving the set of algebraic equations (D). There are in principle two classes of terms GS11 . Class I is made of those coefficients $\eta_{ijk\ell}^{*}$ with $(k,\ell)=\left\{(xx),(xy),(yx),(yy),(zz)\right\}$. The complementary class II is constituted by coefficients with $(k,\ell)=\left\{(xz),(yz),(zx),(zy)\right\}$. Of course, class II (as well as the elements $\eta_{ijzz}^{*}$ of class I) is meaningless in the two-dimensional case ($d=2$).

A careful analysis of the set of algebraic equations shows that the coefficients of the form $\eta_{xzk\ell}^{*}$ and $\eta_{yzk\ell}^{*}$ vanish in class I. In addition, the coefficients of the form $\eta_{xxk\ell}^{*}$ of class II include the first-order contribution to the cooling rate $\zeta_{u,ij}$. However, they obey a set of homogeneous algebraic equations whose solution is the trivial one for arbitrary values of $a^{*}$. A similar behavior is expected for the coefficients of the form $\eta_{xyk\ell}^{*}$, $\eta_{yyk\ell}^{*}$, and $\eta_{zzk\ell}^{*}$. Thus, one can conclude that all the above elements of class II vanish.

The remaining elements of class II are independent of the derivatives $\partial_{a^{*}}P_{ij}^{*}$ and $\partial_{\gamma^{*}}P_{ij}^{*}$. Some of them are given by

[TABLE]

[TABLE]

where the nonzero elements of the pressure tensor $P_{yy}^{*}$ and $P_{xy}^{*}$ are defined by Eqs. (32) and (33), respectively, and the nonlinear shear viscosity $\eta^{*}$ is defined by Eq. (38). The expressions of the remaining elements of class II can be obtained from Eqs.  (D) and (D). Their forms are very long and will be omitted here. Figure 4 shows the dependence of two elements of class I ($\eta_{xzxz}^{*}$ and $\eta_{yzzx}^{*}$) for $\alpha=1$ and 0.8. These two coefficients measure the presence of non-zero values of $P_{xz}$ and $P_{yz}$ due to perturbations of the form $\partial\delta U_{x}/\partial z$ and $\partial\delta U_{z}/\partial x$, respectively. It is quite apparent that, at a given value of $\alpha$, the largest impact of the shear rate on momentum transport occurs on $P_{xz}$. We also observe that $\eta_{xzxz}^{*}$ exhibits a shear-thinning effect more pronounced than that of the nonlinear shear viscosity $\eta^{*}$, as expected from Eq. (77). In addition, the influence of collisional dissipation is very tiny in both generalized transport coefficients.

Finally, the expressions for the non-zero elements of class I contain the derivatives $\partial_{a^{*}}P_{ij}^{*}$ and $\partial_{\gamma^{*}}P_{ij}^{*}$. Those expressions are much more involved than those of class II. In order to illustrate their shear-rate dependence, we consider here the set of coefficients $\left\{\eta_{xxxy}^{*},\eta_{xyxy}^{*},\eta_{yyxy}^{*},\eta_{zzxy}^{*}\right\}$. Note that $\eta_{xxxy}^{*}=-(\eta_{yyxy}^{*}+\eta_{zzxy}^{*})$. In addition, the algebraic equations defining those coefficients show that $\eta_{yyxy}^{*}=\eta_{zzxy}^{*}$. This result is a consequence of the linear version of Grad’s moment method that yields $P_{yy}^{*}=P_{zz}^{*}$. As said before, recent Monte Carlo simulations of granular suspensions ChVG15 have shown that the second normal stress difference is different from zero although its value is very small. The shear-rate dependence of the elements $\eta_{ijxy}^{*}$ is plotted in Fig. 5. The coefficients $\eta_{xyxy}^{*}$ and $\eta_{yyxy}^{*}$ measure the deviations of $P_{xy}$ and $P_{yy}$, respectively, from their unperturbed USF values due to perturbations of the form $\partial\delta U_{x}/\partial y$. While the coefficient $\eta_{xyxy}^{*}$ decreases in general with $a^{*}$ (except for high shear rates), the coefficient $\eta_{yyxy}^{*}$ exhibits clearly a non-monotonic shear-rate dependence regardless the value of the coefficient of restitution. It is also interesting to note that $\eta_{xxxy}^{*}$ is always negative.

VI.2 Thermal conductivity and Dufour-like tensors

The evaluation of the heat flux transport coefficients $\Delta_{ij}\equiv\left\{\kappa_{ij}^{*},\mu_{ij}^{*}\right\}$ is much more involved than that of the shear viscosity tensor $\eta_{ijk\ell}^{*}$. As Eqs. (139) and (140) show, in the steady state the set of transport coefficients $\Delta_{ij}$ also depends on the derivatives of the fourth-degree moments of the USF with respect to $\gamma^{*}$ and $a^{*}$. The evaluation of these derivatives is in general a quite tedious task that can be accomplished by following the steps devised in the Appendix A code .

As mentioned before, we have $\Delta_{xz}=\Delta_{zx}=\Delta_{yz}=\Delta_{zy}=0$ according to the linear shear flow (19). Therefore, there are five nonzero elements of the (scaled) tensors $\Delta_{ij}$: the three diagonal ($\Delta_{xx}$, $\Delta_{yy}$, and $\Delta_{zz}$) and the two off-diagonal elements ($\Delta_{xy}$ and $\Delta_{yx}$). The algebraic equations (D) and (D) also show that the anisotropy induced by the shear flow yields the properties $\Delta_{xx}\neq\Delta_{yy}\neq\Delta_{zz}$ and $\Delta_{xy}\neq\Delta_{yx}$.

To illustrate the shear-rate dependence of the coefficients $\Delta_{ij}$, we consider here the elements $\Delta_{zz}\equiv\left\{\kappa_{zz}^{*},\mu_{zz}^{*}\right\}$ and $\Delta_{xy}\equiv\left\{\kappa_{xy}^{*},\mu_{xy}^{*}\right\}$. The first set of coefficients measures the heat flux along the direction orthogonal to the shearing plane. The second set of coefficients provides information on cross-effects in the thermal conduction since $\kappa_{xy}^{*}$ and $\mu_{xy}^{*}$ measure the transport of energy parallel to the flow direction due to a thermal gradient along the velocity gradient. Figures 6–9 show the generalized coefficients $\kappa_{yy}^{*}$, $\mu_{yy}^{*}$, $\kappa_{xy}^{*}$, and $\mu_{xy}^{*}$ versus $a^{*}$ for $\alpha=1$ and 0.8. We observe first that the deviations of these coefficients with respect to their equilibrium values is significant, regardless of the collisional dissipation. This means that the impact of shear flow on heat transport is in general significant in a region of shear rates where shear thinning is quite important (see Fig. 2). Regarding the diagonal element $\kappa_{yy}^{*}$, it is quite apparent from Fig. 6 that this coefficient decreases with $a^{*}$ in the region of shear rates considered. A similar behavior is found in Fig. 7 for $\mu_{zz}^{*}$ when the collisions are inelastic ($\alpha\neq 1$). On the other hand, for elastic collisions, $\mu_{zz}^{*}$ first increases with $a^{*}$ for small shear rates and then it decreases with the shear rate. In any case, for elastic collisions, the magnitude of $\mu_{zz}^{*}$ is much smaller than that of $\kappa_{zz}^{*}$. Thus, for practical purposes, one can neglect the contribution to the heat flux coming from the term proportional to the density gradient when the collisions are elastic. In accordance with the above results, we conclude that in general the shear flow inhibits the transport of energy along the direction orthogonal to the velocity gradient (vorticity direction). With respect to the influence of $\alpha$ on both generalized coefficients, it appears that the effect of inelasticity is more important in the case of $\mu_{zz}^{*}$ than in the case of $\kappa_{zz}^{*}$.

The absolute values of the off-diagonal elements $\kappa_{xy}^{*}$ and $\mu_{xy}^{*}$ are plotted in Figs. 8 and 9, respectively. As said before, these coefficients measure cross effects in the energy transport. This cross coupling does not appear in the linear regime since the first-order contribution to the heat flux $q_{x}^{(1)}$ is at least of Burnett order (i.e., proportional to $a^{*}\partial_{x}T$). It is quite apparent that the element $\kappa_{xy}$ is negative and its magnitude presents a non-monotonic dependence with $a^{*}$ since it increases first with the shear rate (in the region of small shear rates), reaches a maximum and then decreases with increasing $a^{*}$. This behavior is much more evident in the case of elastic collisions. Regarding the coefficient $\mu_{xy}^{*}$, we observe that it is always negative for granular suspensions ($\alpha\neq 1$) and its magnitude is very small for elastic collisions. Recall that the coefficient $\mu_{ij}^{*}$ vanishes when $\alpha=1$ for vanishing shear rates. As in the case of the diagonal elements, the effect of inelasticity on heat transport is more noticeable for $\mu_{xy}^{*}$ than for $\kappa_{xy}^{*}$.

Finally, it is important to remark that the qualitative shear-rate dependence of $\kappa_{zz}^{*}$ and $\kappa_{xy}^{*}$ obtained here for elastic collisions (ordinary fluids) agrees with the one observed years ago by Daivis and Evans DE93 in molecular dynamics simulations of a thermostatted shear-flow state.

VII Concluding remarks

The influence of gas phase on the transport properties of solid particles under USF has been studied in this paper. In the low-density regime, a viscous drag force term for the interstitial fluid has been incorporated into the Boltzmann kinetic equation to account for the effect of the former on the dynamics of grains. The physical situation is such that the granular suspension is in a state that deviates from the USF by small spatial gradients. Since the system is subjected to a strong shear flow and is not restricted to nearly elastic spheres, the corresponding transport coefficients characterizing momentum and heat transport are nonlinear functions of both the shear rate and the coefficient of restitution. The explicit determination of the above coefficients has been the main objective of the present contribution. The search for such expressions has been prompted by previous results L06 ; G06 obtained for dry granular gases (i.e., in the absence of the viscous drag force). Here, the problem is revisited by considering the effect of the gas phase on transport properties.

Assuming that the USF state is slightly perturbed, the Boltzmann equation (7) has been solved by means of a Chapman-Enskog-like expansion. The new feature of this expansion is that the (local) shear flow distribution is employed as the reference state instead of the usual (local) equilibrium distribution CC70 or the (local) homogeneous cooling state BDKS98 ; GD99 . As already noted in previous works L06 ; G06 ; GMT13 ; GChV13 , since the zeroth-order derivative $\partial_{t}^{(0)}T$ is in general different from zero, the reference base state is not stationary. This fact introduces technical difficulties in the implementation of the perturbation scheme. Thus, in order to get explicit results, the steady-state condition (21) is considered at the end of the calculations. In this state, the (reduced) shear rate $a^{*}$ and the coefficient of restitution $\alpha$ are coupled to the (scaled) friction coefficient $\gamma^{*}$, so that the former two are the relevant parameters of the problem.

To first order of the expansion, the momentum and heat fluxes are given by Eqs. (62) and (63), respectively, where the generalized transport coefficients $\eta_{ijk\ell}$, $\kappa_{ij}$, and $\mu_{ij}$ are defined in terms of the solutions of the set of coupled integral equations(IV.3)–(IV.3). However, since the solution of the above integral equations is in general quite a complex problem, the BGK-like kinetic model (74) has been employed to obtain the explicit shear-rate dependence of the above set of transport coefficients. Although the kinetic model (74) can be considered as a crude representation of the true Boltzmann equation, it gives the same results for the rheological properties as those derived from the Boltzmann equation by means of Grad’s moment method. Given that those theoretical predictions compare quite well with Monte Carlo simulations ChVG15 , it is expected that the results provided by the kinetic model are accurate even for conditions of practical interest, such as strong dissipation and/or large shear rates.

As expected, there are many new transport coefficients in comparison to the case of states close to equilibrium (for ordinary gases) or states near the homogeneous cooling state (for dry granular gases). Here, for the sake of illustration, the shear-rate dependence of some relevant elements of the viscosity tensor, the thermal conductivity tensor, and the Dufour-like tensor have been studied. More specifically, Figs. 4 and 5 show the (reduced) elements $\eta_{xzxz}^{*}$, $\eta_{yzzx}^{*}$, and $\eta_{ijxy}^{*}$, respectively, Figs. 6 and 7 show the diagonal elements $\kappa_{zz}^{*}$ and $\mu_{zz}^{*}$, respectively, and Figs. 8 and 9 show the off-diagonal elements $\kappa_{xy}^{*}$ and $\mu_{xy}^{*}$, respectively. It is apparent that in general the deviation of these coefficients from their equilibrium values (i.e., for $a^{*}=0$ and $\alpha=1$) is quite significant. In addition, the influence of collisional dissipation on transport is much more significant for the heat flux transport coefficients than for the coefficients associated with the pressure tensor.

As said in the Introduction, for ordinary fluids ($\alpha=1$), the thermal conductivity tensor of a thermostatted shear-flow state was determined years ago from the BGK G93 and Boltzmann G95 kinetic equations. The physical situation corresponds to a perturbed steady USF state with $\delta\mathbf{U}=\mathbf{0}$, $p=nT\equiv\text{const.}$ and $\nabla T\neq 0$. Under these conditions, one needs to add an external field that exactly compensates for the increase or decrease of momentum due to the term $\nabla\cdot\mathsf{P}$ GS03 . The addition of this external field affects the value of the thermal conductivity tensor and hence, the situation studied in Refs. G93 ; G95 slightly differs from the one analyzed in the present paper. On the other hand, in order to make a comparison with these previous results G93 ; G95 , one considers particular perturbations such that $\nabla p=0$ and so, $\nabla\ln n=-\nabla\ln T$. Therefore, the heat flux (63) obeys the generalized Fourier’s law

[TABLE]

Figure 10 shows the transport coefficient $\widetilde{\kappa}_{zz}\equiv\kappa_{zz}^{*}-\mu_{zz}^{*}$ versus $a^{*}$ for $\alpha=1$. We observe that the previous predictions made for thermostatted shear flow states from the BGK G93 and Boltzmann G95 equations compare qualitatively well with the results obtained here for arbitrary perturbations. However, at a more quantitative level, it seems that the impact of shear flow on energy transport is more significant in the situation analyzed in this paper than those studied in Refs. G93 ; G95 .

In the case of dry granular gases ($\gamma^{*}=0$ but $\alpha\neq 1$), Saha and Alam SA14 have determined the heat flux of a two-dimensional granular gas under USF. The results were obtained by solving the Boltzmann equation by means of a perturbation expansion around an anisotropic Gaussian distribution. This distribution was employed years ago by Jenkins and Richman JR88 to obtain the rheological properties of USF via Grad’s moment method. The corresponding constitutive relation for the heat flux derived in Ref. SA14 can be written as

[TABLE]

where $\Pi_{ij}$ is the deviatoric or traceless part of the pressure tensor defined by Eq. (27). In Eq. (80), $\kappa_{ij}$ is identified as the thermal conductivity tensor and $\Xi_{ij}$ is a tensor quantifying the contribution to the heat flux coming from the gradient of the deviatoric stress $\Pi_{ij}$. As expected, the tensors $\kappa_{ij}$ and $\Xi_{ij}$ are nonlinear functions of the coefficient of restitution $\alpha$. It appears first that Eq. (80) disagrees with the constitutive relation (63) derived here for the heat flux. On the other hand, in an attempt to make a comparison with the theoretical results obtained in Ref. SA14 for the thermal conductivity tensor, Fig. 11 shows the dependence of the (reduced) coefficients $\kappa_{yy}^{*}$ and $\kappa_{xy}^{*}$ on $\alpha$ for a dry two-dimensional granular gas. Notice that, in order to get analytical expressions, the theoretical results of Ref. SA14 plotted in Fig. 11 were derived by considering terms up to super-Burnett order (i.e., third order in the shear rate). We observe that the $\alpha$-dependence of the diagonal element $\kappa_{yy}^{*}$ is qualitatively different from the one predicted in Ref. SA14 since while in the latter theory $\kappa_{yy}^{*}$ decreases with increasing inelasticity, the opposite happens here. Although a more qualitative agreement is found for the magnitude of $\kappa_{xy}^{*}$, both theoretical results exhibit significant quantitative discrepancies for strong inelasticity. The differences between both theories at the level of the thermal conductivity tensor $\kappa_{ij}$ could be in part due to the different form of the constitutive relation for the heat flux derived in Ref. SA14 . In addition, while the results obtained in the latter work were obtained by solving the Boltzmann equation up to super-Burnett order, the theoretical predictions made in the present paper are based on an exact solution of the BGK-like kinetic model. It would be convenient to perform computer simulations for $\kappa_{ij}$ to check the reliability of the above theories for strong inelasticities.

The explicit results reported in this paper can be useful for studying different problems First, as done in Ref. G06 , an important application is to perform an stability analysis of the hydrodynamic equations with respect to the USF state. This analysis will allow us to identify the conditions for stability in terms of both the shear rate and the coefficient of restitution. Another interesting and challenging problem is to extend the present results by considering the general Langevin-like model proposed in Ref. GTSH12 . This will allow us to provide additional refinements of the predictions obtained here so that, a closer comparison with direct numerical simulations of granular suspensions could be performed. Finally, it would be also relevant to extend the analysis made here for a monodisperse granular suspension to the intriguing and important subject of polydisperse suspensions. A good starting point for this achievement could be the suspension model introduced in Ref. KG13 . Work along the above lines will be carried out in the near future.

Acknowledgements.

I am grateful to Dr. Mariano López de Haro for a critical reading of the manuscript. The present research has been supported by the Spanish Government through grant No. FIS2016-76359-P, partially financed by FEDER funds and by the Junta de Extremadura (Spain) through Grant No. GR15104.

Appendix A Derivatives of the zeroth-order velocity moments with respect to $\gamma^{*}$ and $a^{*}$ in the steady state

The derivatives of the zeroth-order velocity moments with respect to $\gamma^{*}$ and $a^{*}$ in the steady state are determined in this appendix. We start with the pressure tensor $P_{ij}^{(0)}$, whose elements obey Eq. (53). In dimensionless form, Eq. (53) is given by

[TABLE]

where

[TABLE]

and upon deriving Eq. (A) use has been made of the relation (56). Let us consider the elements $P_{yy}^{*}=P_{zz}^{*}$ and $P_{xy}^{*}$. From Eq. (A), one gets

[TABLE]

[TABLE]

The goal here is to evaluate the derivatives $P_{yy,\gamma}^{*}$, $P_{yy,a}^{*}$, $P_{xy,\gamma}^{*}$ and $P_{xy,a}^{*}$ at the steady state. This state is defined by the condition (21). To get these derivatives, we differentiate first Eqs. (A) and (A) with respect to $a^{*}$ and take then the steady-state limit. The result is

[TABLE]

[TABLE]

where here it is understood that all the terms are evaluated at the steady state. To close the problem, we differentiate then Eqs. (A) and (A) with respect to $\gamma^{*}$ and take the steady-state limit with the result

[TABLE]

[TABLE]

The set of nonlinear algebraic equations (A)–(A) can be numerically solved for given values of $a^{*}$ and $\alpha$. In the dry limit case ($\gamma^{*}=0$), the solution to Eqs. (A) and (A) can be written as

[TABLE]

where $P_{xy,a}^{*}$ is the real root of the cubic equation

[TABLE]

Equations (89) and (A) agree with previous results L06 ; G06 derived for a dry granular gas of inelastic hard spheres.

The corresponding derivatives of the fourth-degree velocity moments of the distribution $f^{(0)}$ with respect to $\gamma^{*}$ and $a^{*}$ in the steady state are also needed to determine the generalized coefficients $\kappa_{ij}$ and $\mu_{ij}$ associated with the first-order contribution to the heat flux. To evaluate these derivatives, the BGK kinetic model (74) is considered. The velocity moments of the distribution $f^{(0)}$ are defined as

[TABLE]

These moments verify the equation

[TABLE]

where $k\equiv k_{1}+k_{2}+k_{3}$, $\widetilde{a}\equiv a^{*}/\beta$, $\widetilde{\zeta}_{0}\equiv\zeta_{0}^{*}/\beta$, $\widetilde{\gamma}\equiv\gamma^{*}/\beta$, $\widetilde{\xi}=\widetilde{\gamma}+\widetilde{\zeta}_{0}/2$, and $N_{k_{1},k_{2},k_{3}}$ are the velocity moments of $f_{\text{M}}$. In the steady state, $\widetilde{\xi}=\chi$ where $\chi$ is given by Eq. (35). As in the case of the pressure tensor, the derivative $T\partial_{T}M_{k_{1},k_{2},k_{3}}^{(0)}$ can be written as

[TABLE]

where we have introduced the shorthand notation

[TABLE]

[TABLE]

In dimensionless form, Eq. (A) reads

[TABLE]

where $M_{k_{1},k_{2},k_{3}}^{*}\equiv n^{-1}(m/2T)^{k/2}M_{k_{1},k_{2},k_{3}}$, and

[TABLE]

if $k_{1}$, $k_{2}$, and $k_{3}$ are even, being zero otherwise. Equation (A) provides the expressions of the reduced moments $M_{k_{1},k_{2},k_{3},s}^{*}$ in the steady state (e.g., when $\frac{2}{d}\widetilde{a}P_{xy}^{*}+2\widetilde{\gamma}+\widetilde{\zeta}_{0}=0$).

In order to evaluate the derivatives $M_{k_{1},k_{2},k_{3},\gamma}^{*}$ and $M_{k_{1},k_{2},k_{3},a}^{*}$ in the steady state, we differentiate with respect to $\gamma^{*}$ and $a^{*}$, respectively, both sides of Eq. (A) and then take the steady state condition (21). As an illustration, let us consider the moment $M_{040}^{*}$ which obeys the equation

[TABLE]

From Eq. (A), in the steady state, one gets the identities

[TABLE]

[TABLE]

The solution to the set of linear algebraic equations (A) and (A) gives the derivatives $M_{040,\gamma}^{*}$ and $M_{040,a}^{*}$ in terms of $a^{*}$ and $\alpha$. Proceeding in a similar way, all the derivatives of the fourth-degree velocity moments with respect to both $a^{*}$ and $\gamma^{*}$ can be analytically computed in the steady state.

Appendix B First-order approximation

The kinetic equation for the first-order distribution $f^{(1)}$ is

[TABLE]

The velocity dependence on the right side of Eq. (B) can be obtained from the macroscopic balance equations to first order in the gradients. They are given by

[TABLE]

[TABLE]

[TABLE]

where $\rho=mn$ is the mass density,

[TABLE]

and

[TABLE]

Use of Eqs. (102)–(104) in Eq. (B) yields

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

According to the symmetry properties of $f^{(1)}$, the only nonzero contribution to $\zeta^{(1)}$ comes from the term proportional to the tensor $\nabla_{i}\delta U_{j}$. Thus,

[TABLE]

An estimation of $\zeta_{u,ij}$ has been made in Ref. G06 for a three-dimensional system ($d=3$). The result is

[TABLE]

where $\Pi_{ij}^{*}\equiv\Pi_{ij}/nT$. Of course, when $\alpha=1$, then $\zeta_{u,ij}=0$.

The solution to Eq. (B) has the form

[TABLE]

Note that in Eq. (110) the coefficients $\eta_{ijk\ell}$ are defined through Eq. (64). The coefficients $X_{n,i}$, $X_{T,i}$, and $X_{u,ij}$ are functions of the peculiar velocity ${\bf c}$ and the hydrodynamic fields. In addition, there are contributions from the time derivative $\partial_{t}^{(0)}$ acting on the temperature and velocity gradients given by

[TABLE]

[TABLE]

Substituting Eqs. (112), (B), and (115) into Eq. (B) and identifying coefficients of independent gradients, one finally gets the set of coupled linear integral equations

[TABLE]

[TABLE]

[TABLE]

Upon writing Eqs. (B)–(B), use has been made of the property

[TABLE]

where in the last step we have taken into account that $X$ depends on $\delta{\bf U}$ through ${\bf c}={\bf V}-\delta{\bf U}$.

Appendix C Kinetic model results in the steady USF state

In this appendix, the steady state solution to the BGK-like kinetic model (74) in the steady (unperturbed) USF is briefly analyzed. In this case, $\delta\mathbf{U}=\mathbf{0}$ and so $\mathbf{c}=\mathbf{V}$. In the steady state, $f_{s}(\mathbf{V})$ verifies the kinetic equation

[TABLE]

where here $\zeta$ has been approximated by its Maxwellian approximation $\zeta^{(0)}$ given by Eq. (30). Let us introduce the velocity moments of $f_{s}$ as

[TABLE]

According to the symmetry of the USF distribution $f_{s}$, the only nonvanishing moments correspond to even values of $k_{1}+k_{2}$ and $k_{3}$. In this case, after some algebra, one gets

[TABLE]

where the reduced moments $M_{k_{1},k_{2},k_{3}}^{*}$ are given by

[TABLE]

It is easy to see that the second-degree velocity moments of the BGK model coincide with those obtained from the Boltzmann equation by using Grad’s method, Eqs. (32)–(34).

Appendix D Generalized transport coefficients

The results derived from the BGK-like kinetic model (74) considered to determine the generalized transport coefficients $\eta_{ijk\ell}$, $\kappa_{ij}$ and $\mu_{ij}$ are provided in this appendix. The equations defining the generalized transport coefficients in the BGK model can be obtained from Eqs. (IV.3)–(IV.3) with the replacements (75) and (76):

[TABLE]

[TABLE]

[TABLE]

In order to get the transport coefficients $\kappa_{ij}$, $\mu_{ij}$, and $\eta_{ijk\ell}$, it is convenient to introduce the general velocity moments

[TABLE]

[TABLE]

[TABLE]

These moments provide the explicit forms of the generalized transport coefficients of the perturbed USF problem. To determine them, Eqs. (D)–(D) are multiplied by $c_{x}^{k_{1}}c_{y}^{k_{2}}c_{z}^{k_{3}}$ and integrated over velocity. After some algebra, one achieves

[TABLE]

[TABLE]

[TABLE]

where $M_{k_{1},k_{2},k_{3}}^{(0)}$ are the moments of the zeroth-order distribution function $f^{(0)}$,

[TABLE]

[TABLE]

[TABLE]

and we have introduced the quantities

[TABLE]

[TABLE]

[TABLE]

The integrals (136)–(138) can be computed with the result

[TABLE]

[TABLE]

[TABLE]

Here, $M_{k_{1},k_{2},k_{3}}^{*}$ are the reduced moments of the distribution $f^{(0)}$ defined by Eq. (C). They depend on $n$ and $T$ through their dependence on $\gamma^{*}$ and $a^{*}$. In the steady state, $M_{k_{1},k_{2},k_{3}}^{*}$ is given by Eq. (C) while the derivatives $M_{k_{1},k_{2},k_{3},\gamma}^{*}$ and $M_{k_{1},k_{2},k_{3},a}^{*}$ can be obtained by following the procedure described in appendix A.

The solution to Eqs. (D)–(D) can be written as

[TABLE]

[TABLE]

[TABLE]

where $\widetilde{\omega}_{T}\equiv\omega_{T}/(\nu\beta)$. In Eqs. (D)–(D), we recall that in the steady state the parameter $\widetilde{\xi}=\chi$ is given by Eq. (35), and $\gamma^{*}=\beta\chi-\frac{1}{2}\zeta_{0}^{*}$. The expressions of the generalized transport coefficients $\kappa_{ij}$, $\mu_{ij}$ and $\eta_{ijk\ell}$ can be obtained from Eqs. (D)–(D).