Synchronization of an Exciter and Three Cylindrical Rollers with Different Dry Friction via Dynamic Coupling

Abstract

In this paper, the dynamics of a mechanical exciter and three cylindrical rollers (CRs) with the non-identical friction coefficients interacting through a rigid platform is considered. Sufficient conditions for the existence and stability of synchronous solutions in the coupled system are derived by using the average method of modified small parameters and Routh-Hurwitz principle. The obtained theoretical results are illustrated and analysed based on numerical calculations. In the analysis, the numerical results are presented for simple one-parameter variation, as well as for a group of varied parameters, when the influence of the coupling structure’s parameters on synchronization and stability is studied. An appropriate selection of the key parameters will eventually lead to desired synchronization performance. Finally, the theoretical and numerical results are supported by computer simulations. The stable synchronized states can be observed in the simulations even when there are unavoidably small differences in the three friction coefficients. If we mismatch the friction coefficients of the CRs, they are seen to synchronize with a constant phase difference. The key feature of the proposed coupled system is the dynamic coupling torque, which serves as the vehicle for transferring energy from an induction motor to three CRs without the direct driving sources and the synchronization controller for maintaining the originally synchronous and stable states against the disturbance in the simulations.

1. Introduction

Synchronization is the coordination of events to operate a system in unison, which occurs in coupled dynamical systems. The synchronization phenomena emerge when two or more biological, physical, chemical, or social systems are allowed to interact with each other and regulate their individual rhythms to a uniform rhythm through a common medium named coupling. Interesting examples of the synchronization phenomena are omnipresent: it can be observed in pedestrians walking through a bridge [1,2], in a pack of mighty animals hunting for their prey in a cooperative manner [3]. The fascinating phenomena can also be found in fireflies of Southeast Asia congregated in trees using anticipatory time-measuring and flashing in rhythmic synchrony [4], and in pancreatic islet cells generating insulin [5].

Synchronization of coupled sub-systems in both natural and engineered systems is a commonplace occurrence, and its existence and analysis in mechanical systems has received wide attention and publicity. In the 17th century, the Dutch scientist Ch. Huygens observed synchronized behavior of two pendulum clocks hanging from a common imperceptibly movable support in his serendipitous discovery [6]. Over the span of three centuries after the starting point constituted by Ch. Huygens, there has been an increasing research interest in synchronized pendulums, and many profound theoretical results and synchronization criteria have been contributed towards a better understanding of the charming, insightful phenomenon. For example, in Ref. [7], an improved model for the original Huygens experimental setup of pendulum clocks on synchronization has been presented and new kinds of limit behavior have been observed. Jovanovic and Koshkin [8] applied the Poincare method to the Huygens system and studied two models of connected pendulum clocks synchronizing their oscillations. The studies in Ref. [9] showed that two double pendula self-excited by the escapement mechanism hanging from a horizontally movable beam can reach synchronization. In Ref. [10], the authors studied the phenomena of synchronization in a system of two coupled oscillators and developed a suitable controller that can not only synchronize and stabilize two coupled pendula, but also conserve the mechanical energy when the system evolves towards the synchronous states. More recently, controlled synchronization in asymmetrically coupled oscillators has been studied in Ref. [11] by using a suitable feedback ‘symmetrizing’ controller; consequently, the system exhibits complete in-phase and anti-phase synchronization.

Research and analysis on synchronization in mechanical systems can be naturally associated with the area of engineering sciences, and in particular with vibration utilization engineering. Starting with the pioneering work of I.I. Blekhman [12,13,14,15,16], in which the analytical method of direct motion separation was applied to theoretically explain the synchronization mechanism of exciters, there has been an increasing research interest in this type of synchronization concerning engineering problems. Wen et al. [17,18,19,20,21] simplified Blekhman’s method by selecting the phase differences among exciters as variables and established conditions of implementing synchronization and stability of synchronization for exciters in a vibrating system. On this basis, Zhao et al. [22,23,24,25,26] applied the average method of modified small parameters and introduced variable perturbation parameters to average angular velocity of exciters and their phase differences, so as to provide significant investigations on the feature of frequency capture and the dynamic characteristics of identical or non-identical exciters in a vibrating system of plane or spatial motion. In Refs. [27,28,29,30], the authors studied the synchronization regime of multiple identical or non-identical exciters in a single rigid base, and concluded that using multiple exciters instead of two exciters in a single rigid base cannot improve the effective power of a far-resonant vibrating system, since the exciting forces are mutually offset due to strong inner coupling. Now that the self-synchronization of multi-exciters cannot satisfy the requirements of more energy supply and highly efficient productive processing, the controlled and composite synchronization methods presented in Refs. [30,31,32,33,34] are feasible choices to meet these demands. Lately, the synchronization of two DC motors supported by a rectangular plate and the plate vibration control in synchronous states has been analyzed in Ref. [35], while Kong et al. [36,37] analyzed the Sommerfeld effect and self-synchronization of two or three non-ideal induction motors installed on a simply supported beam structure.

In this paper, we consider the synchronization of an exciter and three CRs with non-identical friction coefficients in a far-resonant vibrating system driven by an exciter. The synchronization mechanism can be realized if an exciter and three CRs rotate synchronously around their rotating centers. Based on this synchronization mechanism, a new type of vibrating crusher or vibrating mill with multiple working chambers can be developed. The materials to break, such as ores, minerals, slags, and concretes, can be located in the clearance area between three CRs and their corresponding inner wall. During the synchronous operation process, three CRs rotate along the inner wall with the effective extrusion, and the goal of crushing the materials is achieved. In practical applications, the friction coefficient $f_{Ri}$ between the surface of CR $i$ and the surface of its cavity is inevitably different from the other two friction coefficients due to the variability in their manufacturing process. This results in an asymmetric coupling: the influence of the rigid platform on each CR will be different from the influence that the rigid platform exerts on the other two CRs. As a result of the lack of symmetry in the coupled system, the synchronization of rotating elements cannot be implemented in some coupled dynamical systems [11,38,39,40,41,42]. Therefore, it is urgent to consider what may happen in the proposed vibrating system and whether an exciter and three CRs can synchronize or not as a result of the non-identical dry friction coefficients. We endeavor to provide an answer to these questions in our studies.

In addition, the significance and necessity of undertaking such a study are stated as follows:

(1) In a sense, the proposed coupled system in this paper can be treated as an improved model of that in Ref. [43]. In Ref. [43], the in-phase synchronization of two exciters is a fundamental precondition for the synchronous rotation of the CR. However, the $r_l$ ($r_l$ is the length ratio of the distance between the rotary center of an exciter and the mass center of the rigid frame to the equivalent rotating radius of the vibrating system about the centroid of the rigid frame) should be large enough to synchronize two exciters in in-phase. In this case, with a high value of $r_l$, the vibrating system will most certainly come in an extra-large size, which is not actually desired in practice. The structural weakness confines the design and applications of vibrating machines. This paper introduces a new coupling structure, which generates an ever-present and stable synchronization regime under variation of the key parameter $r_l$, resulting in a small-sized and compact structure very easy to design and implement in various engineering practices.

(2) To reach the synchronous operation of the exciter and CRs, the induction motor would have to offer extra energy to the CRs indirectly since the direct driving sources of the CRs are absent. It is thus much more difficult to perform the synchronous motion of the exciter and CRs, compared to that of the multiple exciters. Accordingly, seeing whether the energy supplied by the induction motor is sufficient to initiate and sustain this type of synchronization is essential. In this paper, the total energy balance is considered, and the concept of the necessary electromagnetic torque is introduced in the synchronization criterion. If an induction motor’s rated electromagnetic torque surpasses the necessary electromagnetic torque for the vibrating system, then an induction motor would have the capability to achieve the synchronous motion. Additionally, since there are some parameters in the expression of the necessary electromagnetic torque that dominantly influence the occurrence of the synchronization, a detailed investigation of the key parameters regarding the onset of the synchronous motion is conducted. In light of this investigation, when using the synchronization mechanism in this study to develop certain types of vibratory crushers or vibratory mills, the chosen values of the key parameters can ensure that the induction motor in the vibrating system provides adequate energy.

