Macroscopically equivalent granular systems with different numbers of particles

Granular flows often show a remarkable similarity with those of molecular fluids [1, 2]. The success of granular hydrodynamic theories in predicting many complex granular behaviours indicates that such a relation is not only superficial [3–9]. But despite continued development, the defining properties of granular materials, such as the dissipative nature of the particles' interactions, still present a challenge for continuum theories, specially for high packing densities or strong dissipations [10–12]. An additional fundamental difficulty stems from the enormous difference in the total number of constituents between granular and molecular systems; while in molecular media the microscopic relevant length-scale is orders of magnitude smaller than the macroscopic one, in granular media macroscopic fields may vary in distances of the order of a few particle diameters. This possibly big influence of a few particles on macroscopic quantities implies the existence of inherently large fluctuations, which can drastically modify the global dynamics, especially near transitions [13–15]. Deepening our understanding of the role played by these fluctuations is thus of fundamental importance for the development of a successful continuum description of granular media.

In the following paper we first analyse the influence that finite-number-driven fluctuations have on the macroscopic states of flowing granular matter. For this, we study the possibility of constructing macroscopically identical granular systems with significantly different number of particles. Starting from a given macroscopic hydrodynamic state, and using physical arguments and expressions for the transport coefficients from granular hydrodynamic models, we derive the dependency on particle size of all terms of the conservation laws. As in systems with the same macroscopic density the particle size is directly related to the total number of particles, we essentially see the dependency of the macroscopic states on the total number of particles present in the system. We demonstrate that in general the granular hydrodynamic equations are not particle-size invariant. Nevertheless, we show that, by properly scaling the restitution coefficient, the limit of vanishing particle size becomes well defined and leads to invariant conservation equations. Although superficially our scaling analysis is similar to earlier studies of the granular hydrodynamic equations, our purpose is a very different one, namely the investigation of whether it is possible to construct a hydrodynamic formulation that, given the macroscopic fields, remains invariant when changing the size of the particles, and therewith with the total particle number.

Subsequently, to study the relevance of the total number of particles on the strength of fluctuations in flowing granular matter we focus on the granular Leidenfrost state [7, 16, 17]. Using hard-sphere simulations and the derived scalings we create systems with up to four orders of magnitude difference in the particle number density with identical average hydrodynamic fields. Although in average the systems are identical, the amplitude of the oscillations previously observed in the granular Leidenfrost state are seen to decrease with the number of particles while, surprisingly, the frequency of these oscillations remains invariant.

In the following we study the particle-size dependency of the three-dimensional granular hydrodynamic equations. Our goal is to recreate equivalent hydrodynamic states with significantly different total numbers of particles N. Two hydrodynamic states are considered identical if the time-averaged invariant macroscopic hydrodynamic fields are the same in all space. As invariant hydrodynamic fields, functions of the spatial coordinates ${\boldsymbol{r}}=(x,y,z)$ and time t, we consider the packing fraction, $\phi ({\boldsymbol{r}},t)$, the velocity ${\boldsymbol{u}}({\boldsymbol{r}},t)$, and the fluctuations in velocity ${T}^{* }({\boldsymbol{r}},t)$. All these are obtained from the position and velocities of the particles, ${{\boldsymbol{r}}}_{i}$ and ${{\boldsymbol{v}}}_{i}$, through a coarse-graining procedure. The packing fraction is given by $\phi ({\boldsymbol{r}},t)={mn}({\boldsymbol{r}},t)/{\rho }_{p}$, with m, d and ${\rho }_{p}=3\,m/4\pi {\left(d/2\right)}^{3}$ the mass, diameter and density of the particles. The number density field is a convolution of the particles' positions, $n({\boldsymbol{r}},t)={\sum }_{i\,=\,1}^{N}{ \mathcal W }({\boldsymbol{r}}-{{\boldsymbol{r}}}_{i}(t))$, where ${ \mathcal W }$ is a smoothing or coarse-graining function, i.e., a continuous and integrable replacement of the Dirac delta function which allows us to obtain well defined macroscopic quantities [18, 19]. The velocity field is then ${\boldsymbol{u}}({\boldsymbol{r}},t)={\sum }_{i}{{\boldsymbol{v}}}_{i}(t){ \mathcal W }({\boldsymbol{r}}-{{\boldsymbol{r}}}_{i}(t))/n({\boldsymbol{r}},t)$. Finally the fluctuations in velocity ${T}^{* }({\boldsymbol{r}},t)\equiv {k}_{B}T({\boldsymbol{r}},t)/m=\tfrac{1}{3}{\sum }_{i}{\left({{\boldsymbol{v}}}_{i}(t)-{\boldsymbol{u}}({\boldsymbol{r}},t)\right)}^{2}{ \mathcal W }({\boldsymbol{r}}-{{\boldsymbol{r}}}_{i})/n({\boldsymbol{r}},t)$, with $T({\boldsymbol{r}},t)$ the granular temperature field, and k_B the granular equivalent of the Boltzmann constant⁵ . Strong invariance would imply the same packing fraction, velocity and granular temperature distributions in space and time for systems with different number of particles and initial conditions. Naturally, this condition is too strong for what are essentially fluctuating systems, and we therefore focus on invariance of the fields averaged over a timescale τ larger than fluctuations driven by the microscopic scale, in which case we say that the systems are equivalent.

