Synchronization and Stability of a Nonlinear Vibrating Mechanical System Characterized by Asymmetrical Piecewise Linearity

In previous studies about the synchronization of vibrators, the restoring forces of springs are mainly treated as linear directly, whereas the nonlinear features are rarely considered in vibrating systems. To make up this drawback, a dynamical model of a nonlinear vibrating mechanical system with double rigid frames (RFs), driven by two vibrators, is proposed to explore the synchronization and stability of the system. In this paper, the nonlinearity is reflected in nonlinear restoring forces of springs characterized by asymmetrical piecewise linear, where the nonlinear stiffness of springs is linearized equivalently using the asymptotic method. Based on the average method and Hamilton’s principle, the theory conditions to achieve synchronization and stability of two vibrators are deduced. After the theory analyses, some numerical qualitative analyses are given to reveal the coupling dynamical characteristics of the system and the relative motion properties between two RFs. Besides, some experiments are carried out to examine the validity of the theoretical results and the correctness of the numerical analyses results. Based on the comparisons of the theory, numeric and experiment, the ideal working regions of the system are suggested. Based on the present work, some new types of vibrating equipment, such as vibrating discharging centrifugal dehydrators/conveyers/screens, can be designed.


Introduction
When it comes to vibration, in most cases, it is harmful and can damage the machine sometimes. Hence, some research works [1][2][3] stressed the theme of vibration control and suppression enough. In engineering, however, some machines are designed and operated by utilizing vibrations, such as vibrating feeders, vibrating conveyors, probability screen, most of which are based on the selfsynchronization theory of vibrators (unbalanced rotors driven by induction motors separately).
As a typical nonlinear phenomenon, synchronization expresses great potential in plenty of application areas.
Refs. [4] and [5] researched the synchronization of robot manipulators and the predictive control of the permanent magnet synchronous motor, respectively. Huygens [6] and Karmazyn et al. [7] studied the synchronization of the pendulums attached to different structures. As a complex dynamical behavior, the synchronization of chaoticsystem has become the focus of chaos research [8,9].
For the synchronization of vibrators, the first definition, as well as the early explanation for the synchronization mechanism of two identical vibrators, can be found in Refs. [10,11]. After that, many scholars studied the synchronization theory of vibrators. The theory of multiple cycle synchronization of the vibrating mechanical system driven by vibrators was given in Refs. [12][13][14]. Wen et al. [15] further expanded the synchronization theory and applied it to engineering successfully. The composite synchronous motion of vibrators was analyzed in Refs. [16,17]. Hou et al. [18] studied the influences of the structural parameters on synchronous characteristics between the rotors in an anti-resonance system. Dudkowski et al. [19,20] investigated the synchronization behavior and stability of classical lightly supported pendula systems. To maintain stability synchronization of vibrators in the multi-motor-pendulum vibration system, a combined synchronous control strategy was proposed in Ref. [21]. The authors in Ref. [22] investigated the synchronization of two vibrators with common rotational axis. Besides, the influence of non-ideal power source and the particular phenomena around the resonant point of the system were also investigated deeply, for example, some short comments on the synchronization of two or four nonideal vibrators and the Sommerfeld effect were presented in Refs. [23][24][25]. Du et al. [26] considered synchronous characteristics of a non-resonant system.
Recently, some novel researches on synchronization problems of vibrators were presented. For instance, by the theoretical and simulation studies, Refs. [27,28] revealed the synchronization mechanism and stability of the synchronous states in a vibrating system with two homodromy vibrators in detail. From the perspective of electromechanical coupling, Zhang et al. [29] gave the explanation for the synchronous phenomenon and stability of two vibrators. Besides, Chen et al. [30] investigated the stability and coupling dynamical characteristics in a vibrating system with a single RF, which is driven by four vibrators.
Several methods have been used to study synchronization of vibrators. The first method, known as the direct separation motion method [10,11], has been shown to be effective. Especially in dealing with the synchronization problems of two identical vibrators, this method exhibited good effects in simplifying the investigation by combining the differential equations of the two rotors into the differential equations of the phase difference. Taking the influences of the damping into account, Wen et al. [15] further developed the synchronization theory by means of the average method.
The aforementioned methods have been widely applied to solve the synchronization problems of vibrators. However, according to the authors' knowledge, in the previous literatures, the restoring forces of springs of the vibrating system were mainly treated as linear directly. This is because that when the nonlinear features of the system are taken into account, investigation of the synchronization and stability of the system is generally difficult. In engineering, however, we should emphasize that many vibrating systems need to utilize its nonlinear effect, i.e., the damping and the restoring forces of springs may be nonlinear in most cases. To solve the synchronization problems with nonlinear characteristics, the key problem is to deduce the motion differential equations of the system, i.e., how to introduce nonlinear features into the differential equations and make the system be manageable by the numerical or simulation methods. According to Ref. [15], the main investigation methods to deal with the nonlinear problems include the asymptotic methods, the method of small parameters, the equivalent linearization method, and so on.
To make up the drawbacks of previous works, in this paper, a dynamical model with nonlinear springs characterized by asymmetrical piecewise linearity, driven by two vibrators, is proposed to explore the synchronization and stability of the corresponding system by theory, numeric and experiment. Our goal is to introduce the nonlinear features into the vibrating system to make the investigate approach more accurate, as well as reveal the dynamical characteristics with nonlinear restoring forces. The approach proposed in this paper is based on the traditional average method but with the nonlinear characteristics, where the nonlinear features are reflected in the nonlinear stiffness of springs obtained by the asymptotic method. The present work, to a certain extent, is an extension and deep investigation of the previous literatures, which can provide a guidance for designing some new types of vibrating machines corresponding to the considered model.
In this paper, according to the ratio between the operating frequency to the natural frequencies, denoted by z, we generally divide the resonant regions of a vibrating system into four types: sub-resonant (z<0.9), near subresonant (0.9<z<1), near super-resonant (1<z<1.1), and super-resonant (z>1.1).
The rest of the paper is organized as follows: In Section 2, the description of the nonlinear vibrating mechanical system is presented, accompanied by deducing the motion differential equations, and the equivalent stiffness is derived. In Section 3, the theoretical conditions to implement synchronization and ensure its stability are deduced by the average method and Hamilton's principle, respectively. Section 4 is devoted to numerical qualitative analyses, followed by experimental verifications given in Section 5. Finally, conclusions are drawn in Section 6.
In this paper: 'RF' is the abbreviation of 'rigid frame' , and 'SPD' is the abbreviation of 'stable phase difference' .