(3) The synchronization and coupling dynamical characteristics of multiple exciters with direct driving sources have been investigated extensively [21,22,23,24,25,26,27,28,29], whereas the synchronization characteristics and stability analysis of the synchronous states between multiple coupled CRs with non-identical dry friction coefficients are less studied in the available literature. Therefore, a detailed investigation into the dynamics and synchronization of three CRs is essential for not only scientific interest, but also for application purposes. For instance, the authors must figure out how variations in system parameters would affect the synchronization ability, stability, and clinging ability of three CRs, so as to tune the synchronization performance to the ideal case in actual projects.

The structure of the paper is organized as follows. Section 2 presents the modelling of the vibrating system. In Section 3, the existence and stability of synchronous solutions are analytically discussed. In particular, the concept of dynamic coupling torque is introduced. Section 4 explores the nonlinear relationships between the structure parameters and the coupling synchronization characteristics, and reveals how synchronization performance is consequently influenced. Then, the analytic results are illustrated by means of computer simulations in Section 5. Cases of identical and nonidentical friction coefficients are performed in order to determine the possible synchronization behaviour. Additionally, the role played by dynamic coupling torques is analyzed. Finally, in Section 6, the main conclusions are given.

2. Description of the Dynamic Model

Consider the mechanical system depicted in Figure 1a, which represents a vibrating setup of an exciter driven by the induction motor and three CRs placed inside the cavities of a rigid platform. The rigid platform is elastically attached to a fixed support by four weightless springs in a symmetrical arrangement. The induction motor is installed in the center of mass or centroid of the rigid platform, and three cavities are uniformly distributed in a circumference of the rigid platform. In the process of operation, the cylindrical rollers absorb vibration energy from the rigid platform, which becomes a rotational motion along the surface of the cavities and thus automatically tracks the rotating exciter. The coordinate system is shown in Figure 1b. It consists of a fixed frame denoted by $oxy$, a non-rotating moving frame denoted by $G{x}^{\prime }{y}^{\prime }$ and a moving frame (attached at the rigid platform) denoted by $G{x}^{″}{y}^{″}$.

The rigid platform can vibrate in horizontal and vertical directions and swing about its center of mass, its movements are described by coordinates $x$, $y$ and $\psi$, respectively. The motion of an exciter is described by angle ${\phi }_1$ and the rotation angle of CR $i$ is denoted by ${\phi }_{Ri}$.

In the moving frame $G{x}^{″}{y}^{″}$, the coordinates of the mass center of an exciter and three CRs, $x_E^{″}$ and $x_{Ri}^{″}$ ($i=1,2,3$), can be written as:

$$x_E^{″}=\left\{\begin{array}{c}r\cos {\phi }_1\\ -r\sin {\phi }_1\end{array}\right\}$$

$$x_{Ri}^{″}=\left\{\begin{array}{c}{l}_0\cos {\beta }_i+r_{Ri}\cos {\phi }_{Ri}\\ l_0\sin {\beta }_i-r_{Ri}\sin {\phi }_{Ri}\end{array}\right\},i=1,2,3$$

(1)

In the fixed frame $oxy$, the coordinates of the mass center of an exciter and three CRs, $x_E$ and $x_{Ri}$ ($i=1,2,3$), can be expressed as:

$$x_i=x+R{x}_i^{″},i=E,R_1,R_2,R_3$$

(2)

The expression of the kinetic energy is given as follows:

$$\begin{array}{l}T=\frac{1}{2}m({\dot{x}}^2+{\dot{y}}^2)+\frac{1}{2}{J}_m{\dot{\psi }}^2+\frac{1}{2}{J}_1{\dot{\varphi }}_1^2+\frac{1}{2}{m}_1{\dot{\mathbf{x}}}_E^{\mathrm{T}}{\dot{\mathbf{x}}}_E\\ +\frac{1}{2}{\displaystyle \sum _{i=1}^3{m}_{Ri}{r}_{Ri}^2{\dot{\varphi }}_{Ri}^2}+\frac{1}{2}{\displaystyle \sum _{i=1}^3{m}_{Ri}{\dot{\mathbf{x}}}_{Ri}^{\mathrm{T}}{\dot{\mathbf{x}}}_{Ri}}+\frac{1}{2}{\displaystyle \sum _{i=1}^3{I}_{Ri}{\dot{\theta }}_{Ri}^2}.\end{array}$$

(3)

During the running process of the vibrating system, the coordinates of the connection points between spring $i$ ($i=1,2,3,4$) and the rigid platform in the fixed frame $oxy$ is described by

$$x_{ki}=x+R{x}_{k0i},i=1,2,3,4$$

(4)

Hence, the expression of the potential energy is given as follows:

$$V=\frac{1}{2}{\displaystyle \sum _{i=1}^4{(x_{ki}-x_{k0i})}^{\mathrm{T}}{K}_i(x_{ki}-x_{k0i})}$$

(5)

The viscous dissipation function is determined by

$$D=\frac{1}{2}{\displaystyle \sum _{i=1}^4{\dot{\mathbf{x}}}_{ki}^{\mathrm{T}}{F}_i{\dot{\mathbf{x}}}_{ki}}$$

(6)

The dynamic differential equations of motion can be formulated with the use of Lagrange’s equation:

$$\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial (T-V)}{\partial {\dot{q}}_i}-\frac{\partial (T-V)}{\partial q_i}+\frac{\partial D}{\partial {\dot{q}}_i}=Q_i$$

(7)

where $q={\left\{x,y,\psi ,{\varphi }_1,{\varphi }_{R1},{\varphi }_{R2},{\varphi }_{R3}\right\}}^{\mathrm{T}}$ is introduced as the system generalized coordinates, and the generalized forces associated with the system generalized coordinates are given as follows:

$$Q_x=Q_y=Q_{\psi }=0,Q_{\varphi 1}=T_{\mathrm{e}1}$$

$$Q_{\phi Ri}=-f_{Ri}{r}_{Ri}{N}_i,i=1,2,3$$

$$N_i=m_{Ri}(\ddot{x}\cos {\varphi }_{Ri}-\ddot{y}\sin {\varphi }_{Ri}+r_{Ri}{\dot{\varphi }}_{Ri}^2+g\sin {\varphi }_{Ri}),i=1,2,3$$

(8)

The aforementioned modelling process yields to the following equations of motion:

$$M\ddot{x}+f_x\dot{x}+k_{x}x=m_{1}r({\dot{\varphi }}_1^2\cos {\varphi }_1+{\ddot{\varphi }}_1\sin {\varphi }_1)+{\displaystyle \sum _{i=1}^3{m}_{Ri}{r}_{Ri}({\dot{\varphi }}_{Ri}^2\cos {\varphi }_{Ri}+{\ddot{\varphi }}_{Ri}\sin {\varphi }_{Ri})}$$

$$M\ddot{y}+f_y\dot{y}+k_{y}y=m_{1}r(-{\dot{\varphi }}_1^2\sin {\varphi }_1+{\ddot{\varphi }}_1\cos {\varphi }_1)+{\displaystyle \sum _{i=1}^3{m}_{Ri}{r}_{Ri}(-{\dot{\varphi }}_{Ri}^2\sin {\varphi }_{Ri}+{\ddot{\varphi }}_{Ri}\cos {\varphi }_{Ri})}$$

$$\begin{array}{c}J\ddot{\psi }+f_{\psi }\dot{\psi }+k_{\psi }\psi ={\displaystyle \sum _{i=1}^3{l}_0{m}_{Ri}{r}_{Ri}[{\dot{\varphi }}_{Ri}^2\sin ({\varphi }_{Ri}+{\beta }_i)-{\ddot{\varphi }}_{Ri}\cos ({\varphi }_{Ri}+{\beta }_i)]}\\ J_1{\ddot{\varphi }}_1+f_1{\dot{\varphi }}_1=T_{\mathrm{e}1}-m_{1}r(-\ddot{x}\sin {\varphi }_1-\ddot{y}\cos {\varphi }_1+r\ddot{\psi })\end{array}$$

$$\begin{array}{c}{J}_{R1}{\ddot{\varphi }}_{R1}=m_{R1}{r}_{R1}(\ddot{x}\sin {\varphi }_{R1}+\ddot{y}\cos {\varphi }_{R1})-l_0{m}_{R1}{r}_{R1}[{\dot{\psi }}^2\sin ({\varphi }_{R1}+{\beta }_1)+\ddot{\psi }\cos ({\varphi }_{R1}+{\beta }_1)]\\ +m_{R1}{r}_{R1}g\cos {\varphi }_{R1}-m_{R1}{r}_{R1}^2\ddot{\psi }-f_{R1}{r}_{R1}{N}_1\end{array}$$

