Continuously heated granular gas of elongated particles

Some years ago, Harth et al. experimentally explored the steady state dynamics of a heated granular gas of rod-like particles in microgravity [K. Harth et al. Phys. Rev. Lett. 110, 144102 (2013)]. Here, we report numerical results that quantitatively reproduce their experimental findings and provide additional insight into the process. A system of sphero-cylinders is heated by the vibration of three flat side walls, resulting in one symmetrically heated direction, one non-symmetrically heated direction, and one non-heated direction. In the non-heated direction, the speed distribution follows a stretched exponential distribution p(v) ∝ exp(−(|v|/C)1.5). In the symmetrically heated direction, the velocity statistics at low speeds is similar but it develops pronounced exponential tails at high speeds. In the non-symmetrically heated direction (not accessed experimentally), the distribution also follows p(v) ∝ exp(−(|v|/C)1.5), but the velocity statistics of rods moving toward the vibrating wall resembles the indirectly excited direction, whereas the velocity statistics of those moving away from the wall resembles the direct excited direction.


Introduction
Granular media surround us in everyday life. For instance, while walking on a sandy beach our weight is supported by a complex granular force network that draws our footprint in the sand. Besides, centuries of human technological revolution have made possible the handling of many types of particles and grains [1]. Thus, theoretical approaches of granular media are always needed, due to its massive scale of use in industrial processes, and the role they play in many natural phenomena.
Most of the time we observe granular media to exhibit a solid or liquid-like behaviour. However, when dilute granular systems are sufficiently excited, they enter into a gas-like state, which is characterized by a very low collision rate between the grains [2,3]. The dissipative nature of the interactions differentiates such granular gases from usual molecular gases. Indeed, this property makes these out-of-equilibrium systems interesting from the view of fundamental questions of statistical physics, and sometimes even tricky. Such lose ensembles of granular matter occur, e.g. in planetary rings, but also in industrial processes wherever granular materials are strongly excited.
In the past, kinetic and hydrodynamic theories have described the macroscopic properties of granular gases with and without external energy input. When considering the most simple case, the free cooling of an ensemble of spherical particles with constant restitution coefficient, the granular temperature diminishes in time following Haff's law resulting in an asymptotic decay T ∝ t −2 [4][5][6][7][8]. For freely cooling gases of visco-elastic particles, the exponent of the asymptotic algebraic decay of the granular temperature reduces to T ∝ t −5/3 [9]. However, in many cases the particles' interactions are over-simplified, and even fewer quantitative validations against actual experimental data exist.
Experiments investigating the dynamics of elongated particles in a continuously heated steady state were carried out in microgravity by Harth et al. [11,13]. They found that the velocity distributions in the indirectly excited direction were not far from a Gaussian behavior, but much better approximated by a stretched exponential distribution p(v y ) ∝ exp − v y /C 1.5 . The same applies for the rotational velocity distributions. In the directly heated direction, the low-velocity parts of p(v x ) agree well with p(v y ), while pronounced exponential high-velocity tails are observed. Regarding the spatial distribution of particles, no significant persistent inhomogeneities were observed.
In our contribution to Powders & Grains 2021, we report numerical results that reproduce quantitatively the experimental results introduced in Refs. [11,13] and provide additional insight into the process. Our numerical data complement the experimental analysis, accessing the dynamic evolution of each particle in detail. The paper is organized as follow: In section 2, the numerical model is described briefly. In section 3, the numerical data is provided, while our main results are compared with the experimental data. Finally, several further work ideas are mentioned.

