Cluster model for wave patterns of a 3D vertically vibrated granular system

Based on the cluster model for wave-like motions of a 2D vertically vibrated granular system we established previously (Chin Phys Lett, 2010, 27: 124501), a generalization of the cluster model for a 3D granular system is presented in this paper. The 3D cluster model proposes frustums of pyramids as clusters of the 3D granular system, and explains wave patterns as the result of the cluster-boundary and cluster-plate collisions. By analyzing the movement of one cluster in a collision period, we derive a basic equation, which relates the internal characteristic parameter to the external driving parameters. The theoretical results reproduce the behaviors of wave patterns as changes in the driving parameters, and the 3D dispersion relation, which agree with the experimental results.

Trapezoid clusters exist in 2D vertically vibrated granular systems where wave-like motions form ( Figure 1) [1]. The recognition of the trapezoid clusters leads us to construct a cluster model, and the comparisons between theoretical and experimental results support the 2D cluster model [2]. This leads to a consideration of whether clusters exist also in the 3D vertically vibrated granular systems where wave patterns form, and if the clusters exist, whether a 3D cluster model can be constructed. In spite of the inability to observe the movements of local areas of 3D granular systems directly, or the absence of any corroborating experimental reports, the generalization of the cluster model from 2D to 3D granular systems may be carried out using the similarities between 3D and 2D vertically vibrated granular systems.
Some wave patterns, such as squares, strips, and hexagons oscillating at f /2 or f /4, where f is the driving frequency, are observed on the surface of a 3D vertically vibrated granular system [3,4]. Although the patterns of 3D granular systems are richer than those of 2D granular systems [5], there are still some similarities between the two types of granular sys-*Corresponding author (email: xiaohui712@gmail.com) tems, e.g. the ranges of f /2 patterns are both approximately 2.5 < Γ < 4.5, where Γ is the amplitude of dimensionless driving acceleration, the ranges of kinks following the f /2 patterns are both approximately 4.5 < Γ < 5.5, and the ranges of f /4 patterns following the kinks are both approximately Γ > 5.5 [3,5]. Furthermore, experiments point out that the dispersion relations of 3D and 2D granular systems are both 1/ f 2 -dependent [6][7][8].
1 Construction of 3D cluster model 1

