# Heat flux of driven granular mixtures at low density: Stability analysis of the homogeneous steady state

## Abstract

The Navier–Stokes order hydrodynamic equations for a low-density driven granular mixture obtained previously [Khalil and Garzó, Phys. Rev. E 88, 052201 (2013)] from the Chapman–Enskog solution to the Boltzmann equation are considered further. The four transport coefficients associated with the heat flux are obtained in terms of the mass ratio, the size ratio, composition, coefficients of restitution, and the driven parameters of the model. Their quantitative variation on the control parameters of the system is demonstrated by considering the leading terms in a Sonine polynomial expansion to solve the exact integral equations. As an application of these results, the stability of the homogeneous steady state is studied. In contrast to the results obtained in undriven granular mixtures, the stability analysis of the linearized Navier–Stokes hydrodynamic equations shows that the transversal and longitudinal modes are (linearly) stable with respect to long enough wavelength excitations. This conclusion agrees with a previous analysis made for single granular gases.

## I Introduction

It is well established that when a granular material is externally excited the motion of grains resembles the random motion of atoms or molecules in an ordinary gas. These conditions are referred to as rapid flow conditions and has been an active area of research in the past several decades (1); (2); (3); (4); (5); (6). On the other hand, since the collisions between grains are inelastic, the energy monotonically decays in time so that external energy must be added to the system to keep it under rapid flow conditions. In real experiments, the external energy is injected into the granular gas from the boundaries (for instance, shearing the system or vibrating its walls (7); (8)), by bulk driving (as in air-fluidized beds (9); (10)) or by the presence of the interstitial fluid (11); (12); (13). When this external energy compensates for the energy dissipated by collisions, then a nonequilibrium steady state is reached.

Nevertheless, when the granular gas is locally driven, strong spatial gradients appear in the bulk domain and consequently the usual Navier–Stokes hydrodynamic equations are not suitable. Thus, in order to avoid the difficulties linked to the theoretical description of far from equilibrium states, it is usual in computer simulations (14); (15); (16) to drive the gas by means of an external force. Following the terminology employed in nonequilibrium molecular dynamics simulations of ordinary gases (17), this type of external forces are usually called *thermostats*.

Although several kind of thermostats have been proposed in the literature (18); (19), an interesting and reliable model was proposed by Puglisi and coworkers (14) to homogeneously fluidize a granular gas by an external force. More specifically, the thermostat is constituted by two different terms: (i) a drag force proportional to the velocity of the particle and (ii) a stochastic force (Langevin model) where the particles are *randomly* accelerated between successive collisions. This latter force has the form of a Gaussian white force with zero mean and finite variance (20). At a kinetic level and in the low-density regime, the corresponding kinetic equation for this thermostat has the structure of a Fokker-Planck equation (21); (22) plus the Boltzmann collision operator. Apart from using these external forces as thermostatic forces, it is worth noting that the above Langevin-like model has the same structure as several kinetic equations proposed in the literature (23); (24); (25) to model granular suspensions for low Reynolds numbers. In this context, the viscous drag force mimics the friction of solid particles on the viscous surrounding fluid while the stochastic force models the transfer of energy from the interstitial fluid to the granular particles.

More recently (26), the Langevin-like model has been extended to the case of granular mixtures. Since the model attempts to incorporate the influence of gas phase into the dynamics of grains, the drag force is defined in terms of the “peculiar” velocity instead of the instantaneous velocity of the solid particles. Here, is the mean velocity of the interstitial gas and is assumed to be a known quantity of the suspension model. Thus, the fluid-solid interaction force is constituted by an additional term (apart from the drag and stochastic forces) proportional to the difference between the mean velocities of gas and solid phases ( being the mean flow velocity of the granular particles).

Aside introducing the model for driven granular binary mixtures at low density, the Navier–Stokes hydrodynamic equations were derived in Ref. (26) by solving the corresponding Boltzmann kinetic equation by means of the Chapman–Enskog method (27). As in the free cooling case (28); (29); (30); (31), the transport coefficients are given in terms of the solution to a set of coupled linear integral equations. These integral equations can be approximately solved by considering the leading terms in a Sonine polynomial expansion. This task was in part carried out in Ref. (26) where a complete study of the dependence of the mass flux (four diffusion coefficients) and the shear viscosity coefficient on the parameters of the mixture (masses, sizes, concentration, and coefficients of restitution) was worked out. Therefore, a primary objective here is to demonstrate the variation of the transport coefficients of the heat flux by using the same Sonine polynomial approximation as was found applicable for undriven mixtures (29). As expected, the results clearly show that the influence of thermostats on heat transport is in general significant since the dependence of inelasticity on the four heat flux transport coefficients is different from the one found before for undriven mixtures (29).

