Energy non-equipartition in vibrofluidized particles

TL;DR

This study uses DEM simulations to explore how particle interaction parameters affect energy distribution between translational and rotational modes in vibrofluidized granular systems.

Contribution

It reveals the impact of realistic spring stiffness ratios and friction on energy non-equipartition in vibrofluidized particles, highlighting the role of particle roughness.

Findings

Energy ratio decreases with friction for certain stiffness ratios.

Non-equipartition persists at high friction and specific stiffness ratios.

Energy distribution behavior varies with particle roughness and interaction parameters.

Abstract

The aim of the present work is to investigate the influence of the realistic model parameters for particle interactions, specifically the spring stiffness coefficient for the tangential force between particles on the energy equipartition in a vibrofluidized system. To achieve this, a three-dimensional vertically vibrated granular system consisting of spherical particles is simulated using the discrete element method (DEM) implemented in the open-source software LAMMPS. Interparticle and wall-particle interactions are determined using the linear-spring dashpot model. Simulations are performed for particles ranging from nearly perfectly smooth to nearly perfectly rough. Two different values for the ratio of the tangential to normal spring stiffness coefficient $\kappa$ ($2/7$ and $3/4$) are chosen for most of the simulations. The ratio of the translational to the rotational kinetic energy ($K$) monotonically decreases with an increase in the friction coefficient, $\mu$ for $\kappa=2/7$; however, for $\kappa = 3/4$, after an initial reduction with $\mu$, $K$ increases and plateaus at $\approx 5$, indicating the absence of equipartition of energy between the translational and rotational modes. Further simulations performed for $0.67 \le \kappa < 1$ confirm non-equipartition of energy for particles with a very high friction coefficient.

Equations

$$\Delta \left(\frac{1}{2}{v}^{\prime 2}\right)=-\frac{1}{4}(1-e_n^2)({\hat{r}}_{ij}\cdot {\vec{v}}_{ij}{)}^2$$

$$-{\eta }_2\left({\hat{r}}_{ij}\times \vec{G}\right)\cdot \left({\hat{r}}_{ij}\times \vec{G}\right)$$

$$+{\eta }_2^2{\left|{\hat{r}}_{ij}\times \vec{G}\right|}^2-{\eta }_2\left({\vec{\omega }}_{ij}\cdot \left({\hat{r}}_{ij}\times \vec{G}\right)\right)$$

$$\Delta \left(\hat{I}{\omega }^{\prime 2}\right)={\eta }_2\left({\vec{\omega }}_{ij}\cdot \left({\hat{r}}_{ij}\times \vec{G}\right)\right)$$

$$+\frac{{\eta }_2^2}{\hat{I}}\left({\hat{r}}_{ij}\times \vec{G}\right)\cdot \left({\hat{r}}_{ij}\times \vec{G}\right)$$

$$\Delta E=\Delta \left(\frac{1}{2}{v}^{\prime 2}\right)+\Delta \left(\hat{I}{\omega }^{\prime 2}\right)=-\frac{1}{4}(1-e_n^2)({\hat{r}}_{ij}\cdot {\vec{v}}_{ij}{)}^2-$$

$$\frac{\hat{I}}{1+\hat{I}}\frac{1-{\beta }^2}{4}\left({\hat{r}}_{ij}\times \vec{G}\right)\cdot \left({\hat{r}}_{ij}\times \vec{G}\right)$$

$$T_o$$

$$T_o=\frac{\sqrt{2}{U}_o^2}{\pi N{d}^2[(1-e_n^2)+\frac{\hat{I}}{1+\hat{I}}(1-{\beta }^2)+\frac{T_{R,o}}{T_o}\frac{1}{1+\hat{I}}(1-{\beta }^2)]}$$

Full text

I Introduction

