## Abstract

The entropy inequality, commonly taken as an axiom of continuum mechanics, is found to be spontaneously violated in macroscopic granular media undergoing collisional dynamics. The result falls within the fluctuation theorem of nonequilibrium thermodynamics, which is known to replace the Second Law for finite systems. This phenomenon amounts to the system stochastically displaying negative increments of entropy. The focus is on granular media in Couette flows, consisting of monosized circular disks (with 10 to 104 disks of diameters 0.01 m to 1 m) with frictional-Hookean contacts simulated by molecular dynamics accounting for micropolar effects. Overall, it is determined that the probability of negative entropy increments diminishes with the Eulerian velocity gradient increasing, while it tends to increase in a sigmoidal fashion with the Young modulus of disks increasing. This behavior is examined for a very wide range of known materials: from the softest polymers to the stiffest (i.e., carbyne). The disks’ Poisson ratio is found to have a weak effect on the probability of occurrence of negative entropy increments.

## 1 Introduction

Recent works [1,2] have shown that Couette flows of molecular fluids exhibit violations of the entropy inequality. In other words, the randomly evolving entropy does sometimes have negative increments (Fig. 1(a)), while it is a positive-valued function on average, and its time-integral is a weakly increasing function on average. Here, “sometimes” means randomly, while “on average” generally means one of these three cases: (i) volume averaging; (ii) ensemble (or statistical) averaging; and (iii) time averaging. If the averaging is done over sufficiently large, respectively, (i) volumes, (ii) ensembles, or (iii) time intervals, then the negative increments of entropy vanish and one is back to conventional thermodynamics. The adjective “sometimes” is another way of saying that the Second Law is spontaneously being violated—the phenomenon of these spontaneous violations has been intensively studied in nonequilibrium thermodynamics since 1993 [3]. A word of caution: the Second Law is commonly understood to apply to systems of infinitely many elements (e.g., atoms), but in this area, it is often used to apply to finite systems.

*t*, quantifying the thermodynamic reversibility of a trajectory taken by a thermodynamic system, and

*A*is the value of $\varphi t$, Fig. 1(b).

We address this question in this article: Under what circumstances do macroscopic granular flows also follow the fluctuation theorem? To this end, we investigate granular media in Couette flows, consisting of monosized circular disks (with 10–104 disks of diameters 0.01 m to 1 m) with frictional-Hookean contacts simulated by molecular dynamics accounting for micropolar effects.

A micropolar random medium model is set up to assess the energy flow and possibility of negative entropy increments. The probability of such events is studied under several influences: Eulerian velocity gradient and grain elasticity. The present study expands the very recent investigation [5], while continuing the line of studies [6,7]. Note that a fluidized granular medium was experimentally found to display positive and negative fluctuations in the power flux according to the fluctuation theorem [8].

## 2 Stochastic Homogenization of a Random Granular Medium

The granular matter evolves according to deterministic laws of Newtonian dynamics. If we are interested in a description in terms of just a few parameters, we need to homogenize it. What continuum description do we adopt? First, consider these two extremes:

Case 1. Number of grains

*N*is (very) small: from a few to a few hundred;Case 2. Number of grains is (almost) infinite:

*N*→ ∞.

In case 1, the medium’s response—e.g., the overall stress under a prescribed Eulerian velocity gradient—displays random fluctuations, which can be best described by a stochastic process. Such fluctuations are typical of mesoscale physics. In principle, albeit with a great difficulty, this can be studied by physical experiments although only computer experiments will be employed here. There also exists a theoretical possibility of prescribing a stress state, but it would be only doable on a computer.

*B*

_{N}(

*ω*) parametrized by sample events

*ω*of the space of elementary events Ω. The parameter

*N*plays the role of a mesoscale $\delta =N1/d$ in a

*d*-dimensional problem. We are working with planar Couette flows, so

*d*= 2.