$$\begin{array}{c}{J}_{R2}{\ddot{\varphi }}_{R2}=m_{R2}{r}_{R2}(\ddot{x}\sin {\varphi }_{R2}+\ddot{y}\cos {\varphi }_{R2})-l_0{m}_{R2}{r}_{R2}[{\dot{\psi }}^2\sin ({\varphi }_{R2}+{\beta }_2)+\ddot{\psi }\cos ({\varphi }_{R2}+{\beta }_2)]\\ +m_{R2}{r}_{R2}g\cos {\varphi }_{R2}-m_{R2}{r}_{R2}^2\ddot{\psi }-f_{R2}{r}_{R2}{N}_2\end{array}$$

$$\begin{array}{c}{J}_{R3}{\ddot{\varphi }}_{R3}=m_{R3}{r}_{R3}(\ddot{x}\sin {\varphi }_{R3}+\ddot{y}\cos {\varphi }_{R3})-l_0{m}_{R3}{r}_{R3}[{\dot{\psi }}^2\sin ({\varphi }_{R3}+{\beta }_3)+\ddot{\psi }\cos ({\varphi }_{R3}+{\beta }_3)]\\ +m_{R3}{r}_{R3}g\cos {\varphi }_{R3}-m_{R3}{r}_{R3}^2\ddot{\psi }-f_{R3}{r}_{R3}{N}_3\end{array}$$

(9)

3. Existence of Synchronous Motion

In what follows, it is assumed that

$$m_{R1}=m_{R2}=m_{R3}=m_R,r_{R1}=r_{R2}=r_{R3}=r_R$$

(10)

Additionally, the average phase among the exciter and three CRs is assumed to be $\phi$, and their phase differences are $2{\alpha }_1$, $2{\alpha }_2$ and $2{\alpha }_3$, i.e.,

$$\phi =({\phi }_1+{\phi }_{R1}+{\phi }_{R2}+{\phi }_{R3})/4,{\phi }_1-{\phi }_{R1}=2{\alpha }_1$$

$${\phi }_{R1}-{\phi }_{R2}=2{\alpha }_2,{\phi }_{R2}-{\phi }_{R3}=2{\alpha }_3$$

(11)

Then, Equation (11) is transformed to

$${\phi }_1=\phi +\frac{3}{2}{\alpha }_1+{\alpha }_2+\frac{1}{2}{\alpha }_3=\phi +{\varsigma }_1$$

$${\phi }_{R1}=\phi -\frac{1}{2}{\alpha }_1+{\alpha }_2+\frac{1}{2}{\alpha }_3=\phi +{\varsigma }_{R1}$$

$${\phi }_{R2}=\phi -\frac{1}{2}{\alpha }_1-{\alpha }_2+\frac{1}{2}{\alpha }_3=\phi +{\varsigma }_{R2}$$

$${\phi }_{R3}=\phi -\frac{1}{2}{\alpha }_1-{\alpha }_2-\frac{3}{2}{\alpha }_3=\phi +{\varsigma }_{R3}$$

(12)

Since a vibrating system’s motion is periodic, the variation in average mechanical angular velocity $\dot{\phi }$ is periodic as well. If the least common multiple period is assumed to be $T_{\mathrm{LCMP}}$, then the following should be satisfied:

$${\omega }_{\mathrm{m}}=\frac{1}{T_{\mathrm{LCMP}}}{\displaystyle {\int }_{t^{\prime }}^{t^{\prime }+T_{\mathrm{LCMP}}}\dot{\phi }(t)\mathrm{d}t}=\mathrm{constant}$$

(13)

According to the average method of modified small parameters [22,23,24,25,26,27], it is assumed that the instantaneous change coefficient of $\dot{\phi }$ around ${\omega }_{\mathrm{m}}$ and that of ${\dot{\alpha }}_i$ to be ${\epsilon }_0$ and ${\epsilon }_i$ (${\epsilon }_0$ and ${\epsilon }_i$ are functions of time $t$, $i=1,2,3$). This yields

$$\dot{\phi }=(1+{\epsilon }_0){\omega }_{\mathrm{m}},{\dot{\alpha }}_i={\epsilon }_i{\omega }_{\mathrm{m}},i=1,2,3$$

(14)

Introducing Equation (14) into the derivative of Equation (12), we get

$${\dot{\phi }}_1=(1+{\epsilon }_0+\frac{3}{2}{\epsilon }_1+{\epsilon }_2+\frac{1}{2}{\epsilon }_3){\omega }_{\mathrm{m}}=(1+{\upsilon }_1){\omega }_{\mathrm{m}}$$

$${\dot{\phi }}_{R1}=(1+{\epsilon }_0-\frac{1}{2}{\epsilon }_1+{\epsilon }_2+\frac{1}{2}{\epsilon }_3){\omega }_{\mathrm{m}}=(1+{\upsilon }_{R1}){\omega }_{\mathrm{m}}$$

$${\dot{\phi }}_{R2}=(1+{\epsilon }_0-\frac{1}{2}{\epsilon }_1-{\epsilon }_2+\frac{1}{2}{\epsilon }_3){\omega }_{\mathrm{m}}=(1+{\upsilon }_{R2}){\omega }_{\mathrm{m}}$$

$${\dot{\phi }}_{R3}=(1+{\epsilon }_0-\frac{1}{2}{\epsilon }_1-{\epsilon }_2-\frac{3}{2}{\epsilon }_3){\omega }_{\mathrm{m}}=(1+{\upsilon }_{R3}){\omega }_{\mathrm{m}}$$

$${\ddot{\phi }}_1={\dot{\upsilon }}_1{\omega }_{\mathrm{m}},{\ddot{\phi }}_{R1}={\dot{\upsilon }}_{R1}{\omega }_{\mathrm{m}},{\ddot{\phi }}_{R2}={\dot{\upsilon }}_{R2}{\omega }_{\mathrm{m}},{\ddot{\phi }}_{R3}={\dot{\upsilon }}_{R3}{\omega }_{\mathrm{m}}$$

(15)

At the steady-state, the ${\ddot{\phi }}_1$ and ${\ddot{\phi }}_{Ri}$ can be neglected in the first three formulae of Equation (9), and the responses are given as:

$$x=-\frac{r_{m}r}{{\mu }_x}\cos (\phi +{\varsigma }_1+{\gamma }_x)-{\displaystyle \sum _{i=1}^3\frac{{\eta }_m{\eta }_r{r}_{m}r}{{\mu }_x}\cos (\phi +{\varsigma }_{Ri}+{\gamma }_x)}$$

$$y=\frac{r_{m}r}{{\mu }_y}\sin (\phi +{\varsigma }_1+{\gamma }_y)+{\displaystyle \sum _{i=1}^3\frac{{\eta }_m{\eta }_r{r}_{m}r}{{\mu }_y}\sin (\phi +{\varsigma }_{Ri}+{\gamma }_y)}$$

$$\psi =-{\displaystyle \sum _{i=1}^3\frac{{\eta }_m{\eta }_r{r}_{m}r{r}_l}{{\mu }_{\psi }{l}_{\mathrm{e}}}\sin (\phi +{\varsigma }_{Ri}+{\beta }_i+{\gamma }_{\psi })}$$

(16)

The electromagnetic torque of the induction motor can be written as the following expression [25,26] in the situation of operating in the vicinity of ${\omega }_{\mathrm{m}}$:

$$T_{\mathrm{e}1}=T_{\mathrm{e}01}-k_{\mathrm{e}01}{\upsilon }_1$$

(17)

3.1. Synchronization Criterion

By differentiating Equation (16), the $\ddot{x}$, $\ddot{y}$, $\dot{\psi }$ and $\ddot{\psi }$ can be determined. Introducing these derivatives into the last four formulae of Equation (9), and integrating these formulae over $\phi =0~2\mathsf{\pi }$, the averaged balance equations of an exciter and three CRs are given by

$$T_{\mathrm{o}1}=T_{\mathrm{e}01}-f_1{\omega }_{\mathrm{m}}=T_{\mathrm{u}}({\chi }_{\mathrm{f}1}+{\chi }_{\mathrm{a}1})$$

(18)

$$T_{\mathrm{o}2}=T_{\mathrm{e}02}-f_{R1}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01})=T_{\mathrm{u}}({\chi }_{\mathrm{f}2}+{\chi }_{\mathrm{a}2})$$