Needless to say, the knowledge of the complete set of the Navier–Stokes transport coefficients of the mixture opens up the possibility of studying several problems. Among them, an interesting application is to analyze the stability of the homogeneous steady state (HSS). This state plays a similar role to the homogeneous cooling state (HCS) for freely cooling granular flows. It is well known (32); (33) that the HCS becomes unstable when the linear size of the system is larger than a certain critical length . The dependence of on the control parameters of the system can be obtained from a linear stability analysis of the Navier–Stokes hydrodynamic equations. The calculation of is interesting by itself and also as a way of assessing the reliability of kinetic theory calculations via a comparison against computer simulations. In fact, previous comparisons made for between theory and simulations for undriven granular gases (34); (35); (36); (37); (38) have shown a good agreement even for strong inelasticity and/or disparate masses or sizes in the case of multicomponent systems.

A natural question is whether the HSS may be unstable with respect to long enough wavelength excitations as occurs in the HCS. The results derived here by including the complete dependence of the transport coefficients on the coefficients of restitution show that the HSS is linearly stable with respect to long wavelength perturbations. This conclusion agrees with a previous stability analysis carried out for single driven granular gases (39). On the other hand, the quantitative forms for the dispersion relations obtained in this paper are quite different from the ones derived for monocomponent inelastic gases.

The plan of the paper is as follows. In Sec. II the model for driven granular mixtures is introduced and the corresponding hydrodynamic equations to Navier–Stokes order are recalled. Next, the four transport coefficients associated with the heat flux are obtained in Sec. III in terms of the parameters of the mixture and the driven parameters of the model. The results for these coefficients are also illustrated for a common coefficient of restitution and same size ratio as functions of at a concentration for several values of the mass ratio. The results clearly show a significant deviation of these coefficients from their values for ordinary (elastic) gases, specially for strong inelasticity as expected. Section IV addresses the linear stability analysis around the HSS and present perhaps the main relevant findings of the paper. The paper is closed in Sec. V with some concluding remarks.

## Ii Hydrodynamics from Boltzmann kinetic theory

We consider a granular gas modeled as a binary mixture of inelastic hard spheres in dimensions with masses and diameters (). The inelasticity of collisions among all pairs is characterized by three independent constant coefficients of normal restitution , , and , where . Here, is the coefficient of restitution for collisions between particles of species and . We also assume that the particles interact with an external bath. The influence of the bath on the dynamics of grains is encoded in two different terms: (i) a stochastic force assumed to have the form of a Gaussian white noise and (ii) a drag force proportional to the velocity of the particle. Under these conditions, in the low-density regime, the set of coupled nonlinear Boltzmann equations for the one-particle distribution function of each species reads (26)

(1) | |||||

where is the Boltzmann collision operator (2). In Eq. (1), is the drag (or friction) coefficient, is the strength of the correlation in the Gaussian white noise, and and are arbitrary constants of the driven model. In addition, and , where

(2) |

is the mean flow velocity of solid particles. In addition, as said before, can be interpreted as the mean velocity of the fluid surrounding the grains. In Eq.(2), is the total mass density and

(3) |

is the local number density of species . In the case of monodisperse granular gases and for and , the Boltzmann equation (1) is similar to the one proposed in Ref. (25) to model the effects of the interstitial fluid on grains in monodisperse gas-solid dense suspensions. The only difference between both descriptions is that in the latter case the parameters and are functions of the Reynolds number, the solid volume fraction, and the difference . In this context, the results derived here can be useful for instance to understand the stability of homogeneous bidisperse suspensions.

The parameters and can be seen as free parameters of the model. Thus, when , the Boltzmann equation (1) describes the time evolution of a granular mixture driven by the stochastic thermostat employed in several previous works (40); (41). On the other hand, the case and reduces to the Fokker–Planck model for ordinary (elastic) mixtures (14). In this context, our model can be understood as a generalization of previous driven models. Moreover, dimensional analysis clearly shows that the dependence of and on the masses and depends on the specific values of both and considered in each particular situation.

