Linearization and parametric vibration analysis of some applied problems in multibody systems

Multibody System Dynamics, Apr 2009

The paper deals with the application of the Runge–Kutta method for calculating steady-state periodic vibrations of the parametric vibration systems governed by linearized differential equations. The numerical calculation is also demonstrated by two models of multibody systems and measurements on real objects. Good agreement is obtained between the numerical and experimental results. Consequently, the obtained results can also be applicable to investigate other complicated models of multibody systems which perform the steady-state motions.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs11044-009-9156-4.pdf

Linearization and parametric vibration analysis of some applied problems in multibody systems

Van Khang Nguyen 0 Phong Dien Nguyen 0 Manh Cuong Hoang 0 0 M.C. Hoang Maritime University , Haiphong City, Vietnam The paper deals with the application of the Runge-Kutta method for calculating steady-state periodic vibrations of the parametric vibration systems governed by linearized differential equations. The numerical calculation is also demonstrated by two models of multibody systems and measurements on real objects. Good agreement is obtained between the numerical and experimental results. Consequently, the obtained results can also be applicable to investigate other complicated models of multibody systems which perform the steady-state motions. The differential equations of motion of a multibody system are fully nonlinear in both the independent or dependent coordinates [1-6]: - M(q, t )q + k(q , q, t ) = h(q , q, t ). The solution of these nonlinear equations is required in order to simulate the dynamic behavior of multibody systems that undergo large displacements and rotations. It is very difficult or impossible to find the analytical solution. The efficient way to solve the problem is by the numerical methods [711]. Many technical systems work mostly on the proximity of an equilibrium position or, especially, in the neighborhood of a desired motion which is usually called programmed motion, desired motion, fundamental motion, inputoutput motion, etc., according to specific problems. In this study, we use the term desired fundamental motion for this object. The desired fundamental motion of a robotic system, for instance, is usually described through state variables determined by prescribed motions of the end-effector. For a mechanical driver system, the desired fundamental motion can be the motion of working components of the driver system, in which the driver output rotates uniformly and all components are assumed to be rigid. It is very convenient to linearize the equations of motion of this configuration in order to take advantage of the linear analysis tools [1216]. In other words, linearization makes it possible to use tools for studying linear systems to analyze the behavior of multibody systems in the vicinity of a desired fundamental motion. For this reason, the linearization of the equations of motion is most useful in the study of control [12, 13], machinery vibrations [14, 16] and the stability of motion [1519]. Mathematically, the linearized equations of motion of a multibody system usually form a set of linear differential equations with time-varying coefficients. Considering steady-state motions of the multibody system alone, one obtains a set of linear differential equations having time-periodic coefficients M(t )x + C(t )x + K(t )x = d(t ). In the steady state of a machine, the working components perform stationary motions [14 17]. Matrices M(t ), C(t ), K(t ) and vector d(t ) in (2) are time-periodic with the least period T . Equation (2) can then be expressed in the compact form as where we use the state variable x such that x = P(t )x + f(t ) x = x = y y and the matrix of coefficients P(t ), vector f(t ) are defined by M1C , f(t ) = M1d , where I is the identity matrix. To calculate the periodic vibrations of mechanical systems described by differential equations (1) or (2), the harmonic balance method, the shooting method and the finite difference method are usually used [9, 10, 21, 24]. To find the steadystate solution of elastic mechanism and machine, we could also use the Newmark or the RungeKutta method. Following this introduction, an overview of the numerical algorithms for calculating the dynamic stability conditions and periodic vibrations of multibody systems in the steady state is presented in Sect. 2. Sections 3 and 4 are important parts of the paper, in which the numerical calculation using the RungeKutta method is demonstrated and validated by two models of multibody systems and measurements on real objects. where P(t ) is a continuous T -periodic n n matrix. According to Floquet theory [15, 17, 18], the characteristic equation of (6) is independent of the chosen fundamental set of solutions. Therefore, the characteristic equation can be formulated in the following way. Firstly, we specify a set of n initial conditions xi (0) for i = 1, . . . , n, their elements xi(s)(0) = and [x1(0), x2(0), . . . , xn(0)] = I, where I denotes the n n identity matrix. By implementing numerical integration of (6) within interval [0, T ] for n given initial conditions, we obtain respectively n vectors xi (T ), i = 1, . . . , n. The matrix (t ) defined by (T ) = x1(T ), x2(T ), . . . , xn(T ) 2 An overview of the calculation of dynamic stability conditions and periodic vibrations 2.1 Dynamic stability conditions We shall consider a system of homogeneous differential equations x = P(t )x (k = 1, . . . , n). When the Floquet multipliers or Floquet exponents are known, the stability conditions of solutions of the system of linear differential equations with periodic coefficients can be easily determined using the Floquet theorem [1518]. 2.2 Periodic vibrations of multibody systems Now we consider only the periodic vibration of a multibody system in the steady state, which is governed by a set of linear differential equations with the periodic coefficients. As is called the monodromy matrix of (6) [17]. The characteristic equation of (6) can then be written in the form = 0. Expansion of (9) yields an nth-order algebraic equation n + a1n1 + a2n2 + + an1 + an = 0 where unknowns i (i = 1, . . . , n), called Floquet multipliers, can be determined from (10). Floquet exponents are given by x = P(t )x + f(t ) x = P(t )x. x(0) = x(T ). already mentioned in the previous section, these differential equations can be expressed in the compact matrix form where x is the vector of state variables, matrix P(t ) and vector f(t ) are periodic in time with period T . The system of homogeneous differential equations corresponding to (12) is As is well known from the theory of differential equations, if (13) has only nonperiodic solutions except the trivial solution, then (12) has a unique T -periodic solution. This periodic solution can be obtained by choosing the appropriate initial condition for the vector of variables x and then implementing numerical integration of (12) within interval [0, T ]. An algorithm is developed to find the initial value for the periodic solution [9, 10, 16, 17]. Firstly, the T -periodic solution must satisfy the following condition: The interval [0, T ] is now divided into m equal subintervals with the step-size h = ti ti1 = T /m. At the discrete times ti and ti+1, xi = x(ti ) and xi+1 = x(ti+1) represent the states of the system, respectively. Using the fourth-order RungeKutta method, we get a numerical solution [8] k(1i1) = h P(ti1)xi1 + f(ti1) , k(4i1) = h P(ti ) xi1 + k(3i1) Substituting (16) into (15), we obtain where matrix Ai1 is given by xi = Ai1xi1 + bi1 h P2 ti1 + 2 1 h P(ti1) + 2 P(ti )P2 ti1 + 2 P(ti1) (i = 1, . . . , m) and vector bi1 takes the form h P2 ti1 + 2 f(ti1) . Expansion of (17) for i = 1 to m yields i=m1 xm = I i=m1 Ai x0 = cm where c0 = 0, c1 = A0c0 + b0, c2 = A1c1 + b1, . . . , cm = Am1cm1 + bm1. Using the boundary condition according to (14), the last equation of (20) yields a set of the linear algebraic equations where I is the n n identity matrix. The solution of (21) gives us the initial value for the periodic solution of (12). Finally, the periodic solution of (12) with the corresponding initial value can be calculated without difficulty by using numerical methods. Based on the described above RungeKutta method, a computer program for calculating periodic vibrations of multibody systems has been developed at Hanoi University of Technology. The computer program has been used for the numerical calculation of the following two application examples. 3 Periodic vibration of the transport manipulator of a forging press The most common forging equipment is the mechanical forging press. Mechanical presses function by using a transport manipulator with a cam mechanism to produce a preset at a certain location in the stroke. The kinematic schema of such mechanical adjustment unit is depicted in Fig. 1. The dynamic model of this system is schematically shown in Fig. 2. The mechanical system of the driver shaft, the flexible transmission mechanism and the operating mechanism can be considered as rigid bodies connected by massless springdamping elements with time-invariant stiffness ki and constant damping coefficients ci , i = 1, 2. The rotating components are modeled by two rotating disks with moments of inertia I0 and I1. Let us Fig. 1 Kinematic schema of the transport manipulator of a forging press, 1the first gearbox, 2driving shaft, 3the second gearbox, 4cam mechanism, 5operating mechanism (hammer) Fig. 2 Dynamic model of the transport manipulator introduce into our dynamic model the nonlinear transmission function U (1) of the cam mechanism as a function of the rotating angle 1 of the cam shaft, the driving torque from the motor M (t ) and the external load F (t ) applied on the system. The kinetic energy and the potential energy of the considered system can be expressed in the following form: By using the generalized coordinates 0, 1 and q2, we get the generalized forces T = 21 I002 + 21 I112 + 21 m2x22, 1 1 = 2 k1(1 0)2 + 2 k2(x2 y)2. Q0np = M (t ) + c1(1 0), Q1np = F (t )U c1(1 0), Q2np = F (t ) c2q2 U ( t + q1) = U + U q1 + 21 U q12 + where we used the notations Since the system performs small vibrations, i.e. there are only small vibrating amplitudes q1 and q2, substituting (30) into (26) and (27) and neglecting nonlinear terms, we obtain the linear differential equations of vibration for the system: where the prime represents the derivative with respect to the generalized coordinate 1. Substitution of (22)(24) into the Lagrange equation of the second type yields the differential equations of motion of the system I00 c1(1 0) k1(1 0) = M (t ), I1 + m2U 2 1 + m2U q2 + m2U U 1 + c1(1 0) + k1(1 0) 2 = F (t )U , 2 m2U 1 + m2q2 + m2U 1 + c2q2 + k2q2 = F (t ). When the angular velocity of the driver input is assumed to be constant in the steady state, one leads to the following relation: where q1 is the difference between rotating angles 0 and 1 due to the presence of the spring element k1 and the damping element c1, resulted from the flexible transmission mechanism. If we assume that 1 varies little from its mean value during the steady-state motion, then the transmission function y = U (1) depends essentially on the input angle 0 = t . Using the Taylor series expansion around t , we get d = q = We consider now the function U (), called the first grade of the transmission function U (), where the angle is the rotating angle of the cam shaft. In steady-state motion of the cam mechanism, function U () can be approximately represented by a truncated Fourier series [14, 16] k=1 U () = (ak cos k + bk sin k). The functions U , U , U in (34) can then be calculated using (35) for = t . The following parameters are used for numerical calculations: rotating speed of the driver input n = 50 (rpm) corresponding to = 5.236 (1/s), stiffness k1 = 7692 N m; k2 = 106 N/m, damping coefficients c1 = 18.5 N m s; c2 = 2332 N s/m, mass moment of inertia I1 = 1.11 kg m2 and mass m2 = 136 kg. The Fourier coefficients ak in (35) with K = 12 are given in Table 1 for two different cases and coefficients bk = 0. We consider only periodic vibrations which are a commonly observed phenomenon in the system. The periodic solutions of (34) can be obtained by choosing appropriate initial conditions for the vector of variables q. In most cases, the force F (t ) can be approximately a periodic function of the time or a constant. Thus, (32) and (33) form a set of linear differential equations with periodic coefficients. Finally, the linear vibration equations of the mechanical adjustment unit can be expressed in the compact matrix form as Table 1 The Fourier coefficients of the function U () Table 2 The ||max values of the characteristic equation Fig. 3 Calculating results for q1, (a) curve 1using the proposed numerical method, curve 2using WKB method, (b) frequency spectrum corresponding to curve 1 To verify the dynamic stable condition of the vibration system, the maximum of absolute value ||max of the solutions of the characteristic equation, according to (34), is now calculated. The obtained values corresponding to both cases are given in Table 2. It can be concluded that these values satisfy the stable condition ||max < 1 and the system is dynamically stable for both cases. Some calculating results for periodic vibrations of the mechanical adjustment unit are shown in Figs. 35. The periodic solutions of (34) obtained by the RungeKutta method are compared with the calculating results using the WKB method [20] in Figs. 3a and 4a. Comparing both, a considerable difference in the amplitude can be recognized. In addition, the frequency spectra show harmonic components of the rotating frequency, such as 2 , 4 , 6 (Fig. 3b) or , 3 , 5 (Fig. 4b). These spectra indicate that the considered system performs stationary periodic vibrations only. To validate the calculating results using the RungeKutta method, the dynamic load moment of the mechanical adjustment unit was measured on the driving shaft (see also Fig. 1). A typical record of the measured moment is plotted in Fig. 5, together with the curves calculated from the dynamic model by using the WKB method, the kinesto-static calculation and the RungeKutta method. Comparing the curves displayed in this figure, it can be observed that the results using the RungeKutta method more closely agree with the experimental results than the results obtained by the WKB method agree with the kinesto-static calculation. Fig. 4 Calculating results for q2, (a) curve 1using the proposed numerical method, curve 2using WKB method, (b) frequency spectrum corresponding to curve 1 Fig. 5 Dynamic moment acting on the driving shaft of the mechanical adjustment unit 4 Parametric vibration of a gear-pair system with faulted meshing Dynamic modeling of gear vibrations offers a better understanding of the vibration generation mechanisms as well as the dynamic behavior of the gear transmission in the presence of gear tooth damage. Since the main source of vibration in a geared transmission system is usually the meshing action of the gears, vibration models of the gear-pair in mesh have been developed, taking into consideration the most important dynamic factors such as effects of friction forces at the meshing interface, gear backlash, the time-varying mesh stiffness and the excitation from gear transmission errors [2729]. Fig. 6 Dynamic model of the gear-pair system with faulted meshing From experimental works, it is well known that the most important components in gear vibration spectra are the tooth-meshing frequency and its harmonics, together with sideband structures due to the modulation effect. The increment in the number and amplitude of sidebands may indicate a gear fault condition, and the spacing of the sidebands is related to their source [22, 26]. However, according to our knowledge, there are in the literature only a few theoretical studies concerning the effect of sidebands in gear vibration spectrum, and the calculating results are usually not in agreement with the measurements. Therefore, the main objective of the following study is to unravel modulation effects which are responsible for generating such sidebands. Figure 6 shows a relatively simple dynamic model of a pair of helical gears. This kind of the model is also considered in Refs. [23, 24, 27, 28]. The gear mesh is modeled as a pair of rigid disks connected by a springdamper set along the line of contact. The model takes into account influences of the static transmission error which is simulated by a displacement excitation e(t ) at the mesh. This transmission error arises from several sources, such as tooth deflection under load, nonuniform tooth spacing, tooth profile errors caused by machining errors, as well as pitting, scuffing of teeth flanks. The mesh stiffness kz(t ) is expressed as a time-varying function. The gear-pair is assumed to operate under high torque condition with zero backlash and the effect of friction forces at the meshing interface is neglected. The viscous damping coefficient of the gear mesh cz is assumed to be constant. The differential equations of motion for this system can be expressed in the form where i , i , i (i = 1, 2) are rotation angle, angular velocity, angular acceleration of the input pinion and the output wheel, respectively. J1 and J2 are the mass moments of inertia of the gears. M1(t ) and M2(t ) denote the external torque loads applied on the system, rb1 and rb2 represent the base radii of the gears. By introducing the composite coordinate q = rb11 + rb22, (36) and (37) yield a single differential equation in the following form: mredq + kz(t )q + czq = F (t ) kz(t )e(t ) cze(t ) F (t ) = mred Note that the rigid-body rotation from the original mathematical model in (36) and (37) is eliminated by introducing in (39) a new coordinate q(t ). Variable q(t ) is called the dynamic transmission error of the gear-pair system [28]. Upon assuming that when 1 = 1 = const, 2 = 2 = const, cz = 0, kz(t ) = k0, the dynamic transmission error of the gear-pair system q is equal to the static tooth deflection under constant load q0 as q = rb11 + rb22 = q0. Equation (39) yields the following relation: F (t ) F0(t ) = k0q0 + k0e(t ). Equation (39) can then be rewritten in the form mredq + kz(t )q + czq f (t ) = 0 where f (t ) = k0q0 [kz(t ) k0]e(t ) cze(t ). In steady-state motion of the gear system, the mesh stiffness kz(t ) can be approximately represented by a truncated Fourier series [29] kz(t ) = k0 + where z is the gear meshing angular frequency which is equal to the number of gear teeth, times the shaft angular frequency, and N is the number of terms of the series. In general, the error components are not identical for each gear tooth and will produce displacement excitation that is periodic with the gear rotation (i.e. repeated each time the tooth is in contact). The excitation function e(t ) can then be expressed in a Fourier series with the fundamental frequency corresponding to the rotation speed of the faulted gear. For instance, when the errors are situated at the teeth of the pinion, e(t ) may be taken in the form e(t ) = Therefore, the vibration equation of gear-pair system according to (42) is a differential equation with the periodic coefficients. In this example, ten dominant coefficients in the Fourier series of the mesh stiffness expressed in (43) are taken into account: n=1 n=1 i=1 where z = z11 and the excitation function e(t ) is expressed by the six first terms of its Fourier series as follows: Table 4 Fourier coefficients and phase angles of the mesh stiffness 0 1 2 3 4 5 6 7 8 9 10 According to the experimental setup which will be described later, the following parameters of the model are used for numerical simulation: J1 = 9.3 102 (kg m2); J2 = 0.272 (kg m2). Other parameters of the gears are shown in Table 3. A nominal pinion speed of 1800 rpm (1 = 60 or f1 = 30 Hz) is chosen. The mesh stiffness of the test gear-pair at particular meshing position was obtained by means of a FEM software [25]. The static tooth deflection is estimated to be q0 = 1.2 105 (m). The values of Fourier coefficients of the mesh stiffness with corresponding phase angles are given in Table 4. So, the mean 0 value of the undamped natural frequency can be determined as 0 = k0/mred 5462 (1/s), corresponding to f0 = 0/2 869 (Hz). Based on the experimental work [30], the mean value of the Lehr damping ratio = 0.024 is used for the dynamic model. The damping coefficient cz can then be determined by cz = 2 0 mred. The numerical calculation for finding the periodic solutions of (42), based on the Runge Kutta method, is realized with the computing program MATLAB. The calculated dynamic transmission errors are shown in Figs. 7 and 8, corresponding to different excitation functions e(t ). The spectra in Fig. 8, (a) and (b), show clearly the meshing frequency and its harmonics with sideband structures. As expected, the sidebands are spaced by the rotational frequency f1 of the pinion. Comparing amplitudes of these sidebands in both spectra, it can be concluded that the excitation function e(t ) caused by the tooth errors is responsible for generating sidebands. Fig. 7 Modeling result: dynamic transmission error q(t) Fig. 8 Modeling result: frequency spectrum of dq/dt , (a) excitation function e(t) with coefficients e1 = 0.0015, e2 = 0.0035, e3 = 0.0027, e4 = 0.0011, e5 = 0.0005, e6 = 0.0013 (mm) and phase angles 1 = 0.049, 2 = 1.7661, 3 = 0.7286, 4 = 0.5763, 5 = 0.7810, 6 = 1.8172 (radian), (b) excitation function e(t) with larger coefficients e1 = 0.01, e2 = 0.003, e3 = 0.0018, e4 = 0.0011, e5 = 0.0009, e6 = 0.0003 (mm) and phase angles 1 = 1.047, 2 = 1.4521, 3 = 0.5233, 4 = 1.457, 5 = 0.8622, 6 = 1.1966 (radian) The experiment was carried out at an ordinary back-to-back test rig (Fig. 9). The major parameters of the test gear-pair are given in Table 3. The load torque was provided by a hydraulic rotary torque actuator which remains the constant external torque for any motor speed. The test gearbox operates at a nominal pinion speed of 1800 rpm (30 Hz), thus the meshing frequency fz is 420 Hz. A Laser Doppler Vibrometer was used for measuring the Fig. 9 Gearbox test rig Fig. 10 Test rig and measuring configuration (schema) oscillating parts of the angular speed of the gear shafts (i.e. oscillating parts of 1 and 2) in order to experimentally determine the dynamic transmission error (Fig. 10). The measurement was taken with two noncontacting transducers mounted in proximity to the shafts, positioned at the closest place to the test gears. The vibration signals were sampled at 10 kHz. The signal used in this study was recorded at the end of 12hour total test time. At that time a surface fatigue failure occurred on some teeth of the pinion. Figure 11 shows the frequency spectrum of the first derivative of the dynamic transmission error q (t ) determined from the experimental data. The spectrum presents sidebands at the meshing frequency and its harmonics. In particular, the dominant sidebands are spaced by the rotational frequency of the pinion and characterized by high amplitude. This gives a clear indication of the presence of the faults on the pinion. By comparing the spectra displayed in Figs. 12 and 13, it can be observed that the vibration spectrum calculated by the RungeKutta method (Fig. 12) and the spectrum of the measured vibration signal (Fig. 13) show the same sideband structures. 178 Fig. 11 Experimental results: frequency spectrum of dq/dt Fig. 12 Calculating results: frequency spectrum of dq/dt Fig. 13 Experimental results: zoomed frequency spectrum of dq/dt from Fig. 11 Some early calculating results for the same problem of gear vibrations could also be found in the previous study [23], in which the harmonic balance method is employed for solving (42). In study [23], however, only two first terms of the excitation function e(t ) in (46) and three coefficients k1, k2, k3 of the Fourier series of the mesh stiffness in (45) were taken into account since the computational cost increases highly with the number of harmonic terms. On the contrary, by using the RungeKutta method, six terms of the excitation function and up to ten coefficients of the mesh stiffness could be taken in the calculation with an acceptable computational cost. Consequently, the calculating results with the RungeKutta method (Fig. 12) more closely agree with the experimental data than the results calculated by the harmonic balance method in [23]. The obtained results seem promising and an extension to the more complicated geared systems. 5 Conclusions The calculation of dynamic stable conditions and periodic vibrations of elastic mechanisms and machines present an important problem in mechanical engineering. In the this paper, the RungeKutta method is used for calculating the steady-state vibration of the parametric vibration systems governed by linearized differential equations. This numerical method is also demonstrated and tested with two dynamic models of multibody systems and measurements on real objects. Good agreement is obtained between the numerical solutions and experimental results. In the first example, the agreement between the experimental and calculated results with RungeKutta method is closer than the results calculated by WKB method. In the second example, the concordance of the results is obtained for both the amplitude and the frequency content of the parametric vibration of the gear-pair system. The obtained results can also be applicable to investigate other complicated multibody systems which perform the steady-state motions. The application of the RungeKutta method for calculating nonlinear periodic vibrations and bifurcations of multibody systems will be the subject of our future investigation. Acknowledgements This paper was completed with a research grant from the Flemish Inter-university Council for University Development Cooperation (VLIR UOS). An additional support was given by the National Foundation for Science and Technology Development of Vietnam. The authors gratefully acknowledge the generous support of these sponsors. Open Access This article is distributed under the terms of the Creative Commons Attribution Noncommercial License which permits any noncommercial use, distribution, and reproduction in any medium, provided the original author(s) and source are credited.


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs11044-009-9156-4.pdf

Van Khang Nguyen, Phong Dien Nguyen, Manh Cuong Hoang. Linearization and parametric vibration analysis of some applied problems in multibody systems, Multibody System Dynamics, 2009, 163-180, DOI: 10.1007/s11044-009-9156-4