An efficient numerical procedure for calculating periodic vibrations of elastic mechanisms

Vietnam Journal of Mechanics, VAST, Vol. 38, No. 1 (2016), pp. 15 – 25 DOI:10.15625/0866-7136/38/1/5928 AN EFFICIENT NUMERICAL PROCEDURE FOR CALCULATING PERIODIC VIBRATIONS OF ELASTIC MECHANISMS Nguyen Van Khang, Nguyen Phong Dien∗, Nguyen Sy Nam 1Hanoi University of Science and Technology, Vietnam ∗E-mail: dien.nguyenphong@hust.edu.vn Received February 23, 2015 Abstract. This paper proposes a numerical procedure based on the well-known Newmark integration method to determine initial conditions for the periodic solution of a system of linear differential equations with time-periodic coefficients. Based on this, steady-state pe- riodic vibrations of mechanisms with elastic elements governed by linearized differential equations with time-periodic coefficients can be conveniently calculated. The proposed procedure is demonstrated by a dynamic model of a planar four-bar mechanism with the flexible coupler. Keywords: Steady-state vibration, elastic mechanism, Newmark integration method, mode superposition method, dynamic stability. 1. INTRODUCTION In high-speed machines, the motion of the transmission mechanisms is often com- posed of a combination of rigid body motion and elastic deformation [1, 2]. A review on the vibration and stability behavior of mechanisms with elastic links represents an update to earlier literatures surveys on this subject [3–5]. Many researchers have tried to represent the vibration of such mechanisms in a more and more realistic form. Up to now, the following models are used for modeling of flexible links of mechanisms: continuum models [6, 7], lumped parameter models [8, 9], finite element models [10–14]. In general, the mathematical formulation of this vibration problem is quite a com- plicated nonlinear differential equation, for which an exact solution is practically impos- sible. It is possible to calculate the transient solutions by the numerical methods. The linearized equations of motion of an elastic mechanism that performs the steady-state motion can then be expressed approximately by a set of n linear differential equations having time-periodic coefficients M(t)q¨+C(t)q˙+K(t)q = f(t), (1) c© 2016 Vietnam Academy of Science and Technology 16 Nguyen Van Khang, Nguyen Phong Dien, Nguyen Sy Nam where n× n matrices M(t), C(t), K(t) and excitation force vector f(t) in Eq. (2) are time- periodic with the least period T [8–12]. For stability analysis, the homogeneous differen- tial equation of Eq. (1) is considered M(t)q¨+C(t)q˙+K(t)q = 0. (2) This equation can also be represented in the form of Hill’s or Mathieu’s type equa- tion as mentioned in [13, 15, 16]. Because the periodic vibrations are a commonly ob- served phenomenon of mechanisms in the steady-state motion, a number of methods and algorithms were developed for calculating periodic vibrations and dynamic stability analysis [15–29]. Periodic solutions of Eq. (1) can be found directly by other specialized techniques such as the harmonic balance method, the method of conventional oscillator, the WKB method, etc. [24–27]. The T-periodic solution can also be obtained directly and more conveniently by choosing an appropriate set of initial conditions for the vector of variable q, and then solving Eq. (1) within interval [0 , T] under these conditions using a numerical integration methods. For the last approach, an efficient numerical proce- dure was developed to estimate the initial conditions for the T-periodic solution based on the Runge-Kutta method, and tested by a number of applied problems [28, 29]. These studies indicated that the agreement between the experimental and calculating results with Runge-Kutta method is closer than the results calculated by the harmonic balance method and WKB method. To investigate the dynamic stability of elastic mechanisms, we can use Hill method [15] or numerical methods [28, 29]. In this paper, a numerical procedure based on the well-known Newmark integra- tion method is developed to calculate steady-state periodic vibrations of elastic mecha- nisms governed by linearized differential equations with time-periodic coefficients. The proposed procedure is then demonstrated by a dynamic model of a four-bar mechanism with the flexible coupler. 2. NEWMARK PROCEDURE FOR FINDING INITIAL CONDITIONS OF PERIODIC VIBRATION OF LINEAR SYSTEMS The procedure presented below for finding the T-periodic solution of Eq. (1) is based on the Newmark direct integration method. Firstly, the interval [0 , T] is now di- vided into m equal subintervals with the step-size h = ti − ti−1 = T/m. We use notations qi = q(ti) and qi+1 = q(ti+1) to represent the solution of Eq. (1) at discrete times ti and ti+1, respectively. The T-periodic solution must satisfy the following conditions q(0) = q(T), q˙(0) = q˙(T), q¨(0) = q¨(T). (3) Based on the single-step integration method proposed by Newmark, we obtain the following approximation formulas [30, 31] qi+1 = qi + hq˙i + h2 ( 1 2 − β ) q¨i + βh2q¨i+1, (4) q˙i+1 = q˙i + (1− γ) hq¨i + γhq¨i+1, (5) where qi, q˙i, q¨i are approximations to the displacement, velocity and acceleration vec- tors at time ti, β and γ are the constant parameters that define the method. The linear An efficient numerical procedure for calculating periodic vibrations of elastic mechanisms 17 acceleration method, in which a linear variation of the acceleration in the time interval [ti , ti+1] is assumed corresponds to the case γ = 1/4 and β = 1/6. The average accel- eration method is defined by choosing γ = 1/2, β = 1/4. This case corresponds to the assumption that the acceleration is constant over the time interval [ti , ti+1] and equal to 1 2 (q¨i + q¨i+1) [31]. From Eq. (1) we have the following iterative computational scheme at time ti+1 Mi+1q¨i+1 +Ci+1q˙i+1 +Ki+1qi+1 = fi+1, (6) where Mi+1 = M (ti+1) , Ci+1 = C (ti+1) , Ki+1 = K (ti+1) and fi+1 = f (ti+1) . In the next step, substitution of Eqs. (4) and (5) into Eq. (6) yields( Mi+1+γhCi+1+βh2Ki+1 ) q¨i+1= fi+1−Ci+1[q˙i+(1−γ)hq¨i]−Ki+1 [ qi+hq˙i+h2 ( 1 2 −β ) q¨i ] . (7) The use of Eqs. (4) and (5) leads to the prediction formulas for velocities and dis- placements at time ti+1 q∗i+1 = qi + hq˙i + h 2 ( 1 2 − β ) q¨i, q˙∗i+1 = q˙i + (1− γ) hq¨i. (8) Eq. (8) can be expressed in the matrix form as[ q∗i+1 q˙∗i+1 ] = D  qiq˙i q¨i  , (9) with D = [ I hI h2 (0.5− β) I 0 I (1− γ) hI ] , (10) where I denotes the n× n identity matrix, 0 represents the n× n matrix of zeros. Eq. (7) can then be rewritten in the matrix form as q¨i+1 = (Si+1) −1 fi+1 − (Si+1)−1 Hi+1 [ q∗i+1 q˙∗i+1 ] , (11) where matrices Si+1 and Hi+1 are defined by Si+1 = Mi+1 + γhCi+1 + h2βKi+1, (12) Hi+1 = [ Ki+1 Ci+1 ] . (13) By substituting relationships (9) into (11) we find q¨i+1 = (Si+1) −1 fi+1 − (Si+1)−1 Hi+1D  qiq˙i q¨i  . (14) From Eqs. (4), (5) and (8) we get the following matrix relationship qi+1q˙i+1 q¨i+1  = T  q∗i+1q˙∗i+1 q¨i+1  , (15) 18 Nguyen Van Khang, Nguyen Phong Dien, Nguyen Sy Nam where matrix T is expressed in the block matrix form as T =  I 0 Iβh20 I Iγh 0 0 I  . (16) The combination of Eqs. (15), (9) and (14) yields a new computational scheme for determining the solution of Eq. (1) at the time ti+1 in the form qi+1q˙i+1 q¨i+1  = T [ D− (Si+1)−1 Hi+1D ]  qiq˙i q¨i + T  00 (Si+1) −1 fi+1  . (17) In this equation, the iterative computation is eliminated by introducing the direct solution for each time step. Note that the matrices T and D are matrices of constants. By setting xi =  qiq˙i q¨i  , Ai+1 = T [ D− (Si+1)−1 Hi+1D ] , bi+1 = T  00 (Si+1) −1 fi+1  . (18) Eq. (17) can then be rewritten in the following form xi = Aixi−1 + bi (i = 1, 2, . . . , m). (19) Expansion of Eq. (19) for i = 1 to m yields x1 = A1x0 + c1 x2 = A2A1x0 + c2 ................................ xm = ( 1 ∏ i=m Ai ) x0 + cm (20) where c0 = 0, c1 = A1c0 + b1, c2 = A2c1 + b2, . . . , cm = Amcm−1 + bm. Using the condition of periodicity according to Eq. (3), the last equation of Eq. (20) yields a set of the linear algebraic equations( I− 1 ∏ i=m Ai ) x0 = cm. (21) The solution of Eq. (21) gives us the initial value for the periodic solution of Eq. (1). Finally, the periodic solution of Eq. (1) with the obtained initial value can be calculated without difficulties using the computational scheme in Eq. (17). 3. PERIODIC TRANSVERSE VIBRATION OF THE FLEXIBLE COUPLER OF A PLANAR FOUR-BAR MECHANISM One of the most challenging problems in dynamics of machines is the calculation of relative vibrations of elastic members. There are some important cases in which de- formation plays an important role in the dynamic analysis. It happens, for instance, in An efficient numerical procedure for calculating periodic vibrations of elastic mechanisms 19 lightweight high-speed mechanisms [3,4]. The occurrence of large base motion of a mech- anism, in which all bodies are assumed to be rigid, can cause small relative vibrations of flexible links. Conversely, small vibrations of flexible links lead to deviations of its mo- tion from desired base motion, i.e. to dynamic errors. 2 y0 3 x0 O A undeformed rod deformed rod x ( , )w x t E E0 y x  l2 B 4 l4 l1 C Fig. 1. Kinematic schema of a planar 4R-fourbar mechanism with flexible coupler Planar four-bar mechanisms are widely used in reciprocating machines. Steady- state vibrations, dynamic stability analysis and vibration control of the flexible planar four-bar mechanism were the objective of a number of studies, e.g. [6, 7]. In this exam- ple, a planar revolute-jointed 4-link mechanism with the flexible coupler is considered to investigate the problem of relative transverse vibration of the coupler that becomes a se- rious factor at high speed. The kinematic diagram of the mechanism is depicted in Fig. 1. The origin of the ground-fixed coordinate frame {x0, y0} coincides with joint 0 of input link 2. Assuming that the geometrical axis of the coupler (link 3) without deformation is a segment of the straight line that is chosen for x-axis. Neglecting the longitudinal oscilla- tion, the objective of the investigation is to derive analytically the governing equation of the relative transverse vibration of the coupler in the direction of y-axis, and to apply the numerical procedure described in section 2 for finding a periodic solution of the obtained equation. When the angular velocityΩ of input link 2 is assumed to be constant in the steady state, the loop equations can be expressed in the form l2 cosΩt + l3 cos ϕ3 − l4 cos ϕ4 − l1 cos θ = 0, l2 sinΩt + l3 sin ϕ3 − l4 sin ϕ4 − l1 sin θ = 0, (22) where l1, l2, l4 is the length of ground link 1, input link 2 and output link 4 respectively, l3 the undeformed length of coupler 3. Rotation angles ϕ3(t), ϕ4(t) of the coupler and out- put link as well as their time derivatives ϕ˙3, ϕ¨3, ϕ˙4, ϕ¨4 can be determined from Eq. (22) 20 Nguyen Van Khang, Nguyen Phong Dien, Nguyen Sy Nam using a recursive algorithm as Newton-Raphson method. Neglecting higher order non- linear terms, we obtain the partial differential equation of the relative transverse vibration for the coupler in the following simplified form [32] ∂4w ∂x4 + α ∂5w ∂x4∂t − ρ E ∂4w ∂x2∂t2 − [ f4(t) + f3(t)− f2(t) x 2 2 ] ∂2w ∂x2 − [ f3(t)− f2(t)x] ∂w ∂x + µ EI ∂2w ∂t2 + cy EI ∂w ∂t − f2(t)w = − f0(t)− f1(t) (23) where functions fi(t) for i = 0, 1, . . . , 5 are defined by f0(t) = µ EI [ g cos ϕ3 − l2Ω2 sin(Ωt− ϕ3) ] + cyl2 EI Ω cos(Ωt− ϕ3), (24) f1(t) = 1 EI ( cy ϕ˙3 + µϕ¨3 ) , (25) f2(t) = µ EI ϕ˙23, (26) f3(t) = µ EI [ g sin ϕ3 − l2Ω2 cos(Ωt− ϕ3) ] + cxl2 EI Ω sin(Ωt− ϕ3), (27) f4(t)= JC ϕ¨4 + m4gs4 cos(ϕ4 + α4) + Me EI l4 sin(ϕ4 − ϕ3) + + cot (ϕ4 − ϕ3) [ ρ E ϕ¨3 + f0(t) l3 2 + f1(t) l32 3 ] + f2(t) l32 3 − f3(t)l3. (28) The boundary conditions at x = 0 and x = l must be specified for the solution of Eq. (23). These boundary conditions are given by w(0, t) = ∂2 ∂x2 w(0, t) = 0, w(l, t) = ∂2 ∂x2 w(l, t) = 0. (29) Using the mode superposition principle, a solution of Eq. (23) corresponding to the boundary conditions (29) is assumed in the form w(x, t) = n f ∑ i=1 qi(t) sin ipix l3 , (30) where qi(t) are generalized coordinates to be determined. Eq. (23) takes the following form n f ∑ i=1 {[ ρ E ( i pi l )2 + µ EI ] q¨i+ [ α ( i pi l )4 + cy EI ] q˙i+ [( i pi l )4 + ( i pi l )2 f4(t)− f2(t) ] qi } sin ( i pi l x ) + n f ∑ i=1 ( i pi l )2 f3(t)qix sin ( i pi l x ) − n f ∑ i=1 ( i pi l )2 1 2 f2(t)qix2 sin ( i pi l x ) − n f ∑ i=1 ( i pi l )2 f3(t)qi cos ( i pi l x ) + n f ∑ i=1 ( i pi l )2 f2(t)qix cos ( i pi l x ) = − f0(t)− f1(t)x. (31) An efficient numerical procedure for calculating periodic vibrations of elastic mechanisms 21 Multiplying both sides of Eq. (31) by sin(ipix/l) and then integrating the newly obtained equation from 0 to l, we get the system of ordinary differential equations[ ρ E ( i pi l )2 + µ EI ] q¨i + [ α ( i pi l )4 + cy EI ] q˙i + [( i pi l )4 + ( i pi l )2 f4(t) + j2 pi2 2l f3(t)+ − ( 1 4 + j2 pi2 6 ) f2(t) ] qi − n f ∑ i=1,i 6=j i+j=2k+1 αij [ 2 l f3(t)− f2(t) ] qi − n f ∑ i=1,i 6=j i+j=2k+1 αij f2(t)qi = hj(t), (32) where we use indexes j = 0, 1, 2, . . . , n f and k = 1, 2, . . . , and functions hj, αij are hj =  2l jpi f1(t) for j = 2k − 1 jpi [4 f0(t) + 2l f1(t)] for j = 2k + 1 (33) αij = [ 1 (i− j)2 − 1 (i + j)2 ] for i 6= j. (34) Eq. (32) can be rewritten in the compact matrix form as Mq¨(t) + Bq˙(t) +C(Ωt)q(t) = h(Ωt), (35) where M, B, C are n f × n f matrices having coefficients as mij = [ ρ E ( i pi l )2 + µ EI ] δij, bij = [ ρ E ( i pi l )4 + cy EI ] δij, (36) cij =  ( i pi l )4 + ( i pi l )2 f4(t) + j2pi2 2l f3(t)− ( 1 4 + j2pi2 6 ) f2(t) for i = j −αij f2(t) for i 6= j, i + j = 2k −αij [ 2 l f3(t)− f2(t) ] for i 6= j, i + j = 2k + 1 (37) with function δij = { 1 for i = j 0 for i 6= j. The initial value for the periodic solution of Eq. (35) are then determined using Eq. (21), in which parameters γ = 1/4, β = 1/6 and step-size h = 10−4(sec .) were used for the numerical calculation. The calculating parameters are given in Tab. 1. The ob- tained results of the maximal value of coordinates qi in Eq. (35) for a range of rotating speeds of the crank listed in Tab. 2. It can be clearly seen that the transverse vibration re- sponse of the rod can be closely approximated by the first mode since the higher modes are insignificant. The result for transverse vibrations of the connecting rod correspond- ing to different angular velocities of the crank are shown in Fig. 2. In addition, Fig. 3 shows a spectrum calculated by FFT that includes harmonic components of the rotating frequency, such as Ω and 2Ω. The spectrum indicates that the connecting rod performs stationary periodic transverse vibrations only. 22 Nguyen Van Khang, Nguyen Phong Dien, Nguyen Sy Nam Table 1. Calculating parameters Parameters Value Parameters Value l1(m) 1.0 ρ (kg m−3) 7860 l2(m) 0.05 I(m4) 45× 10−10 l3(m) 0.8 cx (kg m−1s−1 ) 0.001 l4(m) 0.8 cy(kg m−1s−1 ) 0.001 s4(m) 0.4 JC(kg m2) 3.35 m3(kg) 3.74 α (s) 10−4 m4(kg) 15 α4(rad.) 0 E(Nm−2) 2.1× 1011 θ(rad.) 0 F(m2) 6× 10−4 Table 2. Calculating results for three modes (n f = 3) Crank speed (rpm) max |q1|(mm) max |q2|(mm) max |q3|(mm) 600 0.2750 0.0090 0.0011 900 0.6309 0.0204 0.0025 1200 1.1583 0.0364 0.0046 1500 1.8704 0.0572 0.0072 Table 2. Calculating results for three modes ( 3)fn  Crank speed (rpm) 1max q (mm) 2max q (mm) 3max q (mm) 600 0.2750 0.0090 0.0011 900 0.6309 0.0204 0.0025 1200 1.1583 0.0364 0.0046 150 1.8704 0.0572 0.0072 Figure 2. Midpoint deflection of the flexible coupler versus crank rotating angle Figure 3. Midpoint deflection of the flexible coupler at 162.8 s (a) versus time and (b) versus angular frequency [s]t 2  1[s ]  ( )b ( )a  radt 62.8 94.2 125.7 Fig. 2. Midpoint deflection of the flexible coupler versus crank rotating angle To verify the correctness of the proposed procedure using Newmark equation, the procedure using the fourth-order Runge-Kutta formula presented in [29] is also applied to solve the same problem. The obtained results with both approaches are identical, but An efficient numerical procedure for calculating periodic vibrations of elastic mechanisms 23 Table 2. Calculating results for three modes ( 3)fn  Crank speed (rpm) 1max q (mm) 2max q (mm) 3max q (mm) 600 0.2750 0.0090 0.0011 900 0.6309 0.0204 0.0025 1200 1.1583 0.0364 0.0046 1500 1.8704 0.0572 0.0072 Figure 2. Midpoint deflection of the flexible coupler versus crank rotating angle Figure 3. Midpoint deflection of the flexible coupler at 162.8 s (a) versus time and (b) versus angular frequency [s]t 2  1[s ]  ( )b ( )a  radt 62.8 94.2 125.7 Fig. 3. M dpoint deflection of the flexible coupler at Ω = 62.8 s−1 (a) versus time and (b) versus angular frequency the computation time with the proposed procedure is greatly reduced, especially in the cases of large number of time steps. 4. CONCLUDING REMARKS In this study, the problem on the calculation of steady-state periodic vibrations of elastic transmission mechanisms is addressed. The differential equations of motion of the mechanism is established and linearized to obtain a system of linear differential equations having time-varying coefficients, known as a parametric vibration system. A numerical algorithm based on Newmark integration method is proposed to determine initial conditions for the periodic solution of a system of linear differential equations with time-periodic coefficients. Using the obtained initial conditions, the pe- riodic solution can be found by a common numerical integration method. For linear systems, this numerical procedure is simpler and easier to implement than one based on the fourth-order Runge-Kutta algorithm which was presented in [29, 32]. Although this approach has been applied to only one example of a flexible four- bar mechanism, but the obtained results are wider applicability to more complicated transmission systems that perform the steady state motions. The problem of vibration control of elastic mechanisms, as addressed in [33, 34], using the periodic solution of the linearized vibration equations will be the subject of our future investigation. ACKNOWLEDGEMENTS This paper was completed with the financial support by the Vietnam National Foundation for Science and Technology Development (NAFOSTED). REFERENCES [1] H. Bremer. Elastic multibody dynamics. Springer, (2008). [2] M. Ge´radin and A. Cardona. Flexible multibody dynamics: A finite element approach. Wiley, (2001). 24 Nguyen Van Khang, Nguyen Phong Dien, Nguyen Sy Nam [3] A. G. Erdman and G. N. Sandor. Kineto-elastodynamics - A review of the state of the art and trends. Mechanism and Machine Theory, 7, (1), (1972), pp. 19–33. [4] G. G. Lowen and C. Chassapis. The elastic behavior of linkages: An update. Mechanism and Machine Theory, 21, (1), (1986), pp. 33–42. [5] S. K. Dwivedy and P. Eberhard. Dynamic analysis of flexible manipulators, a literature re- view. Mechanism and Machine Theory, 41, (7), (2006), pp. 749–777. [6] H. Houben. Untersuchungen u¨ber die stabilita¨t elastischer bewegungen in der koppel eines vierge- lenkgetriebes. PhD thesis, TH Aachen, (1966). [7] H. Horn. Ein beitrag zur dynamischen analyse ebener koppelgetriebe unter beru¨cksichtigung der biegeelastizita¨ten der glieder. PhD thesis, TH Karl-Marx-Stadt, (1974). [8] J. P. Sadler and G. N. Sandor. A lumped parameter approach to vibration and stress analysis of elastic linkages. Journal of Manufacturing Science and Engineering, 95, (2), (1973), pp. 549– 557. [9] G. N. Sandor and X. Zhuang. A linearized lumped parameter approach to vibration and stress analysis of elastic linkages. Mechanism and Machine Theory, 20, (5), (1985), pp. 427–437. [10] B. Bahgat and K. Willmert. Finite element vibrational analysis of planar mechanisms. Mech- anism and Machine Theory, 11, (1), (1976), pp. 47–71. [11] A. Midha, A. G. Erdman, and D. A. Frohrib. Finite element approach to mathematical model- ing of high-speed elastic linkages. Mechanism and Machine Theory, 13, (6), (1978), pp. 603–618. [12] P. K. Nath and A. Ghosh. Kineto-elastodynamic analysis of mechanisms by finite element method. Mechanism and Machine Theory, 15, (3), (1980), pp. 179–197. [13] P. K. Nath and A. Ghosh. Steady state response of mechanisms with elastic links by finite element methods. Mechanism and Machine Theory, 15, (3), (1980), pp. 199–211. [14] W. L. Cleghorn, R. G. Fenton, and B. Tabarrok. Finite element analysis of high-speed flexible mechanisms. Mechanism and Machine Theory, 16, (4), (1981), pp. 407–424. [15] V. V. Bolotin. Kinetische stabilita¨t elastischer systeme. Deutscher Verlag der Wissenschaften Berlin, (1961). [16] M. Badlani and W. Kleinhenz. Dynamic stability of elastic mechanisms. Journal of Mechanical Design, 101, (1), (1979), pp. 149–153. [17] A. Midha, A. G. Erdman, and D. A. Frohrib. A closed-form numerical algorithm for the periodic response of high-speed elastic linkages. Journal of Mechanical Design, 101, (1), (1979), pp. 154–162. [18] W. L. Cleghorn, R. G. Fenton, and B. Tabarrok. Steady-state vibrational response of high- speed flexible mechanisms. Mechanism and Machine Theory, 19, (4), (1984), pp. 417–423. [19] S. D. Yu and W. L. Cleghorn. Dynamic instability analysis of high-speed flexible four-bar mechanisms. Mechanism and Machine Theory, 37, (11), (2002), pp. 1261–1285. [20] V. Masurekar and K. Gupta. Stability analysis of four bar mechanism. Part I: With the as- sumption that damping is absent. Mechanism and machine theory, 23, (5), (1988), pp. 367–375. [21] V. Masurekar and K. N. Gupta. Stability analysis of four bar mechanism. Part II: Taking damping into consideration. Mechanism and machine theory, 23, (5), (1988), pp. 377–382. [22] W. O. Schiehlen and P. C. Mu¨ller. Lineare schwingungen. Akademische Verlagsgesellschaft, Wiesbaden, (1976). [23] J. G. Malkin. Theorie der Stabilitaet einer Bewegung. Akademie Verlag, Berlin, (1959). [24] H. Dresig and I. I. Vul’fson. Dynamik der Mechanismen. Deutscher Verlag der Wissenschaften, Berlin, (1989). [25] J. I. Vul’fson. Analytical investigation of the vibration of mechanisms caused by parametric impulses. Mechanism and Machine Theory, 10, (4), (1975), pp. 305–313. An efficient numerical procedure for calculating periodic vibrations of elastic mechanisms 25 [26] I. I. Vulfson and E. I. Rivin. Vibroactivity of branched and ring structured mechanical drives. Hemi- sphere Pub, (1989). [27] J. Ro¨ssler. Dynamik von mechanismen-antriebssystemen im textil-und verarbeitungsmaschinenbau. PhD thesis, TH Karl-Marx-Stadt, (1985). [28] N. V. Khang. Dynamische stabilita¨t und periodische schwingungen in mechanismen. PhD thesis, TH Karl-Marx-Stadt, (1986). [29] N. V. Khang, N. P. Dien, and H. M. Cuong. Linearization and parametric vibration analysis of some applied problems in multibody systems. Multibody System Dynamics, 22, (2), (2009), pp. 163–180. [30] N. M. Newmark. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, 85, (3), (1959), pp. 67–94. [31] M. Ge´radin and D. Rixen. Mechanical vibrations, Vol. 2. Wiley, Chichester, (1994). [32] V. V. Khiem. Calculation of periodic vibration of mechanisms with flexible link using the Runge- Kutta method. PhD thesis, Hanoi University of Science and Technology, (1998). [33] M. Karkoub and A. S. Yigit. Vibration control of a four-bar mechanism with a flexible coupler link. Journal of Sound and Vibration, 222, (2), (1999), pp. 171–189. [34] D. E. Hill. Dynamics and control of a rigid and flexible four bar coupler. Journal of Vibration and Control, 20, (1), (2014), pp. 131–145.