Apart from the fields and , the other relevant hydrodynamic field of the mixture is the granular temperature . It is defined as

(4) |

where is the total number density. At a kinetic level, it is also convenient to introduce the partial kinetic temperatures for each species defined as

(5) |

The partial temperatures measure the mean kinetic energy of each species. According to Eq. (4), the granular temperature of the mixture can also be written as

(6) |

where is the concentration or mole fraction of species .

Upon deriving the Boltzmann equation (1), it has been assumed that the typical collision frequency for collisions between the solid particles and the bath is much larger than the corresponding frequency for particle collisions (22). In addition, given that the Boltzmann collision operator is not affected by the surrounding fluid, one expects then that the suspension model defined by Eq. (1) will accurately describe situations where the stresses exerted by the interstitial fluid on particles are sufficiently small so they only have a weak influence on the dynamics of grains.

The macroscopic balance equations for the number density of each species , flow velocity , and temperature can be derived from the set of Boltzmann equations (1). They are given by (26)

(7) |

(8) |

(9) |

In the above equations, is the material derivative, is the partial mass density of species , is the mass flux of species relative to the local flow, is the pressure tensor, is the heat flux, and is the cooling rate. In the case of a binary mixture, there are independent hydrodynamic fields, , , , and . To obtain a closed set of hydrodynamic equations, one has to express the fluxes and the cooling rate in terms of the above hydrodynamic fields. These expressions are called “constitutive equations.” Such expressions were derived in Ref. (26) by solving the Boltzmann equation from the Chapman–Enskog method.

It is worth remarking that several physical assumptions are required to obtain the above constitutive equations. First, the hydrodynamics fields , , and are assumed to be the slowest magnitudes of the system and hence, for any initial condition, all the other magnitudes (for instance, the partial temperatures , the fluxes, the cooling rate, ) become functionals of the hydrodynamic fields for times longer than the mean free time. Second, the functional dependence of the distribution functions on the hydrodynamic fields is well approximated by the linear terms (Navier–Stokes approximation) of a series expansion in powers of the spatial gradients. The reference state in this expansion is in general supposed to be the local version of the time-dependent homogeneous state described in Refs. (22) and (26). Third, the difference of the mean velocities is assumed to be at least of first order in the spatial gradients.

Moreover, instead of providing the mass and heat fluxes in terms of the partial densities , it is more convenient to express such fluxes in terms of a different set of experimentally more accessible fields like the mole fraction and the pressure . The corresponding hydrodynamic balance equations for and can be easily derived from Eqs. (7) and (II) by a simple change of variables. They are given by

(10) |

(11) |

where is the temperature ratio for species .

The constitutive equations up to the Navier–Stokes order are

(12) | |||||

(13) |

(14) |

(15) |

The transport coefficients in these equations are

(16) |

As for elastic collisions (27), the Navier–Stokes transport coefficients are given in terms of the solution to a set of coupled linear integral equations. In the steady state, these integral equations can be approximately solved by considering the leading terms in a Sonine polynomial expansion. The evaluation of the diffusion coefficients and the shear viscosity coefficient was accomplished in Ref. (26). The transport coefficients associated with the heat flux will be explicitly determined in Sec. III in terms of the parameters of the mixture (masses, sizes, concentration, and coefficients of restitution) and the driven parameters and of the suspension model.

Once the complete set of transport coefficients is known, the Navier–Stokes constitutive equations (12)–(15) are substituted into the exact balance equations (8)–(II) to obtain the corresponding Navier–Stokes hydrodynamic equations for a driven binary granular mixture. They are given by

(17) |

(18) |

(19) |

(20) |

Here, as mentioned in several previous works (29); (42), the general form of the cooling rate should include second-order gradient contributions in Eqs. (II) and (II). However, as shown for a one-component dilute granular gas (43), these contributions are found to be very small with respect to the remaining contributions and hence they can be neglected in the evaluation of the cooling rate. Apart from this approximation, Eqs. (17)–(II) are exact to second order in the spatial gradients for a low-density driven granular binary mixture.

## Iii Heat flux transport coefficients

The evaluation of the heat flux transport coefficients , , , and requires to consider the second Sonine approximation. The expressions of these transport coefficients are

(21) | |||||

(22) | |||||

(23) | |||||