Description of the System
The considered nonlinear vibrating mechanical system, illustrated in Figure 1, consists of two RFs (i.e., m 1 and m 2 ) linked by the shear rubber springs and the gapactivated compression rubber spring (see k 1y and Δk 1y in Figure 1, respectively), where the former is symmetrically installed, while the latter is asymmetrically arranged on a single side with the average gap μ.
Here m 1 and m 2 are the main vibrating and isolative RFs, respectively, driven by a pair of identical vibrators mounted on both sides of the isolative RF 2. As illustrated in Figure 1, two vibrators can rotate with the same or opposite directions and generate certain exciting forces. The centroid of the RF 1 coincides with that of RF 2, which is also the equilibrium point of the system, marked as o, and the coordinate system oxy is established with o as the origin.
It should be pointed that the shear rubber springs are constrained to generate stiffness only in the horizontal direction and thus the relative motion between two RFs is restricted in y-direction. The motion of the system is considered as the plane motion, and the system has four degrees of freedom: the displacements of two RFs in x-, y-, and ψ-directions, denoted by x, y 1 , y 2 , and ψ, respectively. The angular positions of both vibrators are denoted by φ 01 and φ 02 .
To derive the differential equations of motion for the system, the nonlinear damping force F 1 (y 1 , y 2 ,ẏ 1 ,ẏ 2 ) and the nonlinear restoring force F 2 (y 1 , y 2 ) are firstly given by Eq. (1): where f 1y and Δf 1y are the damping constant of relative motion between two RFs for the shear rubber spring and gap-activated compression rubber spring, respectively; k 1y and Δk 1y are the stiffness of the shear rubber spring (1) and gap-activated compression rubber spring, respectively; μ is the gap between the gap-activated compression rubber spring Δk 1y and RF 1.