(19)

$$T_{\mathrm{o}3}=T_{\mathrm{e}03}-f_{R2}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01})=T_{\mathrm{u}}({\chi }_{\mathrm{f}3}+{\chi }_{\mathrm{a}3})$$

(20)

$$T_{\mathrm{o}4}=T_{\mathrm{e}04}-f_{R3}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01})=T_{\mathrm{u}}({\chi }_{\mathrm{f}4}+{\chi }_{\mathrm{a}4})$$

(21)

With $T_{\mathrm{o}i}$ ($i=1,2,3,4$) denoting the output electromagnetic torque of motor $i$ (motors 2, 3 and 4 are three fictional motors, i.e., $T_{\mathrm{e}02}=T_{\mathrm{e}03}=T_{\mathrm{e}04}=0$). Explicit expressions for the coefficients in Equations (18)–(21) are provided in Appendix A.

From Equations (18)–(21), calculating the difference of $T_{\mathrm{o}1}$ and $T_{\mathrm{o}2}$, $T_{\mathrm{o}2}$ and $T_{\mathrm{o}3}$, $T_{\mathrm{o}3}$ and $T_{\mathrm{o}4}$, we obtain

$$\Delta T_{\mathrm{o}12}=T_{\mathrm{o}1}-T_{\mathrm{o}2}=T_{\mathrm{u}}[({\chi }_{\mathrm{f}1}+{\chi }_{\mathrm{a}1})-({\chi }_{\mathrm{f}2}+{\chi }_{\mathrm{a}2})]$$

(22)

$$\Delta T_{\mathrm{o}23}=T_{\mathrm{o}2}-T_{\mathrm{o}3}=T_{\mathrm{u}}[({\chi }_{\mathrm{f}2}+{\chi }_{\mathrm{a}2})-({\chi }_{\mathrm{f}3}+{\chi }_{\mathrm{a}3})]$$

(23)

$$\Delta T_{\mathrm{o}34}=T_{\mathrm{o}3}-T_{\mathrm{o}4}=T_{\mathrm{u}}[({\chi }_{\mathrm{f}3}+{\chi }_{\mathrm{a}3})-({\chi }_{\mathrm{f}4}+{\chi }_{\mathrm{a}4})]$$

(24)

Additionally, calculating the sum of Equations (22) and (23), Equations (23) and (24), Equations (22)–(24), we can also obtain the differences of the output electromagnetic torque between motors 1 and 3 ($\Delta T_{\mathrm{o}13}$), motors 2 and 4 ($\Delta T_{\mathrm{o}42}$), motors 1 and 4 ($\Delta T_{\mathrm{o}14}$). Consequently, Equations (22)–(24) can describe the coupling characteristics of the exciter and three CRs.

Then, Equations (22)–(24) are put into a new arrangement as follows:

$$\frac{\Delta T_{\mathrm{o}12}}{T_{\mathrm{u}}}-W_{\mathrm{s}01}+{\eta }_m^2{\eta }_r^2{W}_{\mathrm{s}02}=T_{\mathrm{c}12}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)$$

(25)

$$\frac{\Delta T_{\mathrm{o}23}}{T_{\mathrm{u}}}=T_{\mathrm{c}23}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)$$

(26)

$$\frac{\Delta T_{\mathrm{o}34}}{T_{\mathrm{u}}}=T_{\mathrm{c}34}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)$$

(27)

The left side of Equations (25)–(27) stands for the differences of the dimensionless residual torques between two motors, while $T_{\mathrm{c}ij}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)$ from the right side is described as the dimensionless coupling torque. From the expression of $T_{\mathrm{c}ij}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)$ listed in Appendix B, one can obtain that $T_{\mathrm{c}ij}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)$ is a limited function of ${\overline{\alpha }}_1$, ${\overline{\alpha }}_2$, and ${\overline{\alpha }}_3$, i.e.,

$$\left|T_{\mathrm{c}ij}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)\right|\le T_{\mathrm{c}ij\max }(ij=12,23,34)$$

(28)

Therefore, the synchronization criterion is given by

$$\left|\frac{\Delta T_{\mathrm{o}12}}{T_{\mathrm{u}}}-W_{\mathrm{s}01}+{\eta }_m^2{\eta }_r^2{W}_{\mathrm{s}02}\right|\le T_{\mathrm{c}12\max }$$

(29)

$$\left|\frac{\Delta T_{\mathrm{o}23}}{T_{\mathrm{u}}}\right|\le T_{\mathrm{c}23\max }$$

(30)

$$\left|\frac{\Delta T_{\mathrm{o}34}}{T_{\mathrm{u}}}\right|\le T_{\mathrm{c}34\max }$$

(31)

The criterion states that “the absolute values of the differences of the dimensionless residual torques between two motors should not exceed the maximum of their dimensionless coupling torques”. When the criterion is satisfied, the variables ${\overline{\alpha }}_1$, ${\overline{\alpha }}_2$, ${\overline{\alpha }}_3$ and ${\omega }_{\mathrm{m}}$ can be obtained by solving Equations (18)–(21). The solution to Equations (18)–(21) is represented as ${\overline{\alpha }}_{10}$, ${\overline{\alpha }}_{20}$, ${\overline{\alpha }}_{30}$ and ${\omega }_{\mathrm{m}0}$, respectively.

Combining Equations (18)–(21), we have

$${\displaystyle \sum _{i=1}^4{T}_{\mathrm{e}0i}}-f_1{\omega }_{\mathrm{m}}-{\displaystyle \sum _{i=1}^3{f}_{Ri}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01})-{\displaystyle \sum _{i=1}^4{T}_{\mathrm{u}}({\chi }_{\mathrm{f}i}+{\chi }_{\mathrm{a}i})}}=0$$

(32)

Here, Equation (32) is referred to as the averaged balance equation of the total system at the steady-state. Referring to the literature [32,33], for the purpose of reaching synchronization, the motor’s electromagnetic torque should be less than or equal to its rated value. Thus, the condition of implementing synchronous motion is also given by

$$T_{\mathrm{e}01}=f_1{\omega }_{\mathrm{m}}+{\displaystyle \sum _{i=1}^3{f}_{Ri}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01})}+{\displaystyle \sum _{i=1}^4{T}_{\mathrm{u}}({\chi }_{\mathrm{f}i}+{\chi }_{\mathrm{a}i})},\left|T_{\mathrm{e}01}\right|\le T_{\mathrm{eN}1}$$

(33)

where $T_{\mathrm{e}01}$ represents the necessary electromagnetic torque for the vibrating system; $T_{\mathrm{eN}1}$ represents the motor’s rated electromagnetic torque. Next, the expression of the averaged nondimensional loading torque is given as:

$$T_{\mathrm{a}}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)=\frac{1}{4T_{\mathrm{u}}}{\displaystyle \sum _{i=1}^4{T}_{\mathrm{o}i}}$$

(34)

Additionally, its explicit expression is provided in Appendix B. From the expression, it can be seen that $T_{\mathrm{a}}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)$ is a limited function, i.e.,

$$\left|T_{\mathrm{a}}({\overline{\alpha }}_1,{\overline{\alpha }}_2,{\overline{\alpha }}_3)\right|\le T_{\mathrm{amax}}$$

(35)

The synchronization ability coefficient is described by

$${\zeta }_{ij}=\frac{T_{\mathrm{c}ij\max }}{T_{\mathrm{a}\max }},ij=12,23,34$$

(36)

3.2. The Stability of Synchronous Solutions

Equations (18)–(21) after linearization around ${\overline{\alpha }}_{10}$, ${\overline{\alpha }}_{20}$, ${\overline{\alpha }}_{30}$ and ${\omega }_{\mathrm{m}0}$, are respectively,

$$k\left(\frac{3}{2}{\dot{\overline{\alpha }}}_1+{\dot{\overline{\alpha }}}_2+\frac{1}{2}{\dot{\overline{\alpha }}}_3\right)=-{\displaystyle \sum _{i=1}^3{T}_{\mathrm{u}}{\left[\frac{\partial ({\chi }_{\mathrm{f}1}+{\chi }_{\mathrm{a}1})}{\partial {\overline{\alpha }}_i}\right]}_0\Delta {\alpha }_i}$$

(37)

