Exact Results for Non-Newtonian Transport Properties in Sheared Granular Suspensions: Inelastic Maxwell Models and BGK-Type Kinetic Model

The Boltzmann kinetic equation for dilute granular suspensions under simple (or uniform) shear flow (USF) is considered to determine the non-Newtonian transport properties of the system. In contrast to previous attempts based on a coarse-grained description, our suspension model accounts for the real collisions between grains and particles of the surrounding molecular gas. The latter is modeled as a bath (or thermostat) of elastic hard spheres at a given temperature. Two independent but complementary approaches are followed to reach exact expressions for the rheological properties. First, the Boltzmann equation for the so-called inelastic Maxwell models (IMM) is considered. The fact that the collision rate of IMM is independent of the relative velocity of the colliding spheres allows us to exactly compute the collisional moments of the Boltzmann operator without the knowledge of the distribution function. Thanks to this property, the transport properties of the sheared granular suspension can be exactly determined. As a second approach, a Bhatnagar–Gross–Krook (BGK)-type kinetic model adapted to granular suspensions is solved to compute the velocity moments and the velocity distribution function of the system. The theoretical results (which are given in terms of the coefficient of restitution, the reduced shear rate, the reduced background temperature, and the diameter and mass ratios) show, in general, a good agreement with the approximate analytical results derived for inelastic hard spheres (IHS) by means of Grad’s moment method and with computer simulations performed in the Brownian limiting case (m/mg→∞, where mg and m are the masses of the particles of the molecular and granular gases, respectively). In addition, as expected, the IMM and BGK results show that the temperature and non-Newtonian viscosity exhibit an S shape in a plane of stress–strain rate (discontinuous shear thickening, DST). The DST effect becomes more pronounced as the mass ratio m/mg increases.


Introduction
A very usual way of assessing the effect of the surrounding fluid on the dynamics properties of solid particles is through an effective fluid-solid force [1][2][3][4].In some models, this force is simply proportional to the velocity particle (Stokes linear drag law) [5][6][7][8][9][10][11][12].This type of force attempts to mimic the energy dissipated by grains due to their friction on the interstitial viscous gas.A more sophisticated model [13][14][15][16][17][18][19][20][21] incorporates also a Langevin-like stochastic term that accounts for the energy transferred to grains due to their "collisions" with particles of the background gas.However, although this coarse-grained approach has provided reliable results in the past, it would be desirable to consider a suspension model that takes into account the real collisions between grains and particles of the surrounding (molecular) gas.This sort of suspension model (which has been inspired in a previous work of Biben et al. [22]) has been recently proposed [23].In this model, granular particles are assumed to be sufficiently rarefied so that, they do not disturb the state of the molecular (background) gas.As a consequence, the interstitial gas may be treated as a thermostat at the temperature T g .Moreover, although the concentration (mole fraction) of grains is quite small, apart from the elastic collisions between solid and molecular gas particles one has to consider the inelastic collisions between grains themselves.This model can be useful to analyze transport properties in particle-laden suspensions [24] where very dilute particles (like aerosols) are immersed in a fluid (like air).
The rheological properties of a granular suspension under simple (or uniform) shear flow (USF) has been recently determined [25].In contrast to previous attempts [5][6][7][8][9][10][11][12]14,15,19,21], the results obtained in Ref. [25] have been derived from the collisional model proposed in Ref. [23].On the other hand, a limitation of these results is that they have been approximately obtained by employing Grad's moment method [26], namely, a method based on the truncation of a series expansion of the velocity distribution function in (orthogonal) Sonine polynomials.The use of this approximate method is essentially motivated by the form of the collision rate for inelastic hard spheres (IHS) appearing inside the Boltzmann collision operator.The collision rate for IHS is proportional to the magnitude of the normal component of the relative velocity of the two spheres that are about to collide.This velocity dependence of the collision rate for IHS prevents the possibility of deriving exact expressions for the transport properties in the USF problem, even in the case of elastic collisions.
A possible way of overcoming the technical difficulty of the hard-sphere kernel is to consider the so-called inelastic Maxwell models (IMM).As for the conventional Maxwell molecules [27,28], the collision rate of IMM is independent of the relative velocity of the two colliding spheres [29].The use of IMM instead of IHS opens up the possibility of getting exact analytical results of the Boltzmann equation in some specific non-equilibrium situations, like the USF.In particular, the knowledge of the collisional moments of the Boltzmann equation for IMM [30,31] enables a clear exploration of the impact of inelasticity on the non-Newtonian transport properties of the granular suspension without introducing uncontrolled approximations.
Another possible alternative for obtaining exact results is to consider a kinetic model that retains the relevant physical properties of the Boltzmann collision operator but turns out to be more tractable than the true kinetic equation.This route has been widely employed in the past in the case of molecular dilute gases [28] where it has been shown that several exact solutions in far from equilibrium states agree very well with Monte Carlo simulations of the Boltzmann equation.Here, as in previous works [17], we will consider a kinetic model for granular suspensions [32] to complement the theoretical expressions obtained from the Boltzmann equation for IMM.Since this kinetic model is based on the well known Bhatnagar-Gross-Krook (BGK) model [33] for molecular gases, we will referred to it as a BGK-type kinetic model.
The objective of this paper is to determine the rheological properties of granular particles immersed in a bath of elastic hard spheres under USF.At a macroscopic level, the USF is characterized by constant number densities for solid and gas particles, a uniform temperature, and a (common) linear velocity profile U g,x = U x = ay, where a is the constant shear rate.Here, U g and U denote the mean flow velocities of the molecular and granular gases, respectively.Since we are interested here in the steady state where the system admits a non-Newtonian hydrodynamic description, an external thermostat force (proportional to the peculiar velocity) must be introduced to keep constant the temperature T g of the molecular gas.
The use of IMM as well as a BGK-type kinetic model allows us to exactly compute the rheological properties of the granular suspension.These properties are expressed as nonlinear functions of the (reduced) shear rate a * = a/γ (where γ is a drift coefficient characterizing the friction of solid particles on the viscous gas), the coefficient of restitution, the (reduced) background temperature T * g , and the diameter σ/σ g and mass m/m g ratios.Here, σ g and m g are the diameter and mass of the particles of the molecular gas, respectively, while σ and m are the diameter and mass of the solid particles, respectively.As occurs for IHS [25], our results show that the kinetic granular temperature and the non-Newtonian viscosity exhibit a discontinuous shear thickening (DST) effect for sufficiently large values of the mass ratio m/m g .In fact, in the Brownian limiting case (m/m g → ∞) the expressions of the rheological properties derived here reduce to those previously obtained [17] from a coarse-grained description based on the Fokker-Planck operator.This agreement justifies the use of this latter approach to analyze the DST effect in dilute granular suspensions [34].
Apart from the transport properties (which are related with the second-degree velocity moments), the explicit forms of the higher degree velocity moments as well as the velocity distribution function of the granular gas have been also obtained from the BGK model.This is one of the main advantages of using a kinetic model instead of the true Boltzmann equation.Our results show in particular that the fourth-degree moments of the distribution function also exhibit a DST effect.With respect to the velocity distribution function, as expected we find that its distortion from equilibrium is more significant as both the mass m/m g and diameter σ/σ g ratios depart from 1.In addition, a comparison between the BGK results and numerical solutions of the Boltzmann equation from the direct simulation Monte Carlo (DSMC) method [35] for IHS shows a generally good qualitative agreement between both approaches.
The plan of the paper is as follows.The Boltzmann kinetic equation for a granular gas immersed in a bath of elastic hard spheres under USF is presented in Section 2. The balance equations for the temperatures of the molecular and granular gases are also displayed.Section 3 deals with the calculations carried out for IMM of the rheological properties of the granular suspension.While a shear thinning effect is always found for the nonlinear shear viscosity of the molecular gas, the corresponding shear viscosity of the granular gas exhibits a DST effect for sufficiently large values of the mass ratio m/m g .The results derived from the BGK-type kinetic model are provided in Section 4 while a comparison between the theoretical results obtained for IHS, IMM and BGK model is displayed in Section 5 for several systems.Our results highlight a good agreement for the rheology between the three different approaches.Moreover, theoretical results obtained from IMM and BGK model are also compared against computer simulations in the Brownian limit (m/m g → ∞) showing a good agreement.The paper is closed in Section 6 with some concluding remarks.