The equipartition theorem states that the kinetic energy is equally distributed among all degrees of freedom in a fluid [22]; however, the equipartition of energy deviates in gases in round vessels [15], bio-molecules [5], laser-cooled atoms [2], non-spheroidal molecules [6], granular mixtures [17, 30, 8, 19], homogeneously cooling systems [1], and granular gases with rough particles [14]. A seemingly simple system of a vibro-fluidized smooth particles deviates from equipartition of energy [12] and anisotropy in the fluctuating kinetic energy $T_{x,y,z}=\frac{1}{2}\langle(u_{x,y,z}-\langle u\rangle_{x,y,z})^{2}\rangle$ is observed. The anisotropy in $T_{x,y,z}$ in a vertically vibrated granular system is due to the fact that the fluctuating kinetic energy is transferred from the bottom plate to the particles in the vertical direction. The energy is then distributed in the other two directions due to subsequent inter-particle interactions. The isotropic mean squared fluctuating kinetic energy of the particles $(T_{o})$, which is obtained equating the rate of energy input to the system due to bottom-wall particle collision and rate of dissipation due to inelastic inter-particle collisions at the leading order in a moment expansion method [11], scales as $\frac{U_{o}^{2}}{N\,d^{2}\left(1-e_{n}^{2}\right)}$, where $U_{o}$ is the wall velocity, $N$ is the number density, $d$ the particle diameter and $e_{n}$ the normal coefficient of restitution. The difference in the $T_{x},T_{y}$ and $T_{z}$ is maximum near the bottom wall and monotonically decreases along the bed height and the $T_{x,y,z}$ asymptotically approaches $T_{o}$ for $N\,d^{2}\left(1-e_{n}^{2}\right)<<1$ [11, 25]. The behaviour differs when dissipation due to air drag is considered.

In an assembly of realistic granular particles, the partitioning of energy between the translational and rotational modes depends on the surface roughness ([20] and references therein). In the limit of the nearly smooth particles, the rotational and translational kinetic energies are independently balanced, and the kinetic energy is not equally distributed between the rotational and the translational degrees of freedom. In the other limit of nearly perfectly rough particles, the partition of kinetic energy depends on the particle inelasticity and the surface roughness quantified in terms of the normal ($e_{n}$) and the rotational ($\beta$) coefficient of restitution [20]. McNamara and Luding [14] defined a ratio $R=\frac{\Delta E^{o}}{\Delta\overline{E}+\Delta E^{o}}$ to quantify the partition of energy between the translational and rotational modes. Through event-driven simulations, they showed that $R$ weakly depends on the coefficient of restitution for $\beta\approx 0$. In the energy conserving limits, i.e., $\beta\sim-1\,\,\mathrm{or}\,\,+1$, $R$ is independent of $e_{n}$. Grasselli et al [9] studied the anisotropy in a 2-D vibro-fluidized granular bed in micro-gravity. They determined the rotational coefficient of restitution, the ratio of the rotational to the translational kinetic energy ($R_{T}$), and observed that the anisotropy in the fluctuating kinetic energy depends on the area fraction. Castillo et. al [3] discussed the departure from equipartition in the context of a granular system in a magnetically levitated bed.

Although the non-equipartition of kinetic energy in granular systems is widely observed, Nichol and Daniels [16] reported nearly equipartition of energy between the translational and rotational modes for a dense bimodal mixture subject to bidirectional periodic excitation on an air table. In their experiment, the granular assembly is excited by electromagnetic bumpers on three walls, with the fourth wall serving as a constraint.

Potiguar [18] performed numerical simulations for the experimental set-up discussed in [16]. They used a linear spring dashpot model [4] to determine the contact force between the colliding disk-shaped particles, with the spring stiffness constant, $k_{n}=5\times 10^{4}\frac{mg}{d}$, and $\gamma$, the dissipation coefficient, as two parameters. The tangential force is determined from the sliding friction coefficient ($\mu$). Simulations were performed for a wide range of $\gamma$, resulting in $0.2<e_{n}<0.9$ and for $\mu=0.5$. They observed that the ratio of translational to rotational kinetic energy depends on the number density, coefficient of restitution, frequency, and magnitude of the energy injected.

The literature indicates that the non-equipartition of energy in a granular system depends on several factors, including particle properties, dissipation mechanisms, and particle number density. The objective of the present work is to systematically investigate the effect of the friction coefficient and the ratio of the tangential to the normal stiffness coefficients on the partition of fluctuating kinetic energy between the translational and rotational modes of vibro-fluidized particles using the Discrete Element Method (DEM) and analyze the results in the framework proposed in [14, 20]. To that end, simulations are performed using the open-source software LAMMPS with large values of spring stiffness constant, $k_{n}>10^{7}\frac{mg}{d}$ [21, 27] to ensure binary collisions.

