Coarsening dynamics of granular heaplets in tapped granular layers

A semi-continuum model is introduced to study the dynamics of the formation of granular heaplets in tapped granular layers. By taking into account the energy dissipation of collisions and screening effects due to avalanches, this model is able to reproduce qualitatively the pattern of these heaplets. Our simulations show that the granular heaplets are characterized by an effective surface tension which depends on the magnitude of the tapping intensity. Also, we observe that there is a coarsening effect in that the average size of the heaplets, V, grows as the number of taps k increases. The growth law at intermediate times can be fitted by a scaling function V~kz but the range of validity of the power law is limited by size effects. The growth exponent z appears to diverge as the tapping intensity is increased.


81.2
heaplets merge as observed in experiments [3]. The model is then studied numerically and the results indicate some of the relationships between tapping intensity and the effective surface tension of granular layers. Finally analysis of the results shows that the average size of the heaplets V grows as a power law with the number of taps k, V ∼ k z for a limited range of k values. The exponent z appears to diverge as the tapping intensity is increased.

Model
Consider a thin layer of N 1 granular particles spread out over a flat plate. The plate is then tapped repeatedly at a low pace with constant shock intensity. In one complete tapping cycle, there are two different processes, each of which requires a different description. In the tapping phase, the granular layer is perturbed by an external shock and granular particles acquire kinetic energy to allow them to explore the phase space and hop on to neighbouring sites. The hopping range of the granular particles depends on the amount of kinetic energy received by each of the individual particles. The amount of energy supplied to the system is controlled by the dimensionless acceleration Γ = U/g ∆t, where U is the velocity of the plate when it is in motion and ∆t is the time interval over which it is in motion. Note that U is also the vertical velocity gained by an elastic particle sitting on the plate when it is tapped. Our particles are not, however, elastic. Therefore, the vertical velocity gain is αU , where 0 ≤ α ≤ 1 is the departure coefficient analogous to the coefficient of restitution and it characterizes the degree of dissipation of the system. After the tap, the particles undergo a ballistic flight and fall back again onto the static plate where they relax, subject to avalanche dynamics before the next tap. In this relaxation phase, the granular particles may stay immobile or move about, depending on the local slope. If the local slope is less than a critical slope then the profile remains stationary. Conversely, if the local slope is greater than the critical slope, matter moves down the slope collectively as an avalanche, until the slope is less than the critical slope.
Let the area density n(x, k) be the number of granular particles above a unit area of the plate after k taps. Here x = (x, y), where x and y are orthogonal coordinates parallel to the plate. The average diameter D of the granular particles is chosen to be 1, so that n(x, k) is now equal to the local thickness of the granular layer, so long as there is no compactification throughout the tapping process. Of course this is just an idealization: the packing fraction of the granular material will be changed [9] if the granular layer is thick or if it is intentionally prepared in a low-density state, but such cases will not be considered here.
We can approximate the acceleration A p of the plate's movement due to tapping as a sum of M delta functions: where the velocity U gives the intensity of the shocks which are assumed uniform over the whole plate. τ specifies the time interval between taps. Here τ must be greater than any time scale in the relaxation process. This is to ensure that next tap only occurs after the system is fully relaxed. In the dilute limit, a single inelastic grain sitting on the plate experiences a net force