In case 2, the fluctuations tend to diminish as *N* increases and they vanish as *N* → ∞. Determining the response as a function of *N* is one of the key objectives.

*u*

_{i}, $u\u02d9i\u2261vi$, $\phi i$, and $\phi \u02d9i\u2261wi$. For every

*B*

_{N}(

*ω*), in a continuum picture, the fields (3) lead to a micropolar medium in which the linear momentum, angular momentum, and energy conservation equations are as follows:

*I*the inertia of the grain,

*u*the internal energy density (per unit mass),

*q*

_{i}the heat flux,

*d*

_{ij}:=

*v*

_{i},

_{j}Eulerian velocity gradient, $\tau ji$ the Cauchy force-stress tensor,

*κ*

_{ji}:=

*w*

_{i},

_{j}the gradient of spin rate, and

*μ*

_{ji}the couple-stress tensor; ɛ

_{ijk}is the Levi-Civita symbol. Note that the last term mentioned earlier is twice the dual vector of $\tau jk$ scalarly multiplied with the rotation. The stress tensors are defined by linear transformations from the surface unit normal (

*n*

_{j}) into the surface force and moment tractions $ti(n)=\tau jinj$ and $mi(n)=\mu jinj$.

_{2}) must be nonnegative, but its contribution to granular flows under study (i.e., heating of the grains) has a negligible effect on the dynamics of the granular medium and shall, henceforth, be disregarded.

In the third relation, Eq. (6_{3}), we have introduced the dissipation function $\varphi $ (per unit mass) which, in view of Eq. (5_{3}), must be nonnegative. Equation (6_{3}) is the Clausius–Duhem (CD) inequality of interest to us, expressing the Second Law in the setting of micropolar continuum mechanics concepts.

If the size of the granular system (i.e., the number of grains *N*) is small, there are random fluctuations in the system—hence, the dissipation evolves as a stochastic process. But, is it always nonnegative? In other words, we ask: Does the inequality (6_{3}) always hold for any aggregate of *N* particles?

_{3}) needs to be treated as a volume average (denoted by an overbar as $f\xaf:=(1/V)\u222bfdV$), leading us to write:

*d*

_{ki}:=

*v*

_{i},

_{k}, we recall that, under the kinematic uniform boundary condition,

*B*of volume

*V*= |

*B*| under this boundary condition:

*N*is the number of grains in the volume

*V*, while

*N*′ includes periodic image (ghost) grains outside

*V*when periodic boundary conditions are used in the simulations in the following sections. When periodic boundary conditions are not used,

*N*′ =

*N*. Next,

*m*

^{(n)}, $rk(n)$, $vk(n)$, and $fni$ are the mass, position, velocity, and force vectors of grain

*n*, respectively. The body

*B*, which contains

*N*grains, is typically meant to be the hypothetical representative volume element (RVE). In the present study,

*B*signifies, more generally, the randomly evolving statistical volume element (SVE).

*w*

_{i}($\u2261\phi \u02d9i$), equals the volume-averaged one $\kappa ki\xaf$. Next, under the boundary condition (11), consider the volume average of

*μ*

_{ji}

*w*

_{i},

_{j}≡

*μ*

_{ji}

*κ*

_{ji}:

*k*

_{R}stands for the kinetic energy density of rotational motion, while the (generally nonsymmetric) applied curvature–torsion rate tensor $\kappa ki0$ is conjugate to the volume averaged couple-stress tensor

Note that, if *k*_{R} and $\epsilon ijkTjkwi\xaf$ are disregarded in Eq. (12), one obtains the homogenization condition of Hill–Mandel type, $\mu jiwi,j\xaf=\mu ik\xaf\kappa ki0$, which applies to smoothing of slowly deforming inhomogeneous inelastic media of micropolar type.

**v**= (

*v*

_{1},

*v*

_{2}, 0) and $\phi =(0,0,\phi 3)$. The velocity of the top plate is $v1(h)=d210h$,

