Hydrodynamics of a Granular Gas in a Heterogeneous Environment

Abstract

We analyze the transport properties of a low density ensemble of identical macroscopic particles immersed in an active fluid. The particles are modeled as inelastic hard spheres (granular gas). The non-homogeneous active fluid is modeled by means of a non-uniform stochastic thermostat. The theoretical results are validated with a numerical solution of the corresponding the kinetic equation (direct simulation Monte Carlo method). We show a steady flow in the system that is accurately described by Navier-Stokes (NS) hydrodynamics, even for high inelasticity. Surprisingly, we find that the deviations from NS hydrodynamics for this flow are stronger as the inelasticity decreases. The active fluid action is modeled here with a non-uniform fluctuating volume force. This is a relevant result given that hydrodynamics of particles in complex environments, such as biological crowded environments, is still a question under intense debate.

1. Introduction

Biological Systems, and more specifically, systems with active particles, are in general out of equilibrium, which implies that for these systems the results of equilibrium statistical mechanics do not in general apply. The term active matter refers to a system with an ensemble of particles that are able to transform energy—internal or environmental—into movement [1]. A correct description of transport phenomena in this type of systems can be achieved via a specific kinetic theory, that takes into account the peculiar energy processing of these particles, see for example [2,3,4,5] for the well known Vicsek-like models [6].

In recent years, the study of active particles in crowded environments has become the subject of intense study [7]. One of the motivations is the large number of promising applications derived from the use of active particles as micro-machines [7]. In addition, the behaviour of inert particles inside complex active environments remains an open question. For instance, there are experimental studies on inert Brownian particles trapped by a harmonic potential in a fluid composed by bacteria [8]. Also, Sándor et al. [9] study active Brownian particles on a travelling-wave substrate. Other works study particles which can be trapped by holes [10,11]. More realistic active systems such as biological tissues, share important analogies with granular matter, namely glass transition, cage effect, fluidization, solidification and the appearance of giant density fluctuations [12,13,14,15]. There, active complex non-homogeneous forces appear. Also the non-Gaussian behaviour of diffusing non-active particles in heterogeneous media has been studied recently for a spatially varying diffusion coefficient [16].

On the other hand, the transport properties of granular matter are of interest at a fundamental level and also for applications, since granular matter is present in a variety of industry and technology sectors [17,18,19,20]. Granular transport theories usually make a connection with traditional continuum mechanics theories, due to fundamental similarities observed with classical fluids [17,19]. Granular particles loose a fraction of their kinetic energy after collision and for this reason we need to continuously excite the system, in order to sustain the dynamics. This excitation usually comes from the system boundaries, which results in inhomogeneous flows that may eventually become steady [17,19,21].

We propose in this work a study that tries to model a system of macroscopic inert particles that are immersed in an active fluid with inhomogeneous temperature. For this, we model the system as a granular fluid [22] subject to a nonuniform stochastic volume force (in the form of a white noise [23]). The inelastic cooling term, inherent in a granular fluid [19], and a stochastic force allow us in this work to model the energy sink and source terms that are characteristic to active particle systems [1,7]. More specifically, we focus on a very low density granular system (i.e., a granular gas) surrounded by a low density (inhomogeneous) fluid that can be considered as a coarse grained version of an active matter system. Since the system has very low density, particle velocity correlations can be neglected and thus the the inelastic Boltzmann equation applies [19,24]). Moreover, systems of solid particles suspended in a low density interstitial fluid usually fall into the grain-grain collision dominated regime [25], as opposed to the interstitial fluid viscosity-dominated regime, that applies for most cases of interstitial liquids [26].

In previous theoretical works, granular gas fluidization has been studied only in the simplistic case of a homogeneous interstitial fluid [27,28,29,30,31,32]. And yet this theoretical homogeneous state is not quite a realistic situation [33]. Neat examples of inhomogeneous granular gas suspensions are the large sand plumes that can be observed in the Earth’s atmosphere for weeks [34]. A homogeneous energy source is obviously not a realistic situation either in active fluids, at a biological level [1]. For this reason we focus in this work on the more elaborate case of a non-uniform interstitial fluid.

However, the intent of this work is not experimental validation but to extract a more clear picture of the relevant physics of this kind of systems. As a result, we show there are granular gas flows immersed in an inhomogeneous active fluid that can be accurately described with Navier-Stokes (NS) hydrodynamics, even at high inelasticities. Furthermore, we surprisingly found cases where NS hydrodynamics is a better approximation for the granular gas than for a perfectly elastic gas subjected to the same boundary forcing conditions. In particular, we will show that, once the granular gas is fluidized, the intensity of the spatial gradients in response to external excitations (those coming from the boundaries) can actually be weaker for granular gases than for elastic gases.