In terms of these macroscopic fields the granular hydrodynamic equations can be written as:

$$\begin{eqnarray}&&{D}_{t}\phi =-\phi {\rm{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rho }_{p}\phi {D}_{t}{\boldsymbol{u}}=-{\rm{\nabla }}\cdot (p{\boldsymbol{I}})+{\rm{\nabla }}\cdot (\eta \hat{\sigma })+{\rm{\nabla }}(\gamma ({\rm{\nabla }}\cdot {\boldsymbol{u}}){\boldsymbol{I}})+{\rho }_{p}\phi {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{3}{2}{\rho }_{p}\phi {D}_{t}{T}^{* } & = & -p({\rm{\nabla }}\cdot {\boldsymbol{u}})+\gamma {\left({\rm{\nabla }}\cdot {\boldsymbol{u}}\right)}^{2}+\eta (\hat{\sigma }:{\rm{\nabla }}{\boldsymbol{u}})\\ & & +{\rm{\nabla }}\cdot (\kappa {\rm{\nabla }}{T}^{* })+{\rm{\nabla }}\cdot (\mu {\rm{\nabla }}\phi )-I.\end{array}\end{eqnarray} \tag{ 3 }$$

where ${D}_{t}={\partial }_{t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}$ is the material derivative, ${\boldsymbol{I}}$ the identity matrix and the deviatoric strain rate $\hat{\sigma }={\rm{\nabla }}{\boldsymbol{u}}+{\rm{\nabla }}{{\boldsymbol{u}}}^{T}-\tfrac{2}{3}({\rm{\nabla }}\cdot {\boldsymbol{u}}){\boldsymbol{I}}$. These correspond to the mass conservation equation; the momentum conservation equation, with p the pressure, $\eta$ the shear viscosity, γ the second viscosity and ${\boldsymbol{g}}$ the acceleration of gravity; and the granular temperature equation, i.e. energy balance, with κ the coefficient of thermal conductivity (note that Fourier's heat law has been used), $\mu$ an additional transport coefficient present in granular media, and $-I$ the sink of energy-density.

In order for equations (1)–(3) to be d-invariant all terms of any equation should scale equally with d. In what follows we determine these scaling relations, starting with heuristic arguments and, after that, using a specific closure for the transport coefficients and the sink term available in the literature.

First, note that the continuity equation is manifestly invariant, since it only contains hydrodynamic fields and its derivatives, which are by definition identical in equivalent systems. Similarly, the left-hand terms of (2) and (3), as well as the gravitational term in (2), are also d-invariant, if we assume that the particles are made of the same material, i.e., that their density ${\rho }_{p}$ is independent of their size. In principle, one could let ${\rho }_{p}$ scale in any arbitrary way with the particle size, but this would lead to the same conclusions, only complicating the algebra. Additionally, the pressure is proportional to the ideal gas law ($p={\rho }_{p}\phi {T}^{* }$) multiplied by a function of the packing fraction only, and therefore the terms in (2) and (3) involving the pressure are also d-invariant.

The dependency of the transport coefficients on the size of the particles can be deduced from physical arguments. Transport of mass, momentum or energy from one fluid element into a neighbouring one happens in a layer that has a thickness scaling with the mean free path, ${l}_{\mathrm{fp}}$. Since the packing fraction ϕ is invariant, ${l}_{\mathrm{fp}}$ should be proportional to the particle diameter d. Therefore all terms in (2) and (3) involving the transport coefficients $\eta$, γ, κ and $\mu$ are expected to scale as d¹.

In contrast, the dissipation term I can be expressed as the number density squared (scaling as ${d}^{-6}$) times the energy loss per collision (∼d³) times the cross-sectional area (∼d²), leading to $I\sim {d}^{-1}$. Physically this stands to reason: when the particle size decreases, the growth of the number of collisions is faster than the shrinking of the typical energy loss per collision, and therefore the dissipation per unit volume increases.

Using square brackets, ${\left[x\right]}_{d}$, to denote the scaling of a quantity x with d, the above discussion is summarized as

$$\begin{eqnarray}&&{\left[p\right]}_{d}={d}^{0},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\left[\eta \right]}_{d}={\left[\gamma \right]}_{d}={\left[\kappa \right]}_{d}={\left[\mu \right]}_{d}={d}^{1},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\left[I\right]}_{d}={d}^{-1}.\end{eqnarray} \tag{ 6 }$$

The immediate conclusion is that the general granular hydrodynamic model (1)–(3) is not d-invariant, since it contain terms that scale differently with particle size. In fact, when the particle size decreases, all transport terms decrease, whereas the energy sink term increases. In general, it is therefore not possible to create identical hydrodynamic states with (significantly) different numbers of particles.

Even if we consider a steady state (${\partial }_{t}=0$) without macroscopic flow (${\boldsymbol{u}}=0$), in which case equations (1)–(3) reduce to

$$\begin{eqnarray}&&{\rm{\nabla }}p={\rho }_{p}\phi {\boldsymbol{g}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\kappa {\rm{\nabla }}{T}^{* })+{\rm{\nabla }}\cdot (\mu {\rm{\nabla }}\phi )=I,\end{eqnarray} \tag{ 8 }$$

the conservation of momentum (7) is d-invariant, while conservation of energy (8) still isn't. However, for small d (with respect to the typical scale of variations of the fields, L) invariance can be obtained to order ${ \mathcal O }({d}^{2})$. Both I and $\mu$ are proportional to the inelasticity $\varepsilon \equiv (1-{r}^{2})$, where r is the coefficient of normal restitution. If we now let r depend on the particle diameter d such that

$$\begin{eqnarray}&&{\left[\varepsilon \right]}_{d}={d}^{\alpha },\end{eqnarray} \tag{ 9 }$$

for α > 1 the $d\to 0$ limit becomes well defined. In what follows we choose α = 2, not only for being the next natural number that leads to convergence, but also to make the dissipation I scale as the heat-flux κ, which will become useful when studying the equivalent vibrated systems. This also implies ${\left[\mu \right]}_{d}={d}^{3}$, as will be demonstrated further below. Therefore, as the number of particle increases, the influence of the granular transport coefficient diminishes rapidly (as d³) rendering the set of equations (7)–(8) independent of d to second order.

When the condition (9) is used in the general flow case then, for small d, equations (1)–(3) result in

$$\begin{eqnarray}&&{D}_{t}\phi =-\phi {\rm{\nabla }}\cdot {\boldsymbol{u}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rho }_{p}\phi {D}_{t}{\boldsymbol{u}}=-{\rm{\nabla }}p+{\rho }_{p}\phi {\boldsymbol{g}}+{ \mathcal O }(d),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\rho }_{p}\phi {D}_{t}{T}^{* }=-p({\rm{\nabla }}\cdot {\boldsymbol{u}})+{ \mathcal O }(d).\end{eqnarray} \tag{ 12 }$$

