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.

(a) Sample shear stress histories with evident negative excursions and (b) histograms of shear stresses
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.
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.
In the third relation, Eq. (63), we have introduced the dissipation function (per unit mass) which, in view of Eq. (53), must be nonnegative. Equation (63) 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 (63) always hold for any aggregate of N particles?
Note that, if kR and are disregarded in Eq. (12), one obtains the homogenization condition of Hill–Mandel type, , which applies to smoothing of slowly deforming inhomogeneous inelastic media of micropolar type.
3 Planar Couette Flows of Disks

(a) A realization B of a simulated granular medium at area fraction = 60% and (b) a continuum picture of spatially homogenized Couette flow. In both cases, the bottom plate is stationary and the upper one is moving at vplate → 0.
In all the simulated cases, the disks have the same friction coefficient μ = 0.1 and 2d density , while, as will see in Secs. 4 and 5, several other parameters are varied.
Due to a finite number of disks, the system behavior exhibits random fluctuations—a typical case of stochastic mesoscale physics, with , 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.
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
4 Rate Effects
The rate effect in a Couette flow is controlled by the speed of the top plate vplate. 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 h, where 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 107 dt increments with the time-step dt = tcol/50, so as to capture the collisions properly. Simulations indicate that when the velocity of the upper plate is relatively low (), the granular flow starts really slowly, and for this reason, the simulations required to run for a relatively long time (107dt), 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 vplate → 0.
P(−) for three different sizes systems (N = 10, 24, 50) with the same diameter of the disks, d = 1 m
N | P(−) | |
---|---|---|
10 | 0.1 | 0.45 |
1 | 0.33 | |
5 | 0.16 | |
10 | 0.076 | |
25 | 0.011 | |
50 | 0.00 | |
24 | 0.1 | 0.45 |
1 | 0.14 | |
5 | 0.01 | |
10 | 0.00 | |
25 | 0.00 | |
50 | 0.00 | |
50 | 0.1 | 0.1 |
1 | 0.00 | |
5 | 0.00 | |
10 | 0.00 | |
25 | 0.00 | |
50 | 0.00 |
N | P(−) | |
---|---|---|
10 | 0.1 | 0.45 |
1 | 0.33 | |
5 | 0.16 | |
10 | 0.076 | |
25 | 0.011 | |
50 | 0.00 | |
24 | 0.1 | 0.45 |
1 | 0.14 | |
5 | 0.01 | |
10 | 0.00 | |
25 | 0.00 | |
50 | 0.00 | |
50 | 0.1 | 0.1 |
1 | 0.00 | |
5 | 0.00 | |
10 | 0.00 | |
25 | 0.00 | |
50 | 0.00 |
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 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.

The probability of negative entropy increments in Couette flow for a range of Young’s modulus from 0.1 GPa to 32, 100 GPa, for d = 1 m

The probability of negative entropy increments in Couette flow for a range of Young’s modulus from 0.1 GPa to 32, 100 GPa, for d = 0.1 m

The probability of negative entropy increments in Couette flow for a range of Young’s modulus from 0.1 GPa to 32, 100 GPa, for d = 0.01 m
P(−) for Young’s modulus range [0.1; 32, 100] GPa
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 . 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.

The probability of negative entropy increments in Couette flow for a range of Poisson’s ratio from −1 to 0.5, for d = 1 m

The probability of negative entropy increments in Couette flow for a range of Poisson’s ratio from −1 to 0.5, for d = 0.1 m
P(−) for Poisson’s ratio range [0.1; 0.5]
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.
For example, in the case of spontaneous violations of the Second Law occurring in heat conduction [15], one works with , note Eq. (63). 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 () 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.