Numerical model
We used a home-made hybrid GPU-CPU discrete element algorithm (DEM) introduced earlier [16,17]. This numerical tool was adapted to simulate confined systems with moving walls. The model considers an ensemble of N = 256 spherocylinders with length L and sphero-radius r, defining an aspect ratio of ζ = L/(2r). The interaction force between two particles i and j, F i j , is obtained using an algorithm of interacting spheropolygons [20]. The particular case of the spherocylinder is described by two vertices and the spheroradius; thus, its surface is delineated by all the points at distance r from the edge defined by the two vertices. The problem of contact detection between the spherocylinders reduces to the determination of the closest point between two edges. Thus, the overlap distance δ results from the simple overlap of two spheres of radius r. For the interaction between the particles and the walls, we used the same model but assuming the interaction of a spherocylinder with an infinite moving plane. Figure 1 displays a snapshot of the system composed of spherocylinders. The walls on the sides and in the back vibrate sinusoidally (colored blue), the others are stationary. The simulation limits are marked by the wire-frame.
The force F i j acting on particle i by the particle j can be decomposed as F i j = F n · n + F t · t, where F n is the component normal to the contact plane and F t acts in tangential direction t. Here, the normal interaction F n is a linear elastic force, depending on the overlap distance δ n between two spherocylinders. The normal dissipation is included using a velocity dependent viscous damping. The total normal force reads F n = −k n δ n − γ n v n rel , where k n is the spring constant in the normal direction, γ n is the normal damping coefficient and v n rel is the normal relative velocity between the contacting particles i and j. The tangential force F t also contains an elastic term and a tangential frictional term accounting for static friction between the grains. Taking into account Coulomb's friction constraint, which reads rel is the tangential component of the relative contact velocity of the overlapping pair. ξ represents the elastic deformation of an imaginary spring with spring constant k t at the contact, which increases as d ξ(t)/dt = v t rel as long as there is an overlap between the interacting particles. The elastic tangential elongation ξ is kept orthogonal to the normal vector (truncated if necessary) [21], µ is the friction coefficient. A velocity Verlet numerical algorithm [22] integrated the 3D translational equations of motion of each particle, while the rotational degrees of freedom were resolved using a modified leap-frog algorithm [23].
Similar to the experimental setup of [11], we de-  25; 4.25] cm. In the experiment described in Ref. [11,13], the dimensions were ∆X = 8.5 cm, ∆Y = 10 cm, and ∆Z = 6.5 cm. The two X-walls, as well as one of the Z walls moved sinusoidally in their normal direction with a fixed amplitude A and frequency f . The other three walls were fixed. The particle aspect ratio was ζ = 9.2, i.e. L = 12 mm, r = 0.65 mm, and mass m = 0.046 g. Moreover, we have used a particle Young modulus of Y = 10 GPa and a restitution coefficient of e n = 0.6. The friction coefficient of the particle-particle contact was µ p = 0.1, while for the particle-wall contacts µ w = 0.2.
In our analysis, we assume that the system is in steady state conditions. Initially, the particles are randomly distributed within the system domain, and the particle speeds are randomly generated in the range [−0.1; 0.1] mm/s. Then, the excitation begins, and after a short transient, the system reaches a steady state characterized by a nonincreasing mean kinetic energy. For each excitation condition (fixing A and f ), the particle velocities are sampled during 20 seconds of simulation time.