II Methodology

II.1 Background theory

In the simplest hard-sphere model the collisions are characterized by two parameters: the normal ($e_{n}=-\frac{v_{n}^{\prime}}{v_{n}}$) and the rotational coefficients of restitution $(\beta=-\frac{\vec{v}_{s}^{\prime}\cdot\hat{k}_{s}}{\vec{v}_{s}\cdot{k}_{s}})$, where, $v_{n}$ is the component of the relative velocity of particles along the line joining the centres of particles, and $\vec{v}_{s}=\vec{v}_{ij}-\vec{v}_{n}+\frac{d}{2}\left(\hat{r}_{ij}\times\vec{\omega}_{ij}\right)$ is the slip velocity of the point of contact. The primed quantities represent the post-collision properties. The particle indices are $i$ and $j$, $\hat{r}_{ij}$ is the unit vector drawn from particle centre of $i$ to the centre of $j$, the subscript $n$ refers to the normal direction, $\vec{v}_{ij}$ is the relative velocity of the particle $i$ with respect to $j$, and $\vec{\omega}_{ij}$ is the sum of angular velocities about their respective centers, $\hat{k}_{s}$ is the unit vector in the slip velocity direction, defined as $\hat{k}_{s}=\frac{\vec{v}_{s}}{|\vec{v}_{s}|}$ [28].

The changes in the translational, rotational, and total kinetic energy, during a collision, are expressed by Equations  II.1 – II.1 [20]:

[TABLE]

[TABLE]

[TABLE]

where $\Delta\left(\frac{1}{2}v^{\prime 2}\right)$ = $\frac{1}{2}(v_{i}^{2}-v_{i}^{\prime 2})+\frac{1}{2}(v_{j}^{2}-v_{j}^{\prime 2})$, $\eta_{2}=\frac{1}{2}(1+\beta)\hat{I}/(1+\hat{I})$, $\vec{G}=\vec{v}_{ij}+\vec{r}_{ij}\times(\vec{\omega}_{i}+\vec{\omega}_{j})$, $\hat{I}=4I/md^{2}$, and $I$ is the moment of inertia. Mass($m$) and diameter ($d$) of the particles are used as scaling parameters. The translational and rotational velocities are normalized with $U_{o}$ and $2U_{o}/d$, where $U_{o}=2\pi Af$ is the maximum velocity of the vibrating base. The kinetic energy is conserved for a perfectly elastic ($e_{n}=1$) collision between two perfectly smooth ($\beta=-1$) particles. For perfectly elastic particles with rough surfaces, the dissipation of energy is solely due to friction between particles and depends only on $\beta$. Otherwise, the change in the kinetic energy is a function of $e_{n}$ and $\beta$. The term $E_{\text{exch}}=\eta_{2}\left|\vec{\omega}_{ij}\cdot\left(\hat{r}_{ij}\times\vec{G}\right)\right|$ accounts for the gain in the rotational energy compensating for the loss in translational energy and is a function of $\beta$. The ratio of the transfer of energy from the translation to the rotational mode to the energy dissipation $\Theta=\frac{E_{\text{exch}}}{\Delta E}$, in general, depends on both $e_{n}$ and $\beta$. For perfectly elastic particles, $\Theta$ is an explicit function of the rotational coefficient of restitution ($\beta$). $\beta$ depends on the friction coefficient ($\mu$), the impact angle ($\gamma$), and the normal coefficient of restitution ($e_{n}$). The collision is said to be sliding if $-1\leq\beta\leq 0$. in this regime, $\beta=-1+\frac{7}{2}\mu(1+e_{n})\cot{\gamma}$ [29]. In the stick-slip regime when $0<\beta\leq 1$, $\beta$ is a complex function of the material properties [10]. Luding and McNamara[14] showed the dependence of the distribution of mean fluctuating kinetic energy in the rotational and the translational model on $\beta$ using event-driven simulations. They have shown that the distribution is independent of the normal coefficient of restitution.

