Investigation on the influence of unbalanced shaft component in gearbox on displacement using the Newmark-β method

TL;DR

This paper presents a new numerical method to analyze how unbalanced shaft components in gearboxes affect displacement and torsional deformation.

Contribution

The integration of a variable time-step algorithm with the Newmark-β method improves accuracy and speed in dynamic analysis.

Findings

The model shows less than 6% deviation in displacement compared to experimental data below critical speed.

The method enables rapid convergence and reliable predictions for system dynamics.

Parametric analysis reveals the impact of rotational speed, trial mass, and phase angle on system behavior.

Abstract

This study introduces an enhanced numerical approach for analyzing the dynamic behavior of a rotor-bearing system subjected to unbalanced excitation from a gearbox drive shaft. The Newmark-β method with the integration of a variable time-step algorithm was used, allowing the system to be solved rapidly and accurately without compromising stability. This technique enables a precise computation of displacement and torsional deformation of the rotating shaft during its operational cycle. The proposed computational model is validated against experimental data, showing deviations of displacement in normal operation below the critical speed of about 6%. A comprehensive parametric analysis is conducted to evaluate the influence of rotational speed, trial mass, and initial phase angle on the system dynamics. The findings confirm that our enhanced numerical approach yields rapid convergence and reliable predictions, making it a valuable tool for dynamic analysis of rotating systems.

Figures

Gear shaft assembly in the lathe gearbox.

Schematic diagram of the rotor–bearing system.

Forces and unbalance positions of rotor-bearing system with eccentric masses.

Modelling of the study shaft segment from nodes 2, 4 and 6.

A flowchart for dynamic response of the rotor system.

A vibration testing machine of rotating shaft.

A schematic plot of a rotor operating at the first critical speed (mode 1).

Displacement of node 4 along X-axis at 2000 rpm.

Dimension of experimental unbalanced specimen.

The measurement results of balancing used to two planes at the angular position 00−00; mtrial = 30 grams.

Horizontal displacement of node 4 for mtrial = 20 g at various speeds.

Horizontal displacement of node 4 under 20 g trial mass and speed of 2000 rpm, for the different angular positions: 00–00, 00–900, and 00–1800.

Motional orbit of the shaft center at 2000 rpm for various angular positions.

Equations

Full text

1. Introduction

The drive shafts in gearboxes are common components with variety of applications in industries, operating at various speeds and load types. Among many factors causing failure of shafts, vibration due to unbalance is a major cause. Single and coupled vibration modes, including torsional, longitudinal and transverse vibrations induce fatigue, fracture and tribological issues on the rotating shaft components [1–3]. These vibrations result in displacement, performance of gear transmission, wear and cracks [4–6].

Many studies have been carried out on dynamic aspects of the of rotor systems. The modelling of the rotor systems often uses the Jeffcott rotor model having a massless axis and a mass disk placed in the middle of the shaft. Modern analytical methods have been utilized to provide a foundational understanding of rotor dynamics and to conduct simple model experiments [7,8]. In addition, [9] and [10] studied a rotor system with a dynamic model, in which a disc is placed in the middle of a massless elastic rotating shaft. The equations of motion are obtained to Lagrangian dynamics for transverse – torsional vibrations. [11] derived the equation of motion by assuming that the diesel engine drive system can be approached as a simple rotor model such as the Jeffcott rotor. A modified version of this rotor model was also used for analysing the coupled torsional–transverse vibrations of a propeller shaft resulting from misalignment induced by shaft rotation [12]. Besides, [13] modelled a flexible rotating shaft system subject to bending and torsion coupled with the shaft and disk moving away from the center point of the shaft. In these studies, it is evident that shaft displacement caused by vibration is a complex problem. Therefore, it is necessary to investigate the shaft’s displacement in detail due to unbalance, using a new approach and model in the study of rotor dynamics.

