Parametric Vibration Analysis of Transmission Mechanisms Using Numerical Methods

3.1. Numerical procedure based on Runge-Kutta method

Now we consider only the periodic vibration of a dynamic system which is governed by a set of linear differential equations with periodic coefficients. As already mentioned in the previous section, these differential equations can be expressed in the compact matrix form

where xis 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 Eq. (12) is

As well known from the theory of differential equations, if Eq. (13) has only non-periodic solutions except the trivial solution, then Eq. (12) has an 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 Eq. (12) within interval[0,T]. An algorithm is developed to find the initial value for the periodic solution [18, 19]. Firstly, the T-periodic solution must satisfy the following condition

The interval [0,T]is now divided into mequal subintervals with the step-size h=ti−ti−1=T/m.At the discrete times tiand ti+1,xi=x(ti)and xi+1=x(ti+1)represent the states of the system, respectively. Using the fourth-order Runge-Kutta method, we get a numerical solution [5]

where

Substituting Eq. (16) into Eq. (15), we obtain

where matrix Ai−1is given by

and vector bi−1takes the form

Expansion of Eq. (17) for i=1  to    myields

where c0=0,c1=A0c0+b0,c2=A1c1+b1,..., cm=Am−1cm−1+bm−1.Using the boundary condition according to Eq. (14), the last equation of Eq. (20) yields a set of the linear algebraic equations

The solution of Eq. (21) gives us the initial value for the periodic solution of Eq. (12). Finally, the periodic solution of Eq. (12) with the corresponding initial value can be calculated using the computational scheme according to Eq. (15).

3.2. Numerical procedure based on Newmark integration method

The procedure presented below for finding the T-periodic solution of Eq. (2) is based on the Newmark direct integration method. Firstly, the interval [0,T]is also divided into mequal 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. (2) at discrete times tiand ti+1respectively. The T-periodic solution must satisfy the following conditions

Based on the single-step integration method proposed by Newmark, we obtain the following approximation formulas [6-7]

Constantsβ,γare parameters associated with the quadrature scheme. Choosing γ=1/4and β=1/6leads to linear interpolation of accelerations in the time interval [ti,ti+1]. In the same way, choosing γ=1/2,β=1/4corresponds to considering the acceleration average value over the time interval [6, 7].

From Eq. (2) we have the following iterative computational scheme at timeti+1

where Mi+1=M(ti+1),Ci+1=C(ti+1),Ki+1=K(ti+1)and di+1=d(ti+1).

In the next step, substitution of Eqs. (23) and (24) into Eq. (25) yields

The use of Eqs. (23) and (24) leads to the prediction formulas for velocities and displacements at time ti+1

Eq. (27) can be expressed in the matrix form as

with

where 0represents the n×nmatrix of zeros. Eq. (26) can then be rewritten in the matrix form as

where matrices Si+1and Hi+1are defined by

By substituting relationships (28) into (30) we find

From Eqs. (23), (24) and (27) we get the following matrix relationship

where matrix Tis expressed in the block matrix form as

The combination of Eqs. (28), (33) and (34) yields a new computational scheme for determining the solution of Eq. (2) at the time ti+1in the form

In this equation, the iterative computation is eliminated by introducing the direct solution for each time step. Note that Tand Dare matrices of constants.

By setting

Eq. (36) can then be rewritten in the following form

Expansion of Eq. (38) for i=1  to  myields the same form as Eq. (20)

where c0=0,c1=A1c0+b1,c2=A2c1+b2,..., cm=Amcm−1+bm.

Using the condition of periodicity according to Eq. (22), the last equation of Eq. (39) yields a set of the linear algebraic equations

The solution of Eq. (40) gives us the initial value for the periodic solution of Eq. (2). Finally, the periodic solution of Eq. (2) with the obtained initial value can be calculated without difficulties using the computational scheme in Eq. (36).

Based on the proposed numerical procedures in this section, a computer program with MATLAB to calculate periodic vibrations of transmission mechanisms has been developed and tested by the following application examples.

4.1. Steady-state parametric vibration of an elastic cam mechanism

Cam mechanisms are frequently used in mechanical transmission systems to convert rotary motion into reciprocating motion (Figure 1). At high speed, the vibration of cam mechanisms causes transmission errors, cam surface fatigue, wear and noise. Because of that, the vibration problem of cam mechanisms has been investigated for a long time, both theoretically and experimentally.

The dynamic model of this system is schematically shown in Figure 2. This kind of model was also considered in a number of studies, e.g. [25-26]. The mechanical system of the elastic cam shaft, the cam with an elastic follower can be considered as rigid bodies connected by massless spring-damping elements with time-invariant stiffness kiand constant damping coefficients cifor i=1,2,3.Among them k1is the torsional stiffness of the cam shaft. Parameter k2is the equivalent stiffness due to the longitudinal stiffness of the follower, the contact stiffness between the cam and the roller, and the cam bearing stiffness. Parameter k3denotes the combined stiffness of the return spring and the support of the output link. The rotating components are modeled by two rotating disks with moments of inertia I0and I1.Let us introduce into our dynamic model the nonlinear transmission function U(φ1)of the cam mechanism as a function of the rotating angle φ1of the cam shaft, the driving torque from the motor M(t) and the external load F(t) applied on the system.

The kinetic energy, the potential energy and the dissipative function of the considered system can be expressed in the following form

The virtual work done by all non-conservative forces is

Using the generalized coordinates φ0,φ1,q2,q3,we obtain the following relations

Substitution of Eq. (45) into Eqs. (41-44) yields

where the prime represents the derivative with respect to the generalized coordinate φ1.The generalized forces of all non-conservative forces are then derived from Eq. (49) as