II.2 Simulation Method

Figure 1 shows the schematic representation of the computational domain. The domain is periodic in the gravity normal direction. The upper wall is placed at a height of ${1000d}$ to mimic a semi-infinite domain. The bottom wall vibrates sinusoidally with the maximum energy of ${U_{o}^{2}=4\pi^{2}A^{2}f^{2}}$, where $f$ and $A$ are the frequency and amplitude of the vibration, respectively. The base frequency is maintained constant at $f^{*}=2.26\sqrt{d_{p}/g}$.

The amplitude is varied between ${0.3d\leq A\leq 0.7d}$, resulting in the non-dimensional acceleration ($\Gamma=4\pi^{2}Af^{2}/g$) in the range $60\leq\Gamma\leq 140$ [7]. Simulations are performed for $0\leq\mu\leq 10$, and a wide range of $\epsilon=Nd^{2}(1-e_{n}^{2})$, where $N$ is the number density of particles (number per unit base-area)

The linear spring dashpot (LSD) model is used to determine the contact force (Equations presented in Sec I of [26]). The normal spring stiffness constant ($k_{n}=10^{8}\frac{mg}{d}$) is selected to ensure that $\frac{t_{c}}{t_{f}}\ll 1$, where $t_{c}$ is the contact time and $t_{f}$ is the average time between two successive collisions [24, 21, 27]. Simulations are performed for two different values of $\kappa=\frac{k_{t}}{k_{n}}$, where, $\kappa=$ $\frac{2}{7},\frac{3}{4}$ [24, 27]. The viscous dissipation coefficient $\gamma_{n}$ is set such that $e_{n}$ varies between 0.85 - 1 [23]. A wide range of friction coefficient $\mu$ ($0\leq\mu\leq 10$) is used in the present study. A very large value of $\mu$ is included in the simulation to mimic a nearly perfectly rough case [20].

Simulations are performed in two stages. First, the simulation is performed for $7\times 10^{7}t_{c}$, where $t_{c}=3.1\times 10^{-4}\sqrt{\frac{d}{g}}$, is the time spent at contact. The time step of $\Delta t=\frac{t_{c}}{10}$ is used at this stage. Once the fluidized particles’ total kinetic energy reached a steady state, simulations were run with a smaller time step, $\Delta t=t_{c}/100$. Instantaneous linear and angular velocities of particles obtained from the DEM simulations are further analyzed to determine the profiles of the mean fluctuating translational $\left(T_{(x,y,z)}=\frac{1}{2}\langle(v_{(x,y,z)}-\langle v_{(x,y,z)}\rangle)^{2}\rangle\right)$ and rotational $\left(T_{R,(x,y,z)}=\frac{\hat{I}}{2}\langle(\Omega_{(x,y,z)}-\langle\Omega_{(x,y,z)}\rangle)^{2}\right)$ kinetic energy. In a vibro-fluidized bed, there is no mean linear and angular velocity, and the simulation results confirm that both $\langle v_{(x,y,z)}\rangle\rangle$and $\langle\Omega_{(x,y,z)}\rangle)$ are $\approx 0$. The linear and angular velocities are normalized as $\vec{v}=\vec{V}/{U_{o}}$ and $\vec{\Omega}=d\vec{\omega}/2U_{o}$. Here, $U_{o}=2\pi Af$ is the maximum velocity of the vibrating base. The particle’s vertical position is measured relative to the vibrating base.

The bed height averaged mean translational and rotational fluctuating kinetic energy ($\mathrm{KE_{T}}=\frac{1}{2\,h}\int_{0}^{h}\left\langle\left\|\vec{v}\right\|^{2}\right\rangle dz$, $\mathrm{KE_{R}}=\frac{\hat{I}}{2\,h}\int_{0}^{h}\left\langle\left\|\vec{\Omega}\right\|^{2}\right\rangle dz$) are also determined. More than $10^{3}$ configurations are used to determine the ensemble averages represented within $\langle\,\rangle$.

III Results

III.1 Profiles of the mean fluctuating kinetic energy