*h*being the channel height), while the applied gradient of spin is $\kappa ki0=0$. Then, the

*volume averaged dissipation function*simplifies to

## 3 Planar Couette Flows of Disks

*i*due to particle

*j*are

*k*

_{n}(

*k*

_{t}) is the elastic constant in normal (tangential) contact, $\gamma n$ ($\gamma t$) is the viscoelastic damping constants in normal (tangential) contact, $\delta ij$ is the overlap distance of two particles,

*m*

_{eff}=

*M*

_{i}

*M*

_{j}/(

*M*

_{i}+

*M*

_{j}) is the effective mass of two particles of masses

*M*

_{i}and

*M*

_{j}, $utij$ is the elastic tangential displacement,

**n**

_{ij}is the unit vector along the line connecting the centers of two particles, and $vijn$ ($vijt$) is the normal (tangential) component of the relative velocity of two particles. When the distance

*r*between two particles of radii

*R*

_{i}and

*R*

_{j}is greater than their contact distance

*d*=

*R*

_{i}+

*R*

_{j}, there is no force between the particles.

In all the simulated cases, the disks have the same friction coefficient *μ* = 0.1 and 2d density $\rho =200kg/m2$, while, as will see in Secs. 4 and 5, several other parameters are varied.

*k*

_{t}= 2/7

*k*

_{n}and $\gamma t=1/2\gamma n$, where $\gamma n$ is set such as to make

*e*= 0.5 for the normal restitution coefficient $e=exp(\u2212\gamma n\pi /4kn/(meff\u2212\gamma n2))$. So as to capture the collisions properly, the simulations are run at the time-step d

*t*=

*t*

_{col}/50, where

Due to a finite number of disks, the system behavior exhibits random fluctuations—a typical case of stochastic mesoscale physics, with $\rho \varphi \xaf$, being a stochastic process. A homogeneous, deterministic fluid-like behavior may (but does not need to) be expected in the limit of an inifnite number of particles. Each simulation is run under periodic boundary conditions at the entrance and exit ends of the channel.

*P*

^{(−)}+

*P*

^{(+)}= 1, where $P(+)=\u222b0\u221ep\rho \varphi \xaf(x)dx$ is the probability of no violations, i.e., the probability of no negative entropy increments.

To assess *P*^{(−)}, we run LAMMPS simulations in a Monte Carlo sense. Depending on the case, we find *P*^{(−)} > 0. These results stand in contradistinction to the conventional continuum mechanics view, which holds that the CD inequality is always true, i.e., that one should have *P*^{(+)} = 1 and *P*^{(−)} = 0.

To gain insight into the probability of violations of the CD inequality (7) for granular systems, we distinguish several control parameters:

Disk diameter (

*d*= 0.01, 0.1, 1 m)Number of disks in a given simulation (

*N*= 10, 24, 50, 104)Area fraction of disks

Rate effect (top plate speed)

Disk elasticity effect

*t*

_{R}until that state sets in. Then, we divide the total simulation into blocks in which we determine the average values. The nonequilibrium steady state is achieved when fluctuations in the cumulative average are no longer statistically significant. The system in this state is described by the asymptotic steady-state fluctuation theorem

*A*is any positive value of the dissipation function.

## 4 Rate Effects

The rate effect in a Couette flow is controlled by the speed of the top plate *v*_{plate}. Here, we report what was observed as that speed varies. The probability of negative excursions *P*^{(−)} for three different size systems (*N* = 10, 24, 50) with the same diameter of the disks, *d* = 1 m, have been go through. To investigate the effect of upper plate’s speed, for each system size *N* simulated $vplate=d210$*h*, where $d210$ is set in a range from 0.1 to 50 and *h* is the channel height (see the Appendix, Table 1). Figure 3 shows a summary planar graph of the output stored data. All the data that have been archived after the system had reached the steady state. To attain that state, the simulations have been performed for 10^{7} d*t* increments with the time-step d*t* = *t*_{col}/50, so as to capture the collisions properly. Simulations indicate that when the velocity of the upper plate is relatively low ($d210<1$), the granular flow starts really slowly, and for this reason, the simulations required to run for a relatively long time (10^{7}*dt*), even if systems are very small (*N* = 10). Thus, as the speed of the top plate decreases, the probability of negative entropy increments tends to increase although it does correctly go down to zero as *v*_{plate} → 0.