Boltzmann kinetic equation for sheared granular suspensions
We consider a set of solid particles (granular gas) of mass m and diameter σ which are immersed in a solvent (molecular gas) constituted by particles of mass m g and diameter σ g .As usual, the granular gas is modeled as a gas of hard disks (d = 2) or spheres (d = 3) with inelastic collisions.In the simplest model, the inelasticity of collisions is characterized by a constant (positive) coefficient of normal restitution α ≤ 1, where α = 1 refers to elastic collisions.On the other hand, collisions between solid particles and particles of the molecular gas are elastic.We also assume that the number density of grains is much smaller than that of solvent so that, the state of the latter is not perturbed by the presence of the former.In these conditions, we can treat the molecular gas as a bath or thermostat at the temperature T g (once the parameters of the system, specifically the shear rate, have been set).Moreover, although the granular gas is sufficiently rarefied, we take into account the collisions among grains in its corresponding Boltzmann kinetic equation.
We assume that the system (granular particles plus solvent) is subjected to USF.As said in the Introduction section, this state is characterized by constant densities n g and n, uniform temperatures T g and T, and by a (common) linear profile of the x component of the flow velocities along the y axis: a being the constant shear rate.Here, n g , U g , and T g are the number density, the mean flow velocity, and the temperature, respectively, of the molecular gas.In terms of its one-particle velocity distribution function f g (r, v; t), these hydrodynamic fields are defined as where V = v − U is the peculiar velocity.Note that in Eq. ( 4) the Boltzmann constant k B = 1.We will take this value throughout the paper for the sake of simplicity.In addition, in Eqs. ( 1)-(3), n, U, and T denote the number density, the mean flow velocity, and the (granular) temperature, respectively, of the granular gas.They are defined as {n, nU, dnT} = dv 1, v, mV 2 f (v). ( Since the only spatial gradient present in the USF problem is the shear rate, the pressure tensor of the molecular gas, and the pressure tensor of the granular gas are the relevant fluxes in the problem.They provide information on the transport of momentum across the system.Our main target is to determine P g and P for arbitrary shear rates.One of the main advantages of the USF at a microscopic level is that it becomes a spatially homogeneous state when the velocities of the particles are referred to a Lagrangian frame moving with the linear velocity U i = a ij r j .In this new frame and in the steady state, the distribution functions of the molecular and granular gases adopt the form In addition, as the state of the solvent is not perturbed by the solid particles, the temperature T g in the USF state increases in time due to the viscous heating term −aP xy > 0. Thus, as usual in nonequilibrium molecular dynamics simulations [36], an external nonconservative force (thermostat) is introduced in the molecular gas to achieve a stationary state.Among the different possibilities, for simplicity, a force proportional to the particle velocity (Gaussian thermostat) of the form F g = −m g ξV is considered in this paper.The parameter ξ is chosen to be a function of the shear rate by the condition that T g reaches a constant value in the long time limit.Analogously, the granular gas is also subjected to this kind of Gaussian thermostat (i.e., F = −mξV), where ξ is the same quantity for the solvent and the solid particles.Under the above conditions, in the low-density regime, the distribution function f g (V) of the molecular gas obeys the nonlinear (closed) Boltzmann equation while the distribution function f (V) of the granular gas obeys the kinetic equation Here, J g [ f g , f g ] and J[ f , f ] are the nonlinear Boltzmann collision operators for the molecular and granular gases, respectively, and J BL [ f , f g ] is the linear Boltzmann-Lorentz collision operator [37,38].
The balance equations for T g and T can easily obtained by multiplying both sides of Eqs. ( 9) and ( 10) by m g V 2 and mV 2 , respectively, and integrating over velocity.The results are −aP g,xy = dξ p g , ( 11) where p g = n g T g and p = nT are the hydrostatic pressures of the molecular and granular gases, respectively, and the partial production rates ζ and ζ g are defined, respectively, as The cooling rate ζ gives the rate of kinetic energy loss due to inelastic collisions between particles of the granular gas.It vanishes for elastic collisions (α = 1).The term ζ g gives the transfer of kinetic energy between the particles of the granular gas and the solvent.This term vanishes when the granular and molecular gases are at the same temperature (T g = T).Equation ( 11) implies that in the steady state the viscous heating term (−aP g,xy > 0) is exactly balanced by the heat extracted in the gas by the external thermostat.On the other hand, since ζ g can be positive or negative, Eq. (12) implies that in the steady state the term −aP xy − (d/2)pζ g is exactly compensated for the cooling terms arising from collisional dissipation (ζ p) and the thermostat term (ξ p).
The USF state is in general a non-Newtonian state characterized by shear-rate dependent transport coefficients.In particular, one can define the non-Newtonian shear viscosity of the molecular gas as Analogously, the non-Newtonian shear viscosity of the granular gas is given by In addition, beyond the Navier-Stokes domain, normal stress differences are expected in the USF.This means that P g,xx = P g,yy = P g,zz and P xx = P yy = P zz .
It is quite evident that the evaluation of the rheological properties of the molecular and granular gases requires the knowledge of the pressure tensors P g and P. The non-zero elements of these tensors can be obtained by multiplying by m g VV and mVV both sides of Eqs. ( 9) and (10), respectively, and integrating over V.However, to achieve explicit forms for P g and P, one has to compute the collisional moments In the case of IHS, the collisional moments A g , B, and C cannot be exactly computed.As said in the Introduction section, a good estimate of them for IHS has been made in Ref. [25] by means of Grad's moment method [26].This method is based on the expansion of the distributions f g (V) and f (V) in a complete set of orthogonal polynomials; the coefficients being the corresponding velocity moments of those distributions.The above expansion generates an infinite hierarchy of moment equations that must be truncated at a given order.This truncation allows one to arrive to a closed set of coupled equations for the velocity moments that can be recursively solved.Thus, since the results derived in Ref. [25] for the rheological properties of molecular and granular gases are approximated, it is convenient to revisit the problem and determine the exact expressions for the non-Newtonian transport properties of the granular suspension.To achieve such exact forms, two independent approaches will be considered in this paper: (i) the Boltzmann kinetic equation for IMM and (ii) a BGK-type kinetic model for IHS.This task will be carried out in the next two Sections.

Rheology from inelastic Maxwell models
In this Section, we will consider IMM, namely, a collisional model where the collision rate of the two colliding spheres are independent of their relative velocity.In this case, the Boltzmann collision operator J g [ f g , f g ] 1 of the molecular gas can be written as [28,33] where is the total solid angle in d dimensions and ν M g is an independent-velocity collision frequency.In Eq. ( 18), the primes on the velocities denote the initial values {V ′′ 1 , V ′′ 2 } that lead to {V 1 , V 2 } following a binary collision: The effective collision frequency ν M g can be seen as a free parameter of the model to be chosen to get agreement with the properties of interest of the original Boltzmann equation for IHS.For instance, to correctly capture the velocity dependence of the original IHS collision rate, we can assume that the IMM collision rate is proportional to √ T g .In the context of IMM, the inelastic Boltzmann collision operator while the Boltzmann-Lorentz collision operator The relationship between 20) is while in Eq. ( 21) is where In addition, as in the case of the collision frequency ν M g , the collision frequencies ν M and ν M 0 for granular-granular and granular-molecular collisions, respectively, can be chosen to optimize the agreement with the results derived from IHS.We will chose them later.
As mentioned in previous works on IMM [16,30,39], the main advantage of computing the collisional moments of the Boltzmann operator for Maxwell models (both elastic and inelastic models) is that they can be exactly provided in terms of the velocity moments of the distribution functions without the explicit knowledge of the latter.This property has been exploited to compute the second, third, and fourth-degree collisional moments of IMM for monocomponent [30] and multicomponent 1 This is a simple version of the Boltzmann collision operator for Maxwell molecules [31,40] granular gases.The exact knowledge of the second-degree collisional moments allow us to get exact expressions for the rheological properties of the molecular and granular gases.Let us evaluate separately the rheology of both gases.