Figures  2(a) and  2(b) show the vertical profiles for $T_{(x,y,z)}$ and $T_{R,(x,y,z)}$ for $Nd^{2}(1-e_{n}^{2})=0.39$, $e_{n}=0.95$, $\kappa=2/7$, $\mu=0.01$ and a very high value of friction coefficient $(\mu=10)$ respectively. Figure  2(a) also plots the $T_{(x,y,z)}$ for perfectly smooth sphere ($\mu=0$). $T_{z}$ monotonically decreases with height and reaches the equilibrium temperature ($T_{o}$) beyond $\frac{z}{d}=30$. The velocity fluctuations are anisotropic near the base, $T_{z}>T_{x},T_{y}$, consistent with earlier studies [11]. However, equipartition of translational fluctuating kinetic energy among the three degrees of freedom was achieved at $\frac{z}{d}\approx 10$ for the smooth particles. The temperature $T_{(x,y,z)}$ for nearly perfectly smooth particles ($\mu=0.01$) shows a very similar trend as shown by the perfectly smooth particles with $\mu=0$. There is equipartition for the translational energy at $\frac{z}{d}\approx 25$. Profiles for $T_{R,(x,y,z)}$ show a monotonic decrease to an equilibrium value, with height. Also, $T_{R,(x,y,z)}$ is $\mathcal{O}(10^{-2})$ smaller than $T_{(x,y,z)}$ in this case.

For slightly rough particles ($\mu=0.01$, and $\beta\approx-1$), the translational and rotational energy conservation equations are nearly independently satisfied as $\eta_{2}\ll 1$ (Eqs  II.1 and II.1).

On the other hand, for nearly perfectly rough particles ($\beta\approx 1$), $\eta_{2}\rightarrow\frac{\hat{I}}{1+\hat{I}}$, the maximum possible value, which leads to a strong coupling between the conservation of the translational and rotational fluctuating kinetic energy during a contact. Figure  2(b) shows that $T_{(x,y,z)}$ and $T_{R,(x,y,z)}$ have equal values within numerical accuracy for $\frac{z}{d}>10$. Analysis of the collisions with $\kappa=\frac{2}{7}$ for particles with $\mu>0.5$ confirmed that the rotational coefficient of restitution, $\beta\approx+1$. Figure  2(c) shows that there is a significant difference in the equilibrium values of $T$ and $T_{R}$, this is because for $\kappa=\frac{3}{4}$, $\beta<+1$ [27]. The value of $\beta$ depends on the nature of the particle contact defined in [13]. The fraction of the particle-contacts following the gross sliding determined using the method outlined in [27], for different $\kappa$ is presented in [26]. Figure 2(d) shows the solid volume fraction in log-linear scale. For all cases, the profiles of the solid volume fraction show a maximum followed by an exponential decay with a slope $\frac{-gz}{T_{o}}$, where $T_{o}$ is the equilibrium temperature. The volume fraction profile of particles with very high surface roughness $\mu=10$ with $\kappa=\frac{3}{4}$ shows a steeper decay than the other three cases, which is consistent with the results in Figure  2(c).

Figure 3 plot $\frac{T_{(x,y,z)}}{T_{R,(x,y,z)}}$ and $\frac{T(z)}{T_{R}(z)}$, where $T_{(R)}(z)=\frac{1}{3}\left(T_{x(R)}(z)+T_{y(R)}(z)+T_{z(R)}(z)\right)$. For $\kappa=\frac{2}{7}$, $\mu=0.01$, the ratio of the $T$ and $T_{R}$ at equilibrium $\frac{T_{o}}{T_{Ro}}=\mathcal{O}(10^{2})$ and for $\mu=10$ the ratio is close to unity. Figure 3(c) shows that $\frac{T_{o}}{T_{Ro}}$ plateaus at $\approx 5$ for $\kappa=\frac{3}{4}$ and $\mu=10$.

The equilibrium “granular temperature” $(T_{o})$ is obtained by equating the rate of energy input due to the collision with the vibrating wall and the rate of dissipation due to inelastic and frictional contacts truncated at the leading order. For perfectly smooth and nearly perfectly rough particles, $T_{o}$ can be expressed by Eqs 4 and 5 respectively [11, 20]:

[TABLE]