The Newmark-β method, that is a widely used numerical integration technique in finite element analysis, is particularly effective for simulating dynamic systems. A study by [14] used this method to analyze the dynamic characteristics of spur and helical gear systems, while a subsequent study focused on the influence of housing flexibility on gear transmission dynamics [15]. An improved Newmark-β method was applied in [16] to determine nonlinear dynamics and reduce crankshaft torsional vibration. An enhanced version of the method was proposed for long-term simulations [17], offering improved convergence in evaluating the effects of nonlinearities on engine crankshaft torsional behavior. Collectively, these studies confirm that the Newmark-β method is a suitable and reliable tool for addressing time-varying oscillation problems.

In general, the current studies lack experimental evidence to support the theoretical computation. This study represents a computational dynamic model for a rotor-bearing system under unbalanced excitation, using the Newmark-β method with a variable time-step algorithm to determine the unbalance amount and shaft displacement. Various key experimental parameters, including rotational speed, trial mass, and initial phase angle, are put into account. The displacement and the amount of unbalance, determined from experiment, are compared with the numerically computational values to verify the validity of the computational model and thus to evaluate their effects on system dynamics. The motional orbit, representing shaft displacement, is also determined to evaluate the influence of the unbalance on the deflection of the shaft. This will allow to predict the fatigue strength, thus to accurately predict performance during operation.

2. Dynamic model of the rotor system

2.1 The proposed model

This issue arises from the eccentricity and imbalance of rotating components. If the gearbox operates under long-term conditions, it may lead to wear and eventual failure due to fatigue fracture. Fig 1 shows a shaft assembly in the gearbox of a lathing machine, comprising the drive shaft 3, bearings 1 and 5, and gears 2 and 4. Fig 2 represents the schematic diagram of a generalized model comprising a shaft and various assembled disks.

Gear shaft assembly in the lathe gearbox.

Schematic diagram of the rotor–bearing system.

2.2 Parameters of the model

The proposed model of transmission shaft with two disks is shown in Fig 3. The specimens for evaluating the fatigue strength of metallic materials were used, according to ISO 1143:2010. This specimen type is also well-suited for evaluating location to failure [18]. Seven nodes, from node 1 to node 7, on the rotating shaft are selected to be investigated. Nodes 1 and 7 are the bearing positions. Nodes 2 and 6 are the disk placements and nodes 3 and 5 correspond to the positions with the maximum cross-section on the shaft. Node 4 has the smallest cross-section, where the fracture occurs.

Forces and unbalance positions of rotor-bearing system with eccentric masses.

The displacement vector q_i_ for the shaft nodes i (i = 1–7) is:

$$[eqn]$$

where the displacements x_i_ and y_i_ along the X- and Y-axes and torsional angle α_i_ at the investigated positions are

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Since the loads are applied onto the gears at nodes 2 and 6, the rotation angle on the shaft segment between nodes 2 and 6 are thoroughly investigated. The torsional angles at nodes 3 and 5 are assumed to be small and negligible (α_3_ = α_5_ = 0). Thus, the torsional angles at nodes 2, 4, and 6 are $[eqn]$ where i = 2, 4, 6. The governing equation for the rotation angle is given as:

$$[eqn]$$

where $[eqn]$ (t) is the total rotation angle at node i as time t; $[eqn]$ is the initial phase angle at node i; $[eqn]$ is the torsional angle at node i; $[eqn]$ is the angular velocity at node i (i = 2, 4, 6). The kinetic energy of the system is represented as:

$$[eqn]$$

where T t is translational kinetic energy; T_r_ is rotational kinetic energy; $[eqn]$ , $[eqn]$ are the vectors of translational velocities in the X and Y directions; $[eqn]$ is the vector of angular velocities of the components; $[eqn]$ is the mass matrix corresponding to the system from node 1 to node 7; [ $[eqn]$ ] is the static moments of inertia concerning the study nodes 2, 4, 6. The elastic potential energy of the system is:

$$[eqn]$$

where [Kx} and [Ky] are respectively the stiffness matrices corresponding to translational stiffness in the X- and Y-axes, [Kt} is the torsional stiffness matrix related to angular displacements. According to Fig 3, an examination of the shaft segment from node 2 to node 6 represented by l3 to l6 was performed. Fig 4 shows the spring modeling of the shaft segment.

Modelling of the study shaft segment from nodes 2, 4 and 6.