81.3
where α is the departure coefficient mentioned earlier. During a small interval [t k − ∆t/2, t k + ∆t/2], where t k is shorthand for t k = (k − 1)τ , the time of the kth tap and the maximum velocity transferred to the particle is In order for the particle to hop, ∆v z must be positive, or αU > v 0 ≡ g ∆t so that external acceleration overcomes the gravitational force. In what follows this is always the case.
In the continuum limit one has to make two crucial modifications. α must be replaced by an effective departure coefficient of the granular bulk. α is expected to be a monotonic decreasing function of n(x, k). This is because, when the number density is high, inter-grain collisions occur more often and more energy is dissipated, so that grains depart with a smaller departure velocity. While the precise form of α(n) is unknown, in this paper it is taken to be for reasons of simplicity and to make a comparison with Duran's model [3]. Here,n is the average density of the system and α 0 is the effective departure coefficient for average densityn. Also, the second term in equation (2) needs to be changed. In the continuum limit, according to Duran [4] there is an effect dependent on the position of a grain in the heap. If the grain is not at the top of the heap it supports a fraction of the weight of the grains above it; hence its effective mass is increased. As a result, in order for a granular particle sitting on the inclined side of the granular layer to hop, it requires a larger velocity kick than a particle sitting on the flat region of the layer. This screening effect can be incorporated into equation (2) by replacing m in the second term with an effective mass Here D is the average grain diameter and it is set to 1 henceforth. n T is the altitude of the nearest peak in the corrugated granular layer. p is a parameter of unknown value which is set equal to 5 in reference [4] and this value appears to give a match to experimental results. θ c is the angle of repose and is close to the value of π/6 [4]. After these modifications, equation (3) can now be written as where A = α 0n U v 0 is the tapping intensity. There are close similarities between our model and Duran's [3]. In both models, the resulting pattern formation is due to competition between the upward hopping motion of particles and the downward screening effect of avalanches. However, the mechanism causing upward motion is different in the two models. Duran [3] and Shinbrot [10] conjectured that the hopping of granular particles is due to the upcoming air flux through the porous bed of the granular layer. The velocity of the air flux at height n can be approximated by Darcy's law v air = Kp/n, where K is the permeability of the granular bed, p is the pressure differential across the granular bed. The velocity of the air flux v air must be greater than the free fall velocity of a granular particle v f in order to eject this particle. However, in our model we adopt a simpler picture where kinetic 81.4 energy is transferred to the granular particles via direct collisions between plate and particle and between particle and particle. Energy is dissipated via the effective departure coefficient. The form of the effective departure coefficient is chosen as in equation (4) so that the velocity transfer from plate to particle has the same 1/n dependency (n ≥ 1) as v air in Darcy's law. Obviously, we have the freedom to choose other functional forms of α(n): the simple choice here is just to achieve comparison with Duran's model. As observed from the simulation, however, it turned out that the exact form of α(n) is unimportant so long as it is a monotonic decreasing function.
Now it remains to specify the hopping process in the tapping phase. When the granular layer is tapped, at an unstable site where ∆v z > 0, particles can hop to any neighbouring site within a circle defined by a maximum horizontal hopping range, R max (∆v). Here we assume on the unstable site that all the granular particles are re-distributed. Each grain hops with equal probability in any direction and with a random hopping range, R < R max . We can estimate the maximum hopping range R max (∆v) by the following simple arguments. The maximum velocity received by a particle during one tap is ∆v z vertical to the plate. In this case R max (∆v) is zero, for the motion is strictly vertical. It is unlikely that particles will always hop in the vertical direction: an inter-particle collision can change the hopping direction to any angle. Assume that a particle is ejected with a speed ∆v z cos φ with φ denoting the angle from the vertical axis of the plate. Then, the horizontal hopping range is where t is the flight time of the particle in vacuum. R is maximum when φ = π/6, so that the estimated value for R max is given by The estimated R max is very crude and does not take into account the details of the air-particle and particle-particle interactions except phenomenologically. However, one can consider ξ in R max as an adjustable parameter, which varies from one experiment to another experiment depending on the roughness of the material used in the experiment and the humidity of the surrounding environment.
There are several ways to describe the avalanche process during the relaxation phase. Here we use a slope-dependent diffusion equation. We start from the equation of continuity ∂ t n = −∇ · J with the constitutive current J given by where 1/ √ β = tan θ c is the critical slope and η sets the diffusion rate. The current is chosen in such a way that, when the gradient is greater than the critical slope, mass current flows downhill to smooth out density fluctuations, but no mass is transferred if the gradient is less than the critical slope.