(24) | |||||

where . In Eqs. (21)–(24), the explicit forms of the dimensionless Sonine coefficients , , , and have been determined in the Appendix. Moreover, the (reduced) diffusion coefficients , , , and are defined by the relations

(25) |

where

(26) |

is the reduced mass and

(27) |

is an effective collision frequency. Here, is a thermal speed for a binary mixture and . The above diffusion coefficients were obtained in Ref. (26) in the first Sonine approximation.

Before considering a binary mixture, it is interesting to check the consistency between the present expression of with the one derived for a monocomponent granular gas in Ref. (39) when . Therefore, for mechanically equivalent particles (, , ), the Dufour coefficient vanishes as expected and the heat flux (13) can be written as

(28) |

where

(29) |

Note that the relation has been used in writing the heat flux in the form (28). In the limit of mechanically equivalent particles, the results derived in the Appendix yield

(30) |

(31) |

where and the remaining quantities are defined in the Appendix. Equations (30) and (31) agree with the expressions obtained in Ref. (39) when and in the case that one neglects non-Gaussian corrections to the zeroth-order distribution function (namely, by taking in the expressions provided in Ref. (39)). This confirms the relevant known limiting case for the granular mixture results described here.

Furthermore, in order to compare the present results with those obtained for *undriven* granular mixtures (28); (29), it is convenient to consider an equivalent representation in terms of the heat flow defined as (44)

(32) |

where use has been made of the requirement in the second equality of Eq. (32). The difference between and is the contribution to the heat flux coming from the diffusion flux . In particular, for ordinary mixtures (elastic collisions) and in the context of Onsager’s reciprocal relations (44), is the flux conjugate to the temperature gradient in the form of the entropy production. Moreover, the thermal conductivity coefficient in a mixture is measured in the absence of diffusion (i.e., when ). To identify this coefficient one has to express in terms of , , , and . In this representation, the corresponding coefficient of defines the thermal conductivity coefficient (44). To express in terms of , one may use Eq. (12) to write the gradient of mole fraction as

(33) | |||||

The heat flow is obtained by substituting Eq. (13) into Eq. (32) and eliminating by using the relation (33). The final form of is

(34) |

where

(35) |

(36) |

(37) |

(38) |

As for elastic collisions, the coefficient is the thermal conductivity while is called the thermal diffusion (Soret) factor. The coefficient is also present for ordinary mixtures (if both species are mechanically different) while there is a new contribution to proportional to that vanishes for elastic collisions. This latter contribution defines the transport coefficient .

### iii.1 Some illustrative systems

It is quite apparent that the dimensionless forms of the heat flux transport coefficients depend on many parameters: . Moreover, the driven parameters ( and ) along with the parameters and characterizing the class of model considered must be also specified. A complete exploration of the full parameters space is simple but beyond the goal of this paper. Thus, to illustrate the differences between ordinary and granular mixtures the transport coefficients are scaled with respect to their values in the elastic limit. In addition, for the sake of simplicity, we consider a common coefficient of restitution (), a common size , a mole fraction , and four different values of the mass ratio: and . These are the same systems as those studied before for undriven granular mixtures (29). We choose a three-dimensional system () with , , , and .

Figures 1–4 show the heat flux transport coefficients versus the (common) coefficient of restitution for the systems mentioned before. It is understood that all the coefficients have been evaluated with respect to their elastic values, except the case of since this coefficient vanishes for elastic collisions. In this latter case, we have considered the (reduced) coefficient

(39) |

Figure 1 shows the thermal diffusion factor . Note that when . It is seen that while exhibits a non-monotonic dependence on inelasticity when the defect component is lighter than the excess component , this transport coefficient decreases monotonically with decreasing when . Moreover, the impact of inelasticity on the functional form of thermal diffusion is more important for than in the opposite case. The thermal conductivity is shown in Fig. 2. We observe that the (scaled) thermal conductivity decreases with increasing inelasticity, except for mechanically equivalent particles (). This behavior contrasts clearly with the results found for undriven mixtures (29) where the ratio increases with increasing dissipation, regardless of the value of the mass ratio. It is also apparent that, at a given value of , there is a significant mass dependence of the thermal conductivity, especially at strong inelasticity. The coefficient is illustrated in Fig. 3 where it is found that its magnitude is larger than the one obtained before for both the thermal diffusion and thermal conductivity coefficients. As for undriven mixtures (29), is negative for the cases studied here. Finally, the (scaled) transport coefficient is plotted in Fig. 4. This coefficient vanishes for mechanically equivalent particles. In contrast to the thermal conductivity coefficient, increases with inelasticity when while the opposite happens when . Furthermore, the impact of collisional dissipation on this coefficient is quite significant for small mass ratios.