In the torsional shaft problem, x_2_, x_4_, and x_6_ are equated to the torsional angles α_2_, α_4_ and α_6_ of the shaft at nodes 2, 4 and 6, respectively. From the modelling of the shaft segment in Fig 4, the elastic potential energy is written as:

$$[eqn]$$

The dissipated energy of the system is given by:

$$[eqn]$$

where D is the total energy dissipated due to damping in the system; [Cx], [Cy}, and [Ct} are respectively the damping matrices corresponding to translational damping in the X, Y directions and axial torsional damping. The dissipated energy due to damping is written as:

$$[eqn]$$

where $[eqn]$ are the torsional damping at nodes 2, 4, and 6. The Lagrange’s equations of the second kind can be written in terms of the system as follows [19]:

$$[eqn]$$ $$[eqn]$$

where T is the total kinetic energy, V is the total potential energy of the system; $[eqn]$ is force due to the mass of the node; $[eqn]$ is centrifugal force due to unbalance. The differential equation of motion for the rotor at this time is:

$$[eqn]$$

where $[eqn]$ is the mass matrix of the system; $[eqn]$ is damping matrix; $[eqn]$ is a stiffness matrix. The displacement vector $[eqn]$ includes the displacements in the X-direction, Y-direction, as well as angular rotation, $[eqn]$ is excitation force vector:

$$[eqn]$$ $$[eqn]$$

By solving equation (10) for i = 2, 4, 6; we obtain:

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

2.3 Solution method

Using the Newmark time-stepping analysis method to solve Equation (9), the following variables and matrices are involved: $[eqn]$

Take the derivation of Equation (9), we have:

$$[eqn]$$

where

$$[eqn]$$

Transforming Eq. (12), we obtain:

$$[eqn]$$

Equation (13) can be generalized as:

$$[eqn]$$

where $[eqn]$ is the effective stiffness matrix, determined as:

$$[eqn]$$

$[eqn]$ is the effective force vector, determined as:

$$[eqn]$$

In this study, the implicit Newmark–β method was employed, since it does not require adherence to the time step condition as $[eqn]$ . The parameters $[eqn]$ were adjusted during the time-stepping loop and set to $[eqn]$ A constant time step of Δt = 0.01s was selected for the simulations to ensure stability and computational efficiency [20]. Fig 5 shows the algorithm flowchart, used for computing the displacement $[eqn]$ , the velocity $[eqn]$ , and the acceleration $[eqn]$ .

A flowchart for dynamic response of the rotor system.

3. Experimental procedure

An AISI 1045 steel transmission shaft as shown in Fig 3 with the unbalances m2 and m6 was prepared to investigate the behavior of nodes 1–7 during rotational operation. Table 1 shows the geometry dimensions and mechanical properties of the shaft [12,21,22].

Table 1: Symbols, parameters, and unit of power transmission shaft system.

Fig 6 shows the vibration testing machine, designed and fabricated for study the vibration behavior of the rotating shaft. The shaft was driven by an AC motor and an encoder was attached to the motor spindle to record the actual speed of the rotating shaft. A laser sensor measured the displacement of the rotor. The accelerometer and processor recorded the amount and positions of unbalance on the rotor. The operating tests were performed in controlled modes.

A vibration testing machine of rotating shaft.

4. Results and discussions

4.1 The effect of rotational speed on shaft displacement

In Fig 7, the rotor oscillates when operated arround the first critical speed (mode 1). Under the effect of centrifugal force, the shaft is deflected and the shaft center will oscillate with displacement at position 4 with amplitude y_max_.

A schematic plot of a rotor operating at the first critical speed (mode 1).

Neglecting the damping coefficient C, the first critical speed of the system can be determined according to [23] using the following expression:

$$[eqn]$$

In this case, the critical speed is:

$$[eqn]$$

To ensure the stability of the rotor in the model and corresponding to the actual working speeds of the equipment, the experimental speed ranges of the rotor were selected as follows: n1* *= 800 rpm, n2 = 1500 rpm and n3 = 2000 rpm, in which n3 is selected to be less than 70% of critical speed of 2132 rpm [24]. Fig 8 shows the displacement of node 4 in the X direction at speed of n3 = 2000 rpm, determined from the experiment and the Newmark–β numerical simulation. Generally, the amplitude measured from the experiment is larger than those using Newmark–β numerical computation because various factors in the manufacturing process may affect the unbalanced conditions.

