## 5 Simulation Verification
The reason is that an increased difference value between two CRs' friction coefficients can lead to a bigger difference in the synchronization abilities between an exciter and a CR, which can cause bigger fluctuations in the stable phase differences among three CRs, as shown in Figures 4c and 5c. Simulation studies (see Figures 11-13) have revealed that the three coupled CRs can show synchronized behavior even when there are unavoidably small differences between their friction coefficients. The next set of simulation results looks at the vibrating system with \(f_{R1}=0.04\), \(f_{R2}=0.04\), \(f_{R3}=0.04\). As can be seen in Figure 14, the stable and synchronous states become settled at \(t=2.5\) s and have a synchronous rotation velocity of \(992.6\) r/min (stabilizing at or oscillating around it), and phase differences of \(2\alpha_{1}=184^{\circ}\sim 197^{\circ}\) (average value \(190.5^{\circ}\)), \(2\alpha_{2}=0^{\circ}\), and \(2\alpha_{3}=0^{\circ}\). Similar to the previous simulation analysis of identical friction coefficients, the vibrating system exhibits partial synchronization in an exciter and a CR, and in-phase synchronization in three CRs. By comparing the results of the four simulation studies already conducted, it is observed that the greater the sum of the three friction coefficients, the smaller the synchronous rotation velocity. Additionally, the stable values of the phase differences from the simulation results are fairly consistent with those shown in Figure 8.  
For a better clarity, we conduct quantitative comparisons between numerical results and simulation results with changes in \(f_{Ri}\) (\(i=1,2,3\)). Table 1 lists the average phase differences of an exciter and three CRs in steady state of numerical analyses and computer simulations. By comparison, the stable phase differences obtained from the numerical results are in accord with the stable phase differences obtained from the simulation results. Therefore, the simulation results can verify the theoretical analysis and the numerical results above. for \(\eta_{r}\) and \(\eta_{m}\), which are fixed to \(\eta_{r}=0.18\), \(\eta_{m}=1.1\). The results of this simulation can be seen in Figure 15. Instead of reaching synchronous operation among an exciter and three CRs as in the previous results, at the steady-state, no definite relationship exists between an exciter and a CR and they move away from a synchronized state, whereas the in-phase swing motion is observed in the three CRs. In this case, the simulation results can illustrate the effectiveness of the theoretical and numerical analysis. Besides this case, the three identical friction coefficients are changed to three non-identical ones and the computer simulations concerned have been carried out. In these simulations, however, even the three CRs are not able to perform the partial-synchronous swing, not to mention the synchronized motion of an exciter and a CR. Figure 16 shows the time evolution of the dynamic coupling torques in the first four simulations. From the results presented in Figure 16, it should be noted that in whatever simulation, the total dynamic coupling torque \(T_{\text{DC}}\) is kept at zero throughout the entire process. This finding is consistent with a previous discussion in coupling dynamic characteristics. That is to say, the dynamic coupling torques serve as a medium of energy exchange and transmission among an exciter and three CRs; they never generate or consume energy in the vibrating system. At the synchronous and stable operational stage of the system, the dynamic coupling torque \(T_{\text{DC1}}\) arrives at a positive range and the dynamic coupling torques \(T_{\text{DC2}}\), \(T_{\text{DC3}}\), and \(T_{\text{DC4}}\) fall into the negative range. Hence, \(T_{\text{DC1}}\) acts as the coupled load torque of the exciter and \(T_{\text{DC2}}\), \(T_{\text{DC3}}\), and \(T_{\text{DC4}}\) play the part of the coupled driving torques for the three CRs (see Equations (49)-(52)).