In a previous experimental work [35], local balance between the total energy input (fluctuating force plus boundary heating) and the energy sink was found. We use now this kind of balance condition. As we will show below, this local balance results, in the specific set-up considered in the present work, in a steady flow with uniform heat flux throughout the system. This balance condition is in fact analogous to the one occurring in well known non-Newtonian granular flows like the uniform shear flow [36,37].

The structure of the paper is as follows: in Section 2 we present a description of the system and the corresponding steady state equations. In Section 3 we describe the flows where there is a balance between inelastic cooling coming from particle collisions and volume energy input (and with no shear). Section 4 is devoted to a brief discussion on the steady flows resulting from adding a weak shear Section 5 presents a discussion on the transport coefficients and rheology properties of the studied flows and the paper conclusions are drawn in Section 6. The simulation method and transport coefficient equations are discussed in the Appendices Appendix A and Appendix B respectively.

2. System and Steady Base State Equations

2.1. Description of the System

Our system consists of a set of identical inert smooth hard spheres (or disks) immersed in an active fluid. The system is limited by two infinite parallel hard walls, located at planes $y=\pm L/2$ respectively. The walls act like two distinct kinetic energy sources, characterized with temperatures $T_{\pm }$. They may have relative movement (wall velocities $U_{\pm }$ respectively), eventually inducing the particles to continuously flow along the channel between them.

Collisions between the inert particles do not preserve energy (i.e., particle collisions are inelastic). The particles have a diameter (radius) that we denote as $\sigma$ and their mass is m. For inelastic smooth hard spheres/disks, inelasticity may be characterized accurately, in a range of experimental situations, by a constant coefficient of normal restitution $\alpha$ [38,39]. This coefficient ranges from 1 for purely elastic collisions to 0 for purely inelastic collisions. In addition, we model the interstitial active fluid injection of thermal energy onto the granular gas particles as a stochastic force ${\mathbf{F}}^{\mathrm{st}}$. The equation of motion for a particle i can be written

$$m{\dot{\mathbf{v}}}_i={\mathbf{F}}_i^{\mathrm{st}}(\mathbf{r},t)+{\mathbf{F}}_i^{\mathrm{coll}},$$

(1)

where ${\mathbf{F}}_i^{\mathrm{coll}}$ is the force due to inelastic collisions, ${\mathbf{F}}_i^{\mathrm{st}}(\mathbf{r},t)$ is the force exerted by the heterogeneous medium and m the mass of particles (set to 1 for simplicity).

We can model this interaction by means of a fluctuating volume force [40], that we denote as ${\mathbf{F}}^{\mathrm{st}}\left(\mathbf{r}\right)$ and fulfills the conditions [32,35,40,41]

$$\begin{array}{c}\langle {\mathbf{F}}_i^{\mathrm{st}}(\mathbf{r},t)\rangle =\mathbf{0},\hfill \\ \langle {\mathbf{F}}_i^{\mathrm{st}}(\mathbf{r},t){\mathbf{F}}_j^{\mathrm{st}}(\mathbf{r},t^{\prime })\rangle =\mathbf{1}{m}^2\xi {\left(\mathbf{r}\right)}^2{\delta }_{ij}\delta (t-t^{\prime }),\hfill \end{array}$$

(2)

being ${\xi }^2\left(\mathbf{r}\right)$ the fluctuating force intensity [32,41] (that has dimensions of velocity squared times time), $\mathbf{1}$ the unit matrix in d dimensional space, ${\delta }_{ij}$ is the Kronecker delta function and $\langle A\rangle$ indicates average of magnitude A over a sufficiently long time interval (long compared with the characteristic microscopic time scale). It has been shown in the case of the large mass limit and homogeous energy injection, that Equation (1) for a particle can be written as a Langevin equation corresponding to a granular brownian particle [42].

On the other hand, collisions between two grains follow the inelastic smooth hard sphere collisional model

$$\begin{array}{c}{\mathbf{v}}_1^{\prime }={\mathbf{v}}_1-\frac{1+\alpha }{2}(\sigma ·\mathbf{g})\sigma \hfill \end{array}$$

(3)

$$\begin{array}{c}{\mathbf{v}}_2^{\prime }={\mathbf{v}}_2-\frac{1+\alpha }{2}(\sigma ·\mathbf{g})\sigma \hfill \end{array}$$

(4)

with $\sigma =({\mathbf{r}}_1-{\mathbf{r}}_2)/|{\mathbf{r}}_1-{\mathbf{r}}_2|$ and $\mathbf{g}={\mathbf{v}}_1-{\mathbf{v}}_2$. The primes denote pre-collisional velocities.