.1 Physical configuration
Due to the similarities between 3D and 2D granular systems described above, we assume that there are also cluster structures in the 3D vertically vibrated granular systems where wave patterns form. Note that the wave-like motions of 2D granular systems are the result of the cluster-plate and cluster-cluster collisions, and the periodicity of the wave-like motions is determined by the spatial periodicity of the arrangement of the clusters and the temporal periodicity of the cluster-plate and cluster-cluster collisions [2]. The wave pat- terns of 3D granular systems are also periodic phenomena, so the surface wave patterns may imply the internal cluster structures. One may imagine that the clusters of square patterns are frustums of square pyramids, the clusters of strip patterns are frustums of narrow rectangular pyramids, and the clusters of hexagon patterns are frustums of regular hexagonal pyramids. To understand the characteristics of the 3D clusters, we assume that the granules are identical ideal smooth balls [9], and construct the frustums of pyramids from the bottom to the top as shown in Figure 2.
During the construction of frustums of pyramids, we find that for a cluster of the square pattern ( Figure 2(a)), the particle numbers of each layer from the bottom to the top are N 2 , (N − 1) 2 , (N − 2) 2 , · · · , (N − H + 1) 2 , where N is the particle number of each side of the bottom square, and H the num-ber of layers. That is, the layers from the bottom to the top are squares with the particle numbers of each layer decreasing quadratically. For a cluster of the strip pattern ( Figure  2(b)), the particle numbers of each layer from the bottom to the top are where N 1 is the particle number of each long side of the bottom narrow rectangle, and N 2 the particle number of each short side. That is, the layers from the bottom to the top are narrow rectangles with the particle numbers of each layer also decreasing quadratically. Finally, for a cluster of the hexagon pattern ( Figure 2(c)), the particle numbers of each layer from the bottom to the top are where N is the particle number of each side of the bottom regular hexagon. In this case, every two layers from the bot- tom to the top are regular hexagons with the particle numbers of each layer decreasing quadratically.
It should be emphasized that in a 2D granular system, the alternating arrangement between positive and negative trapezoids form wave-like motions, as shown in Figure 1, but there are no negative frustums of pyramids coupling with positive ones in 3D granular systems. The particles filling in the space of an array of positive frustums of pyramids integrate as a boundary of the wave pattern (see the pink areas of Figure  2(a), (b) and (c)). So, the processes of exchanging base angles, positive frustums of pyramids converting to negative ones, wave valleys becoming wave peaks, and the reverse processes, which occur in 2D granular systems [2], may not occur in 3D granular systems. According to experiments on wave patterns, we propose that for the square and strip patterns (in these two cases, the boundaries of wave patterns do not change [4], as shown in Figure 2(a) and (b)), the physical configurations of the 3D cluster model are like those in sequence (take one positive frustum of pyramid for example): First, the positive frustum of pyramid collides with the plate of the container (cluster-plate collision), while the boundary flies freely in the air (Figure 3(a)). During the cluster-plate collision, first the bottom layer of the cluster collides with the plate, the second layer then collides with the bottom layer, the third layer then collides with the second layer, and so on, until the top layer collides with the layer below it, thus, a complete chain-type collision is accomplished. In fact, the collisions occurring in the cluster are complicated, and multiple chain-type collisions occur in sequence, which result in the exhaustion of the energy received from the plate, and the cluster moving together with the plate for some time after the cluster-plate collision (Figure 3(b)), until it leaves the plate with a ballistic motion. The processes described above are the same as those in 2D granular systems [2]. In the air, the upward cluster collides with the downward boundary (clusterboundary collision). However, the exchange of base angles between the cluster (the positive frustum of pyramid) and the boundary does not occur, possibly since the mass of all clusters is much larger than that of the boundary. This is especially the case when the size of the cluster is very large and the number of layers is very small, and the cluster-boundary collision just makes the boundary rebound without exchanging base angles between the cluster and boundary ( Figure  3(c)). Thus, there is no peak-valley alternation, the pattern remains the same, and the boundary is the previous one, except for a periodic slight change of displacement in the vertical direction. What is different for the hexagon pattern, is that there is a peak-valley alternation (Figure 2(c) and (d)), so this case is similar with those of the 2D granular systems, except for the difference between the 3D cluster-boundary and 2D cluster-cluster collisions, i.e. the cluster-boundary collisions result in a periodic change of pattern [4] (Figure 3(d)). In the following analysis, we do not discuss the cluster-boundary collision in detail, except for making a rough estimate of the influence on the movement of the cluster, while we place emphasis on the cluster-plate collision and the parametric relationship derived from it.

Characteristic parameter C H
In the 2D cluster model, we characterized a cluster using the coefficient C H [2]. Here we try to do this in the 3D case. Consider that a positive frustum of square pyramid collides with the plate. Assume that the precollisional velocities of the cluster and plate are v c and v p , respectively. Using the inelastic collision model to calculate the velocity v H of the top layer of the cluster after the first complete chain-type collision, similar to what we did in the 2D granular system, we with Figure 3 Physical configurations of the 3D cluster model (vertical sections). In each figure, the pink represents peaks, the white represents valleys, and the arrows represent vertical movements. (a) The cluster-plate collision; (b) the cluster moves together with the plate for some time; (c) the cluster-boundary collision makes the boundary rebound with no peak-valley alternation; (d) the cluster-boundary collision converts the positive frustum of pyramid to a negative one, and a peak-valley alternation occurs.
where e is the coefficient of restitution which has been assumed to be a constant for all collisions in our model [10,11], and C H,square is the coefficient, C H , of the square pattern. Here we have assumed that N H according to the experiments [4], and eq. (2) is an approximation. Likewise, the coefficient C H of the strip pattern is In this case, we have assumed that N 1 N 2 and N 1 H [4]. Lastly, the coefficient C H of the hexagon pattern is Take a system of steel spheres with a coefficient of restitution e ≈ 0.9 [10] for example, and the dependence of C H − 1 on N for different H's is shown in Figure 4.
It is obvious that for different numbers of layers and different wave patterns, all C H − 1's decrease monotonically with N, and the relative positions of wave patterns in Figure 4 are almost the same. Therefore, the coefficient C H represents the 3D cluster (or wave pattern), that is, different C H 's correspond to different wave patterns. Furthermore, it could also be said that C H represents the wavelength of the wave pattern.