Rheological properties of the molecular gas
The pressure tensor of the molecular gas is defined by Eq. ( 6).To obtain the non-zero elements of this tensor, one multiplies both sides of the Boltzmann equation ( 9) by m g V i V j and integrates over velocity.The result is where p g = n g T g is the hydrostatic pressure of the molecular gas.Upon obtaining Eq. ( 25), use has been made of the result [30] The (reduced) elements of the pressure tensor P * g = P g /(n g T g ) can be easily obtained from Eq. ( 25).They are given by Here, we have introduced the quantities The constraint P * g,xx + (d − 1)P * g,yy = d leads to a cubic equation relating ξ and a: The real root of Eq. ( 29) gives the shear-rate dependence of ξ( a).It is given by [28,41] Comparison between the results derived here for P * g,ij with those recently [25] obtained for IHS by means of Grad's moment method [26] shows that both results are identical if the effective collision frequency ν M g is given by Henceforth, we will take the choice (31) for ν M g .From Eqs. (27), one can identify the (dimensionless) non-Newtonian shear viscosity η * g = ν M g η g /p g = −P * g,xy / a and the (dimensionless) normal stress difference Ψ * g = P * g,xx − P * g,yy as Note that the results derived here for Maxwell molecules yield P g,yy = P g,zz .This result contrasts with the one obtained for hard spheres by numerically solving the Boltzmann equation by means of the DSMC method [35] where it has been shown that P g,yy = P g,zz .However, the difference P g,yy − P g,zz found in the simulations is in general quite small [9].It is also important to remark that in the case of Maxwell molecules there is an exact equivalence between the description with and without the drag force −mξV.Nevertheless, for non-Maxwell molecules this type of force does not play a neutral role in the transport properties of the system [42].The shear-rate dependence of η * g and Ψ * g is plotted in Fig. 1 for a three-dimensional system (d = 3).As expected, the nonlinear viscosity η * g decreases with increasing the (reduced) shear rate a (shear thinning effect).The opposite effect is observed for the normal stress difference function Ψ * g since it increases with the shear rate.Figure 1 also highlights the excellent agreement found between the theoretical results for IMM with those obtained by numerically solving the Boltzmann equation for hard spheres from the DSMC method [35].