$$k_{\mathrm{R}1}\left(-\frac{1}{2}{\dot{\overline{\alpha }}}_1+{\dot{\overline{\alpha }}}_2+\frac{1}{2}{\dot{\overline{\alpha }}}_3\right)=-{\displaystyle \sum _{i=1}^3{T}_{\mathrm{u}}{\left[\frac{\partial ({\chi }_{\mathrm{f}2}+{\chi }_{\mathrm{a}2})}{\partial {\overline{\alpha }}_i}\right]}_0\Delta {\alpha }_i}$$

(38)

$$k_{\mathrm{R}2}\left(-\frac{1}{2}{\dot{\overline{\alpha }}}_1-{\dot{\overline{\alpha }}}_2+\frac{1}{2}{\dot{\overline{\alpha }}}_3\right)=-{\displaystyle \sum _{i=1}^3{T}_{\mathrm{u}}{\left[\frac{\partial ({\chi }_{\mathrm{f}3}+{\chi }_{\mathrm{a}3})}{\partial {\overline{\alpha }}_i}\right]}_0\Delta {\alpha }_i}$$

(39)

$$k_{\mathrm{R}3}\left(-\frac{1}{2}{\dot{\overline{\alpha }}}_1-{\dot{\overline{\alpha }}}_2-\frac{3}{2}{\dot{\overline{\alpha }}}_3\right)=-{\displaystyle \sum _{i=1}^3{T}_{\mathrm{u}}{\left[\frac{\partial ({\chi }_{\mathrm{f}4}+{\chi }_{\mathrm{a}4})}{\partial {\overline{\alpha }}_i}\right]}_0\Delta {\alpha }_i}$$

(40)

where ${\left[•\right]}_0$ denotes the values for ${\overline{\alpha }}_1={\overline{\alpha }}_{10}$, ${\overline{\alpha }}_2={\overline{\alpha }}_{20}$, ${\overline{\alpha }}_3={\overline{\alpha }}_{30}$; $\Delta {\alpha }_i={\overline{\alpha }}_i-{\overline{\alpha }}_{i0}$, $i=1,2,3$; $k=k_{\mathrm{e}01}+f_1$; $k_{Ri}={\left[\frac{\mathrm{d}(f_{Ri}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01}))}{\mathrm{d}{\omega }_{\mathrm{m}}}\right]}_{{\omega }_{\mathrm{m}}={\omega }_{\mathrm{m}0}}$, $i=1,2,3$.

Calculating the differences of Equations (37) and (38), Equations (38) and (39), Equations (39) and (40), and writing the results into a matrix form, then we have

$$E\Delta \dot{\alpha }=D\Delta \alpha$$

(41)

with $E={[e_{ij}]}_{3\times 3}$ and $D={[d_{ij}]}_{3\times 3}$. The explicit expressions for $e_{ij}$ and $d_{ij}$ are provided in Appendix C.

Rewriting Equation (41), we obtain

$$\Delta \dot{\alpha }=C\Delta \alpha ,C=E^{-1}D$$

(42)

Solving the determinant equation $\det (\lambda I-C)=0$, we can deduce the characteristic equation for the eigenvalue $\lambda$ as follows:

$${\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0$$

(43)

By means of Routh-Hurwitz principle, it can be determined that when Equation (44) is satisfied, the zero solutions of Equation (42) are stable.

$$a_1>0,a_3>0,a_1{a}_2>a_3$$

(44)

That is to say, Equation (44) is the stability criterion of synchronous operation. Simplifying Equation (44), we have

$$H_1>0,H_2>0,H_3>0$$

(45)

Here, $H_i$ is referred to as the stability coefficient. The greater the stability coefficient $H_i$ is, the stronger the stability is. Explicit expression for $H_i$ is provided in Appendix C.

3.3. Clinging Ability Coefficients of Three CRs

When the vibrating system operates synchronously at the steady-state, the three CRs should perform rotational motion along the surfaces of the cavities and cling to them tightly, for the purpose of successfully crushing the materials. Consequently, it should be clear that the reaction force on the three CRs exerted by the cavities’ surfaces is more than zero, i.e.,

$$N_i=m_R(\ddot{x}\cos {\varphi }_{Ri}-\ddot{y}\sin {\varphi }_{Ri}+r_R{\dot{\varphi }}_{Ri}^2+g\sin {\varphi }_{Ri})>0,i=1,2,3$$

(46)

Simplifying Equation (46), we obtain

$$r_R{\omega }_{\mathrm{m}}^2>-\ddot{x}\cos {\varphi }_{Ri}+\ddot{y}\sin {\varphi }_{Ri}-g\sin {\varphi }_{Ri},i=1,2,3$$

(47)

Since in the far-resonant vibrating system, the operating frequency is much greater than the natural frequency and the damping constants are small [19,20], it is considered that ${\gamma }_x\approx {\gamma }_y={\gamma }_0$, ${\mu }_x\approx {\mu }_y={\mu }_0$. Introducing Equations (12) and (16) into Equation (47), we have

$$C_{Ri}<1,i=1,2,3.$$

(48)

Here, the symbol $C_{Ri}$ is referred to as the clinging ability coefficient, and the corresponding detailed expressions are provided in Appendix D. Only if $C_{Ri}$ is much less than one do the CRs possess the strong clinging ability, which is indispensable for crushing the materials.

3.4. Analysis of the Coupling Dynamic Characteristics of an Exciter and Three CRs

Equations (18)–(21) can also be written as:

$$T_{\mathrm{e}01}-f_1{\omega }_{\mathrm{m}}=T_{\mathrm{u}}({\chi }_{\mathrm{f}1}+{\chi }_{\mathrm{a}1})$$

(49)

$$-f_{R1}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01})=T_{\mathrm{u}}({\chi }_{\mathrm{f}2}+{\chi }_{\mathrm{a}2})$$

(50)

$$-f_{R2}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01})=T_{\mathrm{u}}({\chi }_{\mathrm{f}3}+{\chi }_{\mathrm{a}3})$$

(51)

$$-f_{R3}{\eta }_m{\eta }_r^2{T}_{\mathrm{u}}(2+{\eta }_m{W}_{\mathrm{c}01})=T_{\mathrm{u}}({\chi }_{\mathrm{f}4}+{\chi }_{\mathrm{a}4})$$

(52)

On the right sides of Equations (49)–(52), $T_{\mathrm{u}}{\chi }_{\mathrm{f}i}$ and $T_{\mathrm{u}}{\chi }_{\mathrm{a}i}$ ($i=1,2,3,4$) are the load torques that a vibrating system acts upon four motors. In the expression of $T_{\mathrm{u}}{\chi }_{\mathrm{f}i}$, the components independent of the phase differences represent the load torques originating from the vibration excited by its own exciter or CR (referred to as, simply, $T_{Ai}$, $i=1,2,3,4$); the components dependent on f_Ri cosine functions of phase differences describe the coupled load torques originating from the vibrations excited by others besides its own exciter or CR (referred to as, simply, $T_{Bi}$, $i=1,2,3,4$); the components dependent on $f_{Ri}$ and sine functions of phase differences exist only in ${\chi }_{\mathrm{f}2}$, ${\chi }_{\mathrm{f}3}$ and ${\chi }_{\mathrm{f}4}$, which depict the exclusive coupled load torques in CR $i$ originating from the vibrations excited by the other CRs and an exciter (referred to as, simply, $T_{Ci}$, $i=1,2,3$).

In the expression of $T_{\mathrm{u}}{\chi }_{\mathrm{a}i}$, the components dependent on sine functions of phase differences are to act as the coupled load torques of a phase leading exciter or CR and the coupled driving torques of a phase lagging exciter or CR, which are defined as the dynamic coupling torques among an exciter and three CRs (referred to as, simply, $T_{DCi}$, $i=1,2,3,4$). In time of synchronous and steady states, the total dynamic coupling torques (referred to as, simply, $T_{DC}$, i.e., ${\displaystyle \sum _{i=1}^4{T}_{DCi}}=T_{DC}$) are equal to zero. In other words, the total dynamic coupling torques do not generate or consume energy but transmit energy between an exciter and three CRs through the vibrating system. The other components dependent on $f_{Ri}$ and cosine functions of phase differences express the exclusive coupled load torques in CR $i$ originating from the vibrations excited by the other CRs and an exciter (referred to as, simply, $T_{Di}$, $i=1,2,3$). Since an exciter and three CRs associate with each other by means of sine and cosine functions of phase differences, the phase differences are referred to as the coupled phases of an exciter and three CRs.