Notice that in the limit $d\to 0$, the above equations converge to those of a perfect (non-dissipative) fluid. Nevertheless, it is important to remark that solutions to these equations are not expected to coincide with solutions of equations (1)–(3) in the $d\to 0$ limit (i.e. the order of the limits is crucial). For example, it is clear that the absence of dissipation in equations (10)–(12) will fail to capture the steady state of any driven granular system, as energy would continue to increase indefinitely. Convergence at d = 0 can thus be expected only for non-driven system, where the steady state corresponds to all material being motionless.

We will now explicitly derive the dependency of the transport coefficients on d considering the expressions obtained by Garzó and Dufty, who used the Chapman-Enskog method to solve the kinetc equation of the Revised Enskog Theory [20, 21]. When expanded using Sonine polynomials, the transport coefficients χ take the general form $\chi ={\chi }_{0}\tilde{\chi }(\phi ,\varepsilon )$, where χ₀ are the values for the low-density and elastic limit [21]. The corrections for excluded-volume and finite dissipation $\tilde{\chi }(\phi ,\varepsilon )$ are quite involved; here we are just interested in the fact that they depend only on the packing fraction ϕ and the inelasticity ε.

The low-density and elastic limit expressions of all transport coefficients can be expressed as functions of ${\eta }_{0}$, as ${\kappa }_{0}=15{\eta }_{0}/4$, ${\mu }_{0}=15{T}^{* }{\eta }_{0}/4\phi$ and ${\gamma }_{0}={\eta }_{0}$ [21]. Their explicit form is given by ${\eta }_{0}=(5/96){\rho }_{p}d\sqrt{\pi {T}^{* }}$. The energy-density dissipation term follows a similar expression,