Rheological properties of the granular gas
As in the case of the molecular gas, the rheology of the granular gas can be also determined by multiplying both sides of Eq. ( 10) by mV i V j and integrating over V.After some algebra, one achieves the result where use has been made of the results [30,31,40] In Eq. (33), is the (reduced) cooling rate for the granular gas, χ = T/T g is the temperature ratio, and The partial cooling rate ζ g can be exactly obtained from Eq. ( 34) as where is the ratio of the mean square velocities of granular and molecular gas particles.The forms ( 36) and ( 38) can be employed to fix the values of the free parameters ν M and ν M 0 .They are chosen under the criterion that ζ and ζ g of IMM are the same as that of IHS of diameters σ and σ 0 .In this latter case, the above cooling rates are estimated by using Grad's approximation [25].In this approximation, where σ = (σ + σ g )/2.Equations ( 38), (40), and (41) yield the identities To compare with the rheological properties of IHS [25], it is convenient at this level of the description to identify the friction (or drift) coefficient γ appearing in the Brownian limiting case (m/m g → ∞) when the molecular gas is at equilibrium.In fact, this limiting case is the situation considered when one employs a coarse-grained approach [1][2][3][4] to assess the impact of the interstitial gas on the dynamics properties of grains.In this limiting case, the expression (35) of the collisional moment C ij reduces to where P * ij = P ij /(nT g ) and we have taken into account that in the Brownian limit µ g → m g /m and (1 + θ)/θ → 1 in the expression (42) of ν M 0 .The form of C ij derived in Ref. [17] by replacing the Boltzmann-Lorentz collisional operator (21) by the Fokker-Planck operator is Comparison between Eqs. ( 43) and ( 45) allows us to identify γ for IMM as The expression (46) for the friction coefficient γ for IMM is the same as the one obtained for IHS [23,25].
We are now in conditions to obtain the nonzero elements of the (reduced) pressure tensor P * ij .From Eqs. ( 33) and (46), one gets the equation where and Here, is the solid volume fraction of the granular gas, T * g = T g /mσ 2 γ 2 , and upon deriving Eq. ( 48) use has been made of the identity As occurs for the rheology of the molecular gas, Eq. (47) shows that the diagonal elements of the pressure tensor P * ij orthogonal to the shear plane xy are equal to P * yy (i.e., P * yy = P * zz = • • • = P * dd ).This implies that the xx element is given by P * xx = dχ − (d − 1)P * yy .The yy and xy elements of the (reduced) pressure tensor can be written as where Note that the elements of the pressure tensor P * g,yy and P * g,xy of the molecular gas must be expressed in terms of the (reduced) shear rate a * and the (reduced) thermostat parameter ξ * .For doing it, one has to take into account the relationships between a and ξ with a * and ξ * , respectively.They are given by a = (γ/ ν M g )a * and ξ = (γ/ ν M g )ξ * , where The equation defining the temperature ratio χ can be easily derived from Eq. ( 47) as From Eqs. ( 52) and (57), one gets finally a * in terms of the parameter space of the system: As happens in the case of IHS [25], the temperature ratio χ cannot be expressed in Eq. ( 58) as an explicit function of the (reduced) shear rate and the remaining parameters of the system.On the other hand, for given values of the parameter space Ξ ≡ χ, α, σ/σ g , m/m g , φ, T * g , the temperature ratio can be implicitly determined from the physical solution to Eq. (58).

Brownian limit
Before illustrating the shear-rate dependence of the rheological properties of the molecular gas for arbitrary values of the mass ratio m/m g , it is convenient to check the consistency of the present results with those derived in Ref. [17] for IMM by using the Fokker-Planck operator (44).This consistency is expected to apply in the Brownian limit m/m g → ∞.In this limiting case, at a given value of the (reduced) shear rate a * , θ → ∞, γ/ν M g ∝ m g /m → 0, a ∝ m g /m → 0, ξ ∝ a 2 ∝ (m g /m) 2 → 0, and ξ * = ξ( ν M g /γ) ∝ m g /m → 0. Consequently, P * g,ij = δ ij and the expressions of P * yy , P * xy , and a * are Equations ( 59) and (60) are consistent with Eqs. ( 32), (33), and (35) of Ref. [17].It is important to note that, to assess consistency with the Fokker-Planck results, the size ratio has been kept constant or proportional to the mass ratio so that ξ * → 0.

