Rigidflexible coupling dynamics simulation of planetary gear transmission based on MFBD
Chiyu Hao^{1} , Guangbin Feng^{2} , Huagang Sun^{3} , Haiping Li^{4}
^{1, 4}Mechanical Engineering College, Shijiazhuang, 050003, China
^{2, 3}New Technology Application Institute, Shijiazhuang, 050003, China
^{1}Corresponding author
Journal of Vibroengineering, Vol. 19, Issue 8, 2017, p. 56685678.
https://doi.org/10.21595/jve.2017.18208
Received 26 January 2017; received in revised form 12 June 2017; accepted 20 June 2017; published 31 December 2017
JVE Conferences
This paper deals with the problem of getting dynamic characteristics in complex mechanical multibody system. Based on the MFBD (MultiFlexibleBody Dynamics) technology, the rigidflexible coupling dynamic simulation method is proposed, and then the method is applied to the planetary gear transmission. The results show that the dynamic stress distribution of planetary gear can be obtained to determine the dangerous location, and the dynamic response characteristics are more obvious. Then the simulation model of planetary gear transmission with broken tooth is established, the fault feature extraction in timefrequency domain is carried out using the acceleration signal. In addition, industrial data is also used to validate the effectiveness of the proposed method.
Keywords: MFBD, rigidflexible coupling simulation, planetary gear transmission, RecurDyn.
1. Introduction
Acquisition of state information is a key aspect in PHM (Prognostics and Health Management), but it is difficult to acquire the dynamical characteristics of key components in the largescale, complex equipment [1, 2]. Development of virtual prototype technology has greatly improved the status. Compared with the traditional multirigid body model, the description of component status in multiflexible body system can be more accurate, but the high nonlinear and complex coupling characteristics also lead to difficulties of the dynamic modeling and solving, therefore it remains need to be further studied on the theory and application [35].
The traditional research of rigidflexible coupling modeling is mainly on the modal reduction method which the stiffness of flexible structure is expressed by the modal by finite element calculation in advance, but the method is limited because it can’t solve the large deformation and the nonlinear problems [68]. Based on the MFBD technology, the modal reduction method is extended by using the node method which put the finite element analysis and multibody dynamics analysis in the same solver, therefore it can directly obtain structure characteristics and effectively solve large deformation and nonlinear problems [911].
The remaining sections of this paper are organized as follow. In Section 2, the theory of rigidflexible coupling dynamics was illustrated. In Section 3, the framework of simulation method was proposed and the rigidflexible coupling model of the planetary gear transmission is established. In Section 4, the simulation results were analyzed. Section 5 validated the proposed method using industrial data. Finally, the conclusions were drawn in Section 6.
2. Rigidflexible coupling dynamics theory [12, 13]
2.1. Multibody dynamics theory overview
The object speed in absolute coordinate system can be defined as Eq. (1):
wherein, ${n}_{c}$ is the number of coordinates in the absolute coordinate system.
The object speed in relative coordinate system can be defined as Eq. (2):
wherein, ${n}_{r}$ is the number of coordinates in the relative coordinate system.
Thus, the Cartesian speed of all objects can be denoted as Eq. (3):
The coefficient matrix $B$ can be denoted by relative coordinates.
The system motion equations can be derived by the first kind Lagrange equation as Eq. (4):
wherein, $M$ is mass matrix, $\ddot{q}$ is acceleration vector, $\dot{B}$ is derivative of coefficient matrix $B$, $\mathrm{\Phi}$ is cutting joints constraint, $\lambda $ is Lagrange multipliers, $Q$ is force vector.
2.2. Motion equations of flexible bodies
As shown in Fig. 1, position vector of any point P can be defined as Eq. 5:
wherein, $r$ is $P$ point vector in inertial coordinate system; ${r}_{0}$ is the origin vector of the relative coordinate in the inertial coordinate system; $A$ is direction cosine matrix; ${\mathrm{s}}_{p}$ is $P$ point vector in relative coordinate system while flexible body is undeformed; ${u}_{P}$ is relative deformation vector, the following equation shows that ${u}_{P}$ can be denoted by the modal coordinates as Eq. (6):
wherein, ${\mathrm{\Phi}}_{p}$ is deformation modal matrix, ${q}_{f}$ is generalized coordinates of deformation. Therefore, the velocity vector and acceleration vector of any point at Flexible body can be denoted as Eq. (78):
Fig. 1. Deformation model of flexible body
2.3. Dynamics equations of flexible body
Considering position, direction and modal of the node $P$ with the deformation, generalized coordinates of flexible body can be defined as Eq. (9):
Dynamics equation of flexible body is derived by the Lagrange equation as Eq. 10:
wherein, $\mathrm{\Psi}$ is constraint equation, $\lambda $ is Lagrange multiplier corresponding to the constraint equation, $\xi $ is generalized coordinates as Eq. (9) defined, $Q$ is generalized force projected to $\xi $, $L$ is defined as $L=TW$, $T$ is kinetic energy and $W$ is potential energy, $\mathrm{\Gamma}$ is energy loss function, the differential equations of motion can be denoted as Eq. (11):
wherein, $\xi $, $\dot{\xi}$, $\ddot{\xi}$ is generalized coordinates time derivatives of flexible body, $M$, $\dot{M}$ is the mass matrix of the flexible body and its time derivative. $K$ is stiffness Matrix, $D$ is damping matrix. The differential Eq. (11) is solved and imported into the Eq. (78), all generalized velocity and acceleration can be obtained.
3. Method and case simulation
3.1. Method and application object
In this paper, based on 3D modeling software SolidWorks, finite element analysis software ANSYS and multibody dynamics software RecurDyn, a method to build rigidflexible coupling simulation model is proposed, the framework is constructed as shown in Fig. 2, the model could be effectively built with interface between these software programs.
Fig. 2. Framework of rigidflexible coupling simulation
Fig. 3. Working principle of integrated transmission device
As shown in Fig. 3, the integrated transmission device is dualpower flow form, the power is input from the engine, the directflow power flow is transmitted to the ring gear of planetary gear transmission through the gear shift mechanism; while the steering power flow is transmitted to the sun wheel of planetary gear transmission through the steering control mechanism, and the two power streams converge and output at the planetary gear transmission to drive the vehicle to realize directdriving and steering working condition.
Planetary gear transmission is a key subsystem in the integrated transmission device and its operation reliability directly affect the whole system performance. But in actual using process, due to its complex structure, heavy workload, closed and bad operating environment, it is found of highly failure such as wear, ablation and other faults. Therefore, the dynamic characteristics, fault mechanism, fault diagnosis and feature extraction of the planetary gear transmission should be paid more attention [1416].
In order to get the dynamic response characteristics considering the effect of elastic deformation. In the case, based on the MFBD technology, the rigidflexible coupling dynamic model of the planetary gear transmission is studied by using the method proposed.
3.2. Calculation of planetary gear transmission
In the directdriving conditions, the DOF of planetary gear transmission system is 1, and in the steering conditions, the DOF of planetary gear transmission system is 2. In this paper, the dynamic characteristics of the steering status are simulated. The ring gear and sun gear are used as input and the planet carrier is output.
Fig. 4. Epicyclic gear structure
In essence, the planetary gear transmission studied is an epicyclic gear structure as shown in Fig. 4, while the 1th component is the sun gear, the 2th component is the planetary gear, the 3th component is the ring gear, the component $H$ is the planet carrier.
According to actual working conditions, while the engine input speed is 162.2 rad/s, the speed of the directdriving output (the ring gear rotation speed) is 80.6 rad/s, and in order to make the results more convincing, the input rotation speed (the sun wheel rotation speed) is set in two cases, the first case is 10 rad/s, the second case is 30 rad/s. According to the transmission principle and the gear parameters, the follow results can be calculated as Table 1 shows.
Table 1. Calculation values of different case
Case 1