Notice that in (2) we have $\xi =\xi \left(\mathbf{r}\right)$ (a space-dependent noise intensity). The noisy term appears in the inelastic Boltzmann equation as a Fokker-Planck-like operator [43,44]

$$\left(\frac{\partial }{\partial t}+\mathbf{v}·\nabla -\frac{1}{2}\xi {\left(\mathbf{r}\right)}^2\frac{{\partial }^2}{\partial v^2}\right)f(\mathbf{r},\mathbf{v};t)=J\left[\mathbf{v}\right|f,f],$$

(5)

where J is the collisional integral for inelastic hard spheres and whose expression may be found elsewhere [19,45], $\mathbf{v}$ is particle velocity, and v its modulus.

Taking into account the system geometries, we consider only $\xi \left(\mathbf{r}\right)=\xi \left(y\right)$, that is the simplest space-dependence our geometry configuration allows for an inhomogeneous active fluid. The steady base states may be deduced from the reduction of the kinetic Equation (5)

$$v_y\frac{\partial f}{\partial y}-\frac{1}{2}\xi {\left(y\right)}^2\frac{{\partial }^{2}f}{\partial v^2}=J[f,f].$$

(6)

We have numerically solved the kinetic Equation (6) by means of the Direct Simulation Monte Carlo method [18,46], applying in each case the corresponding boundary conditions and heterogeneous fluctuating force properties. For more details on the simulation algorithm, see the corresponding section in the Appendix A.

2.2. Steady Base State Equations

Multiplying by velocity momenta the kinetic Equation (6) and performing velocity integrals, we may easily obtain from (6) the mass, momentum and energy balance equations [43,47]

$$\begin{array}{c}\frac{\partial P_{yy}}{\partial y}=0,0=-\frac{1}{mn}\frac{\partial P_{xy}}{\partial y}\hfill \end{array}$$

(7)

$$\begin{array}{c}T\zeta \left(\alpha \right)-m{\xi }^2\left(y\right)=-\frac{2}{dn}\left(P_{xy}\frac{\partial u_x}{\partial y}+\frac{\partial q_y}{\partial y}\right),\hfill \end{array}$$

(8)

where $\mathbf{u}$ is the flow velocity field (with $u_y$ its y-component), n is the particle density, and P and q are momentum flux tensor (also called stress tensor) and heat flux vector respectively, and $\zeta \left(\alpha \right)$ is the cooling rate, a magnitude that emerges from the energy balance due to inelasticity in the collisions [18,19]. Equations (7) and (8) contain no approximations; i.e., they are exact for this type of geometry, even if the steady state were not hydrodynamic. Notice also that the stress tensor element $P_{yy}$ is homogeneous and is not a function of time. This condition breaks if the flow is not laminar, but we will only consider laminar flows in this work.

We take now into account the characteristics of NS hydrodynamics in (7) and (8). In (7) we use that all diagonal elements of the stress tensor are equal and we also use that the horizontal heat flux $q_x$ is null (they both raise from fluxes terms that are of second order in the gradients [48]) and thus, from the first equation in (7), we obtain that $P_{yy}=p=\mathrm{constant}$. We also use the fluxes expressions at NS order [45]: $P_{xy}=\eta \left(\alpha \right)\partial u_x/\partial y$, $q_y=-\lambda \partial T/\partial y-\mu \partial n/\partial y$, where $\eta$ is the viscosity and $\lambda ,\mu$ are the heat flux transport coefficients. With this, the steady state forms of (7) and (8) are

$$\begin{array}{c}p=\mathrm{constant},P_{xy}=-\eta \left(\alpha \right)\frac{\partial u_x}{\partial y}=\mathrm{constant}\hfill \end{array}$$

(9)

$$\begin{array}{c}T\zeta \left(\alpha \right)-m{\xi }^2\left(y\right)=\hfill \\ \frac{2}{dn}\left(\eta \left(\alpha \right){\left(\frac{\partial u_x}{\partial y}\right)}^2+\kappa \left(\alpha \right)\frac{\partial T}{\partial y}+\mu \left(\alpha \right)\frac{\partial n}{\partial y}\right).\hfill \end{array}$$

(10)