$$\begin{eqnarray}&&I=\displaystyle \frac{144}{5\sqrt{\pi }}\displaystyle \frac{{\rho }_{p}}{d}{\phi }^{2}{T}^{* 3/2}\tilde{I}(\phi ,{T}^{* },d,\varepsilon ),\end{eqnarray} \tag{ 13 }$$

To analyse the influence of finite inelasticity, the functions $\tilde{\chi }(\phi ,\varepsilon )$ can be expanded in terms of ε, which yields

$$\begin{eqnarray}&&\tilde{\eta }={A}_{0}^{\eta }(\phi )+{A}_{1}^{\eta }(\phi )\varepsilon +O({\varepsilon }^{2})\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\tilde{\gamma }={A}_{0}^{\gamma }(\phi )+{A}_{1}^{\gamma }(\phi )\varepsilon +O({\varepsilon }^{2})\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\tilde{\kappa }={A}_{0}^{\kappa }(\phi )+{A}_{1}^{\kappa }(\phi )\varepsilon +O({\varepsilon }^{2}),\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\tilde{\mu }={A}_{1}^{\mu }(\phi )\varepsilon +O({\varepsilon }^{2}),\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\tilde{I}={A}_{0}^{I}(\phi ,{T}^{* },d)\varepsilon +O({\varepsilon }^{2}).\end{eqnarray} \tag{ 18 }$$

with the A coefficients depending only on the shown quantities. The dependency on d of the dissipation constant has the form ${A}_{0}^{I}(\phi ,{T}^{* },d)={A}_{0,0}^{I}+{A}_{0,1}^{I}(\phi ,{T}^{* })d$. We then see that for small dissipations, that is, to leading order in eqs.(14)–(18), ${\left[\eta \right]}_{d}$, ${\left[\gamma \right]}_{d}$ and ${\left[\kappa \right]}_{d}$ remain unchanged, while ${\left[\mu \right]}_{d}=d{\left[\varepsilon \right]}_{d}$ and ${\left[I\right]}_{d}={d}^{-1}{\left[\varepsilon \right]}_{d}+{\left[\varepsilon \right]}_{d}$. If then the dissipation is scaled as equation (9), I does no longer diverge when $d\to 0$, as ${\left[I\right]}_{d}=d+O({d}^{2})$, and $\mu$ decreases faster than the other transport coefficients, ${\left[\mu \right]}_{d}={d}^{3}$. Therefore, we see that the derived scalings agree with our physical argumentation based on the mean free path, which leads to (5)–(6).

In summary, we have derived how the transport coefficients scale explicitly on d. We have seen that a general hydrodynamic flow is not invariant on d, as is to be expected. Nevertheless, invariance can be obtained up to order d³ in steady states by assuming a particular scaling of the restitution coefficient, such that $\varepsilon \propto {d}^{2}$. We now consider a specific granular system and generate equivalent average steady states, which allows us to study the driving mechanism behind an observed collective oscillating behavior.

We consider the granular Leidenfrost state, which consists of a density inverted particle arrangement where a high temperature, gaseous region near a vibrating bottom wall sustains a denser, colder bed of grains on top [16, 17]. As the granular Leidenfrost state corresponds to a steady state with no flow, it is expected to be described by equations (7)–(8) (with appropriate excluded volume corrections) [7, 16]. If the energy input is increased, the bed of grains goes through a transition to a buoyancy-driven convective state [15, 22]. Here we avoid this transition by considering a quasi-one-dimensional geometry with base dimensions ${l}_{x}={l}_{y}=l\ll h$, with h the height of the granular column, thus preventing through geometrical constraints the development of convection [23].