Case 2


Ring gear

${\omega}_{31}={\omega}_{32}=$ 80.6 rad/s


Sun wheel

${\omega}_{11}=$ 10 rad/s

${\omega}_{12}=$ 30 rad/s

Planet carrier

${\omega}_{H1}=$ 52.8 rad/s

${\omega}_{H2}=$ 46.8 rad/s

${f}_{H1}=$ 8.42 Hz

${f}_{H2}=$ 7.45 Hz


Planetary gear

${\omega}_{21}=$ 151.7 rad/s

${\omega}_{22}=$ 167.5 rad/s

${f}_{21}^{H}={f}_{21}{f}_{H1}=$ 15.74 Hz

${f}_{22}^{H}={f}_{22}{f}_{H2}=$ 19.22 Hz


Planetary system

${f}_{m1}=$ 220.3 Hz

${f}_{m2}=$ 269.1 Hz

3.3. Multirigidbody model
After adding constrains and contacts of the components, the driving speed of ring gear and sun gear is set as Table 1 shows, the load torque is set to the output planet carrier as 6000 N·m, the simulation time is set as 1 second, the simulation step is set as 2000.
Fig. 5. Multirigidbody model
Fig. 6. Speed values of planetary gear transmission
a)
b)
The speed values of simulation can be seen in Fig. 6, both in the different case, each component gradually increases to stable condition at about 0.02 s and the rotational speed of the planet carrier and the planetary gear fluctuate in a small range between the calculated value. Compared with the calculation values, the accuracy of the multirigidbody model is proved.
3.4. Rigidflexible coupling model
The finite element model of the planetary gear is built in ANSYS as Fig. 7 shows, then the corresponding model file is generated and exported to RecurDyn.
After importing the finite element model of planetary gear to replace original member, the contact relationship and the original constraints disappear so that the finite element model should be modified [9], as shown in Fig. 8.
(1) A set of nodes and rigid FDR unit is set to replace rotating pair between the planetary gear and planet carrier.
(2) The surface of planetary gear tooth surface is defined as “patch” unit, then the contacts between planetary wheel and contacted parts are built using the finite element flexible body contact module (Surface to Fsurface) in FFLex.
Then contact force parameters are added and then the rigidflexible coupling model is established as shown in Fig. 9.
Fig. 7. planetary gear finite element model
Fig. 8. Modification in FFLex
Fig. 9. Rigidflexible coupling model of planetary gear transmission
4. Simulation results analysis
4.1. Analysis in normal state
The dynamic stress distribution of planetary gear can be obtained in the postprocessing module in RecurDyn. As shown in Fig. 10, the maximum stress of the planetary gear occurs in the red area which is in the end of the planetary gear wheel tooth.
Fig. 10. Stress distribution of planetary gear with max stress
The contact force between planetary gear and sun gear in time domain is shown in Fig. 11. The dynamic contact force in the rigidflexible coupling model is larger than in the rigid body model, but the calculation values show the average dynamic contact force in the rigidflexible coupling model is less than in the rigid body model. Fig. 12 is comparison of the contact force in frequency domain, both the meshing frequency of planetary system (${f}_{m1}=$ 220.3 Hz, ${f}_{m2}=$ 269.1 Hz) and the frequency multiplier appear in the figure, and compared with the rigid model, the response of the flexible body model in time domain and frequency domain is more obvious.
Compared with the rigid model, the results show that the dynamic stress distribution and dangerous position of the component can be obtained in the rigidflexible coupling model, and the response characteristics in time domain and frequency domain are more obvious. Therefore, the proposed method can reflect the dynamic response characteristics better.
Fig. 11. Comparison of contact force in time domain
a)
b)
Fig. 12. Comparison of contact force in frequency domain
a)
b)
4.2. Analysis in fault state
Fig. 13 shows the rigidflexible coupling model with broken tooth. According to the previous analysis, the maximum stress occurs in the end of the planetary gear wheel tooth. In actual using process, the gear tooth is prone to fracture after repeated fatigue deformation.
Fig. 13. Establishment of planetary gear with broken tooth
By using the proposed method, the rigidflexible coupling model with broken fault is analyzed. In engineering practice, it is difficult to test the contact force between the components so that the acceleration signal is used to analysis. As shown in Fig. 14, the acceleration timedomain signal on planetary gear is compared. As Fig. 1516 show is the frequencydomain signal comparison after FFT transform.
Fig. 14. Comparison of acceleration in time domain
a)
b)
Fig. 15. Comparison of acceleration in frequency domain
a)
b)
Fig. 16. Comparison of acceleration in partly frequency domain
a)
b)
In the timedomain comparison (Fig. 14), there are larger loads after the broken tooth failure compared with the normal state, but its periodicity is not obvious. In the frequency domain comparison (Fig. 16), the meshing frequency of planetary system (${f}_{m1}=$ 220.3 Hz, ${f}_{m2}=$ 269.1 Hz) can be seen in normal state, and the sidebands (${f}_{21}^{H}=$ 15.74 Hz, ${f}_{22}^{H}=$ 19.22 Hz) are irregularly distributed on the sides. However, due to the aliasing of the frequency, the frequency multiplier and sidebands are not obvious. The amplitude of main frequency and its frequency multiplier increase while the occurrence of broken teeth, and the frequency components and amplitude also increase significantly on both sides of meshing frequency.
5. Experimental validation
On the viewpoints described above, an industrial data is analyzed to validate the proposed method, and the data is acquired from wind turbines. The type of wind turbine is FL600 and it has been reaching to examine and repair deadline. The detail data of can be seen in document [17]. The wind turbine and its internal structure can be seen in Fig. 1718. In the unpacking process, the planetary gear is found of damage as the Fig. 19 shows. Then the acceleration signals under fault state are compared with the normal state to validate the effective of the proposed method.
According the parameters and test conditions, the gearmesh frequency of planetary system is 203.9 Hz, and the gearmesh frequency of is 40.6 Hz. The data in normal and fault state is analyzed in time and frequency domain, the results are as Fig. 20 show. In time domain, there are larger value of the shocktype vibrations under the fault state; and in the frequency domain, the frequency composition and amplitude increase obviously. The results prove the effectiveness of the method this paper proposed.
Fig. 17. FL600 wind turbine gearbox
Fig. 18. The internal structure
Fig. 19. The damaged planetary gear
6. Conclusions
In this paper, the rigidflexible coupling dynamics simulation method was proposed based on the MFBD technique, and then the planetary gear transmission was analyzed by using the method. The simulation results showed that that the dynamic stress distribution could be obtained to determine the dangerous location and the dynamic response characteristics of rigidflexible coupling model were more obvious than the rigid model.
Fig. 20. Comparison in time and frequency domain
a) Comparison of acceleration in time domain
b) Comparison of acceleration in frequency domain
c) Comparison of acceleration in partly frequency domain
Then a fault model with the broken gear was established, the feature of the broken tooth of the planetary transmission gear system was extracted by using the timefrequency domain analysis of acceleration. In time domain, there were larger value of the shocktype vibrations; in the frequency domain, the amplitude modulation and frequency modulation were generated near the meshing frequency, the frequency components and amplitude increased obviously on both sides. Additionally, industrial data was used to validate the applicability of the proposed method to the practical environment. The results showed that the method this paper proposed still has good performance.
References
 Chen Z. S., Yang Y. M., Hu Z. A technical framework and roadmap of embedded diagnostics and prognostics for complex mechanical systems in prognostics and health management systems. IEEE Transactions on Reliability, Vol. 61, Issue 2, 2012, p. 314322. [Search CrossRef]
 Vogl G. W., Weiss B. A., Helu M. A review of diagnostic and prognostic capabilities and best practices for manufacturing. Journal of Intelligent Manufacturing, 2016, p. 117. [Search CrossRef]
 Zheng E., Zhu R., Zhu S., et al. A study on dynamics of flexible multilink mechanism including joints with clearance and lubrication for ultraprecision presses. Nonlinear Dynamics, Vol. 83, Issues 12, 2016, p. 137159. [Search CrossRef]
 Ren Z. S., Yang G., Wang S. S., et al. Analysis of vibration and frequency transmission of high speed EMU with flexible model. Acta Mechanica Sinica, Vol. 30, Issue 6, 2014, p. 876883. [Search CrossRef]
 Valembois R. E., Fisette P., Samin J. C. Comparison of various techniques for modelling flexible beams in multibody dynamics. Nonlinear Dynamics, Vol. 12, Issue 4, 1997, p. 367397. [Search CrossRef]
 Ki K. O., Kang M. K. Convergence acceleration of iterative modal reduction methods. AIAA Journal, Vol. 39, Issue 1, 2015, p. 134140. [Search CrossRef]
 Barbone P. E., Dan G., Patlashenko I. Optimal modal reduction of vibrating substructures. International Journal for Numerical Methods in Engineering, Vol. 57, Issue 3, 2003, p. 341369. [Search CrossRef]
 Rumpler R., Legay A., Deü J. F. A substructuring FE model for structuralacoustic problems with modalbased reduction of poroelastic interface. Computers and Structures, Vol. 89, Issues 2324, 2011, p. 22332248. [Search CrossRef]
 Moon I. D., Yoon H. S., Oh C. Y. A flexible multibody dynamic model for analyzing the hysteretic characteristics and the dynamic stress of a taper leaf spring. Journal of Mechanical Science and Technology, Vol. 20, Issue 10, 2006, p. 16381645. [Search CrossRef]
 Yang Q. J., Li D. N., Kong L. L., et al. Research on MFBD modeling method for a gantry machining center beam components. International Conference on High Speed Machining, 2011, p. 487492. [Publisher]
 Türtscher A., Faßbender F., Kelichhaus D. I. T. Simulation of dynamic stresses including flexible contacts using MFBD technology. SAE Technical Papers, 2006. [Search CrossRef]
 Deng F., He X., Liang L., et al. Dynamics modeling for a rigidflexible coupling system with nonlinear deformation field. Multibody System Dynamics, Vol. 18, Issue 4, 2007, p. 559578. [Search CrossRef]
 Cardona A., Geradin M. A beam finite element nonlinear theory with finite rotations. International Journal for Numerical Methods in Engineering, Vol. 26, Issue 11, 1988, p. 24032438. [Search CrossRef]
 Chen Z., Zhu Z., Shao Y. Fault feature analysis of planetary gear system with tooth root crack and flexible ring gear rim. Engineering Failure Analysis, Vol. 49, 2015, p. 92103. [Search CrossRef]
 Wang J., Li R., Peng X. Survey of nonlinear vibration of gear transmission systems. Applied Mechanics Reviews, Vol. 56, Issue 3, 2003, p. 222227. [Search CrossRef]
 Chen X. H., Cheng G., Shan X. L., et al. Research of weak fault feature information extraction of planetary gear based on ensemble empirical mode decomposition and adaptive stochastic resonance. Measurement, Vol. 73, 2015, p. 5567. [Search CrossRef]
 Li H., Zhao J., Zhang X., et al. Gear fault diagnosis and damage level identification based on Hilbert transform and Euclidean distance technique. Journal of Vibroengineering, Vol. 16, Issue 8, 2014, p. 41374151. [Search CrossRef]