Parametric relationship
We now establish the relationship between the internal characteristic parameter and the external driving parameters according to the movement of a cluster in a collision period, as in the 2D cluster model [2]. However, unlike the 2D case, we take into account the influence of the delay time of the cluster-plate collision on the velocities of layers of the cluster in the 3D cluster model. The important moments of a collision period are described in Figure 5: At the moment −t 0 , the cluster leaves the plate with a ballistic motion; at the moment t 2 , the cluster-boundary collision occurs, which results in the boundary rebounding; at the moment t 1 , the cluster-plate collision occurs, and the process lasts until the moment t 1 + τ, then the cluster moves together with the plate for some time;   Figure 5 A collision period of the f /2 wave pattern. The solid cosine curve represents the movement of the plate, and the dashed curve the movement of the cluster. Here −t 0 , t 2 , t 1 , t 1 +τ, and 2T −t 0 , where T is the driving period of the plate, are the moments at which the cluster-plate separation, the clusterboundary collision, the cluster-plate collision, the cluster-plate comovement, and the second cluster-plate separation occur, respectively. The thick curve represents the process of the cluster-plate collision.
at the moment 2T − t 0 , the cluster leaves the plate once more, and a new collision period begins.
Taking the time into account, eq. (1) is rewritten as where the additional term −gτ is the influence of the delay time τ of the cluster-plate collision on the velocity of the top layer of the cluster. The postcollisional velocity of the plate is v p (t 1 + τ). Similar to the operation performed in the 2D cluster model, let the postcollisional velocity of the top layer of the cluster equal that of the plate for the sake of the periodic movement of the cluster [2], that is We write the precollisional velocity of the cluster as where the term −gΔt is the influence of the cluster-boundary collision on the velocity of the cluster. The precollisional velocity of the plate is and the postcollisional velocity of the plate is We estimate the typical value of the delay time τ to be one driving period T , and assume that τ = T + Δτ, where Δτ is very small compared with τ, then eq. (9) can be rewritten as Inserting eqs. (7), (8) and (10) into eq. (6), and then reducing, we obtain where Γ = 4π 2 f 2 A/g is the amplitude of dimensionless driving acceleration. Eq. (11) is a basic equation of the 3D cluster model. It points out that the internal characteristic parameter C H of a cluster is dependent on the external driving parameters f and Γ, that is, both the driving frequency and driving acceleration could change wave patterns, which is consistent with experiments [3].

Patterns vs. Γ
We now analyze the behavior of wave patterns as a change in the driving acceleration Γ using eq. (11). First, we estimate the dependence of the collision phase 2π f t 1 on Γ. According to the completely inelastic ball collision model [12], the 2π f t 1 -Γ relationship can be linearized as Considering the influence of the cluster-boundary collision on the movement of the cluster, except for the additional term −gΔt in eq. (7), we assume that the cluster-boundary collision also makes the collision phase of the cluster-plate collision smaller, and then modify eq. (12) to 2π f t 1 ≈ Γ + 0.7 (2.5 < Γ < 4.5), 1.2Γ + 1.6 (4.5 < Γ < 7). (13) Inserting eq. (13) into eq. (11), letting 2π f Δτ = π/20 and 2π f Δt = 2π, we obtain the dependence of C H − 1 on Γ, as shown in Figure 6. It should be noted that modifying the constant terms of eq. (13), in acceptable ranges, is to keep Figure 4(c) and Figure  6 consistent in the following analysis.
Before analyzing Figure 6, we analyze the behavior of f /2 wave patterns as a change in Γ. Combining Figure 4(c) with the part of f /2 wave patterns of Figure 6, we get Figure 7. From Figure 7, the behavior of the f /2 wave patterns can be understood [3]: When Γ is small, the square (or strip) pattern occurs first; As Γ increases (see the arrow of Figure 7(a)), the squares (or strips) enlarge, viz. the wavelength increases (see the arrow  (or arrow ) of Figure 7(b)); When Γ is increased to about 4 (the value is smaller for strips), the squares (or strips) are too large to form a pattern on the surface of the 3D granular system, so the squares (or strips) disappear (see the horizontal dashed line of Figure 7), and the system turns to the hexagon pattern with smaller C H ; Continuing the increase in Γ, the hexagons enlarge, viz. the wavelength increases (see the arrow  of Figure 7(b)); When Γ is increased to about 4.5, the hexagons also disappear.
Returning to Figure 6: After the value of Γ exceeds 4.5, C H − 1 suddenly increases as Γ increases initially, then decreases once more. In the range of 4.5 < Γ < 6, there are no wave patterns, and a flat surface appears. Still increasing Γ, the value of C H − 1 returns to the range of Figure 7, and the f /4 wave patterns occur. The behavior of the f /4 wave patterns as a change in Γ is similar to the f /2 case [3].