Next, we take into account that the NS transport coefficients and cooling rate scale in the following forms with temperature [45]: $\eta \left(\alpha \right)={\eta }^{\ast }\left(\alpha \right){\eta }_0^{\prime }{T}^{1/2}$, $\kappa \left(\alpha \right)={\kappa }^{\ast }\left(\alpha \right){\kappa }_0^{\prime }{T}^{1/2}$, $\mu ={\mu }^{\ast }\left(\alpha \right){\kappa }_0^{\prime }{T}^{3/2}/n$, $\zeta ={\zeta }^{\ast }\left(\alpha \right)p/\left({\eta }_0^{\prime }{T}^{1/2}\right)$. Here, the coefficients ${\eta }_0^{\prime }$, ${\kappa }_0^{\prime }$ are related to the elastic gas [49] NS viscosity ${\eta }_0={\eta }_0^{\prime }{T}^{1/2}$ and thermal conductivity ${\kappa }_0={\kappa }_0^{\prime }{T}^{1/2}$, and thus

$$\begin{array}{c}{\eta }_0^{\prime }=\frac{d+2}{8}\mathrm{\Gamma }(d/2){\pi }^{-(d-1)/2}{m}^{1/2}{\sigma }^{-(d-1)}/8,\hfill \end{array}$$

(11)

$$\begin{array}{c}{\kappa }_0^{\prime }=\frac{d{(d+2)}^2}{16(d-1)}\mathrm{\Gamma }(d/2){\pi }^{-(d-1)/2}{m}^{-1/2}{\sigma }^{-(d-1)}.\hfill \end{array}$$

(12)

where $\mathrm{\Gamma }$ is the gamma function [50] and $d=2$ for disks and $d=3$ for spheres. Thus, the steady state balance Equations (9) and (10) yield

$$\begin{array}{c}-\eta {\left(\alpha \right)}^{\ast }{\eta }_0^{\prime }\sqrt{\frac{T}{T_r}}\frac{\partial u_x}{\partial y}=\mathrm{constant},\hfill \end{array}$$

(13)

$$\begin{array}{c}\sqrt{\frac{T}{T_r}}\frac{\partial }{\partial y}\left(\sqrt{\frac{T}{T_r}}\frac{\partial T}{\partial y}\right)=\hfill \\ \frac{d}{2}\frac{{\zeta }_0^{\ast }}{{\beta }_0^{\ast }}\frac{p^2}{T_r}-\frac{d}{2}\frac{m{\xi }^2\left(y\right)n{T}^{1/2}}{{\beta }_0^{\ast }{T}_r}-\frac{{\eta }_0^{\ast }}{{\beta }_0^{\ast }}\frac{T}{T_r}{\left(\frac{\partial u_x}{\partial y}\right)}^2,\hfill \end{array}$$

(14)

where ${\beta }^{\ast }\left(\alpha \right)\equiv {\kappa }^{\ast }\left(\alpha \right)-{\mu }^{\ast }\left(\alpha \right)$, and we call it effective thermal conductivity and $T_r$ is the reference unit of temperature, that we choose as $T_r=T_{-}$ (temperature at the lower, and colder, wall). Here, ${\zeta }^{\ast }\left(\alpha \right)$ is the dimensionless cooling rate [18,51] and the expressions of the reduced transport coefficients viscosity ${\eta }^{\ast }\left(\alpha \right)$, thermal conductivity ${\kappa }^{\ast }\left(\alpha \right)$, and ${\mu }^{\ast }\left(\alpha \right)$ that we used are the ones obtained by Garzó and Montanero [52] for a granular gas heated by a white noise. Their expressions are displayed in the transport coefficients section in the Appendix B.

We use a scaled space variable l that allows us to find an analytical solution for the flow velocity (13) and the temperature (14) profiles. This variable is defined by the relation $\sqrt{T/T_r}\partial /\partial y=\partial /\partial l$. When expressed in terms of this scaled variable, Equations (9) and (10) read

$$\begin{array}{c}\frac{\partial {\widehat{u}}_y}{\partial \widehat{l}}=a,\hfill \end{array}$$

(15)

$$\begin{array}{c}\frac{{\partial }^2\widehat{T}}{\partial {\widehat{l}}^2}=-\gamma (\alpha ,l),\hfill \end{array}$$

(16)

where a and $-\gamma$ are dimensionless magnitudes in the system and we call them local shear rate and temperature curvature, respectively. In (15), (16) and henceforth we use dimensionless variables indicated with a hat, where we have chosen as set of units: $T_r$ for temperature (already defined), m for mass, p for pressure (the steady state hydrostatic pressure, that as we said is a global constant for the system). As a unit for length we use

$${\lambda }_r=\sqrt{2}(2+d)\mathrm{\Gamma }(d/2){\pi }^{-(d-1)/2}{\left(\sqrt{2}{n}_r{\sigma }^{(d-1)}\right)}^{-1},$$

(17)

that is the reference mean free path [49] and for time we use

$${\nu }_r^{-1}={\lambda }_r\left(\sqrt{m/T_r}\right)$$

(18)