## 5 Grain Elasticity Effects

Taking into account that the disk’s Young’s modulus *E* determines the interparticle contact stiffness, we now explore the effect of *E* over the entire range of known material properties from 0.01 GPa (rubber) to 32,100 GPa (carbyne).

Results for three different diameters of disks (*d* = 1 m, *d* = 0.1, and *d* = 0.01 m), and several different system sizes *N* have been archived (see the Appendix, Table 2). Figures 4–6 show graph of *P*^{(−)}. In all the cases, no negative increments $\rho \varphi \xaf$ have been observed for *E* = 0.01 GPa, so our plots start at 0.1 GPa. We observe that, overall, there is a sigmoidal dependence on *E*, with a decrease in Young’s modulus, i.e., the softening of disks, resulting in no violations of the entropy inequality.

d (m) | N | E (GPa) | P^{(−)} |
---|---|---|---|

1 | 10 | 32,100 | 0.0806 |

10 | 1,050 | 0.083 | |

10 | 50 | 0.076 | |

10 | 50 | 0.076 | |

10 | 10 | 0.049 | |

10 | 1 | 0.001 | |

10 | 0.1 | 0 | |

0.1 | 10 | 32,100 | 0.179 |

10 | 1,050 | 0.212 | |

10 | 50 | 0.186 | |

10 | 10 | 0.01 | |

10 | 1 | 0.003 | |

10 | 0.1 | 0 | |

24 | 32,100 | 0.0036 | |

24 | 1,050 | 0.022 | |

24 | 50 | 0.02 | |

24 | 10 | 0 | |

0.01 | 10 | 32,100 | 0.239 |

10 | 1,050 | 0.245 | |

10 | 50 | 0.24 | |

10 | 10 | 0.135 | |

10 | 1 | 0.018 | |

10 | 0.5 | 0.005 | |

10 | 0.1 | 0 | |

24 | 32,100 | 0.09 | |

24 | 1,050 | 0.093 | |

24 | 50 | 0.058 | |

24 | 10 | 0.005 | |

24 | 5 | 0.0015 | |

24 | 1 | 0 | |

50 | 32,100 | 0.0105 | |

50 | 1,050 | 0.0132 | |

50 | 50 | 0.012 | |

50 | 10 | 0.001 | |

50 | 5 | 0 | |

104 | 32,100 | 0.00135 | |

104 | 1,050 | 0.0015 | |

104 | 50 | 0.0011 | |

104 | 10 | 0 |

d (m) | N | E (GPa) | P^{(−)} |
---|---|---|---|

1 | 10 | 32,100 | 0.0806 |

10 | 1,050 | 0.083 | |

10 | 50 | 0.076 | |

10 | 50 | 0.076 | |

10 | 10 | 0.049 | |

10 | 1 | 0.001 | |

10 | 0.1 | 0 | |

0.1 | 10 | 32,100 | 0.179 |

10 | 1,050 | 0.212 | |

10 | 50 | 0.186 | |

10 | 10 | 0.01 | |

10 | 1 | 0.003 | |

10 | 0.1 | 0 | |

24 | 32,100 | 0.0036 | |

24 | 1,050 | 0.022 | |

24 | 50 | 0.02 | |

24 | 10 | 0 | |

0.01 | 10 | 32,100 | 0.239 |

10 | 1,050 | 0.245 | |

10 | 50 | 0.24 | |

10 | 10 | 0.135 | |

