## 5 Simulation Verification  
In this section, the analytical results are illustrated and supported by means of computer simulations. Consider system (9) and apply a Runge-Kutta routine with adaptive stepsize control [44] to it with the following parameter values: \(m_{1}=6\) kg, \(m_{R}=4\) kg, \(r_{R}=0.005\) m, \(l_{0}=0.3\) m. Note that for the parameter values considered in Sections 4 and 5, the synchronization criterion and stability criterion of synchronous states are satisfied and therefore synchronous motions are expected to occur. The analysis starts by assuming identical friction coefficients, i.e., it is assumed that \(f_{R1}=f_{R2}=f_{R3}\). In a first simulation, the friction coefficients are set to \(f_{R1}=0.08,\ f_{R2}=0.08\), \(f_{R3}=0.08\). Note that an interference of \(\pi/4\) angular displacement is applied to CR 1 at \(t=11\) s. The obtained results are shown in Figure 11. At the initial stage, the induction motor of an exciter starts working and its angular velocity rapidly reaches the operating value. Simultaneously, three CRs begin to accelerate and oscillate wildly in their angular velocities due to the influence, exerted via the coupling structure, of the exciter. As time goes on, with the aid of the dynamic coupling torques, the exciter and three CRs operate synchronously and steadily at \(t=2.6\) s. At the steady-state, the rotational velocity of the synchronization for the exciter is 991.3 r/min. As for the three CRs, under the influence of gravity, their rotational velocities of synchronous operation tend to rise and fall equally and periodically around 991.3 r/min. As a result of the regular variations on the synchronous rotational velocities of the three CRs, the stable value of \(2\alpha_{1}\) changes regularly and periodically around a constant phase difference of \(204^{\circ}\) (referred to as the partial synchronization between the exciter and the CR 1). As can be seen in Figure 11c,d, for identical friction coefficients the three CRs are observed to synchronize in-phase (phase difference of \(0^{\circ}\)), which was predicted in the numerical analyses. Furthermore, from the results presented in Figure 11e,f, it is clear that the steady-state vibrations in \(x\)- and \(y\)-directions are equal to each other, whereas the stable value in \(y\)-direction is zero, which indicates a circular motion of the coupling structure. It is to be noted that the movement in circle is an effective and optimal way to transmit energy from an exciter to three CRs through the rigid platform. Figures 11g, 12g, 13g and 14g show the time series of the reaction force on three CRs exerted by the cavities' surfaces. It is worth noting that even the minimum value of \(N_{i}\) at the steady-state reaches more than 100 N. That is to say, the CRs in the four simulation cases can cling to the cavities' surfaces quite tenaciously when rotating around the cavities' centroids synchronously and steadily. Note that the strong clinging ability of the CRs is essential for crushing the materials. In a second simulation, the friction coefficients are changed to \(f_{R1}=0.08\), \(f_{R2}=0.04\), \(f_{R3}=0.08\). In this simulation, the interference of \(\pi/4\) phase is applied to CR 1 at \(t=11\) s. The obtained results are shown in Figure 12. Since the three friction coefficients are non-identical, the angular accelerations of the three CRs are different during the starting process (Figures 12a and 13a). When time reaches approximately 4 s, the vibrating system comes into the synchronous and steady state. At this moment, the rotational velocity for the exciter and phase difference between the exciter and CR 1 of the synchronous operation are stabilised at 991.7 r/min and around a constant value of \(205^{\circ}\) (partial synchronization), on account of the CR's non-negligible gravity effect.