(inverse of reference collision frequency). These two definitions imply that the unit for velocity is $\sqrt{T_r/m}$. We indicate the spatial coordinate y in dimensionless form as $\widehat{y}$, and analogously for all dimensionless magnitudes.

The reader may infer the physical meaning of a and $\gamma$ with the aid of (7) and (8): a (reduced shear rate) measures how rapidly varies the flow velocity field, in units of the mean free path; ${|-\gamma |}^{1/2}$ measures the degree of imbalance between energy volume sinks and sources per unit time (as we explain above in Equation (18), this unit being the inverse of the collision frequency). The differential Equation (16) has been observed for the granular temperature in similar experimental configurations [33].

With the presence of the stochastic volume force in the energy balance Equation (8), it is straightforward to deduce from (14)–(16) and the definitions above that $-\gamma (\alpha ,l)$ has the form

$$\begin{array}{cc}\hfill -\gamma \left(\alpha \right)\equiv & \frac{2(d-1)}{d(d+2){\beta }^{\ast }\left(\alpha \right)}\times \hfill \\ & \left(\frac{d}{2}({\zeta }^{\ast }\left(\alpha \right)-{\widehat{\xi }}^2\left(\widehat{l}\right)\widehat{T}{\left(\widehat{l}\right)}^{-1/2})-{\eta }^{\ast }\left(\alpha \right)a^2\right),\hfill \end{array}$$

(19)

that, as it can be noticed, depends in general on position through $\widehat{l}$.

3. Steady Base States with Energy Balance and no Shear

We now want to look up for steady states in a granular gas fluidized by an heterogeneous source that may be accurately described by hydrodynamics at NS order. We will denote them as ’reference states’. As in rarefied gases, a close to unity Knudsen number ($\mathrm{Kn}$) signals hydrodynamics failure. Moreover, hydrodynamics at first order in the gradients (NS equations) are usually only valid if $\mathrm{Kn}$ takes low values [19,46]. The Knudsen number $\mathrm{Kn}$ is defined as a representative microscopic time or length scale over a macroscopic time or length scale. However, the choice of reference scales is not unique and picking an appropriate Knudsen number for each specific flow problem is a subtle question [53].

It was shown in a previous work [47] that Couette-Fourier flows (i.e., those occurring in our system, under the appropriate circumstances) are characterized by three representative and independent microscopic vs. macroscopic relative scales (or Knudsen numbers) that happen to be constant throughout the system: a, ${\left|\gamma \right|}^{1/2}$ (look at Equations (15) and (16)) and $\Delta T/L$, with $\Delta T\equiv T_{+}-T_{-}$, where $\Delta T/L$ measures the size of the temperature gradient imposed by the boundaries, this gradient being referred to mean free path units. These numbers arise from the specific form of the relevant differential equations for steady states, and the boundary geometry. Since the geometry of the differential Equations (15) and (16) and boundary conditions are the same in this work, we may safely use the same set of reference Knudsen numbers (see previous work [47] for more details on this issue). Furthermore, another previous analysis [51] for granular gases without an interstitial fluid showed that $\Delta T/L$ has a much weaker effect on the emergence non-Newtonian behavior.

Thus, let us analyze the class of solutions for steady flows where this $\Delta T/L$ is the only active source of spatial gradients; i.e., $a=0,\gamma =0$ and $\Delta T/L\lesssim 1$

$$\begin{array}{c}\frac{\partial {\widehat{u}}_y}{\partial \widehat{l}}=0,\hfill \end{array}$$

(20)

$$\begin{array}{c}\frac{{\partial }^2\widehat{T}}{\partial {\widehat{l}}^2}=0.\hfill \end{array}$$

(21)

Moreover, in the $\alpha =1$ limit this state corresponds to the classic Fourier flow of a molecular gas. In order to obtain $\gamma =0$ (with $a=0$), and taking into account (19) we find

$$\widehat{\xi }{\left(\widehat{l}\right)}^2\widehat{T}{\left(\widehat{l}\right)}^{-1/2}={\zeta }^{\ast }\left(\alpha \right)\Rightarrow \widehat{\xi }{\left(\widehat{l}\right)}^2={\zeta }^{\ast }\left(\alpha \right)\widehat{T}{\left(\widehat{l}\right)}^{1/2}.$$

(22)