Rheology from a BGK-type kinetic model of the Boltzmann equation
We consider in this Section the results derived for the USF from a BGK-type kinetic model of the Boltzmann equation [32,43].In the problem for the granular suspension considered here, one has to replace the true Boltzmann operators J g [ f g , f g ], J[ f , f ], and J BL [ f , f g ] by simpler relaxation terms that retain the relevant physical properties of those operators but are more tractable than the true kinetic equations.As in the case of IMM, let us determine separately the rheological properties of the molecular and granular gases by starting from these kinetic models.

Rheological properties of the molecular gas
In the case of the molecular gas, the Boltzmann collision operator J g [ f g , f g ] is replaced by the conventional BGK kinetic model [33,44]: where ν g is an effective velocity-independent collision frequency and f M g is the Maxwellian distribution Thus, according to Eq. ( 9), the velocity distribution function f g (V) obeys the BGK kinetic equation The nonzero elements of the pressure tensor P g,ij can be easily obtained from Eq. ( 63) by multiplying both sides of this equation by m g V i V j and integrating over V.The BGK expressions of the (reduced) elements of the pressure tensor P * g,ij are given by Eqs. ( 27) with the replacement ν M g → ν g .As a consequence, the results derived from the BGK equation for the rheological properties agree with those obtained from the Boltzmann equation for IHS when

Rheological properties of the granular gas
In the case of the molecular gas, we consider the kinetic model proposed by Vega Reyes et al. [32] for granular mixtures where the relaxation terms for the collision operators J[ f , f ] and J BL [ f , f g ] are defined, respectively, as In Eqs. ( 65) and (66), is the Maxwellian distribution of the granular gas, and the reference distribution function f g (V) is taken to be the same as in the well-known Gross and Krook (GK) model for molecular (elastic) mixtures [45], i.e., In Eqs. ( 65) and (66), the quantities ν ′ , ν, and T are chosen to optimize the agreement with some properties of interest of the Boltzmann equation for IHS.The usual way of obtaining the above parameters is to ensure that the kinetic model reproduces the collisional transfer equations of momentum and energy for elastic collisions (α = 1).However, since U = U g in the USF, we only have one constraint (the one associated with the transfer of energy) instead of two, so that T and ν admit more than one form.Here, we propose a choice (see Appendix A for more technical details) that leads to an excellent agreement with the results obtained for IHS from Grad's moment method [25].More specifically, we take the following values of T and ν: It still remains to completely define the model to chose the effective collision frequency ν ′ .It is defined here to reproduce the collisional moment of the Boltzmann equation for IHS when one takes Grad's trial distribution for f [25].This leads to the expression Therefore, the BGK kinetic equation for the sheared granular gas is given by where ν and ν ′ are defined by Eqs. ( 70) and (72), respectively, and the Maxwellian distribution f g is given by Eq. ( 69) with T = T g .The possibility of determining all the velocity moments of the distribution function is likely one of the main advantages of employing a kinetic model instead of the true Boltzmann equation.In the USF problem, it is convenient to define the general velocity moments Although we are here essentially interested in the three-dimensional case, we will compute the velocity moments for d = 3 and d = 2.Note that for hard disks (d = 2), k 3 = 0 in Eq. ( 74) since the z-axis is meaningless.The hierarchy of moment equations can be obtained by multiplying Eq. ( 73) by V k 1 x V k 2 y V k 3 z and integrating over V.The result is where In Eq. ( 75), for hard spheres (d = 3) if k 1 , k 2 , and k 3 are even, being zero otherwise.For hard disks (d = 2), if k 1 and k 2 are even, being zero otherwise.The solution to Eq. ( 75) can be written as (see Appendix B of Ref. [16] for some details) (79) The nonzero elements of the pressure tensor P ij can be easily obtained from Eq. ( 79).In dimensionless form, the BGK-expressions of the elements of P * ij = P ij /nT g are The (steady) temperature ratio χ = T/T g can be obtained from the constraint P * xx + (d − 1)P * yy = dχ.This yields the implicit equation where a * = a/γ, Here, ξ is given by Eq. ( 30) and we recall that T * g = T g /mσ 2 γ 2 .

Brownian limit
As in the case of IMM, it is quite interesting to consider the Brownian limiting case m/m g → ∞.In this case, P * g,ij = δ ij , µ → 1, θ → ∞, and ξ * → 0. Thus, following similar steps as those made for IMM, one gets the expressions

Velocity distribution of the granular gas
Apart from obtaining all the velocity moments, the use of kinetic models allow us in some cases to explicitly determine the velocity distribution function f (V).The BGK-type equation (73) reads where a = a/(ν ′ + ν), λ = λ/(ν ′ + ν), and The hydrodynamic solution to Eq. ( 88) can be formally written as The action of the velocity operators in Eq. ( 90) on an arbitrary function g(V) is [16] The explicit form of the one-particle velocity distribution function can be explicitly obtained when one takes into account in Eq. ( 90) the action of the velocity operators given by Eqs. ( 91) and ( 92).The result can be written as where Here, we have introduced the tensor a ij = aδ ix δ jy .