Results and discussion
As mentioned before, here we focus on the system kinetics in steady state conditions, where the statistics of the particle linear speeds is described, in detail. The distributions are sampled, including particles away from the walls, at a distance δ > L/2, where L is the longer side of the particles. Figure 2 shows the translational velocity distribution p(v y ) in the non-directly excited direction y (no moving wall), obtained for an oscillation amplitude of A = 1 mm and frequency f = 30 Hz. Note that our numerical data agree very well with the fit to the experimental data [11] shown in blue, p(v y ) = 0.012 exp − v y /C y 1.5 , with C y = 51.77 mm/s. The standard deviation of the data set is σ y = 42.2 mm/s, resembling the experimental value σ exp y = 43.9 mm/s [11]. Note that no significant fat tail is observed on the semi-log scale, which is consistent with the fact that there is practically no missing area in the low velocity range.
This result is in contrast to previous numerical observations of the cooling process of ellipsoids [15] and spherocylinders [17], i.e. without excitation and using periodic boundary conditions. In those cases, Gaussian velocity distributions without fat tails were obtained. However, another experimental study, carried out with similar rod-like particles (in confined conditions), observed non-Gaussian velocity distributions even in the scaling regime of granular cooling [6]. Our validated simulations may aid to reveal the origin of this discrepancy in future work. Figure 3 illustrates the velocity distribution p(v x ), obtained in the directly excited direction x (two moving walls), for oscillation amplitude of A = 1 mm and frequency f = 30 Hz. Once again with reasonable accuracy, we find the statistics agree with a stretched exponential p(v x ) = 0.008 exp − (|v x | /C x ) 1.5 , with C x = 65.5 mm/s and standard deviation σ x = 75.6 mm/s. Note, the same functionality and fitting parameters were used, when describing the experimental data [11], and result in a standard deviation of σ exp x = 70.8 mm/s [11]. Moreover, both the numerical and experimental data exhibit fat tails in the high velocity regime, which are not captured by the stretched exponential fit (see inset of Fig. 3) and Ref. [11]). Here, we can argue that the small deviations come from the slightly different container geometry and from possible artifacts due to the perspective view in the experiment. We also observe a slight accumulation of rods near the nonmoving walls. So, when analysing the whole data (without the δ > L/2 condition), the distributions show a slightly larger probability of slow particles, while the shape of the tails do not change significantly.
Additionally, the simulations also allow accessing the particle velocity in the direction where the system is nonsymmetrically heated (only one moving wall). Note, this data was not accessible in experiments due to technical reasons [11,13]. The statistics corresponding to this degree of freedom is shown in Fig. 4 together with the data of the previously analyzed directions. Remarkably, the velocity distribution p(v z ) becomes asymmetric, reflecting . Speed distribution p(v x ) in the x direction (two moving walls), obtained for A = 1 mm and f = 30 Hz. The inset shows the same but in semi-log scale. In both cases, the best fit of the experimental data Ref [11] is also shown for comparison. the single-sided excitation. The part of the distribution corresponding to negative velocities (towards the moving wall) agrees with p(v y ), whereas the part at positive velocities (away from the moving wall) agrees with p(v x ). Namely, the heated direction exhibits a fat tail, while the non-heated one does not. This result might suggest that the fixed walls behave as an energy sink, removing any trace of the direct excitation once grains have reversed their zdirection of the motion.
Our numerical simulations allow us to perform a systematic study of the system's dynamics when changing the heating intensity. Figure 5  In all cases, the velocity values have been scaled with v c = 2π f A, i.e. the maximum wall velocity. Interestingly, all the scaled speed distributions collapsed on the same curve. The latter indicates that the maximum velocity of the moving walls v c is the relevant characteristic value, denoting that the steady state energy of the system is controlled by the walls. These outcomes would corroborate the guess made in Ref. [11], that the kinetic energy in the system scales with v 2 c ∝ A 2 f 2 , and not the peak acceleration of the walls. It is interesting to note that the energy gain of a single rod bouncing on a plate in this range of excitation parameters under normal gravity rather scales ∝ A 1.5 f instead [24]. We may assume that either microgravity conditions or collective effects alter the scaling of the energy input into our granular gases of rod-like grains.
Outlook: The velocity distribution functions and scaling laws described above require additional investigation and explanation from the point of view of individual rod collision dynamics (with the walls and among rods), whose accurate modeling can still be improved by detailed comparison of experimental and numerical singlecollision data. A second important aspect is the identification of the role of the walls for the detailed shape of the velocity distributions. An interesting challenge will be the comparison of simulations and experiments at higher packing fractions. Here we can expect more influence of particle-particle collisions relative to collisions with the walls, highly pronounced orientational inhomogeneity as well as clustering effects. Taking into account the recent advances in the analysis of experimental data in similar systems [25], one expects more possibilities to fine-tune the simulations and get more insight into the dynamics of granular gases. Further investigations can be envisioned in the direction of ensembles of particles of more elaborate shapes, as well as mixtures of different particles.