Equation (20) implies that flow velocity $\widehat{u}$ is constant respect to $\widehat{l}$, which, from the definition of l leads to $\widehat{u}\left(\widehat{y}\right)$ also constant (and to all physical effects, flow velocity may be considered as null). Moreover, Equation (21) implies that temperature $\widehat{T}$ is linear in $\widehat{l}$, which (since $\sqrt{T/T_r}\partial /\partial y=\partial /\partial l$) yields $\widehat{y}\propto {\widehat{l}}^{3/2}$ and thus $\widehat{T}\left(\widehat{y}\right)\propto {\widehat{y}}^{2/3}$; i.e., the stochastic force generates a Fourier-like temperature profile [37]. At a physical level, this makes sense if we consider that the granular gas is fluidized by a low viscosity interstitial fluid that is heated from the boundaries. It is easy to deduce from (8), (10) and (19) that $\gamma =0$ implies that inelastic cooling is balanced by volume energy input and thus, that heat flux is uniform throughout the system.

Therefore our base steady states are (not previously reported) granular flows with uniform heat flux. Summarizing, these NS steady states have the following properties: (i) $P_{xy}=\eta \left(\alpha \right)\partial u_x/\partial y$, $q_y=-\lambda \partial T/\partial y-\mu \partial n/\partial y$, (ii) constant hydrostatic pressure $\widehat{p}$, (iii) null normal stress differences (${\widehat{P}}_{ii}={\widehat{P}}_{jj}=p$) and null ${\widehat{q}}_x$ [48], (iv) constant ${\widehat{q}}_y$, and (v) linear $\widehat{T}(\widehat{l})$ profiles, or equivalently, $\widehat{T}\left(\widehat{y}\right)\propto {\widehat{y}}^{3/2}$ [37]. Properties (i) to (iii) are fulfilled by all Couette-Fourier steady flows at NS order [47] and conditions (iv) and (v) are fulfilled by all Couette-Fourier steady flows with uniform heat flux [37].

We need to remark that the relation in (22) between fluctuating force intensity and temperature is not a hypothesis but a consequence of the experimental fact that there is local energy balance. On the other hand, there is an essential difference of this new flow respect to previously studied granular flows with uniform heat flux (the simple shear flow included [36,54]).

Very surprisingly, we can show that contrary to other steady granular flows, the kind of state defined by (22) can be accurately described in the context of Navier-Stokes equations. In order to prove this feature, we have simulated the system described by Equations (20)–(22) by means of the DSMC method (more details on simulations are given in the Appendix A). We have checked that the temperature profiles obtained under steady state in effect follow the behaviour $\widehat{T}\left(\widehat{y}\right)\propto y^{2/3}$, inherent to gas steady flows with uniform heat flux [37]. Results can be seen in Figure 1, for which the agreement with the NS theory prediction is excellent, independently of the inelasticity value.

4. Weakly Sheared Steady States

We discuss here weakly sheared states, in order to capture the gradual the eventual appearance of non-Newtonian behavior. Thus, starting out of the reference base states in the previous subsection; i.e., still with condition (22) being fulfilled, we analyze now sheared states ($a\ne 0$), for which we have from (19),

$$-\gamma =-\frac{(d-1)}{(d+2){\beta }^{\ast }\left(\alpha \right)}{\eta }^{\ast }\left(\alpha \right)a^2\ne 0.$$

(23)

For this type of flow we have also performed DSMC simulations. Relation (23) helps us understand that an increase in shear rate a (one of the relevant Knudsen numbers) implies an increase in $\gamma$ (one of the other Knudsen numbers). Thus, for increasing shear we should gradually depart from NS hydrodynamics as both a and ${\left|\gamma \right|}^{1/2}$ increase. However, relation (23) is valid only as long as NS hydrodynamics applies. For non-Newtonian regime, the relation between the temperature curvature $\gamma$, local shear rate a and coefficient of restitution $\alpha$ would be more intricate and in general it must be solved numerically [51]. A non-linear theory is beyond the scope of the present work, but, fortunately, we may use results from DSMC to analyze deviations from NS hydrodynamics.