Comparison between IMM and BGK results
In Sections 3 and 4, we made use of the Boltzmann equation for IMM and the BGK-type kinetic model to investigate the shear-rate dependence of rheological properties in a sheared granular suspension.These properties are expressed in terms of the coefficient of restitution α, the reduced background temperature T * g , and the diameter σ/σ g and mass m/m g ratios.Additionally, there exists a residual dependence on density through the volume fraction φ.To avoid that, one could for instance reduce the shear rate using the effective collision frequencies ν M (T) or ν(T).However, for consistency with simulations and considering the background temperature T g as a known quantity, we opted to employ γ(T g ) as the reference frequency.
Given that in this Section the second-degree moments of the distribution function are compared with molecular dynamics (MD) simulations for IHS in the Brownian limiting case [14], we set fixed values of T * g = 1 and φ = 0.0052 for subsequent analysis.The selection of T * g as a free parameter imposes a constraint between the diameter σ/σ g and mass m/m g ratios [25]: This relation ensures convergence of results to those obtained via the Fokker-Planck equation as m/m g → ∞ since ξ * → 0. Furthermore, since we want to recover the results obtained in Ref. [25] derived from Grad's method, we take n/n g = 10 −3 (rarefied granular gas).
The second-degree moments expressed through the reduced temperature χ, non-Newtonian shear viscosity η * , and the normal stress difference Ψ * are plotted in Fig. 2 for α = 0.9 and 1.Here, Ψ * = P * xx − P * yy = dχ − dP * yy .Equations ( 52) and (80) provide analytical expressions for IMM and BGK-type kinetic model, respectively, from which rheological properties are illustrated.Notably, there is nearly perfect agreement between Grad's solution for IHS as obtained in Ref. [25] and both IMM and BGK-type results for any mass ratio, highlighting the ability of relatively simple models to capture essential properties of granular suspensions.
In particular, a DST transition characterized by S-shaped curves becomes more pronounced as the mass ratio m/m g increases.Specifically, the non-Newtonian shear viscosity η * exhibits a discontinuous transition (at a certain value of a * ) which intensifies as the particles of the granular gas become heavier.The theoretical results are validated with MD simulations [14] in the Brownian limiting case (m/m g → 0), showing generally good agreement despite slight discrepancies in the transition zone.Simulations suggest a sharper transition, likely due to the absence of molecular chaos in highly non-equilibrium situations.To address this, DSMC simulations for IHS are performed in the same limit, showing good agreement with theoretical results and further emphasizing a more pronounced transition.This phenomenon is likely attributable to a sudden growth of higher-order moments causing the distribution function to strongly deviate from the reference Maxwellian approximation (67).Some technical details of the application of the DSMC method are available in the supplementary material of Ref. [23] [https://doi.org/10.1017/jfm.2022.410].
The simplicity of the BGK and IMM models enables exploration beyond second-degree moments.Accordingly, we utilize the BGK-type kinetic equation to compute the fourth-degree moments.Although similar calculations could be performed in the case of IMM, we opt to omit them due to their extensive analytical effort.Additionally, drawing insights from the Fokker-Planck model [17] and dry granular gases [39], we anticipate potential divergences of the moments derived from IMM m/m g =10 3 10 4 and the normal stress difference Ψ * (c) as functions of the (reduced) shear rate a * for two different values of the coefficient of restitution α: 0.9 (left panel) and 1 (right panel).The graphs represent four distinct mass ratio values m/m g : 10 3 (yellow lines), 10 4 (blue lines), 10 6 (red lines), and the Brownian limit (black lines).Here, T * g = 1, d = 3, and φ = 0.0052.The solid lines correspond to the IMM results, the dashed lines are the BGK-like results, and the dotted lines refer to Grad's solution for IHS.Symbols denote computer simulation results performed in the Brownian limit: circles refer to the DSMC data obtained in this paper for IHS while squares are MD results obtained in Ref. [14] for IHS.
under certain shear rate conditions.We focus our efforts on calculating the following fourth-degree moments Thus, in terms of the generic moments M k 1 ,k 2 ,k 3 and according to the expression (79), the moments M * 4|0 and M * 2|xy are given by m/m g =10 3  4|0 (a * )/M * 4|0 (0) as a function of a * for α = 0.9 and 1 and three values of the mass ratio.We observe that variations in the mass ratio do not significantly alter the trends observed in the Brownian limiting case [17].An abrupt transition in the higher-order moments is evident within a small region of a * .Specifically, the kurtosis M * 4|0 increases with the mass ratio m/m g until it converges to the value obtained in the Brownian limit.Consistent with the conclusions drawn in Ref. [25], an increase in the mass of the particles of the granular gas results in an elevation of the granular temperature.Consequently, non-equipartition accentuates and moves the suspension away from equilibrium, leading to an increase in kurtosis as the distribution function deviates from its Maxwellian form.Regarding the influence of collisional dissipation, we observe that the effect of α on M * 4|0 remains relatively discrete.Figure 3(b) illustrates the shear-rate dependence of the (reduced) moment M * 2|xy .This moment vanishes in the absence of shear rate (a * = 0).Similar conclusions to those made for the moment M * 4|0 can be drawn.Theoretical predictions for the fourth-degree moments are compared against DSMC simulations for IHS conducted in this paper in the Brownian limiting case.A qualitative agreement is observed, although simulations suggest a sharper transition.Some quantitative discrepancies are noticeable, which are mainly disguised by the scale.To assess the reliability of the BGK-type results, we focus on the region 0 < a * < 1, where all the fourth-degree velocity moments of IMM are well-defined functions of the shear rate.In addition, non-Newtonian effects are still significant in the range of values of the (reduced) shear rate a * ≤ 1.To this purpose, Fig. 4 shows the (reduced) fourth-degree moments M * 4|0 (a * )/M * 4|0 (0) and M * 2|xy (a * ) for α = 0.7 and 1.These moments are also illustrated as obtained for IMM in the Brownian limit [17].It is worth noting that the results derived in Ref. [17] stem from considering an effective force modeling the interstitial gas, diverging from the limit of a Boltzmann-Lorentz operator modeled by a BGK-type equation as m/m g → ∞.Consequently, since DSMC simulations employ the exact Fokker-Planck operator, they perfectly align with the IMM results, while discrepancies emerge when compared with the BGK-type results.It is noteworthy that the BGK-type model slightly overestimates the deviation from the Newtonian situation (a * = 0), a phenomenon also observed for molecular gases [46].Moreover, non-Newtonian effects become apparent even at low values of a * .
0) (a) and −M * 2|xy (a * ) (b) as functions of the (reduced) shear rate a * obtained from the BGK-type equation for two different values of the coefficient of restitution α: 0.7 (solid lines) and 1 (dashed lines).The graphs represent four distinct mass ratio values m/m g : 5 (yellow lines), 10 (blue lines), 50 (red lines), and the Brownian limit (black lines).Here, T * g = 1, d = 3, and φ = 0.0052.Symbols refer to the DSMC results for IHS in the Brownian limit (squares for α = 1 and circles for α = 0.7).The green lines are the IMM results as obtained in Ref. [16] in the Brownian limit.
Finally, in Fig. 5, the ratio R x ) is plotted for a * = 1 and four different values of the mass ratio.Here, the marginal distribution ϕ x (c x ) is given by Eq. (95).It is evident that the deviation from equilibrium (R x = 1) becomes more significant as the mass ratio m/m g increases.Moreover, a comparison between theory and DSMC simulations reveals some disagreement in the BGK-type solution.Although the relative difference of these discrepancies is relatively small (it is about 8%), this contradicts what has been observed in Ref. [9], where good agreement between the BGK solution and DSMC data is shown in the region of thermal velocities.
x ) for a * = 1 as a function of the (reduced) velocity c x for α = 0.9 and five different values of the mass ratio m/m g : 1 (yellow line), 5 (blue line), 10 (red line), 50 (green line), and the Brownian limit (black lines).Here, T * g = 1, d = 3, and φ = 0.0052.The dashed line refer to the DSMC results for IHS in the Brownian limit.