Absolute Motion Differential Equation
The equivalent stiffness and the equivalent damping constant of the relative motion between two RFs in y-direction are assumed to be k ′ 1y and f ′ 1y , respectively. The absolute motion differential equations for the considered system, formulated by using the Lagrange's principle, are as follows: where M 1 is the mass of RF 1, while M 2 is that of RF 2 (including two vibrators); M is the mass of the total vibrating system; m 01 and m 02 are the mass of the vibrator 1 and 2, respectively; m 0 is the mass of the standard vibrator; r is the eccentric radius of each vibrator; J 0i is the moment of inertia of the vibrator i (i=1, 2); J is the moment of inertia of the total vibrating system about its mass center; l 0 is the distance between the rotation center of each vibrator and the mass center of the total vibrating system; l e is the equivalent rotary radius of the total vibrating system about its mass center; f 0i is the damping coefficient of axis of the induction motor i (i=1, 2); f x , f 2y , and f ψ are the damping constant of the isolative RF in x-, y-, and ψ-directions, while k x , k 2y , and k ψ are the stiffness parameters; β i is the angle between the line from the rotation center of the vibrator i to the mass center of the total vibrating system and x-axis; T ei is electromagnetic torque of the motor i.

Relative Motion Differential Equation
Since the vibrating system considered in Figure 1 consists of two RFs, it is necessary to derive the relative motion differential equation of the system in y-direction.
To facilitate the study, some parameters and omissions are set as follows: (i) The average phase and the phase difference of the two vibrators are set as φ and 2α, respectively, i.e., φ 01 =φ+α, φ 02 =φ−α; (ii) When the system is operated synchronously, the average angular velocity of the two vibrators is denoted by ω m0 ; (iii) Based on Ref. [15], in the steady state, the relationships between acceleration and displacement satisfy the fact of ÿ i = −ω 2 m0 y i (i = 1, 2) , and the coil isolative springs are relatively very soft, generally satisfy the fact of k x (or k 2y ) << k 1y in engineering; (iv) φ 01 and φ 02 in the second formula of Eq. (2) can be neglected when the system is operated synchronously.
Based on the above parameter setting and omissions, the formulae with respect to y 1 and y 2 in Eq. (2), can be rearranged as follows: To obtain the solution of the relative motion between two RFs in y-direction, by the procedure of Eq.
, the motion differential equation described by the relative displacement, relative velocity, and relative acceleration, are obtained, see Eq. (5): where m is called the induced mass of the vibrating system.

Relative Motion Responses
Based on Eq. (5), the natural frequency of the relative motion between two RFs in y-direction, can be deduced as Eq. (6): Besides, based on the method of solving responses in Ref. [31], the relative motion response, i.e., the solution of y 12 , can also be obtained, see Eq. (7): , where F 12 is the coefficient that the single harmonic excitation contributes to the effective vibration amplitude of the relative motion between two RFs; γ 12 is the relative phase lag angle; z 0 is the ratio between the operating frequency and the natural frequency; ξ ′ 1y is the equivalent critical damping ratio of the relative motion between two RFs in y-direction.
According to Eq. (7), the relative vibration amplitude, denoted by δ 12 here, is calculated by δ 12 = 2|F 12 | cos α , where δ 12 is the absolute value of the sum of amplitude vectors.