The results we obtained are quite surprising (see

Table 1. DSMC measurements of Knudsen numbers a and ${\left|\gamma \right|}^{-1/2}$, for steady states described by (20)–(22), with $T_{+}/T_{-}=5$, $\Delta L=15\overline{\lambda }$ (i.e., $\Delta T/\Delta L=(4/15)((T_{+}-T_{-})/\overline{\lambda })$, and $\Delta U$ defined as $\Delta U\equiv (U_{+}-U_{-})/\sqrt{T_{-}/m}$.

	$\Delta \mathit{U}=5$	$\Delta \mathit{U}=5$	$\Delta \mathit{U}=3$	$\Delta \mathit{U}=3$
$\mathbf{\alpha }$	$\mathit{a}$	${\left|\mathit{\gamma }\right|}^{-\mathbf{1}/\mathbf{2}}$	$\mathit{a}$	${\left|\mathit{\gamma }\right|}^{-\mathbf{1}/\mathbf{2}}$
0.99	0.173	0.058	0.115	0.036
0.9	0.191	0.034	0.120	0.018
0.8	0.193	0.027	0.119	0.017
0.7	0.193	0.024	0.119	≲$1\times {10}^{-2}$
0.6	0.188	≲$1\times {10}^{-2}$	0.115	≲$1\times {10}^{-2}$
0.5	0.187	≲$1\times {10}^{-2}$	0.113	≲$1\times {10}^{-2}$

). In effect, we can see that similar values of local shear rate a tend to yield greater values of the other Knudsen number $\left|\gamma \right|$ for more elastic gases than for more inelastic gases. This means that, contrary to what would be expected, the Knudsen numbers in this type of flow tend globally to be smaller for higher inelasticities, which could indicate that departures from NS hydrodynamics are more important, for similar shear rates, for more elastic gases. However, we will confirm this point in the following Section 5, by measuring from DSMC simulations the set of transport coefficients and rheological properties and comparing with NS theory predictions.

We also performed an additional series with varying $\Delta T/\Delta L$ (not shown in figures) and checked that the measured transport coefficients do not depend on its value, analogously to the result in a previous work on granular flows with uniform heat flux [37].

5. Transport Coefficients and Rheology

Results in the preceding section suggest that: (a) our steady states are well described by NS hydrodynamics, even for strong inelasticities, as long as shearing from the boundaries is not too strong; (b) deviations from NS hydrodynamics are stronger for more elastic gases. Both observations would be contrary to what occurs in steady granular flows (like the simple shear flow [36], just to put one example).

Let us analyze in more detail to what extent the reference (no-shear) steady states described in Section 3 and the sheared states described in Section 4 are strictly NS flows or not, by analyzing the relevant transport properties. This analysis will be done as a function of the degree of inelasticity in grains collisions and intensity of shearing. In order to analyze non-Newtonian behavior, we need to define the cross thermal ${\varphi }^{\ast }$ conductivity transport coefficient ${\varphi }^{\ast }$ and the reduced normal stresses, defined by

$$q_x={\varphi }^{\ast }{\lambda }_0\partial T/\partial y,{\theta }_i=P_{ii}/p,$$

(24)

where ${\lambda }_0$ is the thermal conductivity for a gas with elastic particle collisions [49]. For hydrodynamics at NS order, $q_x,{\varphi }^{\ast }=0$ and ${\theta }_i=P_{ii}/p=1$. The degree of divergence from these numerical values quantifies eventual departures from NS hydrodynamics. The appearance of $q_x\ne 0$ and ${\theta }_i\ne 1$ occurs with the emergence in the fluxes of terms of second order in the gradients (Burnett order). For our geometry, the only remaining second order terms are of form $(\partial T/\partial y)(\partial u_x/\partial y)$, ${\partial }^2{u}_x/\partial y^2$ in the heat flux and of the form ${\partial }^{2}T/\partial y^2$, ${(\partial T/\partial y)}^2$, ${(\partial u_x/\partial y)}^2$ in the stress tensor (moment flux) respectively (the rest of the gradient second order terms in the fluxes are null for the geometry of our system) [48,49].

In Figure 2, we present the theory-simulation comparison for reduced viscosity ${\eta }^{\ast }\left(\alpha \right)$ and effective thermal conductivity ${\beta }^{\ast }\left(\alpha \right)$, for different shearing levels. For sufficiently low shearing, the agreement between theory and simulation is very good for the range of inelasticities represented here. In the limit of no shear, the agreement for the thermal conductivity is excellent from $\alpha =1$ down to $\alpha =0.7$. The viscosity, however, cannot be measured on the reference (null shear) steady states. Of course, this limitation is not intrinsic of the granular gas since the same happens for classic fluids. Nevertheless, a relatively good viscosity measurement can be obtained in the series with lowest shear used in the simulations ($a\simeq 0.019$), for which the agreement at high inelasticities between theory and simulation is actually slightly better than for the thermal conductivity measured at $a=0$. For both coefficients, and at high inelasticities ($\alpha \simeq 0.5$), the agreement is improved when the zeroth order distribution function in the Chapman-Enskog method [49] is not the Maxwellian but an approximation to the solution of the homogeneous state for the heated granular gas (slightly different from the Maxwellian [55]). This approximation improve (represented with dashed lines in Figure 2) consists in taking into account the first term in an expansion in Sonine polynomials around the Maxwellian (for more references addressing the Sonine polynomials expansion method in kinetic theory, see for instance [48,49,56]).

We report in Figure 3a the reduced normal stresses ${\theta }_x\equiv P_{xx}/p,{\theta }_y\equiv P_{yy}/p$ measured in the simulations, as a function of the coefficient of restitution $\alpha$ and for several series with different shearing intensities. As we see, deviations from the NS behavior (${\theta }_i\ne 1$) are small (less than $5\%$) even at moderately high values of shear (close to $a=0.2$). In this case, the coefficient of restitution seems not to have an important impact on the behavior of both ${\theta }_x$ and ${\theta }_y$. In Figure 3b we represent the value of the cross conductivity coefficient ${\varphi }^{\ast }$. As we can see, the deviations from linear behavior (i.e., from ${\varphi }^{\ast }\left(\alpha \right)=0$) are more significant here: up to a maximum of approximately $40\%$ for the highest shear rate series represented here. However, as expected, the NS behavior is recovered for small shear rate (always less than $6.5\%$ for the lowest non-null shear rate series). It is remarkable that, again, departures from NS behavior are stronger in the quasi-elastic limit than for more inelastic gases (i.e., more inelastic flows show less significantly non-Newtonian behavior here). This is consistent with the non-trivial observation in the previous section that $\left|\gamma \right|$ increases as we approach the quasi-elastic limit. As we already noticed, a higher $\left|\gamma \right|$ implies higher Knudsen number, and thus, our observations are self-consistent and not an artifact of the theory. The NS behavior (${\theta }_i=1$, ${\varphi }^{\ast }=0$) is strictly recovered for the reference states ($a=0$), for all values of inelasticity.

6. Conclusions

In summary, our results on the transport properties of our system strongly indicate that NS hydrodynamics may apply for steady flows of inert particles in a heterogeneous medium, as long as the relevant Knudsen number scales imposed by the boundary conditions in the problem do keep small. This is in no way different to the condition that is required to a molecular gas in order to obey hydrodynamics at NS order [53].

More generally, the fact that this result is independent of the degree of inelasticity supports the hypothesis that inelasticity by itself does not cause necessarily the breakdown of scale separation in granular gases such as it appears in the Homogeneous Cooling State. A recent work on the aging time to hydrodynamics for the homogeneous cooling state of the granular gas also supports this hypothesis [57]. The dynamics is actually a bit more complex and the energy source necessary to fluidize the system may cancel out the gradients inherently produced by inelasticity, bringing the granular gas back to the realm of NS hydrodynamics. In fact, theory-simulation comparisons seem to indicate that a part of the disagreement between theory and simulation can be ruled out when we consider an improve in the approximation to the reference distribution function in the Chapman-Enskog method (when a Sonine polynomial expansion to first order [48,56] for the homogeneous cooling state [58] is used instead of the Maxwellian [55]) and thus, the relatively small disagreement would not be due to an eventual failure of hydrodynamics but to the inherent weakness of the perturbative solution used in the Chapman-Enskog procedure [18,49] for calculation of the transport coefficients.

Therefore, there is no a priori reason for concern respect to the validity of hydrodynamics for granular gases. This would depend essentially on the specific properties and geometry of the flow of interest, not just the inelasticity.

Considering a granular gas fluidized by an active fluid should be regarded, not as an artifact, but as a natural situation since, in the first place, rapid granular flows are possible just because an energy input is able to compensate for inelastic cooling (otherwise, clustering instabilities and total freezing of the dynamics occur [19]). However, an interstitial fluid is a very frequent configuration for granular dynamics problems present in nature [59]. Applications to biology problems are therefore promising. Moreover, it has been experimentally proved that this type of fluctuating force like the one we use to model the active fluid reaches local equilibrium with a granular particle [35]. Based on these observations, our model for energy input from the active fluid should be closer, in most cases, to experimental situations [60] where, interestingly, particle clustering instabilities also may occur.

In addition, the results in the present work have helped us to produce further results, since we have already applied NS hydrodynamics with a non-uniform stochastic force also to problems and applications like granular segregation [61], where excellent agreement is found between granular impurity segregation criteria resulting from NS theory and computer simulations. Finally, we plan to extend this work by developing the corresponding non-Newtonian hydrodynamic theory and also by calculating the linear stability criteria of theses base flows. On the other hand, this work could be relevant for the study of transport in crowded active enviroments [7] where the time scales of the active fluids is faster than the inert particles relaxation time. The existence of a well defined NS hydrodynamics in heterogeneous active bath [8], is also an important result in order to explain biological systems. Another possible future work is to study the effect of sorting of active swimmers with inert particles in heterogeneous active media. Also the study of the behaviour of particles, mean square displacement, mobility, diffusion coefficient depending on the media. A study of the hydrodynamics of more complex forms of the volume force, for example with time dependency, such as for a traveling wave [9], is also a promising line of research.