4.1. Boundary conditions

In order for the macroscopic system to be completely defined only the set of boundary conditions rests to be determined. At steady states the injected energy must be equal to the total dissipated energy. To account for the energy injection through particle-bottom wall collisions, equation (8) is integrated over the whole domain, resulting in

$$\begin{eqnarray}&&{J}_{\mathrm{in}}={\int }_{0}^{h}I{\rm{d}}z,\end{eqnarray} \tag{ 19 }$$

where ${J}_{\mathrm{in}}$ is the energy-density flux injected to the system through particle collisions with the bottom boundary. We have considered that at the free top boundary $J(z)\to 0$ as $z\to \infty$. In order to simplify our system we consider periodic boundary conditions in the lateral directions, and thus only the bottom wall injects energy. ${J}_{\mathrm{in}}$ can be analytically determined in the dilute limit for a sinusoidal driving of amplitude a and angular frequency ω, when ${u}_{b}\equiv a\omega \ll \sqrt{T(z=0)}$ [24, 25], a condition expected to be valid in the Leidenfrost state, leading to

$$\begin{eqnarray}&&{J}_{\mathrm{in}}={\rho }_{p}\phi (z=0)\,{u}_{b}^{3}{\left(\displaystyle \frac{2{T}^{* }(z=0)}{{u}_{b}^{2}}\right)}^{1/2}.\end{eqnarray} \tag{ 20 }$$

Therefore, using our previously derived scalings, and imposing that ${\left[h\right]}_{d}={d}^{0}$, it is straightforward to see that, by (19), ${\left[{u}_{b}^{2}\right]}_{d}={d}^{1}$, and thus ${\left[a\omega \right]}_{d}={d}^{1/2}$. In our case we use ${\left[a\right]}_{d}={d}^{1}$, to make sure that the amplitude remains smaller than the particle size in the limit $d\to 0$. It then follows that

$$\begin{eqnarray}&&{\left[\omega \right]}_{d}={d}^{-1/2}.\end{eqnarray} \tag{ 21 }$$

The balance of energy (19) also indicates that the relevant number of particles must be considered per unit base area, $N/{l}^{2}$, as in equilibrium the particles in a vertical section dissipate the energy injected by the vibrations of the column's base area. As $N\equiv {\int }_{{\rm{\Omega }}}n({\boldsymbol{r}})\,{\rm{d}}{\boldsymbol{r}}$, with Ω the domain, from the definition and d-invariance of ϕ we see that

$$\begin{eqnarray}&&{\left[N/{l}^{2}\right]}_{d}={d}^{-3}.\end{eqnarray} \tag{ 22 }$$

In order to reduce the total number of particles in the system and thus decrease the computation time, the dimensions of the container are scaled such that ${\left[{l}^{2}\right]}_{d}={d}^{2}$.

4.2. Dimensionless quantities

Numerical simulations are performed using the event-driven discrete particle method [26]. For details about the algorithm's characteristics we refer the reader to [23]. As in the theoretical model, we assume collisions are determined by a single coefficient of restitution r. Different systems are referred to as ${{\bf{S}}}_{d}\equiv \left\{d;{N}_{d},{l}_{d},{r}_{d},{\omega }_{d},{a}_{d}\right\}$, with the variables' subscripts denoting the specific d. In order to produce equivalent systems a reference one must be defined, which we take to be ${{\bf{S}}}_{1}$, and use d₁ = 1 as length-scale for all systems. Time is measured with respect to the acceleration of gravity, in units of $\sqrt{{d}_{1}/g}$.

Following the theoretical analysis, we scale the amplitude proportional to d, $a=0.1d$, such that ${a}_{1}=0.1$. The small prefactor is chosen to minimize the spatial inhomogeneities induced by higher oscillation amplitudes, and thus approach the limit of an effective fixed temperature boundary condition [27]. The frequency of oscillation is taken such that the system is well within the Leidenfrost state, ${\omega }_{1}=30\sqrt{g/{d}_{1}}$. Finally, two reference particle-particle coefficients of restitution will be considered, r₁ = 0.9 and ${r}_{1}^{e}=0.99$, referred to as dissipative and quasi-elastic systems, respectively (the superscript denotes quantities for the quasi-elastic case). As we want to compare systems with similar packing fractions, in the quasi-elastic case we consider a larger number of particles, so that F₁ = 12 and ${F}_{1}^{e}=32$.

5.1. Results