Patterns vs. f
In this subsection, we analyze the behavior of wave patterns as a change in the driving frequency f using eq. (11). To enable a comparison with the experiments [3], for the f /2 wave patterns, we let Γ = 3.5, where only the f /2 square or strip pattern appears. Inserting the first case of eq. (13) into eq. (11), and assuming that Δτ ≈ 0.002 s, we obtain the dependence of C H − 1 on f for the f /2 wave patterns (The value of Δτ is determined as follows: according to the critical value of C H − 1 where only the strip pattern appears, 2π f Δτ ≈ π/5, we then take the critical frequency f ≈ 43 Hz where only the f /2 strip pattern appears [3] into it, then get Δτ ≈ 0.002 s): For the f /4 wave patterns, we let Γ = 6.5, where only the f /4 square or strip pattern appears. Inserting the second case of eq. (13) into eq. (11), and assuming that Δτ ≈ 0.00025 s, we obtain the dependence of C H − 1 on f for the f /4 wave patterns (Here 2π f Δτ ≈ π/22, we then take the critical frequency f ≈ 90 Hz where only the f /4 strip pattern appears [3] into it, then get Δτ ≈ 0.00025 s): Combining eqs. (14) and (15) with Figure 4(c), we can analyze the behavior of the wave patterns as a change in f [3], as shown in Figure 8: When f increases, C H − 1 increases, while N decreases. This means that the granular system converts wave patterns from both squares and strips to only strips (see the arrows  and  of Figure 8(b), under the horizontal dashed line, both the squares and strips can appear, while above the line only the strips appear, since the value of N is so small that the squares cannot exist in a stable manner). For the f /2 patterns, the critical frequency is about 43 Hz; for the f /4 patterns, it is about 90 Hz (see the vertical dashed lines of Figure 8(a)).

Dispersion relation
Finally, we analyze the dispersion relation of 3D granular systems [6]. Letting Γ = 3.5, in the range of 10 < f < 40 Hz, the f /2 square pattern can appear [3]. For the square pattern, the wavelength λ = Nd, where d is the particle diameter. Thus, eq. (2) can be rewritten as the dependence of C H,square − 1 on λ/d: On the other hand, to emphasize the influence of f on λ, we regard the delay time τ (= T + Δτ) as a constant. In fact, τ is f -dependent, and τ is always adjusted so that 2π f τ ≈ 2π. Here, we only consider the average effect of τ. For a granular system of seven layers of steel spheres, taking τ ≈ 0.053 s (Letting f = 20 Hz in 2π f τ = 2π + π/20, get τ ≈ 0.053 s) and the first case of eq. (13) into eq. (11), we obtain the dependence of C H,square − 1 on f : Combining eq. (16) with eq. (17), we obtain the dependence of λ/d on f , viz. the 3D dispersion relation, as shown in Figure 9. From Figure 9, we see that the theoretical dispersion relations are in agreement with the experimental ones, and the larger the particle diameter, the smaller the relative error. For the case of d = 0.4 mm, the theoretical result is almost identical to the experimental one. It is worth noting that, during the construction of the 2D cluster model, we did not take into account the influence of the delay time of the cluster-plate collision on the velocities of layers of the cluster, so the derived 2D dispersion relation is Γ-independent, but there, a correlation coefficient k was introduced to relate λ to Γ [2]. Analyzing reveals that the effect of k and the influence of Γ on λ are equivalent.

Conclusions
In summary, the 3D cluster model proposes frustums of pyramids as the clusters of 3D granular systems, and considers the wave patterns as the result of the cluster-plate and clusterboundary collisions. By analyzing the movement of a cluster in a collision period, a basic equation which relates the internal characteristic parameter to the external driving parameters is obtained. Using the basic equation, we analyze the behaviors of wave patterns as changes in the driving acceleration and driving frequency, and the 3D dispersion relation. In certain ranges of parameters, the theoretical results are in agreement with the experimental results. Of course, some aspects of the model need to be improved, such as, the influence of the cluster-boundary collision on the movement of a cluster, where we have no idea how the influence depends on parameters. With regard to the dependence of the delay time of the cluster-plate collision on parameters and its influence on the results, we are still uncertain.