10 | 1 | 0.018 | |

10 | 0.5 | 0.005 | |

10 | 0.1 | 0 | |

24 | 32,100 | 0.09 | |

24 | 1,050 | 0.093 | |

24 | 50 | 0.058 | |

24 | 10 | 0.005 | |

24 | 5 | 0.0015 | |

24 | 1 | 0 | |

50 | 32,100 | 0.0105 | |

50 | 1,050 | 0.0132 | |

50 | 50 | 0.012 | |

50 | 10 | 0.001 | |

50 | 5 | 0 | |

104 | 32,100 | 0.00135 | |

104 | 1,050 | 0.0015 | |

104 | 50 | 0.0011 | |

104 | 10 | 0 |

Finally, the investigation of grain elasticity effects is accompanied by a study of the influence of Poisson’s ratio, varying between −1 and 0.5, on the probabilities of negative $\rho \varphi \xaf$. Overall, Poisson’s ratio has been found to have a weak effect on *P*^{(−)} over its entire range (see the Appendix, Table 3). These results are plotted in Figs. 7 and 8 in which the value of the probability of negative entropy increments is almost constant for the Poisson’s ratio examined range.

d (m) | N | ν | P^{(−)} |
---|---|---|---|

1 | 10 | 0.1 | 0.079 |

10 | 0.2 | 0.08 | |

10 | 0.3 | 0.076 | |

10 | 0.4 | 0.082 | |

10 | 0.5 | 0.084 | |

0.1 | 10 | 0.1 | 0.181 |

10 | 0.2 | 0.181 | |

10 | 0.3 | 0.186 | |

10 | 0.4 | 0.185 | |

10 | 0.5 | 0.189 | |

24 | 0.1 | 0.0189 | |

24 | 0.2 | 0.0182 | |

24 | 0.3 | 0.019 | |

24 | 0.4 | 0.0203 | |

24 | 0.5 | 0.0205 |

d (m) | N | ν | P^{(−)} |
---|---|---|---|

1 | 10 | 0.1 | 0.079 |

10 | 0.2 | 0.08 | |

10 | 0.3 | 0.076 | |

10 | 0.4 | 0.082 | |

10 | 0.5 | 0.084 | |

0.1 | 10 | 0.1 | 0.181 |

10 | 0.2 | 0.181 | |

10 | 0.3 | 0.186 | |

10 | 0.4 | 0.185 | |

10 | 0.5 | 0.189 | |

24 | 0.1 | 0.0189 | |

24 | 0.2 | 0.0182 | |

24 | 0.3 | 0.019 | |

24 | 0.4 | 0.0203 | |

24 | 0.5 | 0.0205 |

## 6 Stochastic Thermomechanics

### 6.1 Random Fields.

*What continuum description do we adopt?*To answer it, note that the fluctuation theorem directly involves the dissipation function $\varphi $. Thus, we employ the thermomechanics with internal variables (TIVs) to generalize the continuum mechanics to account for spontaneous negative entropy increments. TIV is based on the internal (or free) energy

*u*and the dissipation function $\varphi $, which need to be taken as random fields [14] over the material ($B$) and time (

*T*= (−∞,

*t*]) domains:

_{3}), we have lumped the mass density $\rho $ and $\varphi $ together.

For example, in the case of spontaneous violations of the Second Law occurring in heat conduction [15], one works with $\rho s\u02d9th(i)\xaf$, note Eq. (6_{3}). When a tensorial approximation of the constitutive response is possible, the heat conductivity becomes a tensor-valued random field [13].

The time dependence of stochastic fluctuations of the volume average $\rho \varphi \xaf$ ($=\rho s\u02d9m(i)\xaf$) over the channel of Couette flow has recently been investigated in Ref. [5]. With extensive sampling of LAMMPS-generated space-time trajectories (i.e., realizations of this random process), the covariance could be very well approximated by either the Cauchy or Dagum models [16–18]. These trajectories have a fractal dimension *D* ≃ 1.7 and are antipersistent: the Hurst exponent *H* < 0.5.