Absolute Motion Response
Using the transfer function method and the superposition theorem [31], the absolute responses in x-, y-, and ψ-directions are presented by: with

Derivation of the Equivalent Stiffness Using the Asymptotic Method
Since the damping of vibrating machines is relatively very small, the nonlinear damping force F 1 (y 1 , y 2 ,ẏ 1 ,ẏ 2 ) in Eq.
Based on the asymptotic method [15], the first approximate solution of Eq. (5) is obtained as: with where θ is the phase lag angle of the response for the relative motion.

Theory Condition of Realizing Synchronization of Two Vibrators
To obtain the synchronization criterion of the two vibrators, based on the chain rule, the responses x, y 2 , and ψ listed in Eq. (8) are differentiated for time t. The obtained solutions ẍ , ÿ 2 , and ψ are inserted into the last two formulae in Eq. (2), followed by integrating them over φ = 0-2π and rearranging it, then we can obtain the average balanced equations of two vibrators as follows: 2) represents the electromagnetic torque of an induction motor operating steadily with ω m0 ; T u is the kinetic energy of the standard vibrator.
In the above integration process, the change of 2α with respect to time t is very small compared with that of φ, which can be considered as a slow-changing parameter [10], denoted by its integration middle value 2α.
The dimensionless residual torques of induction motors 1 and 2, denoted by T R1 and T R2 , respectively, are presented by Eq. (15): Rearranging Eqs. (13) and (14) by the procedures of addition and subtraction, respectively, the following expressions are yielded: T where T Load is the load torque of the vibrating system acting on the two motors; T Difference denotes the difference between the dimensionless residual torques of the two motors; and T Capture is called as the torque of frequency capture.
Rearranging Eq. (17), then we have: where D a is defined as the synchronization index of the system. According to Eq. (18), we can learn that only the absolute value of D a be equal or greater than 1, can Eqs. (13) and (14) be solved. Under the precondition of |D a | ≥ 1 , therefore, the theory condition to achieve synchronization can be obtained, see Eq. (20), which requires that the torque of frequency capture be greater than or equal to the absolute value of the difference of the dimensionless residual torques of two motors, i.e., Based on Eqs. (18) and (19), when the value of D a is fixed, there exist two corresponding 2α . Taking the condition of D a = 2 , for example, we can deduce the facts of 2α + θ c = π/6 or 2α + θ c = 5π/6 in this case. Besides, when two motors are completely identical, we have T Difference = 0 while D a → ∞ , which leads to that 2α + θ c = 0 or 2α + θ c = π.
According to Eqs. (13)- (20), the theory condition to achieve synchronization and the final phase difference of the synchronous states are greatly influenced by the rotation directions and the mounting angles of two motors. The necessary condition for achieving the frequency capture and reach the synchronous operation of two vibrators is that the torque of frequency capture should be large enough to overcome the difference in the residual torques of the two motors. The system can implement θ c = arctan(D/C), D/C > 0, π + arctan(D/C), D/C < 0, synchronization operation, due to the contribution of frequency capture torque.

Stability of the Synchronous States
As mentioned in Section 3.1, generally there exist two solutions for 2α , here we should point out that, one of these two solutions is stable, while the other is not. To reveal the stability of the solutions, it is necessary to further deduce the stability criterion of the system.
The stability criterion can be obtained by using Hamilton's principle, where the kinetic energy T and potential energy V are given respectively by: The expression of Hamilton's average action amplitude over one periodic, denoted by I, is presented by: From Ref. [15], the stable synchronous solution of 2α are supposed to correspond to be the minimum of Hamilton's average action amplitude, i.e., the second-order derivative of I should be greater than zero, see Eq. (24): Substituting Eq. (23) into Eq. (24), the stability criterion of the system is given by: with When the rotation directions and mounting angles of two motors are fixed, the stability criterion described in Eq. (25) can be further simplified. For example, for two co-rotating vibrators with mounting angles β 1 = π and β 2 = 0, the stability criterion can be simplified as: where H is defined as the coefficient of the ability of stability of the system. According to Eq. (26), for two corotating vibrators with mounting angles β 1 =π and β 2 =0, the stability criterion is described as that the product between the stability coefficient and the cosine of the phase difference should be greater than zero. According to Eq. (26), we can note that H>0 are in compliance with the stable phase difference (SPD) 2α * 0 ∈ (−π/2, π/2) , , while 2α * 0 ∈ (π/2, 3π/2) holds for the fact of H<0.
It should be noted that the above stability criterion is derived by the principle of minimum of Hamilton's average action amplitude, and we here only discuss the stability of synchronous states for two vibrators, so the type of the stability belongs to the local stability.