The left sides of Equations (50)–(52) serve as a kind of friction torque from three CRs, which will dissipate vibration energy of the rigid platform and kinetic energy of three CRs. Due to the coupled torques $T_{Bi}$, $T_{Ci}$, $T_{DCi}$ and $T_{Di}$, the energy can be transferred from an exciter to three CRs through the moving rigid platform and compensates for three CRs’ energy consumption.

4. The Influence of Structural Parameters on Synchronization and Stability

In Section 3, the theoretical results regarding the synchronization and stability criterions of three CRs and an exciter were presented in their simplified form. Section 4 would serve to give numerical results on these points. This section focuses on the influence of the structure parameters on the synchronization and stability of the proposed vibrating system. For the vibrating system, the parameters are: $m=360 \mathrm{kg}$, $r=0.05 \mathrm{m}$, $J_m=38.5 {\mathrm{kgm}}^2$, $k_x=254.5 \mathrm{kN}/\mathrm{m}$, $f_x=1.366 \mathrm{kN}\cdot \mathrm{s}/\mathrm{m}$, $k_y=203.2 \mathrm{kN}/\mathrm{m}$, $f_y=1.22 \mathrm{kN}\cdot \mathrm{s}/\mathrm{m}$, $k_{\psi }=183 \mathrm{kN}\cdot \mathrm{m}/\mathrm{rad}$, $f_{\psi }=0.375 \mathrm{kN}\cdot \mathrm{s}/\mathrm{rad}$; for the induction motor ($50 \mathrm{Hz}$, $380 \mathrm{V}$, $6-\mathrm{pole}$, $0.2 \mathrm{kW}$, rated speed $980 \mathrm{rpm}$), the parameters are: $R_s=40.4 \mathrm{\Omega }$, $L_s=1.212 \mathrm{H}$, $R_r=12.5 \mathrm{\Omega }$, $L_r=1.222 \mathrm{H}$, $L_m=1.116 \mathrm{H}$, $T_{\mathrm{eN}1}=2.01 \mathrm{N}\cdot \mathrm{m}$, $f_1=0.005$.

4.1. Parameters Influence on the Synchronization

The critical structural parameters are changed to investigate the influence on the synchronization characteristics. The physical meaning of $T_{\mathrm{e}01}$ can be seen in Equation (33), i.e., $T_{\mathrm{e}01}$ is the necessary electromagnetic torque of the vibrating system. $T_{\mathrm{e}01}$ is a limited function of the phase differences and must not be more than its rated value ($T_{\mathrm{eN}1}$) at any time. Here, $T_{\mathrm{e}01\max }$ is defined as the maximum value of $T_{\mathrm{e}01}$ when the values of the phase differences vary. At the stage of unstable state during the entire operational process of the system, the dynamic coupling torques are assumed to manage to perform the stable synchronization by regulating the phase differences. This can easily result in the wide fluctuations of the rotational velocities and the phase differences; consequently, $T_{\mathrm{e}01}$ fluctuates as well and even reaches its all time maximum value ($T_{\mathrm{e}01\max }$). To sum up, the inequation $T_{\mathrm{e}01\max }\le T_{\mathrm{eN}1}$ is the necessary condition of achieving synchronization with sufficient energy.

Figure 2a presents maximum values of the necessary electromagnetic torque for the given sets of friction coefficients under different values of $r_m$ when $r_R=0.005 \mathrm{m}$, $m_R=4 \mathrm{kg}$. The value of $r_m$ can be varied by changing $m_1$. It can be seen that the increasing $r_m$ brings the increasing required value of $T_{\mathrm{e}01\max }$. It also should be noted that the larger the sum of a set of friction coefficients, the nearer a curve of $T_{\mathrm{e}01\max }$ comes to the value of $T_{\mathrm{eN}1}$. When the sum of a set of friction coefficients is equal to another, their corresponding curves overlap with each other. By fixing an arbitrary friction coefficient (here, $f_{R1}=0.08$), Figure 2b exhibits a linear growth of $T_{\mathrm{e}01\max }$ with the increasing $f_{R2}$ and $f_{R3}$. It is worth noting that each of the three friction coefficients has the same effect on the variation of $T_{\mathrm{e}01\max }$. Accordingly, based on the necessary condition of achieving synchronization with sufficient energy, the chosen values of $r_m$ and $f_{Ri}$ should be relatively smaller to ensure that the inequation $T_{\mathrm{e}01\max }\le T_{\mathrm{eN}1}$ holds. Figure 2c shows a plot of $T_{\mathrm{e}01\max }$ with changes in ${\eta }_r$ and ${\eta }_m$. The values of ${\eta }_r$ and ${\eta }_m$ can be varied by changing $r_R$ and $m_R$, respectively. In Figure 2c, the parameters of the vibrating system are set to $m_1=6 \mathrm{kg}$ and $f_{R1}=f_{R2}=f_{R3}=0.08$. The area I indicates the regions where $T_{\mathrm{e}01\max }\le T_{\mathrm{eN}1}$, i.e., the regions where an exciter and three CRs can reach synchronization, whereas some other area denotes the regions of unsynchronized behavior.

Next, with the following set of physical parameters $m_R=4 \mathrm{kg}$, $r_R=0.005 \mathrm{m}$, the maximum values of dimensionless coupling torques ($T_{\mathrm{c}ij\max }$) and synchronization ability coefficients (${\zeta }_{ij}$) are calculated and plotted as the functions of $r_l$ for the different $f_{Ri}$ ($i=1,2,3$) and $r_m$ in Figure 3, Figure 4 and Figure 5. At first, consider the simplified case of the identical friction coefficients. Assume that $f_{R1}=f_{R2}=f_{R3}=0.08$. In Figure 3, it can be found that even with variations in $r_m$ and $r_l$, equations $T_{\mathrm{c}12\max }=T_{\mathrm{c}13\max }=T_{\mathrm{c}14\max }$, $T_{\mathrm{c}23\max }=T_{\mathrm{c}42\max }=T_{\mathrm{c}34\max }$, ${\zeta }_{12}={\zeta }_{13}={\zeta }_{14}$, ${\zeta }_{23}={\zeta }_{42}={\zeta }_{34}$ will constantly hold. The reason is that the system parameters are completely symmetrical in this situation. It is noticed that, increasing $r_m$ leads to the decreasing of all synchronization coefficients $T_{\mathrm{c}ij\max }$ and ${\zeta }_{ij}$ (see Figure 3, Figure 4 and Figure 5). For small $r_l$, the maximum value of a dimensionless coupling torque and synchronization ability coefficient between an exciter and a CR remain largely unchanged with variations in $r_l$; for large $r_l$, one can notice a rapid increase of these two synchronization parameters when the value of $r_l$ is increased. It must also be noted that the maximum value of a dimensionless coupling torque and synchronization ability coefficient between two CRs slightly decrease at first, but shortly increase, with the increasing $r_l$. Recall, an increased value of $T_{\mathrm{c}ij\max }$ results in an increase in the allowable residual torque difference between two coupled objects, Equations (29)–(31). That is to say, it is easier to satisfy the synchronization criterion and implement synchronization of an exciter and three CRs. Consequently, in this case, to make the system implement synchronization more easily and improve the ability of synchronization of the system, a smaller value of $r_m$ and a greater value of $r_l$ should be chosen.

Consider now the cases of different friction coefficients. As depicted in Figure 4 and Figure 5, the results show that the variation tendency of the synchronization parameters $T_{\mathrm{c}ij\max }$ and ${\zeta }_{ij}$ with changes in $r_l$ is about in line with that in the case of identical friction coefficients. From Figure 4 and Figure 5, it is worth noting that, for large $r_l$, the maximum values of dimensionless coupling torques and synchronization ability coefficients between two coupled objects are sensitive to replacement of the two coupled objects (for example, the switch to an exciter and CR 3 from an exciter and CR 1 or from CR 1 and CR 2 to CR 3 and CR 1) when $r_m$ is fixed. The reason behind this is the substitution of one friction coefficient for another. In detail, it is worthwhile to mention that, from a comparison between Figure 4 and Figure 5, one can realize that a larger value of the friction coefficient leads to a decrease in the maximum value of a dimensionless coupling torque and synchronization ability coefficient between an exciter and a corresponding CR; an increased difference value between two arbitrary friction coefficients results in a decrease in the maximum value of a dimensionless coupling torque and synchronization ability coefficient between two corresponding CRs.