Displacement of node 4 along X-axis at 2000 rpm.

Table 2 compares the X and Y displacement of node 4, determined from the Newmark-β method and from the experiment. The displacement determined experimentally is higher than those determined from the theoretical simulation because the actual stiffness of shaft and bearing is not uniform and due to the inaccuracy and clearance between shaft and bearing. For the eperation below the first critical speed, the vibrations remain stable, with displacement amplitudes of X = ± 0.046 and Y = ± 0.081. The relative error of displacement in the X- and Y-axes for speeds of 2000 rpm is 6% and 12%, respectively. As the rotor reaches the first critical speed of 3000 rpm, the vibrational amplitude significantly increases, with displacements of X= ± 0.101 and Y= ± 0.119 mm because of the resonance effects of the rotor components. In general, the displacement in the direction Y is higher than in the X direction because of the gravity force of the components, acting downward, resulting in the higher vertical vibration. The relative errors in displacements are higher by 12% and 18%.

Table 2: Displacement in the X and Y axes at node 4.

A comparison between the numerical results obtained using the Newmark–β method and the experimental data indicated that the displacement deviation of approximately 6%, which is acceptable for engineering analysis and validates the accuracy of the computational model.

4.2 The effect of unbalance on shaft displacement

In engineering, for some design circumstance or unavoidable eccentricity in the manufacturing process, the unbalance exists and should be taken into account. The unbalance, eccentricity, and angular position (the angle between the planes of trial masses m_2_ and m_6_) will be analyzed hereafter. The masses m_2_ and m_6_ were mounted at three positions 0^0^–180^0^, 0^0^–90^0^, and 0^0^–0^0^ such that the direction of m_6_ respectively makes angular positions of 180^0^, 90^0^ and 0^0^ to the direction of m_2_, as shown in Fig 9.

Dimension of experimental unbalanced specimen.

4.2.1 Symmetrical trial mass (angular position: 00–00 on both disks).

The structure of the unbalanced rotor is shown in Fig 9, where the trial masses m_2_ and m_6_ are in the same direction. According to the unbalance standard, the corresponding test masses considered in the range follows as G1, G2.5, G6.3 and G16 [25]. For G1, the allowable eccentricity e (n_3_ = 2000 rpm) is:

$$[eqn]$$

At the minimum speed n_1_, the allowable test mass for G16 is calculated as follows:

$$[eqn]$$

Three kinds of trial masses 10 g, 20 g, and 30 g and rotational speeds of 800, 1500 and 2000 rpm were selected to determine the unbalanced, using the commercial Erbessd – Instruments interface, as shown in Fig 10. For the highest trial mass m_trial_ of 30 grams, the unbalanced amounts are respectively 41.18 mm/s and 39.06 mm/s for two nodes 2 and 6.

The measurement results of balancing used to two planes at the angular position 00−00; mtrial = 30 grams.

Table 3 represents the measurement results of rotor unbalance at angular position 0^0^–0^0^ with corresponding to load values and speeds. The corresponding unbalance significantly increases with the increase of trial mass mtrial and the rotational speed n respectively. The unbalanced amount G was highest of approximately 41.1 mm/s at n of 2000 rpm with mtrial of 30 g and got lowest of 0.58 mm/s at 800 rpm and mtrial of 0 g.

Table 3: Unbalance measurement results when attaching trial mass at angular position 00–00.

Fig 11 and Table 4 represent the variation in horizontal X-axis displacement at node 4 under different operating speeds. The results clearly show the displacement amplitude increase with the rotational speed, particularly as the rotor approaches its first critical speed. This behavior reflects a corresponding rise in dynamic excitation, which is characteristic of resonance phenomena in rotor-dynamic systems.

Table 4: Displacement measurement results at node 4 for 20g trial mass.

Horizontal displacement of node 4 for mtrial = 20 g at various speeds.

4.2.2 Asymmetrically placed trial mass (angular positions 00–900 and 00–1800).