Parameters Setting of the Nonlinear Vibrating Mechanical System
Based on the above theoretical analyses, in this section, taking the condition of β 1 = π, β 2 = 0 for example, we numerically give some results on coupling dynamical characteristics of the nonlinear vibrating system, including the stability coefficient, the SPD, and the relative vibration amplitudes in y-direction, as shown in the following For the parameters of rubber springs, including dynamical stiffness k 1y and static stiffness Δk 1y , the former can only be tested in experimental system operating steadily, while the latter is generally obtained using the tensile test bed. According to our experiences in engineering, the dynamical stiffness of each rubber spring can be treated as linear under the relative small deformation interval. Hence, the nonlinear stiffness of rubber springs in this paper, can be treated with piecewise linear. The stiffness of the shear rubber spring and the gap-activated compression rubber spring in present work are set as: k 1y = 9800 kN/m, and Δk 1y = 3300 kN/m.

Numerical Qualitative Results on Coupling Dynamical Characteristics of the System
Based on the stability criterion and considering both two counter-rotating or co-rotating vibrators, the stability coefficient (i.e., H ), as well as the SPD, are plotted in Figures 2, 3, 4. In Figure 2, for the case of two counterrotating vibrators, the coupling dynamical characteristics of the system can be explained according to two courses concerning ω m0 , i.e., the ω m0 increasing course and the ω m0 decreasing course, which are described by the solid arrows and the dotted arrows, respectively. As illustrated in Figure 2, the H-ω m0 plane is divided into two regions: regions I and II, where in region I, H>0 complies with the fact of -π/2<2α<π/2, while in region II, H<0 holds for the fact of π/2<2α<3π/2. Particularly, for two identical vibrators considered in present work, the phase difference is stabilized in the vicinity of zero or π, see Figures 3 and 4.
For the increasing course of ω m0 , in region I, the value of H continues to increase along with the increment of ω m0 until a sharp decrease as the curve approaches point D (i.e., the resonant point of ω 0 ≈279.3 rad/s). Moreover, the curve occurs a break around point D (jump down from points D to D ′ ), where the value of H changes from positive to negative, and in the meantime the corresponding SPD jumps from zero to π for two identical vibrators, see D to D ′ in Figure 3. After the point D ′ , the value of H increases again along the curve. When it approaches point F, however, the curve changes from region II to region I and the value of H is greater than zero again. In this condition, the curve of SPD changes from G ′ to G in Figure 3. For the ω m0 decreasing course, the trend of the curve of H is plotted with the dotted arrow. Similarly, there exists a breakpoint from points E ′ to E, where the SPD jump down from π to zero simultaneously. Different from Figure 3, when two vibrators rotate in the same direction, there only exits one breakpoint during the ω m0 increasing course (see D to D ′ in the solid arrows route) and the ω m0 decreasing course (see E ′ to E in the dotted arrow route) according to Figure 4.
The characteristics of the stability and SPD, have a direct influence on the motion characteristics of the system. By numerical methods, the changing curves between the relative vibration amplitude in y-direction (i.e., δ 12 ) and the operating frequency ω m0 is plotted in Figure 5. Since in region I the stability coefficient is positive and the corresponding SPD is around zero for two identical vibrators, the positive superposition of the exciting forces of the two vibrators, as well as the relative linear motion between two RFs in y-direction can be achieved in this case. As shown in Figure 5, in the subresonant and near sub-resonant state for ω 0 (the resonant point ω 0 ≈279.3 rad/s is around point D in Figure 5), the vibration amplitude of the relative motion is relatively large, which is consistent with the requirements in engineering. Besides, the corresponding breakpoints are also observed. On the one hand, for both the two counterrotating vibrators and two co-rotating vibrators, there exists a breakpoint around D-D ′ for the ω m0 increasing course and that of E ′ -E for the ω m0 decreasing course. On the other hand, when the value of ω m0 is around points G ′ or G, we can observe an additional breakpoint when two vibrators rotate in the opposite direction.
Besides, as can be seen in Figure 5, although the relative vibration amplitude is greater than zero in the right of point G for two counter-rotating vibrators, the value of δ 12 is not large enough for engineering applications in this curve interval.
Based on the above numerical qualitative analyses, the ideal working points of the vibrating machines corresponding to the considered model, should be selected in the sub-resonant or near sub-resonant states of ω 0 . Especially in the latter condition, the phase difference between two identical vibrators is stabilized in the vicinity of zero, the vibration amplitude of the relative motion between the two RFs in y-direction, is relatively larger. In this case, the energy is saved due to the resonant effect and the relative ideal stability of the system is obtained, see points A and B in Figures 2, 3, 4, 5, these features are just the desires in engineering. It should be pointed out that the points A, B, and C correspond to the experimental verification points in Section 5, where the former two points are located in the ideal working region, while the characteristics of point C is useless in engineering.