Next, assume that the three variable friction coefficients are identical, i.e., $f_{R1}=f_{R2}=f_{R3}=f_R$. Figure 6 shows the relation of $T_{\mathrm{c}ij\max }$ and ${\zeta }_{ij}$ to friction coefficient $f_R$. As presented in Figure 6, it is noticeable that as $f_R$ is increased, the synchronization coefficients $T_{\mathrm{c}ij\max }$ ($ij=12,13,14$) and ${\zeta }_{ij}$ decrease gradually, and for $T_{\mathrm{c}ij\max }$ ($ij=23,42,34$), in contrast, there is a correlated rise. In a vibrating system, the residual torque difference between an exciter and a CR is always much higher as compared to that between two CRs. For purpose of making it easier to satisfy the synchronization criterion and implement synchronization, a larger value of $T_{\mathrm{c}ij\max }$ ($ij=12,13,14$) has more significant effect on the implementation of synchronization. Furthermore, the larger the value of ${\zeta }_{ij}$, the stronger the synchronization ability. Consequently, a small $f_R$ is needed in this case. In Figure 7, the contour plots of the maximum values of dimensionless coupling torques and the synchronization ability coefficients versus ${\eta }_r$ and ${\eta }_m$ are pictured. The system parameters are set to $m_1=6 \mathrm{kg}$ and $f_{R1}=f_{R2}=f_{R3}=0.08$. Looking at Figure 7, which gives insight into variation of the synchronization effect, one sees that there is a marked rise in the synchronization coefficients $T_{\mathrm{c}ij\max }$ and ${\zeta }_{ij}$, with the increasing ${\eta }_r$ and ${\eta }_m$. Therefore, the dimensionless parameters ${\eta }_r$ and ${\eta }_m$ are often preferably high to enhance synchronization ability and make it easy to implement synchronization, but should be reduced when the coordinate $\{{\eta }_r,{\eta }_m\}$ is not in the area I (can be seen in Figure 2) to satisfy the synchronization criterion.

In general, to have the best synchronization effect, an appropriate selection of each structural parameter must be implemented, and further understanding of the structural parameters influence on the stability is also needed.

4.2. Parameters Influence on the Stability

In this part, the parameters of the vibrating system are fixed at $m_1=6 \mathrm{kg}$, $m_R=4 \mathrm{kg}$, $r_R=0.005 \mathrm{m}$. Recall, the balance equations of an exciter and three CRs are nonlinear functions of the phase differences, Equations (18)–(21). Once the synchronization criterion and stability criterion are satisfied and the three friction coefficients are defined, the stable phase differences values can be calculated with variations in $r_l$ by adopting a numerical method. Next, the stability coefficients and clinging ability coefficients can be calculated based on their expressions and the stable phase differences values. The obtained numerical results are shown in Figure 8, Figure 9 and Figure 10.

In Figure 8, it is shown that the consistent emerging values of stable phase differences are ‘immune’ to changes in $r_l$. As depicted in Figure 8a, the constant phase differences are: $2{\alpha }_1={204}^{\circ }$ (or $-{156}^{\circ }$), 2α₂ = 0^∘, 2α₃ = 0^∘; in Figure 8b, the constant phase differences are: $2{\alpha }_1={206}^{\circ }$ (or −154^∘), 2α₂ = 344^∘ (or $-{16}^{\circ }$), $2{\alpha }_3={16}^{\circ }$; in Figure 8c, the constant phase differences are: $2{\alpha }_1={206}^{\circ }$ (or $-{154}^{\circ }$), $2{\alpha }_2={352}^{\circ }$ (or $-8^{\circ }$), $2{\alpha }_3={352}^{\circ }$ (or $-8^{\circ }$); in Figure 8d, the constant phase differences are: $2{\alpha }_1={192}^{\circ }$ (or $-{168}^{\circ }$), $2{\alpha }_2=0^{\circ }$, $2{\alpha }_3=0^{\circ }$. From this figure, it is evident that the three CRs synchronize in-phase for three identical friction coefficients and synchronize with a non-zero phase difference for three non-identical friction coefficients. Moreover, the absolute value of the stable phase difference (ranging from $-{180}^{\circ }$ to ${180}^{\circ }$) between an exciter and CR 1 declines as the friction coefficient corresponding to CR 1 is increased ($2{\alpha }_1$ remaining almost unchanged when fixing $f_{R1}$), whereas that between two arbitrary CRs grows as the difference value between two corresponding friction coefficients is increased.

According to Equation (45), only when $H_i>0$ ($i=1,2,3$) can the stability criterion of synchronous operation be satisfied, and as a result, a stable solution of the phase differences be obtained. Figure 9 shows the dependence of the stability coefficients as a function of the dimensionless parameter $r_l$. It is seen from Figure 9 that even with variations in $r_l$, all elements ($H_1$, $H_2$ and $H_3$) are positive, which accounts for the ongoing and ever-present values of stable phase differences during the entire process of varying $r_l$. Moreover, one of the noticeable features in Figure 9 is the significant growth in three stability coefficients when $r_l$ is increased, even as a set of friction coefficients is changed to another one. As mentioned before, the greater the stability coefficient $H_i$, the stronger the stability of synchronous operation. Therefore, a greater value of $r_l$ should be chosen to improve the stability of synchronous operation.