Simulation
The model is studied on a N × N square lattice with a periodic boundary condition. Initially, the granular layer is prepared by assigning a random height between [0, 5] to each site and letting the granular layer relax. At the beginning, the fluctuation in n(x, 0) is considerably smaller than the height of the heaplets formed later. On each tapping cycle, the velocity of particles at each site is calculated and the density number is updated according to rules defined by equations (6) and (8). After this, the equation of continuity and the corresponding constitutive current equation (10) are solved and iterated numerically until there is no more mass transferred at each site.
The following density plots show typical results from the simulation. Figure 1 shows the density plots of n(x, k) at different values of k, the number of taps, for two set of parameters. The top panel shows extended patterns of ridges and corresponds to a smaller tapping intensity (A = 6.0), while the lower panel shows localized heaplets and corresponds to a greater tapping intensity (A = 16.8). Both simulations coarsen when k increases, eventually reaching their final stage where the entire system is a single granular heap.
One can see from figure 1 that, as the tapping intensity is increased, the patterns favour isolated circular heaplets which suggests that an effective surface tension can be defined which increases with the tapping intensity. Since there are no attractive forces between granular particles, the granular layer cannot have a true surface tension. This 'surface tension' effect is entirely due to the system trying to maximize the local energy dissipation by decreasing the value of the effective departure coefficient so that particles are less mobile. In order to maximize the local energy losses, particles have to be as close to each other as possible, for this will increase the number of inelastic collisions between particles. Therefore, with the constraint that the local slope should not exceed the critical slope, the global attractor of the pattern is a single large heap, where local energy loss is the largest. Akiyama et al [12], Clément [11] and Duran [3] have discussed the effects of convection on the heaping process of granular layers. Duran [3] in particular suggests that the 'surface tension' is brought about by the convective forces dragging 81.6 particles from the surroundings into the granular heaplets. However, in our simulation it is due to the maximization of energy losses, since in our model the granular layer does not contain any information about the internal interactions of the sand heap, and therefore convective forces are not taken into account.
Nevertheless, the simulation result does show the Laplace-Young pressure effect suggested by Duran [3]. A typical example is shown in figure 2. (The parameters here are A = 16.0 and ξ = 0.2207, but similar results are found over a wide range of values of these parameters.) In the centre of the figure there are two connected granular heaplets. After several taps the smaller heaplet merges with the larger heaplet due to matter moving along the connecting neck.
We are also interested in the dynamics of the size of the heaplets. The relevant length scale measure l is the average thickness of the site h i weighted by the local volume, V i , and the average volume of the heaplet is V = l 3 . Since we know that h i = n i and the local volume of the site is just V i = n i σ, where σ is the area element of the lattice site, it follows that There is some limited evidence for power law behaviour V ∼ k z in the intermediate region of the log-log plots of V versus k in figure 3. There are three regions in each plot. The early region is where the granular layer is randomly distributed and the number density fluctuation is not large enough to trigger coalescence. The late-time region is where all the heaplets have joined into one big heap, leaving some remaining individual particles hopping randomly about and occasionally encountering the large single heap. The intermediate region is also the fastgrowing region where heaplets grow and merge into each other. The range of the intermediate region grows and remains linear as the system size increases, as can be seen from the two curves with A = 17.6. This shows that system size limits the range of the intermediate power law region. For very large system sizes we would expect on this basis that the intermediate power law region would have a greater range and remain linear in the log-log plot. One can easily estimate the increase in the average volume of the heaplets V by the following simple argument. As the saturated value of V corresponds to a state where a single heap contains approximately all the grains in the system, it should be proportional tonL 2 . Here the average number of grains per siten is a conserved quantity. Then, an increment in V is given by ∆ log 10 V = 2 ∆ log 10 L. In figure 3 the broken curve corresponds to a lattice size L 1 × L 1 = 80 × 80 and the marked broken curve to L 2 ×L 2 = 120×120. The vertical shift in the graphs is expected to be ∆ log 10 V = 0.352 which may be compared to 0.360 from the simulation. All curves correspond to a lattice size L 1 × L 1 = 80 × 80 and are averaged over ten different runs except for the marked broken curve which corresponds to L 2 × L 2 = 120 × 120 and is averaged over five runs. Each curve is saturated after a long time. This can be clearly seen from the A = 18.0 full curve. The marked curve corresponds to the larger lattice size and shows a wider range of the linear region and a larger saturation value of V .
The exponent z in the power law depends on the tapping intensity A as is shown in figure 4. As we can see z is roughly constant, z ∼ 1.0, for low tapping intensity. When the tapping intensity is increased A > 14, z suddenly increases and appears to diverge near A = 18. For these large values of the tapping intensity, we observed that the system appears to be in a flat configuration for some time and suddenly a single heap formed within a few taps. For still greater tapping intensities, the heaplet patterns disappear and the layer remains flat with a small random variation superimposed. This pattern-disorder transition occurs at a not very precisely defined critical tapping intensity of A ∼ 18. (This critical intensity appears to depend slightly on the initial conditions, which is the reason for our statement that it is not precisely defined.) Above this critical tapping intensity, no patterns are observed. This is due to the fact that the perturbation is so strong that no site is stable, and each site constantly undergoes re-distribution on each tapping cycle. Although the power law parametrization is based on limited evidence, it does provide a useful parameter z to distinguish the pattern-forming region (A < 18) from the non-pattern-forming region (A > 18).  < 14), z is roughly constant (z ∼ 1.0), but appears to diverge for A > 14. Each value of z is averaged over ten different runs.

Conclusion
We have introduced a simple model of the heaping process in a tapped granular layer. The model is capable of reproducing the essential morphologies of tapped granular systems. Qualitatively the effective surface tension of the granular heaps is closely related to the tapping intensity, and it is shown that there is no need for convective forces for this effective surface tension to exist. The length scale of the system coarsens according to the power law l ∼ k z in a limited range. The exponent z appears to diverge as tapping intensity is increased and provides a useful parameter for distinguishing the pattern-forming and the non-pattern-forming regions.