Concluding remarks
In our study, we have explored the non-Newtonian transport properties of a dilute granular suspension subjected to USF using the Boltzmann kinetic equation.The particles are represented as d-dimensional hard spheres with mass m and diameter σ, immersed in an interstitial gas acting as a thermostat at temperature T g .Various models for granular suspensions incorporate a gas-solid force to represent the influence of the external fluid.While some models consider only isolated body resistance via a linear drag law [5,6,8,[10][11][12]47,48], others [13,49] include an additional Langevin-type stochastic term.In this paper, we consider a suspension model where the collisions between grains and particles of the interstitial (molecular) gas are taken into account.Thus, based on previous studies [50,51], we discretize the surrounding molecular gas, assigning individual particles with mass m g and diameter σ g , thereby accounting for elastic collisions between grains and background gas particles in the starting kinetic equation.
Under USF conditions, the system is characterized by constant density profiles n and n g , uniform temperatures T and T g , and a (common) flow velocity U x = U g,x = ay, where a denotes the shear rate.In agreement with previous investigations on uniform sheared suspensions, the mean flow velocity U is coupled to that of the gas phase U g .Consequently, the viscous heating term due to shear and the energy transferred by the grains from collisions with the molecular gas are compensated by the cooling terms derived from collisional dissipation, allowing the achievement of a steady state.A distinctive feature of the USF is that the one-particle velocity distribution function f (r, v) depends on space only through its dependence on the peculiar velocity V = v − U. Consequently, the velocity distribution function becomes uniform in the Lagrangian reference frame moving with V.This means that f (r, v) ≡ f (V).Based on symmetry considerations, the heat flux q vanishes, making the pressure tensor P the relevant flux.Therefore, to understand the intricate dynamics of granular suspensions under shear flow, it is imperative to focus on their non-Newtonian properties -derived from the pressure tensor.These include the (reduced) temperature χ = T/T g , the (reduced) nonlinear shear viscosity η * , and the (reduced) normal stress difference Ψ * .
Given that the most challenging aspect of dealing with the Boltzmann equation lies in the collision operator, it is reasonable to explore alternatives that render this operator more analytically tractable than in the case of IHS.Among the most sophisticated techniques in this regard is to consider the Boltzmann equation for IMM.As in the case of elastic collisions [28,41], the collision rate for IMM is independent of the relative velocity of the colliding particles.As a consequence, the collisional moments of degree k of the Boltzmann collisional operator can be expressed as a linear combination of velocity moments of degree less than or equal to k.To complement the results derived for IMM, we have also considered in this paper the use of a BGK-type kinetic model where the true Boltzmann operator is replaced by a simple relaxation term.Here, we employed both approaches to compute the rheological properties of the sheared granular suspension.Thus, our objective was twofold.Firstly, we aimed to assess the reliability and compatibility of the proposed models with previous results [25] obtained for IHS using Grad's moment method.Additionally, DSMC simulations for IHS were performed as an alternative method to validate any potential discrepancies identified.Secondly, taking advantage of the capabilities provided by the BGK model, we endeavored to calculate the velocity distribution function and the higher-order moments that offer insights into its characteristics.
Before proceeding with the computation of rheological properties, it is necessary to understand the response of the molecular fluid to shear stress.This assessment was also conducted using both IMM and BGK-type kinetic model that were later used to model the granular gas.A novelty here is the application of a force (Gaussian thermostat) of the form F = −mξV to compensate for the energy gained through viscous shear stresses.This force, by consistency, also applies to the granular gas, maintaining convergence to a steady state.As anticipated, the results agree well with those obtained through Grad's moment method [25] [see Eqs. ( 32) and ( 64)].Consequently, once the problem conditions (including the shear rate a) are defined, the molecular temperature T g is determined, effectively serving as a thermostat for the granular gas.
After accurately describing the rheology of the molecular gas, we focused on modeling the granular gas.Using both IMM and BGK-type model separately, we have calculated the non-zero elements of the the pressure tensor.The knowledge of these elements allows us to identify the relevant rheological properties of the granular suspension.As shown in Fig. 2, these quantities are represented as functions of the coefficient of restitution α and the mass ratio m/m g .In particular, we find that the theoretical results obtained from the Grad's method for IHS, IMM and BGK-type model show remarkable agreement, with almost indistinguishable curves.This underlines the effectiveness of structurally simple models in capturing the complexities of sheared granular suspensions.We observe a DST-type transition starting at a certain value of a * , which increases with the mass ratio m/m g .Interestingly, similar to the MD simulations performed for IHS [14], the DSMC data suggest a more abrupt transition than predicted by theory.Given that the main divergences between Grad's (for IHS) and DSMC results arise from the form of the distribution function, the significance of investigating higher-order moments to assess the deviation of the distribution function from its Maxwellian reference is then justified.
Based on previous literature where discrepancies in fourth-degree moments have been observed [16,39], and acknowledging the potential lengthiness of calculations, for the sake of simplicity, we decided to employ only the BGK-type model to compute higher-order moments.Specifically, we concentrated on the (symmetric) fourth-degree moments M 4|0 and M 2|xy .The shear-rate dependence of these moments is illustrated in Fig. 3 for the same parameter values of α and m/m g as those employed for the rheological quantities.Initially, we note that the fourth-degree moments also exhibit an abrupt transition at a value of a * that increases with m/m g until reaching the Brownian limit.Furthermore, DSMC simulations in the Brownian limit qualitatively capture the profile of these moments, although some quantitative disparities are apparent.To ascertain the extent of these discrepancies, we narrowed our focus to the interval 0 < a * < 1 where non-Newtonian effects are apparent.Additionally, we included IMM results directly as obtained from an effective Fokker-Planck-type model [17].Figure 4 illustrates that BGK results overestimate the deviation from the moments computed when no shear stress is applied compared to DSMC simulations and the results obtained using an effective force to model the interstitial gas.These disparities are also observed in the marginal distribution function φ x .
The theoretical findings presented here motivate the comparison with computer simulations.Although the observed agreement in the Brownian limit is encouraging, there is scope to extend this