For nearly perfectly smooth particles in a vibro-fluidized bed where there is no mean angular velocity ($\langle\Omega\rangle=0$), $T_{o}$ and $T_{R,o}$ are distinct, the expression for the $T_{o}$ depends on $T_{R,o}$:

[TABLE]

Figure 4(a) and 4(b) plot $\frac{U_{o}^{2}}{T_{o}}$ vs $Nd^{2}(1-e_{n}^{2})$ for simulations with perfectly smooth particles and $Nd^{2}\left[(1-e_{n}^{2})+(1-\beta^{2})\right]$ for vibro-fluidized particles with high surface roughness, respectively. Results are in excellent agreement with the theory (Eq.  4) for $Nd^{2}(1-e_{n}^{2})<1$ [11]. Eq 5 is a direct extension of the leading order theory for rough particles in the dilute vibro-fluidized bed [11, 20]. An excellent agreement between the simulation results and theory for $Nd^{2}\left[(1-e_{n}^{2})+(1-\beta^{2})\right]<1$ is observed. Figure  4(c) plots $\frac{U_{o}^{2}}{T_{o}}$ vs. $Nd^{2}\left[(1-e_{n}^{2})+\frac{\hat{I}}{1+\hat{I}}(1-\beta^{2})\right]$ for particles with $\mu<0.1$. Eq.  6 can be rewritten neglecting the last term on the right-hand side, since $\frac{T_{R,o}}{T_{o}}\ll 1$ for these cases. It is worth noting that a simple scaling for $T_{R,o}$ cannot be obtained in this regime. Due to the lack of a unique scaling for $T_{o}$ for the wide range of roughness, the results are plotted against the dimensional parameter $U_{o}^{2}$ in Figure  5(a). The height averaged fluctuating kinetic energy, $\mathrm{KE_{T}}$ and $\mathrm{KE_{R}}$ are determined for $0\leq\mu\leq 10$, three different values of $e_{n}$, and $\kappa$ and the ratio $k=\frac{\mathrm{KE_{T}}}{\mathrm{KE_{R}}}$ is examined in Section III.2. The effect of the frictional losses is discussed in Section  III.3.

III.2 Bed height averaged mean fluctuating kinetic energy

The ratio, $K=\frac{\mathrm{KE_{T}}}{\mathrm{KE_{R}}}$ is plotted against $U_{o}^{2}$ for $Nd^{2}=4$, $e_{n}$ = $0.85,0.9,0.95$ and $\mu$ ranging from 0.01 to 10 in Figure 5(a)). The plots in the panel suggest that $K$ is independent of base velocity and the normal coefficient of restitution and depends only on the friction coefficient. As the friction coefficient increases, $K$ approaches unity. Figure  5(b) shows the variation of $K$ with $\mu$ for two different number densities, $Nd^{2}=1,4$, and coefficient of restitution $e_{n}=0.95$. The ratio of the mean fluctuating rotational to the translational kinetic energy, $K$ is found to be independent of the number density. The same final $K$ was obtained across simulations with four different initial configurations, each with a distinct value of $K_{o}=K(t=0)$. This initial condition independence is shown in figure 5(c). Figure  5(d) shows $K$ vs $\mu$ for perfectly smooth particles. In this case, the dissipation is purely due to friction. $K$ decreases monotonically with $\mu$ and plateaus at unity as $\mu$ assumes a very high value, for $\kappa=\frac{2}{7}$. In contrast, for $\kappa=\frac{3}{4}$, the behaviour is non-monotonic; $K$ starts to deviate from the value for $\kappa=\frac{2}{7}$ for $\mu\geq 0.1$. To understand the cause of the deviation, the DEM simulation data are analyzed and presented within the hard-sphere framework. As a first step, the data is processed (a) to identify the contact and determine the effective value of $\beta$ and (b) to determine the terms in Equations  II.1 and II.1.

III.3 Energy balance during contact

