Rheology of dilute cohesive granular gases

TL;DR

This paper investigates the rheological behavior of dilute cohesive granular gases, deriving flow curves from kinetic theory and analyzing stability, clustering, and viscosity characteristics through theoretical and numerical methods.

Contribution

It provides a theoretical and numerical analysis of the rheology of dilute cohesive granular gases, including flow curve derivation and stability conditions.

Findings

Stable shear flow exists only above a critical shear rate.

Viscosity in stable flow matches that of hard-core granular particles.

Below critical shear rate, clustering occurs with viscosity approximated by cluster mean diameter.

Abstract

Rheology of a dilute cohesive granular gas is theoretically and numerically studied. The flow curve between the shear viscosity and the shear rate is derived from the inelastic Boltzmann equation for particles having square-well potentials in a simple shear flow. It is found that (i) the stable uniformly sheared state only exists above a critical shear rate and (ii) the viscosity in the uniformly sheared flow is almost identical to that for uniformly sheared flow of hard core granular particles. Below the critical shear rate, clusters grow with time, in which the viscosity can be approximated by that for the hard-core fluids if we replace the diameter of the particle by the mean diameter of clusters.

Equations

$$U(r)=\{\begin{array}{ll}\infty & (r\le d)\\ -\epsilon & (d<r\le \lambda d)\\ 0& (r>\lambda d)\end{array},$$

$${\overleftrightarrow{P}}^k=\frac{1}{L^3}\sum _{i=1}^{N}m{\mathbf{V}}_i{\mathbf{V}}_i.$$

$$T_{\mathrm{B}}^{\ast }$$

$${\eta }_{\mathrm{B}}^{\ast }$$

$$\frac{1}{M_2}={\alpha }_1-{\alpha }_2(t-t_{\mathrm{c}}),$$

$${\eta }_{\mathrm{eff}}$$

$${\dot{\gamma }}_{\mathrm{eff}}$$

$$c_k\propto k^{-\beta }$$

$$\frac{d{c}_k(t)}{dt}=\frac{1}{2}\sum _{i+j=k}{K}_{i,j}{c}_i(t)c_j(t)-c_k(t)\sum _{j=1}^{\infty }{K}_{k,j}{c}_j(t).$$

$$M_2(t)=\frac{M_2(0)}{1-2M_2(0)K_{0}t},$$

$$\left(\frac{\partial }{\partial t}+{\mathbf{v}}_1\cdot \nabla \right)f(\mathbf{r},{\mathbf{v}}_1,t)=J({\mathbf{v}}_1|f),$$

$$J({\mathbf{v}}_1|f)$$

$$\times [\mathcal{J}\sigma (\chi ,v_{12}^{\prime \prime })f(\mathbf{r},{\mathbf{v}}_1^{\prime \prime },t)f(\mathbf{r},{\mathbf{v}}_2^{\prime \prime },t)$$

$$-\sigma (\chi ,v_{12})f(\mathbf{r},{\mathbf{v}}_1,t)f(\mathbf{r},{\mathbf{v}}_2,t)]$$

$$+\int d{\mathbf{v}}_2\int d\hat{\mathbf{k}}\Theta (\tilde{b}-\min (\lambda ,\mathfrak{N}))v_{12}$$

$$\times [\sigma (\chi ,v_{12}^{\prime \prime })f(\mathbf{r},{\mathbf{v}}_1^{\prime \prime },t)f(\mathbf{r},{\mathbf{v}}_2^{\prime \prime },t)$$

$$-\sigma (\chi ,v_{12})f(\mathbf{r},{\mathbf{v}}_1,t)f(\mathbf{r},{\mathbf{v}}_2,t)].$$

$$\{\begin{array}{l}{\mathbf{v}}_1={\mathbf{v}}_1^{\prime \prime }-A({\mathbf{v}}_{12}^{\prime \prime }\cdot \hat{\mathbf{k}})\hat{\mathbf{k}}\\ {\mathbf{v}}_2={\mathbf{v}}_2^{\prime \prime }+A({\mathbf{v}}_{12}^{\prime \prime }\cdot \hat{\mathbf{k}})\hat{\mathbf{k}}\end{array}$$

$$A=\{\begin{array}{ll}1-\frac{1}{2}(1-e){\mathfrak{N}}^2\frac{{\cos }^2{\theta }_{\mathrm{c}}}{{\cos }^2\theta }& (\tilde{b}\le \min (\lambda ,\mathfrak{N}))\\ 1& (\tilde{b}>\min (\lambda ,\mathfrak{N}))\end{array},$$

$$\left({\partial }_t-\dot{\gamma }{V}_{1y}\frac{\partial }{\partial V_{1x}}\right)f({\mathbf{V}}_1,t)=J({\mathbf{V}}_1|f),$$

$${\partial }_t{P}_{\alpha \beta }^k+\dot{\gamma }({\delta }_{\alpha x}{P}_{y\beta }^k+{\delta }_{\beta x}{P}_{\alpha y}^k)=-{\Lambda }_{\alpha \beta },$$

$$\overleftrightarrow{\Lambda }\equiv -m\int d{\mathbf{v}}_1{\mathbf{V}}_1{\mathbf{V}}_{1}J({\mathbf{V}}_1|f).$$

$$f(\mathbf{V})=f_{\mathrm{M}}(\mathbf{V})\left[1+\frac{m}{2T}\left(\frac{P_{\alpha \beta }^k}{p^k}-{\delta }_{\alpha \beta }\right)V_{\alpha }{V}_{\beta }\right],$$

$$f_{\mathrm{M}}(\mathbf{V})=n{\left(\frac{m}{2\pi T}\right)}^{3/2}\exp \left(-\frac{m{V}^2}{2T}\right).$$

$$\overleftrightarrow{\Lambda }$$

$${\partial }_t{p}^k+\frac{2}{3}\dot{\gamma }{P}_{xy}^k=-\zeta p^k,$$

$${\partial }_t\Delta P^k+2\dot{\gamma }{P}_{xy}^k=-\nu \Delta P^k,$$

$${\partial }_t{P}_{xy}^k+\dot{\gamma }\left(p^k-\frac{1}{3}\Delta P^k\right)=-\nu P_{xy}^k.$$

$$\dot{\gamma }=\sqrt{\frac{3}{2}\frac{{\nu }^2\zeta }{\nu -\zeta }}$$

$$=n{d}^2\sqrt{\frac{\pi \epsilon }{m}}{\left(\frac{\epsilon }{T}\right)}^3\left[(1-e{)}^{1/2}{\dot{\gamma }}^{(1)\ast }+\mathcal{O}\left((1-e)\right)\right],$$

$$P_{xy}^k=-\frac{p^k}{\nu }\sqrt{\frac{3}{2}\zeta (\nu -\zeta )}$$

$$=n\epsilon {\left(\frac{T}{\epsilon }\right)}^{3/2}\left[(1-e{)}^{1/2}{P}_{xy}^{k(1)\ast }+\mathcal{O}\left((1-e)\right)\right],$$

$$\Delta P^k=\frac{3\zeta }{\nu }{p}^k$$

$$=3n\epsilon {\left(\frac{T}{\epsilon }\right)}^{3/2}\left[(1-e)\Delta P^{\ast }(T)+\mathcal{O}\left((1-e{)}^2\right)\right],$$

$${\dot{\gamma }}^{(1)\ast }$$

$$P_{xy}^{k(1)\ast }$$

$$\Delta P^{(1)\ast }$$

$$\eta =\frac{(\nu -\zeta )p^k}{{\nu }^2}$$

Full text

I Introduction

Granular materials, having dissipative interactions between particles, are ubiquitous in daily life as unusual solids, liquids and gases Jaeger1996 . It is important to know rheological properties of granular flows to control the granular materials Jop2006 ; Boyer2011 ; MiDi2004 ; Silbert2001 ; Jenkins1985a ; Jenkins1985b ; Lutsko2005 ; Saitoh2007 ; Gnoli2016 . The rheological properties of granular flows strongly depend on their densities ranging from dilute gases to the jammed solids. When we focus on the rheology of granular flows for the density below the volume fraction $\varphi<0.5$, the description in terms of the Boltzmann-Enskog equation gives quantitatively correct results Brey1998 ; Garzo1999 ; Brilliantov ; Mitarai2007 ; Chialvo2013 ; Hayakawa2017 , while the appropriate theory for the description of denser flow is still controversial Suzuki2015 ; DeGiuli2015 .

So far, most of previous studies on dry granular flows assume that the interactions between grains can be described by repulsive and dissipative forces. Attractive interactions, however, are not negligible for fine powders and wet granular particles Castellanos2005a ; Herminghaus2013 ; Herminghaus2005b ; Mitarai2006 . The origins of such cohesive forces are, respectively, van der Waals force for fine powders and capillary force for wet granular particles. Such attractive forces cause the liquid-gas phase transition and the clustering instability as well as the enhancement of the jamming transition Rahbari2010 ; Chaudhuri2012 ; Gu2014 ; Singh2014 ; Takada2014 ; Saitoh2015 ; Irani2014 ; Irani2016 . Therefore, to know the rheology of flows consisting of cohesive granular particles is important not only for engineers but also for physicists.

In our previous paper Takada2016 , we have developed the systematic kinetic theory of freely cooling dilute cohesive granular particles in terms of the inelastic Boltzmann equation for particles having square-well potentials. Nevertheless, we still need to analyze the rheology of cohesive granular particles under a simple shear flow, because (i) we are interested in a nonequilibrium steady state under the balance between an external force such as shear and the energy dissipation due to inelastic collisions, and (ii) the viscosity of granular flows under a simple shear differs from that for freely cooling granular gases Santos2004 ; Santos2007 . Moreover, we have to consider the contribution of clustering caused by the attractive interaction between grains to the rheology systematically.

In this paper, we try to clarify the rheological properties of dilute granular gases having an attractive interaction described by the square-well potential. The organization of this paper is as follows. In the next section, we explain the setup and the results of the event-driven simulation by DynamO Bannerman2011 . In Sec. III, we derive the shear viscosity in terms of the Boltzmann equation to compare the results with those from the simulation. We also briefly explain the result of the linear stability analysis. In Sec. IV, we summarize our results. Technical details are described in Appendices A–H.

II Molecular dynamics simulation under a simple shear

In this section, we explain our model and the setup of our event-driven simulation in terms of DynamO Bannerman2011 for dilute cohesive granular gases under a uniform shear in Sec. II.1. We present the results of our simulation in Sec. II.2.

II.1 Our model

We consider a collection of monodisperse particles in which the mass and the diameter are, respectively, given by $m$ and $d$. We assume that the interaction between particles is described by the square-well potential

[TABLE]

where $r$, $\varepsilon$ and $\lambda$ are the distance between particles, the well depth and the ratio of the whole potential range to the hard-core repulsive range, respectively. We assume that each collision is inelastic when two particles collide at $r=d$, and collisions are elastic otherwise. Here, the inelasticity is characterized by the restitution coefficient $e$, which is the ratio of the post-collisional relative normal speed to the pre-collisional one. Note that the detailed expressions of collision processes by this potential are presented in Ref. Takada2016 . We also note that we are mainly interested in nearly elastic cases, i.e., $e\lesssim 1$ because the applicability of the kinetic theory for cohesive granular gases is limited in this region Takada2016 . A simple shear flow characterized by the shear rate $\dot{\gamma}$ is applied in $x$-direction under the Lees-Edwards boundary condition Lees1972 . (We show the results under the flat boundary condition and to clarify the artifacts caused by the periodic boundary condition in Appendix A.) The time evolutions of the position $\bm{r}_{i}$ and the velocity $\bm{v}_{i}$ of $i$-th particle are updated by the event-driven simulation for hard-core particles. We mainly simulate the systems of $N=1,372$ particles in a cubic box, whose size is $L=41.6d$. We also simulate the system of $N=10,976$ particles in a cubic box corresponding to $L=83.1d$ to check finite size effects. Throughout this paper, we fix the packing fraction as $\varphi=N(\pi d^{3}/6)/L^{3}=0.01\ll 1$, the inelasticity $1-e=0.01,0.1$, and $0.3$, and the ratio characterizing the potential well $\lambda=1.5$. We measure various quantities by changing the dimensionless shear rate $\dot{\gamma}^{*}\equiv\dot{\gamma}\sqrt{md^{2}/\varepsilon}$. We show the results for $1-e=0.01$ in the main text and present the results for $1-e=0.1$ and $0.3$ in Appendix B to clarify the role of inelasticity.

II.2 Results

Let us present the results of our MD. Figures 1 and 2 exhibit the results of the dimensionless kinetic or granular temperature $T^{*}\equiv T/\varepsilon$ and the shear viscosity $\eta^{*}\equiv\eta d^{2}/\sqrt{m\varepsilon}$ against the dimensionless shear rate $\dot{\gamma}^{*}$, respectively, in steady states above the critical shear rate for $e=0.99$. Here, the shear viscosity is defined by $\eta=-P^{k}_{xy}/\dot{\gamma}$, where the time-averaged kinetic part of the stress tensor $\overleftrightarrow{P}$ is expressed as Alam2003 ; Bannerman2011 :

[TABLE]

Here, $\bm{V}_{i}=\bm{v}_{i}-\dot{\gamma}y_{i}\hat{\bm{e}}_{x}$ is the peculiar velocity of $i$-th particle with the unit vector in $x$-direction $\hat{\bm{e}}_{x}$. Note that the stress tensor for dilute gases should be dominated by the kinetic part even if clustering takes place (see Appendix C). To obtain stabilized data, the stress is time-averaged during the dimensionless time interval $100$.

Figures 1 and 2 indicate the existence of the critical shear rate $\dot{\gamma}_{\rm c}^{*}(\equiv\dot{\gamma}_{\rm c}\sqrt{\varepsilon/md^{2}}=0.023)$ above which there exist steady states. We also plot steady Bagnoldian expressions for the kinetic temperature and the viscosity

[TABLE]

of hard-core dilute granular gases in Figs. 1 and 2 Santos2004 . It is remarkable that Eqs. (3) and (4) give precise results for $\dot{\gamma}>\dot{\gamma}_{\rm c}$ except for the region in the vicinity of $\dot{\gamma}_{\rm c}$. These results can be understood as follows: When the shear rate is sufficiently larger than the critical one, the temperature determined by the energy balance is also larger than the well depth, where the attractive force is negligible. This is the reason why the flow curve reduces to Bagnoldian expressions in the high shear regime.

For small shear rate, there does not exist any steady state. Figure 3 is the time evolution of clustering process observed in our MD for $\dot{\gamma}^{*}=0.01$. The time evolution of the mean cluster size is plotted in Fig. 4, in which two particles belongs to a same cluster when the distance between them is less than $\lambda d$. The mean cluster size is almost unity for large shear rate ($\dot{\gamma}>\dot{\gamma}_{\rm c}$), while it drastically increase after $t^{*}=180$ for small shear rate ($\dot{\gamma}<\dot{\gamma}_{\rm c}$), where we have introduced the dimensionless time $t^{*}\equiv t/\sqrt{md^{2}/\varepsilon}$ and the second moment of the cluster size $M_{2}=\sum_{k=1}^{\infty}k^{2}c_{k}$ with the size distribution $c_{k}$ of the size $k$. It should be noted that the growth rate of the mean cluster size becomes smaller after $t^{*}=250$. It should be noted that $M_{2}$ can be regarded as the mean cluster size because $M_{1}=\sum_{k=1}^{\infty}kc_{k}$ is always equal to the unity. Note that this tendency seems to be insensitive to the system size from the comparison of the results of $N=1,372$ with those of $N=10,976$. The inset of Fig. 4 tries to compare the cluster growth with

[TABLE]

for $t>t_{\rm c}=180\sqrt{md^{2}/\varepsilon}$. Here, we have introduced the fitting parameters $\alpha_{1}$ and $\alpha_{2}$, where the fitting range is $0\leq t^{*}-t_{\rm c}^{*}\leq 15$. The justification of this fitting curve will be discussed later.

Let us introduce the effective shear viscosity and the effective shear rate scaled by $M_{2}$. We adopt two assumptions: First, each cluster can be replaced by a sphere which has the identical size. Next, we ignore the size distribution or size fluctuation of clusters. From these assumptions, the diameter and the mass of the clusters are, respectively, given by $d_{\rm cl}=M_{2}^{1/3}d$ and $m_{\rm cl}=(d_{\rm cl}/d)^{3}m=M_{2}m$. Therefore, we can introduce the effective shear viscosity and the effective shear rate

[TABLE]

respectively, where $\eta_{\rm c}$ is the critical shear viscosity at $\dot{\gamma}_{\rm c}$. We, respectively, plot the time evolution of the temperature versus $\dot{\gamma}_{\rm eff}$ introduced in Eq. (7) and the relationship between Eqs. (6) and (7) in Figs. 5(a) and (b), where the dimensionless effective shear viscosity and shear rate are, respectively, introduced by $\eta_{\rm eff}^{*}\equiv\eta_{\rm eff}d^{2}/\sqrt{m\varepsilon}$ and $\dot{\gamma}_{\rm eff}^{*}\equiv\dot{\gamma}_{\rm eff}\sqrt{md^{2}/\varepsilon}$. In our simulation, the initial temperature $T_{0}$ satisfies Eq. (3) at given $\dot{\gamma}^{*}$. It is noteworthy that the effective viscosity is approximately represented by Bagnoldian expressions in the ranges $10^{-3}\lesssim\dot{\gamma}^{*}\lesssim\dot{\gamma}^{*}_{\rm c}$, though the flow curves evolve with time to access the origin. Note that the kinetic temperature $T^{*}$ is not well approximated by such a crude treatment (see Fig. 5(a)).

We also investigate the cluster size distribution when the clustering proceeds. For small clusters, $c_{k}$ can be well fitted by a power law as

[TABLE]

for $1\leq k\leq 10$, where the exponent is $\beta\simeq 2.67$ (see Fig. 6). This broad size distribution, though the cutoff size is not large, is incompatible with the assumption (ii) in which the sized distribution is negligible. It should be noted that $\beta$ is almost independent of time in the range $180\leq t^{*}\leq 250$ and the system size as shown in Fig. 6, when the mean cluster size drastically increases. We also note that this exponent is insensitive to the shear rate within the range $0.003\leq\dot{\gamma}^{*}\leq 0.015$ and the system size.

We discuss the evolution of cluster size described by Eqs. (5) and (8) in the unstable region. This process might be explained by Smoluchowski’s rate equation Ziff1980 ; Ziff1982 ; Hendriks1983

[TABLE]

Judging from Eqs. (5) and (8), we may use the corresponding coagulation kernel $K_{i,j}=K_{0}ij$. It is noted that this kernel can be applied to systems where all the elements are equally reactive in polymerization processes Ziff1980 ; Ziff1982 ; Hendriks1983 . Because the time evolution of $M_{2}(t)$ can be explicitly solved as $M_{2}(t)=M_{2}(0)/(1-2M_{2}(0)t)$ Ziff1980 , the mean cluster size is given by

[TABLE]

which qualitatively agrees with Eq. (5). We also note that the size distribution in the vicinity of the gelation satisfies $c_{k}\sim k^{-5/2}$ which is similar to Eq. (8). At present, the applicability of Smoluchowski’s equation (9) to our clustering process is not clear. Further investigation along this line will be needed.

Before closing this section, let us briefly summarize the results for smaller $e$ such as $e=0.7$ and $0.9$ (see Appendix B). Even if we are interested in moderately dissipative situations, the qualitative behavior is common, i. e. (i) Bagnoldian expressions can be used for highly sheared cases, (ii) there is a critical shear rate that the uniform state is unstable. In particular, we should note that the the Bagnoldian expression for $\eta_{\rm eff}^{*}$ is still valid in clustering regime for $\dot{\gamma}_{\rm eff}^{*}>10^{-3}$.

III Kinetic theory

In the previous section, we numerically found the existence of the critical shear rate $\dot{\gamma}_{\rm c}$, below which there is no steady state. In this section, let us consider the Boltzmann equation for granular gases having the square-well potential Eq. (1) under a simple shear flow. Although we have considered a shear driven by the Lees-Edwards boundary condition in the simulation, we consider a bulk shear in the treatment of the kinetic theory for simplicity. Thus, we evaluate the steady observables in terms of the Boltzmann equation and compare the theoretical results with those obtained by the simulation in the previous section. We also clarify what determines this critical shear rate. Note that such a theoretical analysis is only possible for nearly elastic cases $e\lesssim 1$.

Let us begin with the Boltzmann equation Chapman

[TABLE]

for a dilute gas consisting of particles interacting through the square-well potential in Eq. (1), where $J(\bm{v}_{1}|f)$ is the collision integral

[TABLE]

Here, we have introduced the step function $\Theta(x)=1$ for $x\geq 0$ and $\Theta(x)=0$ otherwise, the refractive index $\mathfrak{N}\equiv(1+4\varepsilon/mv_{12}^{2})^{1/2}$ Landau ; Goldshtein , $\tilde{b}=b/d$, $v_{12}=|\bm{v}_{12}|=|\bm{v}_{1}-\bm{v}_{2}|$, the Jacobian $\mathscr{J}$ of the transformation between pre-collisional velocities ($\bm{v}_{1}^{\prime\prime}$, $\bm{v}_{2}^{\prime\prime}$) and the post-collisional velocities ($\bm{v}_{1}$, $\bm{v}_{2}$), and the collision cross section $\sigma(\chi,v_{12})$ between particles $1$ and $2$ at the scattering angle $\chi$. For the square-well potential, the relationship between $(\bm{v}_{1}^{\prime\prime},\bm{v}_{2}^{\prime\prime})$ and $(\bm{v}_{1},\bm{v}_{2})$ is written as Takada2016

[TABLE]

where

[TABLE]

with the angle $\theta$ between $\bm{v}_{12}$ and $\hat{\bm{k}}$, and $\theta_{\rm c}$ satisfies $\cos\theta_{\rm c}=\{1-b^{2}/(\mathfrak{N}^{2}d^{2})\}^{1/2}$. This process is equivalent to that used in Ref. Takada2016 . It should be noted that the expression (14) is only valid for nearly elastic cases $1-e\ll 1$.

Let us consider a uniformly sheared flow characterized by $u_{x}=\dot{\gamma}y$, $u_{y}=u_{z}=0$ to derive observables in the steady state. Using the peculiar velocity as $V_{x}=v_{x}-\dot{\gamma}y$, $V_{y}=v_{y}$, $V_{z}=v_{z}$, we can rewrite the Boltzmann equation (11) as

[TABLE]

where we have ignored the spatial fluctuations in Eq. (11). Multiplying $mV_{1\alpha}V_{1\beta}$ with Eq. (15) and integrating over $\bm{V}_{1}$, we obtain the time evolution of the kinetic stress tensor

[TABLE]

where $P^{k}_{\alpha\beta}\equiv\int d\bm{v}mV_{\alpha}V_{\beta}f(\bm{V},t)$ is the kinetic stress tensor and $\Lambda_{\alpha\beta}$ is defined by

[TABLE]

We assume that the velocity distribution function is given by Grad’s moment method Grad1949 ; Herdegen1982 ; Jenkins1985a ; Jenkins1985b ; Tij2001 ; Garzo2002 ; Garzo2005 ; Kremer2011 ; Garzo2013 ; Chamorro2015 ; Hayakawa2016_1 ; Hayakawa2016

[TABLE]

where we adopt Einstein’s rule for Greek indices where duplicated indices take summation over $x$, $y$, and $z$. Here, we have introduced the pressure $p^{k}$ defined by $p^{k}\equiv(P^{k}_{xx}+P^{k}_{yy}+P^{k}_{zz})/3$ which satisfies the equation of state for the ideal gas $p^{k}=nT$. We have also introduced the Maxwellian distribution function $f_{\rm M}(\bm{V})$:

[TABLE]

Using the ansatz Eq. (18), $\overleftrightarrow{\Lambda}$ can be divided into two parts: diagonal and non-diagonal parts as

[TABLE]

(the derivation is given in Appendix D), where $1_{\alpha\beta}=1$ for $\alpha=\beta$ and [math] otherwise. Here, $\nu$ and $\zeta$ are, respectively, the frequency given by Eq. (60) and the dissipation rate given by Eq. (56) in Appendix E.

From Eqs. (16) and (20), we obtain the time evolution equation of the pressure, the normal stress difference $\Delta P^{k}\equiv P^{k}_{xx}-P^{k}_{yy}$, and the shear stress $P^{k}_{xy}$ as

[TABLE]

Unfortunately, the observables in Eqs. (21)–(23) cannot be expressed as functions of $\dot{\gamma}$ explicitly. We note that the second normal difference $P^{k}_{yy}-P^{k}_{zz}$ is not included in the above treatment, which is known to exist not only in denser systems Sangani1996 ; Hayakawa2017 but also in dilute systems Tsao1995 . However, as shown in Appendix F, $P^{k}_{yy}-P^{k}_{zz}$ is much smaller than $\Delta P^{k}$ in this system. Moreover, if we adopt the linearized approximation from the Maxwellian in the evaluation of the nonlinear collision integral, the second normal stress difference disappears in the dilute gas. This is the reason why we only consider the set of $p^{k}$, $P^{k}_{xy}$, and $\Delta P^{k}$.

Let us focus on the steady state. From Eqs. (21)–(23), we can express the shear rate, shear stress tensor, and stress difference as a function of $T$ in the steady state as

[TABLE]

[TABLE]

respectively, where the quantities with asterisk are dimensionless variables such as

[TABLE]

Here, $\nu_{1}^{(0)*}$ and $\zeta_{1}^{(1)*}$ are given by Eqs. (63) and (59), respectively. From Eqs. (24) and (25), we obtain the shear viscosity

[TABLE]

where

[TABLE]

Here, $\nu_{1}^{(1)*}$ and $\nu_{2}^{(1)*}$ are given by Eqs. (64) and (65), respectively. We note that Eq. (24) determines the relationship between the shear rate and the temperature, where it is easy to express the shear rate as a function of the temperature, though the actual control parameter is the shear rate. Then, the relationships between the shear rate and the other observables are also parametrically plotted in terms of the temperature in Figs. 7–9. Figure 10 plots the temperature dependence of $\nu^{*}\equiv\nu\sqrt{md^{2}/\varepsilon}$ and $\zeta^{*}\equiv\zeta\sqrt{md^{2}/\varepsilon}$ whose expressions are presented in Appendix E. The steady expressions in the simple shear flow can only exist above the lower bound temperature $T_{\rm LB}$ which is determined by $\nu(T)=\zeta(T)$ (see Eqs. (24) and (25)). For our choice of parameters $e=0.99,\lambda=1.5$, the lower bound temperature is given by $T_{\rm LB}^{*}(\equiv T_{\rm LB}/\varepsilon)=0.0144$.

There are two branches for the theoretical $\eta$ above the critical shear rate, though $\eta$ on the lower branch is linearly unstable as explained in Appendix G. The shear rate and the stable viscosity, respectively, tend to $\dot{\gamma}^{2}/\omega_{\rm HC}^{2}\to(5/4)(1-e)$ and $\eta/\eta_{\rm HC}\to 1-(5/6)(1-e)$ for $T>T_{{\rm c}}$ and $e\to 1$, where $\omega_{\rm HC}=(16/5)nd^{2}\sqrt{\pi T/m}$ and $\eta_{\rm HC}=5/(16d^{2})\sqrt{mT/\pi}$ are the collision frequency and the shear viscosity of dilute hard-core gases, respectively Santos2004 . It is remarkable that the upper branches in Figs. 7 and 8 reproduces well the simulation results. We also note that the temperature difference $\Delta T=\Delta P^{k}/n$ obtained from both the simulation and the kinetic theory of hard-core dilute granular gases

[TABLE]

agree in high sheared regime, though the theoretical prediction does not agree with the simulation in the low shear regime (this relationship can be derived easily by Ref. Santos2004 ), while the deviations become large near the critical shear rate $\dot{\gamma}\approx\dot{\gamma}_{\rm c}$. Then, $\Delta T/T=5(1-e)/(3-e)$ with the aid of Eq. (3) (Fig. 9).

Let us evaluate the critical shear viscosity $\dot{\gamma}_{\rm c}$ below which the steady state does not exist. The critical condition is given by $\partial\eta/\partial\dot{\gamma}=(\partial\eta/\partial T)/(\partial\dot{\gamma}/\partial T)\to\infty$ at $\dot{\gamma}_{\rm c}$, which is reduced to $\partial\dot{\gamma}/\partial T=0$. As shown in Appendix H, this critical temperature $T_{\rm c}(\equiv T_{\rm c}/\varepsilon)=0.910$ corresponds to the critical shear rate $\dot{\gamma}_{\rm c}^{*}=0.0228$. The theoretical critical shear viscosity agrees with the numerical critical shear viscosity.

We perform the linear stability analysis for the sheared uniform state (the detailed explanation is given in Appendix G). For simplicity, we ignore the spatial degree of freedom. From the set of equations (21)–(23), the uniform shear state is stable for $T>T_{\rm CL}(=0.04\varepsilon)$ in our choice of parameters. We plot the linearly unstable region as crosses in Figs. 7–9.

IV Conclusion

In this paper, we have performed the event-driven molecular dynamics simulation for cohesive granular gases under a uniform shear and clarified the rheological properties of these particles. We have found that there exists a steady state when the shear rate is larger than the critical shear rate, while the clustering proceeds for lower shear rate. Even for lower shear rate, introducing the effective shear rate and the shear viscosity, we have found that the flow curve can be approximately expressed as the Bagnoldian expression if we replace the diameter of the particle by the mean diameter of clusters. We have obtained two branches for the steady uniformly sheared state from the analysis of the inelastic Boltzmann equation, one of which is consistent with the simulation, and the other branch is linearly unstable.

Acknowledgments

The authors thank Andrés Santos, Koshiro Suzuki, Kuniyasu Saitoh, Takeshi Kawasaki, and Michio Otsuki for their useful comments. One of the authors (ST) wishes to express his sincere gratitude to Tomohiko G. Sano and Thorsten Pöschel for their helpful comments. Numerical computation in this work was partially carried out at the Yukawa Institute Computer Facility. This work is partially supported by Scientific Grant-in-Aid of JSPS, KAKENHI (Grant No. 16H04025 and No. JP16H06478).

Appendix A Simulation under the flat boundary condition

In this Appendix, we examine the applicability of the Lees-Edwards boundary condition from the comparison of the simulations under the flat boundary condition Takada2014 . We prepare two flat walls at $y=\pm L/2$, moving in $x$-direction with the velocities $\pm\dot{\gamma}L/2$. When a particle having $(v_{x},v_{y},v_{z})$ hits the walls at $y=\pm L/2$, the velocity changes to $(v_{x}\pm\dot{\gamma}L/2,-v_{y},v_{z})$ after the collision, respectively.

Figure 11 plots (a) the shear viscosity against the shear rate (7) for $\dot{\gamma}>\dot{\gamma}_{\rm c}$ and (b) the effective flow curve for $\dot{\gamma}<\dot{\gamma}_{\rm c}$ which gives the relationship between $\eta_{\rm eff}$ and $\dot{\gamma}_{\rm eff}$. The result under the flat boundary condition for $\dot{\gamma}>\dot{\gamma}_{\rm c}$ almost agrees with that under the Lees-Edwards boundary condition and the Bagnoldian expression (4), which is contrary to the previous studies under the bumpy boundary condition Liem1992 ; Santos1992 ; Tij2001 . For $\dot{\gamma}>\dot{\gamma}_{\rm c}$, the density and velocity gradient are almost constant except for the boundary layers as shown in Fig. 11(c). We also note that the effective viscosity obtained from the simulation under the flat boundary condition becomes much smaller than that under the Lees-Edwards boundary condition for $\dot{\gamma}<\dot{\gamma}_{\rm c}$. In conclusion, to adopt the Lees-Edwards boundary condition does not cause any artifact for $\dot{\gamma}>\dot{\gamma}_{\rm c}$ while it is controversial for $\dot{\gamma}<\dot{\gamma}_{\rm c}$.

Appendix B Results for highly dissipative cases

In this Appendix, we present the results for larger inelasticity, especially for $e=0.9$, and $0.7$. Figure 12 shows the results of the shear viscosity against the shear rate and the effective shear viscosity (6) against the effective shear rate (7). Even for moderately inelastic case, Bagnoldian expressions (3) and (4) give precise results if the steady state exists above the critical shear rate.

However, the time evolutions of the temperature and the effective shear viscosity (6) against the effective shear rate (7) in Fig. 12 cannot be used, if the steady state is unstable for $\dot{\gamma}<\dot{\gamma}_{\rm c}$ for larger inelasticity.

Appendix C Comparison of the stress tensor with the kinetic stress tensor

In this Appendix, we verify whether the stress tensor is dominated by the kinetic stress. To validate this, we measure the ratio of the stress tensor (2) to the kinetic part of the stress tensor by using our full scratched simulation code. Here, the stress tensor is defined as

[TABLE]

where $\overleftrightarrow{P}^{c}$ is the collisional contribution of the stress tensor, defined by

[TABLE]

Here, $\Delta\bm{p}_{i}$ is the change of the momentum of $i$-th particle during a collision, and $\Delta\tau$ is set to $100\sqrt{md^{2}/\varepsilon}$. We have confirmed that the results are insensitive to the time interval $\Delta\tau$ in the range $10\sqrt{md^{2}/\varepsilon}$ to $100\sqrt{md^{2}/\varepsilon}$. Figure 13 represents the ratios for the pressure, the shear stress, and the pressure difference. All quantities satisfy

[TABLE]

for $\dot{\gamma}>\dot{\gamma}_{\rm c}$, and

[TABLE]

for $\dot{\gamma}<\dot{\gamma}_{\rm c}$. in our simulations. Although the data are a little scattered for the unstable regime ($\dot{\gamma}<\dot{\gamma}_{\rm c}$), we conclude that the collisional contribution to the stress tensor is negligible even for $\dot{\gamma}<\dot{\gamma}_{\rm c}$.

Appendix D Derivation of Eq. (20)

In this Appendix, let us show the detailed derivation of $\overleftrightarrow{\Lambda}$. From the collision rule Eq. (13), the following relationship is satisfied:

[TABLE]

Here, we introduce the dimensionless velocity $\bm{c}_{i}\equiv\bm{v}_{i}/(2\varepsilon/m)^{1/2}$. Inserting Eq. (38) into Eq. (17) and integrating over $\bm{C}=(\bm{c}_{1}+\bm{c}_{2})/2$, $\Lambda$ is expressed as

[TABLE]

where $\bm{c}_{12}=\bm{c}_{1}-\bm{c}_{2}$, $P^{*}_{\alpha\beta}=P_{\alpha\beta}/(nT)$, and we have introduced $\overleftrightarrow{\Lambda_{i}^{*}}$ ($i=1,2,3,4$) as

[TABLE]

respectively. Here, $\chi$ is the scattering angle and we expand $\chi$ in terms of the small inelasticity $1-e$ as

[TABLE]

The explicit forms of $\chi^{(0)}$ and $\chi^{(1)}$ are, respectively, given by Takada2016

[TABLE]

with the dimensionless collision parameter $\tilde{b}\equiv b/d$ and

[TABLE]

To evaluate $\overleftrightarrow{\Lambda_{i}^{*}}$ ($i=1,2,3,4$), the following relations are useful:

[TABLE]

It should be noted that these relations can be derived if we choose $\bm{c}_{12}$ as $z$-axis and use polar coordinates to $\hat{\bm{k}}$. Inserting Eq. (50) into Eq. (40), we rewrite $\overleftrightarrow{\Lambda_{1}^{*}}$ as

[TABLE]

Similarly, from Eqs. (41) and (51), $\overleftrightarrow{\Lambda_{2}^{*}}$ reduces to

[TABLE]

Similarly, $\overleftrightarrow{\Lambda_{3}^{*}}$ and $\overleftrightarrow{\Lambda_{4}^{*}}$ are, respectively, given by

[TABLE]

Inserting Eqs. (52)–(55) into Eq. (39), we can obtain Eq. (20).

Appendix E Derivation of the energy dissipation rate $\zeta$ and the frequency $\nu$

In this Appendix, we derive the expression for the energy dissipation rate $\zeta$ and the frequency $\nu$ in terms of the Boltzmann equation.

First, let us evaluate the energy dissipation rate $\zeta$. The energy dissipation rate can be expanded in terms of the series of the small inelasticity $1-e$ as

[TABLE]

where $\zeta^{(0)}$ and $\zeta^{(1)}$ are, respectively, given by Takada2016

[TABLE]

with

[TABLE]

where $\tilde{b}_{\rm max}=\min(\mathfrak{N},\lambda)$ with the introduction of a function $\min(x,y)$ to select the smaller one between $x$ and $y$.

Next, let us derive $\nu$. From the Appendix D, the frequency $\nu$ is written as

[TABLE]

where $\nu^{(0)}$ and $\nu^{(1)}$ are, respectively, given by

[TABLE]

with

[TABLE]

Unfortunately, $\zeta(T)$ and $\nu(T)$ cannot be expressed explicitly. Therefore, we adopt the numerical integrals in Eqs. (56)–(65) over $\tilde{b}$ and $c_{12}$ for each $T/\varepsilon$.

Appendix F The normal stress differences

In this Appendix, let us show that the second normal stress difference $P^{k}_{yy}-P^{k}_{zz}$ is much smaller than the first normal stress difference $\Delta P^{k}$ in our simulations. Figure 14 represents the shear rate dependence of the ratio of the second normal difference to the first normal stress difference obtained from the simulations, which shows that the ratio satisfies

[TABLE]

This result also validates our treatment that we only consider the set of $p^{k}$, $P^{k}_{xy}$, and $\Delta P^{k}$ in our theoretical treatment.

Appendix G Linear stability analysis

In this Appendix, we study the linear stability of the uniform shear state without any spatial fluctuation Hayakawa2016 . Let us consider the steady state of the dimensionless temperature $T^{*}$, the dimensionless pressure difference $\Delta P^{*}=\Delta P/(n\varepsilon)$, and the dimensionless shear stress $P_{xy}^{*}\equiv P_{xy}/(n\varepsilon)$. We add a perturbation around the steady state as $\delta\phi\equiv(\delta T^{*},\delta\Delta P^{*},\delta P_{xy}^{*})$. Introducing the dimensionless quantities $\dot{\gamma}^{*}\equiv\dot{\gamma}t_{0}$, $\zeta^{*}=\zeta t_{0}$ and $\nu^{*}\equiv\nu t_{0}$ with $t_{0}\equiv\sqrt{md^{2}/\varepsilon}$, the time evolution of the fluctuation is linearized as

[TABLE]

where the matrix ${\cal M}$ is defined by

[TABLE]

with $\dot{\gamma}^{*}+\dot{\gamma}^{*}_{T}\delta T^{*}$, $\zeta^{*}+\zeta^{*}_{T}\delta T^{*}$, and $\nu^{*}+\nu^{*}_{T}\delta T^{*}$ with $\dot{\gamma}^{*}_{T}\equiv(\partial\dot{\gamma}^{*}/\partial T^{*})_{T}$, $\zeta^{*}_{T}\equiv(\partial\zeta^{*}/\partial T^{*})_{T}$, and $\nu^{*}_{T}\equiv(\partial\nu^{*}/\partial T^{*})_{T}$. The Laplace transform of Eq. (67) is expressed as

[TABLE]

where $\varphi(s)={\cal L}[\delta\phi(t)]$ is the Laplace transform of $\delta\phi(t)$. Let us assume $s_{i}$ ($i=1,2,3$) as the eigenvalues of the matrix ${\cal M}$. This yields that the matrices $(s\overleftrightarrow{1}-{\cal M})$ and $(s\overleftrightarrow{1}-{\cal M})^{-1}$ have eigenvalues $s-s_{i}$ and $1/(s-s_{i})$ ($i=1,2,3$), respectively. The inverse Laplace transform of $1/(s-s_{i})$ is given by $\exp(s_{i}t)\Theta(t)$, where $\Theta(t)$ is the step function, i.e., $\Theta(t)=1$ for $t\geq 0$ and $\Theta(t)=0$ otherwise. The system becomes unstable if any one of the eigenvalues have positive real part. Let us calculate the eigenvalues numerically. The determinant of $s_{i}\overleftrightarrow{1}-{\cal M}$ is given by

[TABLE]

where $A_{1}$, $A_{2}$, and $A_{3}$ are, respectively, given by

[TABLE]

[TABLE]

We also focus on the eigenvalue whose absolute value is smallest. Neglecting the terms proportional to $s_{i}^{2}$ and $s_{i}^{3}$ in Eq. (70), we can obtain the linear approximation solution of Eq. (70) as

[TABLE]

Figure 15 presents the temperature dependence of the real part of each eigenvalue $s_{i}$ ($\Re s_{1}>\Re s_{2}\geq\Re s_{3}$), where $\Re s_{i}$ indicates the real part of $s_{i}$. We also plot the linear approximation solution (74). This linearized eigenvalue in Eq. (74) gives a good description of the full linear stability analysis. Note that the real part of eigenvalue becomes positive for $T<T_{\rm CL}\equiv 0.04\varepsilon$ in which the steady state is unstable. Near the critical temperature, the magnitude of the smallest eigenvalue is approximately $10^{-5}$. We note that approximately $10^{5}$ collisions per particle are needed to reach the steady state in this regime.

Appendix H Critical behavior

In this Appendix, we calculate the critical temperature for the linear stability analysis, where $\partial\dot{\gamma}/\partial T$ becomes zero. From Eq. (24), $\partial\dot{\gamma}/\partial T$ satisfies

[TABLE]

This means that $\nu=0$ or

[TABLE]

should be satisfied at the temperature satisfying $\partial\dot{\gamma}/\partial T=0$. As shown in Fig. 10, $\nu$ becomes zero at $T\simeq 0.00451\varepsilon$. Let us calculate other temperatures where $\partial\dot{\gamma}/\partial T$ becomes zero. Using the dimensionless quantities, the left hand side of Eq. (76) can be rewritten as

[TABLE]

Figure 16 shows that Eq. (77) has only one solution $T_{\rm c}=0.910\varepsilon$. The corresponding shear rate becomes $\dot{\gamma}_{\rm c}^{*}=0.0228$ and any steady state does not exist for $\dot{\gamma}<\dot{\gamma}_{\rm c}$.

Let us expand the shear viscosity around the critical temperature $T_{\rm c}$. First, the quantities $\nu$ and $\zeta$ are expanded as

[TABLE]

respectively, where $\nu_{\rm c}=\nu(T_{\rm c})$, $\nu_{\rm c}^{\prime}=(d\nu/dT)_{T_{\rm c}}$, $\nu_{\rm c}^{\prime\prime}=(d^{2}\nu/dT^{2})_{T_{\rm c}}$, $\zeta_{\rm c}=\zeta(T_{\rm c})$, $\zeta_{\rm c}^{\prime}=(d\zeta/dT)_{T_{\rm c}}$, and $\zeta_{\rm c}^{\prime\prime}=(d^{2}\zeta/dT^{2})_{T_{\rm c}}$, respectively. Inserting into Eq. (24) and Eq. (30), the shear rate and the shear viscosity are, respectively, expressed as

[TABLE]

where $\dot{\gamma}_{\rm c}$, $\eta_{\rm c}$, $C_{1}$ and $C_{2}$ are, respectively, given by

[TABLE]

In the vicinity of the critical temperature $T_{\rm c}$, the relationship between the shear viscosity and the shear rate becomes

[TABLE]

Figure 8 shows this expansion with well agrees with both the results of our simulation and the relationship obtained by Eqs. (24) and (30).