Time-averaged macroscopic fields are seen to converge as $d\to 0$. A snapshot of the system configuration for three different particle sizes is shown in figure 1, where the significant difference in particle numbers can be directly recognized. Vertical profiles of the time-averaged packing fraction $\hat{\phi }(z)\equiv {\left\langle \phi (z)\right\rangle }_{\tau }$ and the squared velocity fluctuations ${\hat{T}}^{* }(z)\equiv {\left\langle {T}^{* }(z)\right\rangle }_{\tau }$ are shown in figure 2 for several different ${{\bf{S}}}_{d}$. The characteristics of the Leidenfrost state can be readily recognized: low density, high temperature regions near the bottom, below high density, low temperature regions higher up [22]. As expected from (7)–(8), only for small d the conserved fields converge, although convergence is fast enough to allow us to fabricate equivalent systems with a difference of more than four orders of magnitude in N/l². The gaseous region (close to the bottom boundary) presents the most significant differences, although the maximum of $\hat{\phi }(z)$ also decreases slightly as $d\to 0$. Variations in the gaseous region are also significant in ${\hat{T}}^{* }(z)$, which presents a threefold increase accompanied by a rising total temperature ${\hat{T}}_{T}^{* }\equiv \int \hat{T}(z){\rm{d}}z$ as $d\to 0$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The shapes of $\hat{\phi }(z)$ and ${\hat{T}}^{* }(z)$ suggest that an additional source of disagreement comes from finite-size effects. The free-volume near the bottom wall, of course not taken into account in (7)–(8), is proportional to d, leading to significant differences for large d (see figure 2(a)). Secondly, variations of the convergent macroscopic fields can become comparable to d for large particles; notice for example that in the convergent $\hat{\phi }(z)$ (in figure 2(a)), the height of the gaseous region ${h}_{g}\approx 10$, which corresponds to only ${h}_{g}=5{d}_{2}$. In other words, the scalings fail when there is no clear separation of scales.

Beyond finite-size effects, the value of the coefficient of restitution is expected to have a significant influence, as we have taken $\varepsilon \propto {d}^{2}$, indirectly affecting all transport coefficients. Indeed, quasi-elastic systems show a much higher agreement as d is varied, as shown in figure 2 by ${\hat{\phi }}^{e}(z)$ and ${\hat{T}}^{* e}(z)$. The most significant differences are again observed near z = 0, further suggesting that these are finite-size, boundary-layer effects.

The sources of disagreement can be traced by measuring each term of equations (7)–(8) separately. For the energy-density dissipation rate, equation (8) states that I scales inversely linear with d, but simulation measurements show a considerable deviation in the range of d studied, as shown in figure 3(a). In dissipative systems variations are more than ∼50%, while in quasi-elastic systems it is only ∼15%. Nevertheless, a convergent behaviour is observed, and deviations between N/l² ≈ 800 and N/l² ≈ 50000 are already less than 3% in both cases. The general improvement for higher r suggests that the differences stem from neglecting higher order ε dependencies of the transport coefficients.

Zoom In Zoom Out Reset image size

Deviations from the expected d-invariant behaviour of T* are even stronger, especially in the dissipative case, as shown in figure 3(b) via ${\hat{T}}_{T}^{* }\equiv \int {\hat{T}}^{* }(z){dz}$. But again a convergent behaviour is observed as $d\to 0$, and the quasi-elastic case shows a significant improvement. Interestingly, the dependency of ${\hat{T}}_{T}^{* }$ with d changes between dissipative and quasi-elastic systems: in the former case it increase monotonically as $d\to 0$, while in the latter there is a slight overall decrease as d diminishes. The observed convergence to a macroscopic state, even when the individual scaling relations are seen to deviate, may be explained by the convergence of the general equations (1)–(3) to the perfect fluid equations.

5.2. Low-frequency oscillations