By changing the initial phase angle, the results of rotor unbalance at positions 0^0^–90^0^ and 0^0^–180^0^ are shown in Table 5. At a constant rotational speed of 2000 rpm, variations in the load angular position 0^0^–0^0^, 0^0^–90^0^, and 0^0^–180^0^ and trial mass m**ₜᵣᵢₐₗ of 10, 20, and 30g result in corresponding changes in the component unbalance. The experimental results indicate that the lowest unbalance of 2.6 mm/s occurs at a trial mass of 10 grams and angular positioned at 0^0^–180^0^. This is because the centrifugal forces of m2 and m6 canceled out. On the contrary, the unbalance G obtained the highest value of 42.3 mm/s at the mass 30g and angular positioned at 0^0^–0^0^ because of the resonance of centrifugal forces of m2 and m6.

Table 5: Unbalance G and displacement at 2000 rpm for various trial masses and angular positions.

In addition, Fig 12 and Table 6 show a decreasing trend of displacement amplitudes by ±0.098 mm, ± 0.079 mm, and ± 0.054 mm as the load position shifts from 0^0^–0^0^ to 0^0^–90^0^ and then to 0^0^–180^0^, respectively.

Table 6: Displacement measurement results under 20g load at 2000 rpm, with different load angular positions: 00–00, 00–900, and 00–1800.

Horizontal displacement of node 4 under 20 g trial mass and speed of 2000 rpm, for the different angular positions: 00–00, 00–900, and 00–1800.

4.3 Motional orbit

To evaluate the influence of the unbalanced mass position on the rotor’s motion trajectory, a series of experiments was conducted at a rotational speed of 2000 rpm. A trial mass 30g was mounted on the two discs of the rotating shaft, with various combinations of phase angular positions between the loads. Three unbalanced mass distribution cases were investigated for three angular positions 0^0^–0^0^, 0^0^–90^0^ and 0^0^–180^0^. The motion orbit of the shaft center at the speed n3 = 2000 rpm is represented in Fig 13. The orbits for load-free rotation and the unbalanced masses in the opposite angular directions 0^0^–180^0^ obtains the smallest value, while the orbit for angular position 0^0^–0^0^ is most unstable. This is because of the cancelation or resonance of centrifugal forces of the masses, well agreeing with the experimental displacements in Table 5. Furthermore, the data indicate that the oscillation observed with a trial mass of 30g for the angular position 0^0^–0^0^ closely resembles the shaft behavior at its first critical speed.

Motional orbit of the shaft center at 2000 rpm for various angular positions.

It can be observed that changes in the angular positions of the trial mass 0°–0°, 0°–90°, and 0°–180° strongly affect the vibration amplitude, thus the shape of the shaft center orbit and well agrees to the expected behavior predicted based on the influence parameters. The greater amount of unbalance and their positions concentration will generate stress accumulation and micro-cracks, resulting in shaft failure.

5. Conclusions

The following conclusions are made:

a. The effects of such key parameters as rotation speed, eccentricity due to unbalanced and initial phase angle during operation were simulated. The experimental equipment was built and experimental results were compared with the simulated data.b. The vibration becomes severely unstable when rotation speed approaches the critical speed, as shown in Tables 3 and 4. The displacement amplitude significantly increases with the amount of unbalance, as shown in Table 5.c. The initial phase angle of the trial mass position has great influence on the vibration of the rotating shaft, thus the fatigue bending strength and performance of the component, as shown in Tables 5 and 6 and Fig 13. By determining the vibration and unbalance conditions of the shaft component, the fatigue limitation can be predicted to ensure the safety and long-term stability of the rotor system.d. For operation under the critical speed, the unbalance analysis using the Newmark-β method well agrees with the experimental result, as represented in Table 2, showing that the Newmark-β method is a reliable and stable approach for analyzing the behavior of rotor-bearing systems. This allows rapid and high computational efficiency to predict the shaft operation with fast convergence time.

Further research may be proceeded on the following issues:

i. Unbalance analysis for double-phase and composite materials.ii. Analysis for materials with thin films, coating layer or the surface layer with residual stress.iii. Unbalance analysis in high temperature conditions.

Supporting information

S1 DataExperimental data for Figs 8, 11, 12 and 13.(XLSX)