The stability coefficients in a synchronous state are largely influenced by $r_l$, whereas the clinging ability coefficients of three CRs are unaffected by $r_l$ and are determined by the three friction coefficients, as shown in Figure 10. It is observed that all the clinging ability coefficients are very tiny–far less than one with four different sets of friction coefficients. This indicates that three CRs will cling to the surface of cavities quite tenaciously when rotating around their centroids synchronously and steadily. Furthermore, it is worthwhile to mention that from a comparison between two different sets of three identical friction coefficients (Figure 10a,d) or among three non-identical friction coefficients, one can realize that a large $f_{Ri}$ makes for better clinging ability.

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\mathrm{kg}$, $m_R=4\mathrm{kg}$, $r_R=0.005\mathrm{m}$, $l_0=0.3\mathrm{m}$. Note that for the parameter values considered in Section 4 and Section 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 π/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\mathrm{s}$. At the steady-state, the rotational velocity of the synchronization for the exciter is $991.3\mathrm{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\mathrm{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\text{-}$ and $y\text{-}$directions are equal to each other, whereas the stable value in $\psi \text{-}$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.

Figure 11g, Figure 12g, Figure 13g and Figure 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\mathrm{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 (Figure 12a and Figure 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\mathrm{r}/\min$ and around a constant value of ${205}^{\circ }$ (partial synchronization), on account of the CR’s non-negligible gravity effect. Despite the non-identical friction coefficients, the rotational velocities of the synchronization for the three CRs (oscillating regularly around $991.7\mathrm{r}/\min$) are equal because of the regulating effect from the dynamic coupling torques. However, the difference in the friction coefficients makes the limit behavior of the coupled CRs transform from in-phase synchronization to partial synchronization (comparing Figure 12c,d with Figure 11c,d), i.e., the stable value of $2{\alpha }_2$ changes from $0^{\circ }$ to ${344.5}^{\circ }$, 2α₃ from 0^∘ to ${16.2}^{\circ }$. Note that there are small fluctuations in the stable phase differences among three CRs. The reason behind this phenomenon is that, as mentioned in Figure 4 and Figure 5, the mismatch in the friction coefficients can cause the difference in the synchronization abilities between an exciter and a CR, i.e., in the second simulation, the intensity level of being synchronized by the exciter through the rigid frame in CR 2 is higher than that in CR 1 or CR 3. In Figure 12e,f (Figure 13e,f), it is observed that apart from identical vibrations in $x\text{-}$ and $y\text{-}$directions, the vibration in $\psi \text{-}$direction is non-zero and obvious, which indicates a mixed form of motion consisting of a circular motion and a minor swing motion at the steady-state. Note that the swing movement fails to effectively transmit energy.

Figure 13 presents the results of computer simulation when $f_{R1}=0.08$, $f_{R2}=0.06$, $f_{R3}=0.04$. When the vibrating system moves into the synchronized state at $t=7 \mathrm{s}$, the relevant steady-state values are as follows: the synchronous rotation velocity nears 991.9 r/min (or fluctuating around it); the phase differences are $2{\alpha }_1={199}^{\circ }~{212}^{\circ }$ (average value ${205.5}^{\circ }$), $2{\alpha }_2={351}^{\circ }~{352}^{\circ }$ (average value ${351.5}^{\circ }$), and $2{\alpha }_3={352}^{\circ }~{353}^{\circ }$ (average value ${352.5}^{\circ }$). It should be stressed that the fluctuations of the stable values of $2{\alpha }_2$ and $2{\alpha }_3$ in Figure 13c,d are smaller than those from Figure 12c,d. 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 Figure 4c and Figure 5c.

Simulation studies (see Figure 11, Figure 12 and Figure 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 \mathrm{s}$ and have a synchronous rotation velocity of $992.6 \mathrm{r}/\min$ (stabilizing at or oscillating around it), and phase differences of $2{\alpha }_1={184}^{\circ }~{197}^{\circ }$ (average value 190.5^∘), $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. Stable values of the phase differences obtained from the numerical analyses and computer simulations.

Friction Coefficients	Stable Phase Difference	Numerical Results	Simulation Results
$\begin{array}{l}{f}_{R1}=0.08\\ f_{R2}=0.08\\ f_{R3}=0.08\end{array}$	$2{\alpha }_1$	${204}^{\circ }$	${204}^{\circ }$
	$2{\alpha }_2$	0∘	$0^{\circ }$
	$2{\alpha }_3$	$0^{\circ }$	$0^{\circ }$
$\begin{array}{l}{f}_{R1}=0.08\\ f_{R2}=0.04\\ f_{R3}=0.08\end{array}$	$2{\alpha }_1$	${206}^{\circ }$	${205}^{\circ }$
	$2{\alpha }_2$	${344}^{\circ }$	${344.5}^{\circ }$
	$2{\alpha }_3$	16∘	${16.2}^{\circ }$
$\begin{array}{l}{f}_{R1}=0.08\\ f_{R2}=0.06\\ f_{R3}=0.04\end{array}$	$2{\alpha }_1$	${206}^{\circ }$	${205.5}^{\circ }$
	$2{\alpha }_2$	${352}^{\circ }$	${351.5}^{\circ }$
	$2{\alpha }_3$	${352}^{\circ }$	${352.5}^{\circ }$
$\begin{array}{l}{f}_{R1}=0.04\\ f_{R2}=0.04\\ f_{R3}=0.04\end{array}$	$2{\alpha }_1$	${192}^{\circ }$	${190.5}^{\circ }$
	$2{\alpha }_2$	$0^{\circ }$	$0^{\circ }$
	$2{\alpha }_3$	$0^{\circ }$	$0^{\circ }$

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.

Finally, to examine the theoretical and numerical results further, the study is extended to the case when $\left\{{\eta }_r,{\eta }_m\right\}$ does not locate in the regions where $T_{\mathrm{e}01\max }\le T_{\mathrm{eN}1}$ (i.e., the area I, can be seen in Figure 2) and the synchronization criterion cannot be satisfied. The parameter values and dry friction coefficients are the same as in the first simulation, except 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_{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_{DC1}$ arrives at a positive range and the dynamic coupling torques $T_{DC2}$, $T_{DC3}$, and $T_{DC4}$ fall into the negative range. Hence, $T_{DC1}$ acts as the coupled load torque of the exciter and $T_{DC2}$, $T_{DC3}$, and $T_{DC4}$ play the part of the coupled driving torques for the three CRs (see Equations (49)–(52)). In this way, with the help of the dynamic coupling torques, the energy can be transferred from the exciter with a power supply to the three CRs without the direct driving sources. Additionally, with the increase of the sum of three friction coefficients, the stable value of $T_{DC1}$ has increased, and the sum of the absolute values of $T_{DC2}$, $T_{DC3}$, and $T_{DC4}$ is no exception. It indicates that the sum of three friction coefficients equals to a larger value, and the three CRs will subsequently demand more energy from the exciter by virtue of the dynamic coupling torques, hence, the synchronous rotation velocity gets smaller. When time reaches 11 s, an interference of $\pi /4$ angular displacement is applied to CR 1 and the system comes into a self-adjusting stage. To cope with the sudden increase in CR 1’s phase, the $T_{DC2}$ (CR 1’s dynamic coupling torque) immediately changes from negative sign to positive sign, and thus serves as the load torque to decelerate CR 1’s rotational velocity rapidly. In this period, the other dynamic coupling torques make adjustments accordingly. As a result, the phase differences tend to go back to their previously stable values. A short time of readjustment and fluctuation later, the newly synchronous and stable states are in accordance with the originally synchronous and stable states. During this transient process, the dynamic coupling torques play a vital role in helping the vibrating system recover to its originally synchronous and stable state.

The time series corresponding to the necessary electromagnetic torques in the simulations are shown in Figure 17. Clearly, the greater the sum of the three friction coefficients, the larger the maximum value and stable value of $T_{\mathrm{e}01}$. The obtained simulation results coincide with the numerical analyses in Figure 2a. From Figure 17e, it can be seen that the necessary electromagnetic torque at the steady-state exceeds the rated electromagnetic torque, and accounts for the failure of synchronous operation among an exciter and three CRs in the corresponding simulation. Without an adequate and strong supply of energy from the induction motor, the aforementioned synchronous motion cannot be initiated and maintained.

6. Conclusions

In the present work, we have focused on the study of synchronization in an exciter and three CRs with non-identical friction coefficients. It has been demonstrated that the three coupled CRs are seen to synchronize in-phase for friction coefficient matched conditions, and show partial synchronization when there is a small difference in their friction coefficients.

There are two essential components of synchronization criterion (refer to Equations (29)–(31) and (33)). The two components play an analogous role in calculating the numerical values of phase differences, stability coefficients, and clinging ability coefficients at the steady-state. Since the CRs are not equipped with direct driving sources, they obtain energy from the induction motor through the vibrating rigid platform (with the help of dynamic coupling torques). Hence, the induction motor offers energy not only to the exciter but also to the CRs, for the purpose of reaching synchronous rotation of the exciter and CRs. The second component of the synchronization criterion (Equation (33)) takes into consideration the total energy balance and introduces the concept of the necessary electromagnetic torque; as a result, it can be used to determine whether the induction motor can offer adequate energy to achieve synchronization.

Interestingly, in the simulation cases, for an appropriate value of $T_{\mathrm{e}01\max }$ ($T_{\mathrm{e}01\max }\le T_{\mathrm{eN}1}$), the synchronous motion is robust against mismatches in friction coefficients. However, for a large value of $T_{\mathrm{e}01\max }$ ($T_{\mathrm{e}01\max }>T_{\mathrm{eN}1}$), the synchronous swing of three CRs seems to be sensitive to differences in friction coefficients. In this interval of $T_{\mathrm{e}01\max }$ three CRs can swing synchronously only when the friction coefficients for each CR are identical.

Results presented in this paper, although numerical, indicate how important the proper choice of system parameters and coupling structure’s properties is in order to obtain desired, regular behavior in real, practical implementations of engineering problems. The friction coefficient $f_{Ri}$ and dimensionless variables $r_m$, $r_l$, ${\eta }_r$, and ${\eta }_m$ can all affect the synchronization and stability and thus can be used to tune the synchronization performance to the ideal case conveniently at the design and optimization stage of vibrating machines. The combination of all these key parameters will eventually lead to a specific type of synchronous motion.

The proposed vibrating system can produce an ideal synchronization regime that is ever-present and stable under variation of the structural parameter $r_l$, leading to a small-sized and compact structure very easy to design and implement in various engineering practices.

The dynamic coupling torques act as a self-adjusting synchronization controller, which can automatically maintain the original rotational velocities and phase differences and, thus, the obtained simulation results show robust motion synchronization in the exciter-CR vibrating system against the disturbance.

The results of this study demonstrate an innovative solution for designing a type of desired vibrating crusher in various engineering practices by employing this synchronization mechanism. Crushing devices can hold materials between a CR’s surface and an inner wall’s surface. During the process of synchronization, a CR rotates along an inner wall and applies sufficient force to bring the surfaces together to generate enough energy within the materials being crushed, so that its molecules separate from each other. Without a doubt, the materials needing breaking will influence the synchronization between CRs and an exciter. This study can be further improved by addressing the issue of materials’ effects on the synchronous and stable operation of the proposed vibrating system, and this issue will be investigated in future work.