Figure 2 .
Figure 2. Plots of the (steady) granular temperature χ (a), the non-Newtonian shear viscosity η * (b), and the normal stress difference Ψ * (c) as functions of the (reduced) shear rate a * for two different values of the coefficient of restitution α: 0.9 (left panel) and 1 (right panel).The graphs represent four distinct mass ratio values m/m g : 10 3 (yellow lines), 10 4 (blue lines), 10 6 (red lines), and the Brownian limit (black lines).Here, T * g = 1, d = 3, and φ = 0.0052.The solid lines correspond to the IMM results, the dashed lines are the BGK-like results, and the dotted lines refer to Grad's solution for IHS.Symbols denote computer simulation results performed in the Brownian limit: circles refer to the DSMC data obtained in this paper for IHS while squares are MD results obtained in Ref.[14] for IHS.

Figure 3 (
Figure 3(a)  illustrates the ratio M * 4|0 (a * )/M * 4|0 (0) as a function of a * for α = 0.9 and 1 and three values of the mass ratio.We observe that variations in the mass ratio do not significantly alter the trends observed in the Brownian limiting case[17].An abrupt transition in the higher-order moments is evident within a small region of a * .Specifically, the kurtosis M * 4|0 increases with the mass ratio m/m g until it converges to the value obtained in the Brownian limit.Consistent with the conclusions drawn in Ref.[25], an increase in the mass of the particles of the granular gas results in an elevation of the granular temperature.Consequently, non-equipartition accentuates and moves the suspension away from equilibrium, leading to an increase in kurtosis as the distribution function deviates from its Maxwellian form.Regarding the influence of collisional dissipation, we observe that the effect of α on M * 4|0 remains relatively discrete.Figure3(b) illustrates the shear-rate dependence of the (reduced) moment M * 2|xy .This moment vanishes in the absence of shear rate (a * = 0).Similar conclusions to those made for the moment M * 4|0 can be drawn.Theoretical predictions for the fourth-degree moments are compared against DSMC simulations for IHS conducted in this paper in the Brownian limiting case.A qualitative agreement is observed, although simulations suggest a sharper transition.Some quantitative discrepancies are noticeable, which are mainly disguised by the scale.To assess the reliability of the BGK-type results, we focus on the region 0 < a * < 1, where all the fourth-degree velocity moments of IMM are well-defined functions of the shear rate.In addition, non-Newtonian effects are still significant in the range of

Figure 5 .
Figure 5. Plot of the ratio R x(c x ) = ϕ x (c x )/(π −1/2 e −c 2x ) for a * = 1 as a function of the (reduced) velocity c x for α = 0.9 and five different values of the mass ratio m/m g : 1 (yellow line), 5 (blue line), 10 (red line), 50 (green line), and the Brownian limit (black lines).Here, T * g = 1, d = 3, and φ = 0.0052.The dashed line refer to the DSMC results for IHS in the Brownian limit.