Experimental Verifications
From the aforementioned theoretical and numerical qualitative analyses, the coupling characteristics of the system are revealed, and the ideal working region has been given. Based on the synchronization condition and stability criteria of the system obtained from the theoretical analysis, the numerical qualitative analyses of the system is carried out. In order to verify the effectiveness of the theoretical methods and the correctness of the numerical results, two sets of experiments are implemented in this section, one of which is selected in the ideal working region mentioned in numerical analyses, whereas the other is the opposite. It should be pointed

Experimental Illustration
The experimental system corresponding to the dynamical model plotted in Figure 1 is shown in Figure 6, where two identical vibrating motors (i.e., vibrators) are installed on both sides of the isolative RF 2 (m 2 ) symmetrically. The eccentric mass of each vibrator is adjusted to m 0 = 0.9 kg and the parameters of the experimental system are the same as those in numerical analyses. Using Tri-axial acceleration sensors and Hall-sensors, the vibration responses of the system in x-, y-, and ψ-directions, and the phases of the two vibrating motors can be measured, respectively. The corresponding experimental results are shown in the following sub-sections. Since two vibrators are selected as identical, we can observe that the system reaches the synchronous and stable operation quickly after the starting procedure.

Experimental Results in the
As shown in Figure 7a and b, the synchronous rotational velocity for two counter-rotating and co-rotating vibrators are around 2400.1 r/min (ω m0 = 251.2 rad/s) and 2395.2 r/min (ω m0 = 250.7 rad/s), respectively, while the SPD is around −41°. These results are roughly in agreement with the numerical analyses of points A (ω m0 = 251.2) and B (ω m0 = 250.7) in Figures 3 and 4. Under the condition of this phase difference, the system exhibits the relative linear motion between two RFs and the vibration amplitude is large enough to be useful in engineering, see Figure 7c-f. Moreover, from Figure 7c, the phase relationship between y 1 and y 2 is almost the opposite phases, the value of the relative vibration amplitude (y 1 −y 2 = δ 12 ) for two counter-rotating and co-rotating   Figure 5. Based on the above analyses, the experimental results illustrated in Figure 7a-f can be well matched with numerical qualitative results in Figures 2, 3, 4, 5.