The instantaneous positions of the particles are analyzed in a manner similar to that described in [27, 26]. Once the contacts are identified and the binary nature of the collisions is verified, the parameter $\beta=-\frac{v_{s}^{\prime}}{v_{s}}$ is determined from the pre- and the post-collision velocities. The energy change $\Delta E$ (Eq. II.1) and $E_{\mathrm{exch}}(=\eta_{2}\left|\vec{\omega}_{ij}\cdot\left(\hat{r}_{ij}\times\vec{G}\right)\right|$ are determined for each contact detected. The median (Q2) of the distribution of $\Delta E$ and $E_{\mathrm{exch}}$ are plotted as a function of $\mu$ in Figure  6. Data is collected over $10^{3}$ configurations from the simulations performed with $\Delta t=\frac{t_{c}}{100}$.

$\Delta E$ is non-monotonic for $\kappa=\frac{2}{7}$. Dissipation during a collision due to friction increases with $\mu$ up to $\mu=0.1$, and reduces thereafter. The sliding and sticking regimes are mutually exclusive for $\kappa=\frac{2}{7}$. In the sticking regime, $\beta\approx 1$. In this limit, the collisions are energy conserving (Eq  II.1. With $\mu$, the fraction of contact in the sticking regime increases, resulting in more energy-conserving contacts. Nearly 75$\%$ of the contact is in the sticking regime for $\mu=1$ (Figure presented in [26]). In case of $\kappa=\frac{3}{4}$, $\beta<1$, and the fraction of contact in the stick-slip regime plateaus at 55$\%$ beyond $\mu=1$. $E_{\mathrm{exch}}$ increases monotonically with $\mu$ before it plateau for $\kappa=\frac{2}{7}$ and reduces for $\kappa=\frac{3}{4}$. The ratio $\Theta=\frac{E_{\mathrm{exch}}}{\Delta E}$ increases sharply with $\mu$ for $\kappa=\frac{2}{7}$, explaining the equipartitioning of fluctuating energy between the translational and the rotational modes at high $\mu$. In contrast, for $\kappa=\frac{3}{4}$, $\Theta<1$ for all values of $\mu$ and reduces to $\sim 0.03$ for very rough particles ($\mu>1$). The equipartition is not observed for $\kappa=\frac{3}{4}$ because dissipation exceeds the energy exchange between the translation and rotational modes.

IV Conclusion

An assembly of rough, inelastic spherical particles subject to vertical vibration was simulated using the open-source code LAMMPS. The linear spring-dashpot model is used to determine the normal and tangential forces between particles in contact. The normal spring constant is selected such that the collisions are predominantly binary. Two values of $\kappa(=\frac{k_{t}}{k_{n}})$ are selected. The time-period of the normal and the tangential contacts are equal for $\kappa=\frac{2}{7}$ and two mutually exclusive regimes of contacts are obtained. The physical interpretation of the stiffness constant as the inverse of the compliance leads to $0.67\leq\kappa<1$ [27]. The rotational coefficient of restitution obtained from simulations with $0.67\leq\kappa<1$ is qualitatively similar to the experimentally observed results ([27] and references therein). Results obtained for $\kappa=\frac{3}{4}$ are presented here for detailed discussion. Simulations were also conducted for $\kappa=0.67$ and $0.95$, and results confirmed that $K=5.1\pm 0.05,4.91\pm 0.04$ respectively.

The observations from the simulations are:

The equipartition of the mean-squared fluctuating kinetic energy is observed in simulations with $\kappa=\frac{2}{7}$ and for particles with very high friction coefficients. For this range of parameters, $>75\%$ of contacts fall into the energy-conserving sticking regime. 2. 2.

For $\mu\geq 0.3$ and $0.67\leq\kappa<1$, the ratio of translation to rotational fluctuating kinetic energy $K$ increases beyond 1.5, and the equipartition of energy is not observed. This is because the stick-slip collisions are not energy-conserving.

Non-equipartition of energy between different degrees of freedom is relevant for granular rheology. The results presented here also suggest that the selection of $\kappa$ may be crucial in predicting macroscopic flow behaviour of realistic particles.

Acknowledgements.

We acknowledge the Indian Institute of Technology Bombay for the licensed version of Grammarly. The software was used to check the manuscript’s English grammar. VK was supported by funding from the MHRD and the Science and Engineering Research Board, Government of India (Grant no. SR/S2/JCB-31/2006).