### 6.2 Axioms of Continuum Thermomechanics.

As outlined in Ref. [7], the axioms (also called *principles*) of thermomechanics need to be modified so as to admit negative entropy production events. In effect, we modify the hierarchy of axioms of continuum (thermo)mechanics theories:

Axiom of Causality: “The future state of the system depends solely on the probabilities of events in the past.” That is, quoting Evans and Searles [4]: “the probability of subsequent events can be predicted from the probabilities of finding initial phases and a knowledge of preceding changes in the applied field and environment of the system.”

Fluctuation theorem is derived from the Axiom of Causality [4].

Upon spatial, statistical, or time averaging, the Second Law in a conventional deterministic form is obtained as a special case of the fluctuation theorem. At that stage, the deterministic continuum mechanics is recovered.

Axiom of Determinism—saying that the present values of the thermokinetic process (stress, heat flux, internal energy, and entropy) depend on the whole history—is dictated by points 1, 2, and 3. The choice of the thermokinetic process depends on the particular physics involved and may take the form of a classical or generalized (e.g., micropolar as in the present study) stochastic theory.

Axiom of Local Action is to be replaced by the choice of continuum model (classical or generalized), the spatial and temporal smoothing (scale dependence). Taken alone, that axiom makes no clear reference to the microstructure of the medium.

Axiom of Equipresence (all the constitutive quantities depend a priori on the same variables) is dropped given that the violation of Second Law may occur in one physical process present in constitutive relations, not all [19].

## 7 Conclusions

Couette flows of macroscopic granular media randomly exhibit negative entropy increments, just like the already studied analogous flows of molecular fluids. This result, tantamount to spontaneous violations of the entropy inequality, is described by the fluctuation theorem of nonequilibrium thermodynamics, which is known to replace the Second Law for finite systems. The phenomenon is understood as a stochastically occurring release of energy stored during collisional dynamics of granular matter. As the number of disks increases, all the fluctuations go down, and hence, the probability of the violations of the entropy inequality decreases, while the conventional Clausius–Duhem inequality is recovered.

Quantitative results are obtained through a range of simulations in the vein of molecular dynamics (discrete element type) accounting for disk rotations and force-and-torque interactions of disks. In each simulated case, the medium consists of monosized circular disks (with 10 up to 104 disks of diameter 0.01–1 m) with frictional-Hookean contacts. The granular matter is homogenized as a micropolar medium in the sense of Hill–Mandel condition. This homogenization delivers the energy flow and, hence, the entropy evolution as function of time.

It is found that the probability of negative entropy increments diminishes when the Eulerian velocity gradient increases.

As the Young modulus of disks increases, the probability of negative entropy increments tends to increase in a sigmoidal fashion for a very wide range of known materials: from the softest polymers to the stiffest (i.e., carbyne). Conversely, the disks’ Poisson ratio is found to have a negligible effect in this study.

Assessment of two-point statistics of any random process requires more extensive sampling of realizations than for the one-point statistics. Such sampling has recently been carried out on LAMMPS-generated realizations of Couette flows [5]. The covariances can be approximated by the so-called Cauchy or Dagum models. In fact, parametric studies reveal that the space-time entropy evolutions possess fractal dimension of *D* ≃ 1.7 and are antipersistent (with the Hurst exponent *H* < 0.5).

## Acknowledgment

Expert comments on LAMMPS by S. J. Plimpton (Sandia National Laboratories) and on Couette flows by R. E. Khayat (University of Western Ontario) are appreciated. R. L. acknowledges support by the University of Messina. The work of MO-S was funded, in part, by the National Institute of Biomedical Imaging and Bioengineering of the NIH under award R01EB029766.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.

### Appendix

## References

**29**, 359, 2017. 10.1007/s00161-015-0451-4