Experimental Results in the Near Super-resonant State
for ω 0 To further examine the validity of the used analytic method and the obtained numerical results, it is also necessary to give the additional experimental verification for the other resonant region. As shown in Figure 8, two co-rotating vibrators are operated in the near super-resonant state for ω 0 , under the power supply frequency of 48 Hz. From Figure 8a and b, in the steady state, the synchronous rotational velocity is in the vicinity of 2850.4 r/min (ω m0 = 2 98.3 rad/s) and the phase difference is stabilized around 160°, which are approximately consistent with the point C (ω m0 = 298.3) in Figure 4. With this SPD, the vibration responses of the system are approximately zero, see Figure 8c-f. As presented in the theoretical and numerical analyses, in this case, the exciting forces of two vibrators are almost completely canceled, which leaves the system almost at rest or without vibration.
In the meantime, since the system is operated in the near super-resonant state in this condition, the obvious resonant phenomena is witnessed, see Figure 8c-f.
For the characteristics of the relative motion between two RFs, as shown in Figure 8(d), the vibration amplitude of (y 1 −y 2 ) is about 0.5 mm, while the numerical result is 0 (see point C in Figure 5). The reason for this slight distinction between numeric and experiment is that the SPD is not exactly 180° but nears it, this is because that the exciting forces of two vibrators cannot be canceled with each other completely, which leads to the fact that the relative linear motion with small vibration amplitude is possible. Therefore, we may conclude that the experimental results can roughly agree with what is shown in the above theoretical and numerical ones.
From the above Figure 8b (or Figure 7b), one can see that the phase difference in time domain analysis is not a fixed value, this is called as the uncertainty of phase difference, which is generally caused by the following reasons: (i) Although the two driving motors are completely identical, their actual electromagnetic output torque are not completely identical due to the difference of motor characteristics caused by processing and manufacturing; (ii) The fluctuations of supply current (or voltage) may influence the output characteristics of motors; (iii) The random fluctuations of loads may also leads to the uncertainty of phase difference, etc.
Additionally, the practical functions of vibrating machines are mainly depend on the stability of system, and such stability is mainly referred to as the stability of phase difference between the two vibrators. Although there are some small fluctuations (i.e., uncertainty) in time domain for the phase difference in the steady state, generally in engineering these fluctuations are so small that they can not influence the practical function of vibrating machine according to our engineering experiences. So in practical utilization of the analytical results, such stability of synchronous state can guarantee the normal technological treatment process in engineering.

Conclusions
(1) The mathematical modeling for a nonlinear vibrating mechanical system with double RFs driven by two identical vibrators is constructed, where the nonlinear stiffness of springs is treated by the process of equivalent linearization using the asymptotic method. (2) The criterions for implementing synchronization and ensuring its stability are deduced. According to the synchronization criterion, the torque of frequency capture should be greater than or equal to the absolute value of the difference of the dimensionless residual torques of two motors. For the stability characteristics of the synchronous state, taking two co-rotating vibrators with the mounting angles of β 1 = π and β 2 = 0 for example, the stability criterion requires that the product between the stability coefficient and the cosine of the phase difference be greater than zero. (3) Based on the theory, numeric and experiment analyses, the ideal working region of vibrating machines similar to the present model should be selected in the sub-resonant or near sub-resonant states for ω 0 . In this condition, the SPD between two vibrators is close to zero and the positive superposition of the exciting forces can be achieved, which satisfies the requirements in engineering. Especially in the near sub-resonant state for ω 0 , the vibration amplitude of the relative motion between two RFs is relatively larger, and the energy is saved due to the resonant effect. (4) The present work can lay a theoretical foundation in designing some new vibrating equipment with double RFs driven by two vibrators, such as vibrating discharging centrifugal dehydrators, vibrating con-