Shaken beds of grains in density inverted states undergo collective semi-periodic oscillations, referred to as low-frequency oscillations (LFOs) [15, 23, 28, 29]. These are clearly identifiable in all simulations, making it possible to study their properties in equivalent macroscopic systems with different N. Their characteristic amplitude, quantified by the standard deviation of the evolution of the vertical centre of mass, ${a}_{{LFO}}\equiv \sigma ({z}_{{cm}}(t))$, is seen to be proportional to ${N}^{-1/2}$ (${d}^{1/2}$), as shown in figure 4(a). On the other hand, the characteristic frequency ${\omega }_{{LFO}}$, obtained from the fast Fourier transform of z_cm(t), converges to a finite value as $N\to \infty$ ($d\to 0$), as shown in figure 4(b). The ${\left[{\omega }_{{LFO}}\right]}_{d}={d}^{0}$ behaviour for large N is in accordance with a previously derived theoretical expression [23], where disregarding higher order effects, the characteristic frequency was found to be given by ${\omega }_{{LFO}}^{t}={\left(g{\rho }_{g}/{m}_{s}\right)}^{1/2}$ with ρ_g the density of the gaseous region and m_s the total mass of the solid region. ${\omega }_{{LFO}}^{t}$ is thus expected to become d-invariant, as $\hat{\phi }(z)$ converges for $d\to 0$ and both ρ_g and m_s are macroscopic quantities determined by $\hat{\phi }(z)$.

Zoom In Zoom Out Reset image size

From the decrease of ${a}_{{LFO}}$ we can conclude that, in the limit of $d\to 0$, LFOs would be unmeasurable, making it an essentially finite-number (granular) phenomena. Moreover, the $1/\sqrt{N}$ law suggests that low-frequency oscillations are driven by intrinsic fluctuations due to the low number of particles in the system. We propose the following interpretation: for small N, the momentum fluctuations given by particles of the gaseous phase hitting the solid/fluid phase is comparable to the overall momentum of the solid phase, and as such the amplitude of the oscillations are large. On the contrary, for larger N, a significant amount of particles in the gaseous phase would have to transfer momentum to the solid phase at the same time to have an equivalent impact, a situation that becomes increasingly improbable as N increases.

It is interesting to notice that even though the amplitude of LFOs becomes negligible as N increases, the mode is still present in the macroscopic system, as the ${\left[{\omega }_{{LFO}}\right]}_{d}\approx {d}^{0}$ behaviour shows. This further suggests that LFOs are an intrinsic characteristic of density inverted states, as argued in [23, 29]. The situation is curious, as the mode is a macroscopic phenomenon, but its amplitude is driven by microscopic effects. That is, the timescale to obtain invariant macroscopic fields, $\tau \propto 1/{\omega }_{{LFO}}$, is independent of the length scale of averaging d. Furthermore, the evolution of z_cm(t) is seen to become less chaotic and closer to a harmonic oscillation with a clearly defined frequency as $d\to 0$, as increasingly steep peaks in the Fourier transforms indicate (not shown) [29]; this is another sign that low-frequency oscillations are driven by finite-number fluctuations.

We have studied the possibility of creating macroscopically equivalent granular systems in the same volume with significantly different numbers of particles N, by varying their size d. Considering the granular hydrodynamic equations, we have demonstrated that it is not possible to obtain equivalent systems in the most general flow case, as different terms scale differently with the particle diameter d. Nevertheless, after appropriately scaling the restitution coefficient, the limit $d\to 0$ becomes well defined and leads to a d-invariant set of conservation laws corresponding to those of a perfect fluid. As a consequence of the proper dissipation scaling, the steady-state, fluxless equations become d-invariant in the low-dissipation limit, and for small particle sizes, up to ${ \mathcal O }({d}^{3})$.

Simulations of perfect hard-spheres allowed us to test the derived scalings for a considerable range of total numbers of particles. As a test case we considered the granular Leidenfrost state, which was seen to converge to a limit time-averaged macroscopic state as $d\to 0$, with the convergence considerably faster for lower energy dissipation. Furthermore, the collective oscillatory behavior (LFO) present in the granular Leidenfrost state was deduced to be driven by the statistical fluctuations in systems with lower numbers of particles. This follows from the decrease of the amplitude of the oscillations with d for macroscopically equivalent systems. Moreover, the frequency of the LFOs remains approximately constant in the range of d studied, in accordance with previous results predicting a dependency only on macroscopic quantities.

As a final comment, we would like to remark that the same framework could be used for the study of other non-equilibrium steady states in granular matter. Macroscopic convergence can be expected for different N, opening the possibility of studying macroscopically equivalent particle systems with significantly different numbers of particles, which could, for example, help in the up-scaling of simulations to regions relevant in industries.

NR and SL acknowledge support from the NWO-STW VICI grant number 10828. NR also acknowledges support from the Millennium Nucleus 'Physics of active matter' of the Millennium Scientific Initiative of the Ministry of Economy, Development and Tourism (Chile).