In summary, the heat flux transport coefficients for driven granular mixtures differ noticeably from those for ordinary mixtures, especially at strong inelasticity. Depending on the value of the mass ratio, in some cases these transport coefficients decrease with decreasing while in others they increase with inelasticity. This non-monotonic behavior with the mass ratio is also present in the free cooling case (29). Since the expressions of the heat flux transport coefficients obtained here are quite complex (in contrast to more phenomenological approaches), it is very difficult to provide a simple explanation of the non-monotonic trend in the mass ratio observed in Figs. 1–4. Finally, with respect to the impact of masses on the heat transport, we also observe that in general the transport coefficients , , , and exhibit a strong influence of the mass ratio, being this influence larger than the one found in undriven granular mixtures (29).

## Iv Linear stability analysis of the homogeneous steady state

It is well known for undriven granular mixtures that the HCS is unstable against long enough wavelength perturbations (29); (45). These instabilities can be well predicted by a linear stability analysis of the Navier–Stokes hydrodynamic equations. In fact, the solution of the linearized hydrodynamic equations provides a critical length beyond which the system becomes unstable. The theoretical predictions of (42); (29); (45) have been shown to compare quite well with computer simulations for monocomponent (34); (35); (36) and multicomponent (37); (38) granular fluids.

On the other hand, a similar analysis for single driven granular fluids (39) has clearly shown that the HSS is always (linearly) stable. This analytical finding agrees with the results obtained from Langevin dynamics simulations (15). An interesting question arises then as to whether, and if so to what extent, the conclusions drawn for monocomponent driven granular gases (15); (39) may be altered when a binary mixture is considered.

In order to analyze the stability of the homogeneous solution, Eqs. (17)–(II) must be linearized around the HSS. In this state, the hydrodynamic fields take the steady values , , , and , where the subscript denotes the steady state. Moreover, in reduced units, the steady state condition determining the temperature ratios is

(40) |

where , is the (reduced) partial cooling rate, and

(41) |

Here, is defined by Eq. (27) with the replacements and . It is interesting to note that for elastic collisions and mechanically different particles energy equipartition is fulfilled () only when (26). This is the expected result in agreement with the fluctuation-dissipation relation (46). The above relation between and will be kept in the remaining part of this section.

We assume that the deviations are small where denotes the deviations of , , , and from their values in the HSS. Moreover, as usual we also suppose that the interstitial fluid is not perturbed and hence, .

In order to compare the present results with those obtained for undriven granular mixtures (29), the following dimensionless space and time variables are introduced:

(42) |

where . Then, as usual, a set of Fourier transform dimensionless variables are defined as

(43) |

(44) |

where is defined as

(45) |

Note that in Eq. (45) the wave vector is dimensionless.

In terms of the above dimensionless variables, as usual, the transverse velocity components (orthogonal to the wave vector ) decouple from the other four modes and hence can be obtained more easily. Their evolution equation is

(46) |

where the eigenvalue is given by

(47) |

where

(48) |

and . Note that although the subscript has been omitted in Eqs. (46)–(48) for the sake of brevity, it is understood that all the quantities are evaluated in the HSS.

The solution to Eq. (46) is simply given by

(49) |

For mechanically equivalent particles, and so, . This means that the transversal shear mode is always stable for monocomponent granular gases. This finding is consistent with the results derived for simple granular fluids (39). In the general case, the (reduced) transport coefficient is (26)

(50) |

where

(51) |

and

(52) | |||||

Here, is defined by Eq. (87). Since , according to Eqs. (50) and (51) it is easy to prove that

(53) |

and so, one infers the general result

(54) |

Therefore, the transversal shear mode is always (linearly) stable.

The remaining (longitudinal) four modes are the concentration field , the temperature field , the pressure field , and the longitudinal component of the velocity field (parallel to ). As in the undriven case (29), they are coupled and obey the time-dependent equation

(55) |

where here denotes the set of four variables . The square matrices in Eq. (55) are