Substitution of Eqs. (46)-(48) and (50) into the Lagrange equation of the second type yields the differential equations of motion of the system in terms of the generalized coordinates φ0,φ1,q2,q3

When the angular velocity Ωof the driver input is assumed to be constant in the steady state

one leads to the following relation

where q1is the difference between rotating angles φ0and φ1due to the presence of the spring element k1and the damping element c1.Assuming that φ1varies little from its mean value during the steady-state motion, the transmission function y1=U(φ1)depends essentially on the input angle φ0=Ωt.Using the Taylor series expansion around Ωt,we get

where we used the notations

Since the system performs small vibrations, i.e. there are only small vibrating amplitudes q1,q2and q3,substituting Eqs. (57)-(59) into Eqs. (52)-(54) and neglecting nonlinear terms, we obtain the linear differential equations of vibration for the system

In most cases, the force F(t)can be approximately a periodic function of the time or a constant. Thus, Eqs. (61)-(63) form a set of linear differential equations with periodic coefficients. Finally, the linearized differential equations of vibration can be expressed in the compact matrix form as

where

We consider now the functionU′(φ), 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 expressed by a truncated Fourier series

The functions U¯′,U¯′′,U¯‴in Eq. (64) can then be calculated using Eq. (65) for φ=Ωt.Parameters used for the numerical calculation are listed in Table 1. Two set of coefficients akin Eq. (46) are given in Table 2 corresponding to two different cases of cam profile, coefficients bk=0.Without loss of generality, the external force Fis assumed to have a constant value of 100 N.

The rotating speed of the driver input nintakes firstly the value of 100 (rpm) corresponding to angular velocity Ω≈10.47(rad/s) for the calculation. The periodic solutions of Eq. (64) are then calculated using the numerical procedures proposed in Section 3. The results of a periodic solution for coordinate q₃, which represents the dynamic transmission errors within the considered system, are shown in Figures 3 and 4. The influence of cam profile to the vibration response of the system can be recognized by a considerable difference in the vibration amplitude of both curves in Figure 3 and the frequency content of spectrums in Figure 4. In addition, the spectrums in Figure 4 shows harmonic components of the rotating frequency, such as ,Ω 3Ω, 5Ω which indicate stationary periodic vibrations.

Figures 5 and 6 show the calculating results with rotating speed nin=600(rpm), corresponding to Ω≈62.8(rad/s). The mechanism has a more serious dynamic transmission error at high speeds. It can be seen clearly from the frequency spectrums that the steady state vibration at high speeds of the considered cam mechanism may include tens harmonics of the rotating frequency as mentioned in [3].

The calculation of the periodic solution of Eq. (64) was implemented by a self-written computer program in MATLAB environment, and a Dell Notebook equipped with CPU Intel® Core 2 Duo T6600 at 2.2 GHz and 3 GB memory. The calculating results obtained by the numerical procedures are identical, but the computation time with Newmark method is greatly reduced in comparison with Runge-Kutta method as shown in Figure 7, especially in the cases of large number of time steps.

4.2. 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 [31-33].

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 [27], [30]. However, according to our knowledge, there are in the literature only a few of 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 investigation is to unravel modulation effects which are responsible for generating such sidebands.

Figure 8 shows a relative simple dynamic model of a pair of helical gears. This kind of the model is also considered in references [24, 28, 32, 33]. The gear mesh is modeled as a pair of rigid disks connected by a spring-damper 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 transmissions error arises from several sources, such as tooth deflection under load, non-uniform 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 czis 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. J₁ and J₂ are the mass moments of inertia of the gears. M₁(t) and M₂(t) denote the external torques load applied on the system. r_b1 and r_b2 represent the base radii of the gears. By introducing the composite coordinate

Eqs. (66) and (67) yield a single differential equation in the following form

where

Note that the rigid-body rotation from the original mathematical model in Eqs. (66) and (67) is eliminated by introducing the new coordinate q(t) in Eq. (69). Variable q(t) is called the dynamic transmission error of the gear-pair system [32]. Upon assuming that when φ˙1=ω1=const,φ˙2=ω2=const,cz=0,kz(t)=k0,the transmission error qis equal to the static tooth deflection under constant load q0as q=rb1φ1+rb2φ2=q0.Eq. (69) yields the following relation

Eq. (69) can then be rewritten in the form

wheref(t)=k0q0−[kz(t)−k0]e(t)−cze˙(t).

In steady state motion of the gear system, the mesh stiffness k_z(t) can be approximately represented by a truncated Fourier series [33]

where ωzis the gear meshing angular frequency which is equal to the number of gear teeth times the shaft angular frequency and Nis the number of terms of the series.

In general, the error components are no 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. When the errors are situated at the teeth of the pinion, e(t) may be taken in the form

Therefore, the vibration equation of gear-pair system according to Eq. (72) is a differential equation with the periodic coefficients.

According to the experimental setup which will be described later, the model parameters includeJ1= 0.093(kgm²), J₂ = 0.272 (kgm²) and nominal pinion speed of 1800 rpm (f₁ = 30 Hz). The mesh stiffness of the test gear pair at particular meshing position was obtained by means of a FEM software [29]. The static tooth deflection is estimated to be q0 = 1.2×10-5(m). The values of Fourier coefficients of the mesh stiffness with corresponding phase angles are given in Table 4. The mean value of the undamped natural frequency ω¯0=k0/mred≈5462s-1,corresponding to f¯0=ω¯0/2π≈869(Hz). Based on the experimental work, the mean value of the Lehr damping ratio ζ¯=0.024is used for the dynamic model. The damping coefficient czcan then be determined by cz=2ω¯0ζ¯mred.