Fatigue Reliability Prediction Method of Large Aviation Planetary System Based on Hierarchical Finite Element

Abstract

The reliability of planetary equipment determines the economic affordability and service safety, to a large extent, for a helicopter transmission system. However, with the continuous improvement of the progressiveness and large-scale degree of new aviation planetary equipment, the contradiction between reliability design indexes and R&D economy is also gradually highlighted. This paper takes the large aviation planetary system as a research object, aims to accurately evaluate the system reliability level formed in design processes, and deeply excavates the inherent characteristics of the planetary system in functional realization and builds a system fatigue reliability evaluation model accordingly. An advanced hierarchical finite element technology is used to calculate dangerous tooth load histories under the influence of system global elastic behavior, and the tooth probability fatigue strength is obtained through the gear low-cycle fatigue test and life distribution transformation method, so as to provide economic load and strength input variables, respectively, for the reliability model. This prediction method can provide targeted structural optimization guidance in the development and design of the large aviation planetary system and significantly reduce the cost of reliability index realization for this kind of large-scale, high-end equipment in design iteration processes.

1. Introduction

The heavy-lift helicopter is military and civilian general-purpose strategic equipment related to core national interests, as well as an important symbol of the aviation science and technology levels and even the comprehensive national strength. The high-power transmission system technology is a core technology field for improving the performance of the heavy-lift helicopter by reducing its noise and vibration levels and controlling its life-cycle cost. Both the United States and Russia have listed the reliability and economic affordability of the high-power transmission system as key technical indicators in their respective R&D plans for the advanced heavy-lift helicopter, and they have put forward specific low maintenance design requirements for the transmission system [1]. Among the largest number of heavy-lift helicopters currently in service, the large aviation planetary mechanism, as the basis and core of their transmission systems, determines the scientific and technological levels of the transmission systems to a great extent and is one of the bottlenecks restricting the development of transmission system technologies in the heavy-lift helicopter [2].

As a deceleration terminal directly connected with the main rotor, the large aviation planetary mechanism is a power transmission link with the worst load environment and highest strength requirements in the heavy-lift helicopter transmission system [3]. With the continuous improvement of performance standards, such as the power density level and dry running capacity of the advanced heavy-lift helicopter main reducer in developed countries, the speed and bearing limits of the large aviation planetary mechanism continue to break through, and the structure and power sharing forms become more complex. These increase the complexity and uncertainty in structural strength design and analysis for the planetary system and lead to a higher risk of structural failure when the new generation of heavy aviation planetary equipment is pursuing lightweight [4]. The transmission system is one of the three key dynamic systems (engine system, transmission system, and rotor system) in a helicopter. Considering weight or space, a redundancy design for the transmission system cannot be realized. Therefore, the reliability of its core link will directly determine the service safety and life-cycle cost of the helicopter [5].

High reliability and long life have become the core technology development directions for large aviation planetary equipment in the future [4,6]. Compared with small and medium-sized helicopters, the planetary mechanism in the heavy-lift helicopter has larger parts and components, higher requirements for structural mechanics indexes, and more stringent design standards for important geometric elements. If only relying on the design level and experience parameters of existing related products to develop a larger-sized isomorphic mechanism, it may directly lead to the failure to guarantee the reliability and durability for large aviation planetary equipment. In order to meet performance matching requirements of the new generation heavy-lift helicopter main reducer, a system reliability index prediction method that adapts to the structure and operation characteristics of the large-scale aviation planetary equipment is urgently needed. This method will enable developers to accurately evaluate the reliability level of products in the early stage of design, and then provide targeted reference data for the structural optimization and reliability improvement of products. This will significantly reduce the cost of expensive large-scale aviation planetary mechanism products for the realization of reliability indicators in development iteration processes.

In a gear transmission system, the failure of any gear or tooth will affect the transmission capacity of the entire system, so it is generally assumed that the gear transmission system is a series system, with the gear or tooth as basic functional units in the reliability analysis [7,8]. The loads on different individual gears in the system, or even different teeth on the same gear, may be different, but they all have a certain mathematical relationship with the system input power. This load correlation and general load randomness make the failure behavior of each unit no longer independent, so it is unreasonable to simply think that the reliability of the series system is equal to the product of each unit’s reliability. Especially for the large aviation planetary gear train with more complex structure, this “unit independent failure assumption” will even lead to serious errors in the task of system reliability analysis [9]. In addition, most relevant research works simplify the gear transmission system into a general series system or series parallel hybrid system, which completely fails to reflect the special attribute that the functional form of this kind of system changes with time [10]. Different from the “chain type” series system, in which each potential failure unit is loaded at the same time, the energy transfer in the gear system is completed by the alternating loading of each tooth. Therefore, the gear system is not only a series system in the sense of spatial structure, but also a series system in the time domain. The structure and motion forms of the large aviation planetary gear train are more complex, and the teeth meshing timing characteristic will be more prominent. Therefore, the important influence of this behavior mechanism needs to be fully considered in a strict and effective reliability evaluation.

Tooth root bending fatigue strength is one of the most important strength check indexes in high-speed and heavy-duty gear systems. For the planetary mechanism in heavy-lift helicopters, its extreme load conditions and harsh reliability requirements will put forward a higher standard for tooth root bending fatigue resistance. Gear researchers are always trying to find a trade-off solution between obtaining accurate results from gear stress analyses and high-efficiency computing. For both factors, accuracy on the one hand, and high-efficiency computing on the other hand, the accuracy and efficiency of tooth root bending stress calculation determines the validity of fatigue reliability assessment of large aviation planetary wheel systems. At present, the main calculation method of tooth root bending stress is to use general finite element tools [11,12,13], which are flexible in technical processing and are not limited by input conditions, such as geometric characteristics and material properties, and the analysis results are relatively comprehensive (depending on software post-processing ability). Several studies have investigated how to better analyze gear stresses using finite element models. A strategy for the mesh refinement of finite element models of gear drives, based on the application of multi-point constraints, has been applied. It allows for a considerable reduction of the computational cost in the analysis of the whole cycle of meshing. In the model, the refined area mesh and the non-refined area mesh were connected by multi-point constraint (MPC) at the same time, in order to save time on the FEM solution, on the premise of ensuring the accuracy of model analysis [14,15,16]. Chen investigated the effects of a gear tooth spalling in a helical gear by finite elements. The results showed that spalling also caused the rapid increase of tooth root stress in the spalling meshing area [17]. Franco Concli analyzed early crack propagation in single tooth bending fatigue by combining finite element analysis with critical plane fatigue criterion. The results showed that the crack propagation direction at the ${\rho }_{fp}$ did not follow the plane of maximum alternating shear stress, but the plane of maximum damage parameter [18]. However, the general finite element method has high computational cost in model setting and solution operations and is usually only suitable for isolated solutions of gear parts or several teeth, but it is difficult to perform global operations at the system-level. Moreover, it is difficult to predict tooth bending conditions or add corresponding boundary conditions to the model, especially for thin-walled rim gears or gears directly installed on bearings (such as the planetary gears with rotation and revolution characteristics). In addition, the mesh density on the tooth surface cannot be refined to the micron level to reflect the micro-geometric state of the tooth surface, so it is difficult to analyze the influence of tooth surface modification on tooth root stress by the general finite element method.

This paper takes the large aviation planetary system as a research object, aims to accurately evaluate the system reliability level formed in design processes, and deeply excavates the inherent characteristics of the planetary system in functional realization and builds a system fatigue reliability evaluation model accordingly. Advanced hierarchical finite element technology is used to calculate tooth dangerous load histories under the influence of system global elastic behavior, and the tooth probability fatigue strength is obtained through the gear low-cycle fatigue test and life distribution transformation method, so as to provide economic load and strength input variables, respectively, for the reliability model. This prediction method can provide targeted structural optimization guidance in the development and design of the large aviation planetary system and significantly reduce the cost of reliability index realization for this kind of large-scale, high-end equipment in the design iteration processes. At the same time, it can provide important reference data for the first renovation period of relevant finalized products and then provide the technical support for their economic guarantee in the whole life-cycle.

2. Tooth Root Stress Calculation Based on Hierarchical Finite Element Technology

2.1. Principle Analysis of Hierarchical Finite Element Technology

Most tooth root stress calculation standards are analyzed and processed according to the geometric characteristics of tooth root, referring to empirical data, and considering working conditions. The biggest disadvantage of these methods is that they are limited to the specified geometries, materials, and working conditions. For the unconventional tooth geometries commonly existing in the aviation field, the calculation results obtained by these standard methods may not be accurate. In particular, the urgent need for lightweight components for the new generation of large aviation planetary equipment makes the components in the system widely adopt lighter and thinner structural design forms. The introduction of such a large number of flexible features means that the transmission shafts, support frames, gears, and other critical transmission parts will undergo significant elastic deformation under heavy loads. The resulting stiffness problem makes it necessary to fully consider the nonlinear mechanical behaviors for the whole system, such as elastic deformation and meshing misalignment, in order to accurately calculate the tooth root stress.

In order to solve the contradiction between the calculation accuracy and calculation cost of the general finite element method, when facing the mechanical simulation analysis task for the large complex gear system, an advanced hierarchical finite element method for root stress analysis of the large aviation planetary mechanism is proposed in this paper. When calculating the tooth root stress, it is directly based on a detailed tooth 3D finite element sub-model, so it does not need to rely on empirical data and industry experts. In the system-level model, the dynamic load line on the tooth surface was obtained by quasi-static static analysis and tooth surface micro-geometric analysis, and then it was loaded on the tooth surface of the finite element sub-model as a load boundary condition. At the same time, system elastic deformation results, including the deformation behavior of the gear rim and tooth, were extracted and loaded into the finite element sub-model as a displacement boundary condition. In this way, the calculation results of tooth root stress could include the influence of the elastic deformation (including deformations, such as rim distortion and tooth bending and so on), meshing misalignment, and tooth surface micro-geometry in the whole system.

The analysis and solution for the system-level model took the analytical algorithms in the international standard as the default calculation basis, and they have a strong auxiliary calculation ability in the tooth surface micro-geometric analysis. At the same time, the general finite element method was used as a backup advanced analysis method, and the introduction of finite element components into the system-level model made the system elastic deformation results more reliable. In the system-level model, there is no need to perform detailed root stress calculation, which greatly saves the calculation cost for system-level analysis. In addition, tooth surface micro-geometric analysis was run in the system-level model, so it is no longer necessary to include micro-geometric parameters in the secondary finite element sub-model. In the quasi-static analysis results of the system-level model, the physical effects, such as rotation, centrifugation, and thermal expansion of transmission components, were considered, so the boundary conditions imposed on the secondary sub-model were relatively simple. In the face of advanced simulation and analysis tasks for a large aviation planetary system, only considering the convenience of modeling and boundary condition setting, the computational efficiency of the hierarchical finite element method will be much higher than that of the general finite element method. In addition, different from the commercial finite element software using a nonlinear equation solver, the hierarchical finite element method uses an improved simplex solver to ensure convergence within a set number of iterations. Although the total number of freedom degrees in the two-level finite element model may be very large, this hierarchical analysis approach keeps the amount of CPU time and memory required for analytical calculation within the capabilities of a personal computer.

It can be seen that the hierarchical finite element method is a more accurate, fast, and easy method for calculating tooth root stress, and it is especially suitable for the strength checking and analysis tasks of large aviation gear systems. In this paper, a high-fidelity mechanical simulation model was constructed by the hierarchical finite element method, based on the structure and material performance parameters of a certain type of large aviation planetary mechanism products. Based on this, the dangerous tooth root stress under system quasi-static elastic mechanical behavior was calculated, and the effective load input variables were provided for the system reliability evaluation model.

2.2. System-Level Elastic Mechanical Behavior Simulation Analysis

2.2.1. Overall Configuration of System Model

A semi-analytical finite element technique was used to construct the system-level elastic mechanical simulation model of the large aviation planetary mechanism and accurately evaluate the elastic deformation of the large thin-walled parts, the meshing misalignment between teeth, and the micro-modification on tooth surface, and then provide detailed load and displacement boundary conditions for the secondary sub-model of the tooth root stress analysis. The semi-analytical finite element method was used in the system-level analysis process. Although fine mesh densities were used for some large structural members, the effective combination of analytical and finite element calculations resulted in significant improvements in the computational efficiency at the system-level. The system model of the planetary mechanism and the structural parameters of the gear train are shown in Figure 1 and Table 1, respectively. The power flow path inside the system starts from the input shaft, passes through the sun gear and planet gear to a planet carrier, and finally, the planet carrier transmits the motion and power to a main rotor shaft, after a deceleration of 3.334 times. In addition, the rated working condition parameters of the planetary mechanism mainly include an input power of 5000 kW, input speed of 500 rpm, and working temperature of 70 °C.

2.2.2. Deformable Planet Carrier

An FE model and the geometry of a planet carrier are illustrated in Figure 2, which shows an example of an FE grid for a seven-planet system, with its supporting conditions simulated by lumped stiffness elements. Since there are no relative displacements between the axis of rotation of the planet pins and the planet nodes (planet centers), a fixed interface component mode synthesis method can be used to reduce the size of the carrier model. The stiffness of the springs is considered to be significantly higher than that of the bearings and pins, thus creating a rigid link in the radial directions between the contour nodes and planet centers. The classic reduction process of Craig and Bampton [19] was employed, in which the displacements were expressed, in terms of static and dynamic modes, as

(1)$\left[\begin{array}{l}{\mathsf{X}}_{\mathrm{PCint}}\\ {\mathsf{X}}_{\mathrm{PCbou}}\end{array}\right]=\left[\begin{array}{c}{\left[{\mathit{\varphi }}_{ \mathsf{D}}\right]}_{}\left[{\mathit{\varphi }}_{ \mathsf{S}}\right]\\ {\left[\mathbf{0}\right]}_{}{{}_{}}_{}\left[I\right]\end{array}\right]\cdot \left[\begin{array}{c}{\mathsf{q}}_{\mathrm{PC}}\\ {\mathsf{X}}_{\mathrm{PCbou}}\end{array}\right]$

with ${\mathsf{X}}_{\mathrm{PCint}}$ as the vector of the internal degrees of freedom, ${\mathsf{X}}_{\mathrm{PCbou}}$ as the vector of the degrees of freedom at the contour nodes, $\left[{\mathit{\varphi }}_{ \mathsf{D}}\right]$ as the truncated modal matrix of the fixed-interface planet-carrier structure, $\left[{\mathit{\varphi }}_{ \mathsf{S}}\right]$ as the static mode matrix, $\left[\mathsf{I}\right]$ as the identity matrix, $\left[\mathbf{0}\right]$ as the nil matrix, and ${\mathsf{q}}_{\mathrm{PC}}$ as the vector of the planet-carrier modal unknowns.

2.2.3. Planet Bearing Element

This connecting part was composed of seven lumped spring elements across the planet bore, in order to connect node O_j of planet j to the corresponding three contour nodes of the planet-carrier substructure (denoted as N₁, N₂, and N₃ in Figure 3). Planets were modeled as rigid disks with 6 degrees of freedom (DOFs), which are the infinitesimal generalized displacements superimposed on rigid-body motions and represented by screws of coordinates.

(2)$\left\{{\tau }_j\right\}\left\{\begin{array}{l}{\mathsf{u}}_{\mathsf{j}}\left(O_j\right)=v_j{\mathsf{S}}_{\mathsf{j}}+w_j{\mathsf{T}}_{\mathsf{j}}+u_j\mathsf{Z}\\ {\mathsf{\omega }}_{\mathsf{j}}={\phi }_j{\mathsf{S}}_{\mathsf{j}}+{\psi }_j{\mathsf{T}}_{\mathsf{j}}+{\theta }_j\mathsf{Z}\end{array}\right\}$

where S_j, T_j, and Z are the unit vectors of the frame fixed to the sun gear/planet j mesh (Figure 4).

Considering two points, P and Q, belonging to planet j, their displacements were expressed, using the shifting property of the moment of screw $\left\{{\tau }_j\right\}$, as

(3)$\begin{array}{l}{\mathsf{u}}_{\mathsf{j}}(P)={\mathsf{u}}_{\mathsf{j}}(O_j)+\mathit{P}{\mathit{O}}_j\times {\mathsf{\omega }}_{\mathsf{j}}\\ {\mathsf{u}}_{\mathsf{j}}(Q)={\mathsf{u}}_{\mathsf{j}}(O_j)+\mathit{Q}{\mathit{O}}_j\times {\mathsf{\omega }}_{\mathsf{j}}\end{array}$

Denoting K_v, K_w, and K_u the stiffness in the S_j, T_j, and Z directions between P and N₁, Oj, and N₂, and Q and N₃, and the strain energy stored in the spring element 1 connecting node N₁ and point P reads

(4)$U_1=\frac{1}{2}{\mathsf{X}}^T\left[K_u{\mathrm{V}}_{\mathrm{u}}\cdot {\mathrm{V}}_{\mathrm{u}}^T+K_v{\mathrm{V}}_{\mathrm{V}}\cdot {\mathrm{V}}_{\mathrm{V}}^T+K_W{\mathrm{V}}_{\mathrm{W}}\cdot {\mathrm{V}}_{\mathrm{W}}^T\right]\mathsf{X}=\frac{1}{2}{\mathsf{X}}^T\left[{\mathsf{K}}_1\right]\mathsf{X}$

where ${\mathsf{X}}^T=(u_j,v_j,w_j,{\phi }_j,{\psi }_j,{\theta }_j,u_{N1},v_{N1},w_{N1})$ is the vector containing the degrees of freedom associated with nodes O_j (6 DOFs) and N₁ (3 DOFs). The structural vectors V_u, V_v, and V_w are expressed, in terms of the planet width and its angular position ${\mathsf{\Phi }}_j$ (Figure 4), as

(5)$\begin{array}{l}{V}_u^T=⟨1,0,0,0,0,0,-1,0,0⟩\\ V_v^T=⟨0,\cos {\mathsf{\Phi }}_j,-\sin {\mathsf{\Phi }}_j,-\frac{b}{2}\sin {\mathsf{\Phi }}_j,-\frac{b}{2}\cos {\mathsf{\Phi }}_j,0,0,-1,0⟩\\ V_w^T=⟨0,\sin {\mathsf{\Phi }}_j,\cos {\mathsf{\Phi }}_j,\frac{b}{2}\cos {\mathsf{\Phi }}_j,-\frac{b}{2}\sin {\mathsf{\Phi }}_j,0,0,0,-1⟩\end{array}$

A similar procedure was used for the other spring elements attached to the same planet, thus leading to the 15 × 15 stiffness matrix of the planet j bearing element, which connected the 6 degrees of freedom at the planet center $(u_j,v_j,w_j,{\phi }_j,{\psi }_j,{\theta }_j)$ and those at nodes N₁(u_N1, v_N1, w_N1), N₂(u_N2, v_N2, w_N2), and N₃(u_N3, v_N3, w_N3).

2.2.4. Deformable Ring Gears

Ring gears and their supports were modeled by using FE analysis, combined with lumped parameter elements, to account for the ring-gear bearing stiffness and for the contribution of the teeth modeled as lumped masses. The ring-gear FE model was reduced using a component mode synthesis method based on the mode shapes of the undamped isolated structure. After separating the internal displacements, ${\mathsf{X}}_{\mathrm{RGint}}$, and the degrees of freedom at the boundary nodes, ${\mathsf{X}}_{\mathrm{RGbou}}$, i.e., on the root cylinder and in correspondence with the planet ring-gear contacts, the following approximation was used:

(6)$\left[\begin{array}{l}{\mathsf{X}}_{\mathrm{RGint}}\\ {\mathsf{X}}_{\mathrm{GDbou}}\end{array}\right]=\left[{\mathsf{\Phi }}_{\mathrm{N}}\right]{\mathsf{q}}_{\mathrm{RG}}$

with ${\mathsf{X}}_{\mathrm{RGint}}$ as the vector of the internal degrees of freedom, ${\mathsf{X}}_{\mathrm{RGbou}}$ as the vector of the degrees of freedom at the contour nodes, $\left[{\mathsf{\Phi }}_{\mathsf{N}}\right]$ as the truncated modal matrix of the undamped ring-gear structure, and ${\mathsf{q}}_{\mathrm{RG}}$ as the vector of the modal unknowns.

2.2.5. System Support Condition

In this analysis task, it was assumed that the reducer box had an absolutely rigid body; the input shaft was supported by two tapered roller bearings (TRB) with O-type layout, and their outer rings were rigidly connected to the box. All planet gears in the system were evenly distributed on the planet carrier along a circumferential direction. The double-row tapered roller bearing (DRTRB), with an X-type layout, was assembled inside each planet gear, and the bearing inner ring was rigidly connected with the planet shaft. A radial ball bearing (RBB) fixed the planet carrier and made it have a certain floating amount (its size and direction were controlled by the bearing clearance) to counteract the unequal load sharing among planet gears, and its outer ring was rigidly connected with the box. The structural forms and parameters of various main bearings in the system are shown in Figure 5 and Table 2, respectively.

2.2.6. Finite Element Component Modeling

The solid model creation of standard parts, such as gears and bearings, can be realized with the help of professional software, such as RomaxDesigner or GearTrax (GearTrax2020, Camnetics Inc., Oregon, WI, USA). The powerful parametric modeling function for the standard parts in the software can enable designers to complete this work easily. In addition, four key finite element components were also included in the planetary mechanism system model, which are the planet carrier, the ring gear, the planet gear, and the input shaft, integrated with the sun gear. Based on the advantages of 3D modeling software, such as NX (UG2020, Siemens PLM Software, Plano, TX, USA) and SolidWorks (SW2021, Dassault Systemes, Concord, MA, USA), and general finite element software, such as HyperMesh (HyperMesh2020, Altair Troy, MI, USA) and ANSYS (ANSYS 19.0, ANSYS Corporation, Pittsburgh, PA, USA), a high-quality parametric modeling function for finite element components can be realized, and the core technology implementation can refer to the literature [20,21]. The final modeling results and mesh quality parameters are shown in Table 3. It should be emphasized that creating complete teeth features in the gear finite element models can make the system-level flexibility analysis results cover the accurate mechanical behavior characteristics, such as the rim distortion and tooth bending deformations, which can provide more effective displacement boundary conditions for the secondary sub-model of the tooth root stress calculation. Additionally, it is now possible to keep CPU requirements down to about 8 s per time step and memory needs down to 128 MB. This has been made possible by accepting programming complexity, in exchange for an increase in speed.

2.2.7. Setting Node Connection

The above-mentioned solid models, such as the gears and bearings, as well as finite element models, such as the input shaft and planet carrier, were imported into the RotationMaster (RM) software platform, respectively. It can seamlessly integrate a variety of CAD and CAE software, and it is an advanced simulation platform for the comprehensive analysis and calculation of large and complex transmission systems. It has been widely used in related engineering technology fields. The imported solid models and finite element models were positioned and assembled in a unified global coordinate system, and the node connection relationships between the two types of models were established. We comprehensively adjusted the control parameters, such as the search criteria and selection modes, and used RM’s tolerance search technology to screen the connection node groups between the finite element models and the solid models. The analytical schematic diagram of two search criteria in the technology is shown in Figure 6.

The two semicircles in the figure show the geometric dimensions of the cross-sections of the hollow and solid cylinders, respectively, which will be used to define the search spaces corresponding to the two search criteria. Within the range of the specified cylinder height, “shell search” will search for all the nodes between the inner and outer walls of the hollow cylinder and on the walls. This search criterion is suitable for the node connection setting between the solid model of a bearing inner ring and the finite element model of a shaft. “Solid search” will search out all the nodes inside the solid cylinder, including its outer surface. This search criterion is suitable for establishing the node connection relationship between finite element teeth and their meshing lines. RM can automatically generate the pitch circle diameter (PCD), according to the geometric characteristics at the node connection location. The tolerance factor (TF) was used to adjust the diameter parameters of the inner and outer walls of the cylinders to optimize the connection attributes, such as the number of connection nodes and the connection stiffness. In addition, RBE2 element type (detailed in the NASTRAN help file) can bundle the selected connection nodes into a rigid node group, which means that all the nodes will have the same displacement parameters, so it can simplify the model and improve calculation efficiency. RBE3 element type distributes loads evenly across the selected connection nodes. It reduces the stiffness of the connection area by increasing the flexibility around the nodes, so that the elastomer deformation behavior obtained from the analysis will be more accurate, but the calculation cost will also increase significantly. Table 4 shows the node connection results and corresponding node parameters between the sun gear teeth on the input shaft and the meshing lines (pitch circle position), between the planet gear and the outer ring of a tapered roller bearing, between the mounting holes on the ring gear and the reducer box, and between the pin holes on the planet carrier and the planet shafts (single side).

2.2.8. System Elastic Behavior Analysis

The finite element models in the system are polycondensed to extract the corresponding mass and stiffness matrices, and the calculation process can be performed on general finite element software, such as ANSYS or ABAQUS (ABAQUS 2021, Dassault SIMULIA, Providence, RI, USA) or the RM software platform. At the same time, the deformation smoothness was used as an evaluation index, and the performance of each stiffness matrix was tested by means of the load transfer behavior among polycondensation nodes. Finally, the loading boundary conditions were applied to the system model, and the global calculation of quasi-static elastic behavior was performed to obtain the node displacement response of each elastic member under rated working conditions. In the system-level simulation analysis, the planet gear transferred the power received from the sun gear to the double-row tapered roller bearing (DRTRB) and the planet shaft, in turn, and finally, acted on the load on a pin hole of the planet carrier. Figure 7 shows the calculation results of the planet carrier elastic deformation and the load distribution on each DRTRB. Under the rated working conditions of the system, the maximum node resultant displacement (NRD) at the pin hole of the planet carrier exceeded 1000 μm, which will lead to significant movement position error for the planet shaft in working state, and seriously affect the tooth root stress response in the planetary gear train.

In the system-level simulation analysis, it is necessary to obtain the node displacement response of the planet gear rim and teeth, so as to provide necessary displacement boundary conditions for the secondary sub-model of tooth root stress calculation. Figure 5 shows the calculation results of the elastic deformation and stress distribution of a planet gear under rated working conditions. The structural element integrity of the system model ensures that the elastic deformation of the planet carrier and planet shaft, as well as the bearing clearance and other factors, can be taken into account in the results. Comparing the original contour of the planet gear with its loaded deformation contour (the deformation proportion has been enlarged), as shown in Figure 8a, it can be seen that the elastic deformation of the planet gear rim presents an elliptical shape. Figure 8b shows the stress distribution nephogram of the planet gear after being loaded, which identifies the two teeth meshing with sun gear and ring gear, respectively. The high stress zone was located at the tooth root of the loaded side of the two teeth, respectively, and the maximum von Mises stress (vMS) reached 298 MPa. It can also be found from the figure that the tooth root stress value of the teeth that were not engaged in meshing was not zero, which is mainly due to the significant elliptical deformation of the rim structure. The imprecise tooth root stress behavior was included in the system-level simulation analysis results, and more accurate tooth root stress calculation results need to be further obtained from the secondary sub-model.

In the system-level simulation analysis, the node displacement response of the sun gear rim and teeth are reflected in the displacement nephogram of the input shaft, and the calculation result is shown in Figure 9a. It can be seen from Figure 9 that the sun gear and the ring gear have the same deformation trend, and there are obvious elastic deformation zones near the teeth meshing with the planet gears. The flexible design of the sun gear rim structure makes each high deformation zone evenly cover multiple teeth in a circumferential direction, and the thin-walled rim and cantilever support of the ring gear make each high deformation zone offset by the same amount in axial direction. These elastic behavior characteristics will have a significant impact on the tooth root stress distribution and stress level for the corresponding gears.

2.2.9. Calculation of Tooth Surface Load Line

The tooth surface micro-geometric analysis was carried out at the system-level, which used the boundary conditions of the whole system model. Based on the calculation results of system elastic deformation and meshing misalignment, and considering tooth surface micro-geometric information, the dynamic load line on tooth surface was calculated, which provided the load boundary conditions for the secondary sub-model of tooth root stress calculation. The micro-geometric analysis modules in the RotationMaster and RomaxDesigner software (RomaxDesigner20, Romax Technologies Ltd., Nottingham, UK) can carry out very accurate mathematical definition for the tooth surface, including micro-geometric parameters, such as involute shape and tooth surface modification, and the nonlinear contact model of the tooth surface can be automatically created by accurately evaluating the geometry and meshing direction on the tooth surface. The tooth surface modification geometric information in the planetary mechanism is digitally expressed, and the visualized modeling results are shown in Figure 10, which will be used as one of the key input conditions for the system-level tooth surface micro-geometry analysis. The contour maps of the modification amount in the figure, respectively, describes the modification state on the working tooth surface for each gear, which, respectively, cover the micro-modification data, such as inclination, drum amount, and top trimming, in both directions of tooth direction and tooth profile.

The tooth surface dynamic load line results of each gear in the system can be obtained through tooth surface micro-geometric analysis. The results consider the influence of the elastic deformation and meshing misalignment of the system and record the information, such as the meshing line movement trajectory on tooth surface and the tooth surface load distribution under each rotation step. The calculation results of the tooth surface load line of the sun gear, planet gear, and ring gear are shown in Figure 11, which show the tooth surface load distribution state, when the load line moved to the center of tooth height. It can be seen from the figure that the load on each tooth surface was basically evenly distributed along the direction of tooth width, and a peak load was located in the center area of tooth width. The load lines (contact spots) in Figure 11c,d were wider, which was the result of internal meshing forming between the ring gear and planet gear.

2.3. Secondary Sub-Model Construction and Tooth Root Stress Analysis

2.3.1. Construction of Secondary Sub-Model

In the system-level model, the dynamic load lines on tooth surface were obtained through quasi-static static analysis and tooth surface micro-geometric analysis (see Figure 11), which could be loaded on the tooth surface of the secondary sub-model as load boundary conditions. At the same time, the system elastic deformation results, including the deformation behavior of the gear rim and teeth (see Figure 7, Figure 8 and Figure 9), could be extracted and loaded into the sub-model as displacement boundary conditions. In this way, the tooth root stress analysis results of the sub-model naturally included the elastic deformation of the entire system (including rim distortion deformation and tooth bending deformation, etc.), meshing misalignment, and tooth surface micro-geometry, as well as other factors. In addition, the dynamic behavior factors, such as the gear rotation and meshing, were also considered in the system model. When these important factors are fully considered in the system model, the overall structure of the sub-model can be significantly simplified. It could only contain a few teeth and the corresponding rim feature; the configuration is shown in Figure 12, and its modeling and computational costs focused more on the tooth root geometric details. In terms of model quality, the tooth root geometric elements of the sub-model were more comprehensive than those of the system model, and its quality indicators, such as accuracy level and mesh quality were also higher.

The secondary sub-model of the tooth root stress analysis with high-fidelity was established by a general finite element method. The parametrically defined detailed gear solid models in the system-level model were auxiliary modeling elements for the 3D finite element sub-models, which solved the problem of difficult modeling for the general finite element method, to a certain extent. In addition, in order to further alleviate the contradiction between calculation accuracy and operation speed, a first-order hexahedron element was adopted in the finite element sub-model. When the meshes are relatively coarse, the tooth top may appear “jagged” in appearance, due to shear self-locking or hourglass effect from the first-order element, but this does not affect the analysis results of tooth root stress. RM software platform can recommend economical and adaptive mesh density here, which can accurately capture the steep stress gradient at tooth root and is enough to make the calculation of the sub-model converge within predetermined number of iterations. In system-level analysis, the meshing force on teeth is applied through condensation nodes and finite element nodes, but in a secondary sub-model, the meshing force is applied by the tooth surface load line obtained from the system-level analysis.

2.3.2. Result Analysis of Tooth Root Stress

If the planetary system is considered an absolutely rigid system, and the misalignment of the gear mesh will be caused mainly by the clearance between the moving parts. This will generate phase-wise meshing interference between the planetary gears and, thus, will significantly increase the stress level in the tooth roots. Introducing the elastic mechanism into the system model can alleviate the meshing interference behavior between gears, to some extent (the more significant the flexible characteristics of the system within a certain range, the more pronounced this mitigation effect), reducing tooth root stress.

Firstly, the tooth root stress analysis program of the sun gear sub-model was run, and the maximum principal stress value at the tooth root was calculated according to the first strength theory. The calculation results of the system deformation and tooth surface load line are applied to the sub-model as input conditions, and then the model is solved according to the set rotation steps. According to the coincidence degree of gear pairs in the system, the meshing in and meshing out processes of a tooth are divided into 16 rotation steps, and the sub-model performed a finite element solution for each rotation step. Although this made the processes of stiffness decomposition and load vector inverse substitution relatively complex (involving multiple recursive traversals at the substructure layer), it was worthwhile to significantly reduce the calculation cost by appropriately increasing the program complexity.

Figure 13a shows the calculation results of the tooth root stress of the sun gear under 16 rotation steps. The stress results, without considering system elastic deformation, were 11.9%–17.3% higher than those considering this factor, which indicates that ignoring the flexible behavior characteristics of the large aviation planetary mechanism may directly lead to an over-conservative design scheme for corresponding structural strength. Figure 12 shows the evolution process of tooth root stress distribution of sun gear teeth during meshing (corresponding to steps 8~16 in Figure 13a). Load lines successively passed through three teeth, A, B, and C, in the sub-model, and the figure indicates the specific value and position of the maximum tooth root stress (the maximum mentioned here refers to the maximum in position) under each rotation step. As can be seen from the action of load lines in the first three steps that tooth A gradually disengaged meshing, and its root stress kept decreasing, while the adjacent tooth B gradually entered meshing, and its root stress kept increasing. When step 11 was reached, the maximum tooth root stress position jumped from tooth A to tooth B. When the 12th step was reached, the sun gear changed from double-tooth meshing to single-tooth meshing, and the root stress at tooth B increased significantly, then reached a stress peak state at the 14th rotation step. It was the maximum value state that tooth root stress could reach during the meshing process, and it was also one of the most dangerous load states that caused the sun gear to suffer from bending fatigue tooth breakage; then, the gear entered double-tooth meshing again and continuously circulated the above process. In addition, the stress peak of tooth root bending fatigue obtained by the calculation model in ISO6336 is also shown in Figure 13 a, and the calculation model is as follows:

(7)${\sigma }_F=\frac{F_{\mathrm{t}}}{b\cdot m}\cdot K_A\cdot K_V\cdot K_{F\beta }\cdot K_{F\alpha }\cdot Y_F\cdot Y_S\cdot Y_{\beta }\cdot Y_B\cdot Y_{DT}$

The simulation result, without considering system elastic deformation, is in good agreement with this calculation result, which proves the validity of the load boundary conditions and other model parameters in the simulation analysis process.

Figure 13b shows the calculation results of the tooth root stress of the planet gear meshing with the sun gear and the ring gear, respectively, in 16 rotation steps, and it also marks the demarcation line of the tooth root stress grade under single-tooth meshing and double-tooth meshing. One side tooth surface of the planet gear tooth was meshed with the sun gear, and the stress peak at the tooth root reached 272 MPa; the other side’s tooth surface meshed with the ring gear, and the stress peak at the tooth root was 276 MPa. In addition, the calculation result of the tooth root stress peak of the ring gear was 280 MPa. These stress results were obtained under system-rated working conditions, and the effects of the elastic deformation, meshing misalignment, and tooth surface micro-geometry of the whole system were fully considered.

2.3.3. Load Input Variable Conversion

Based on the hierarchical finite element method, the system power peak distribution was transformed into the gear stress peak distribution, which provided direct load input variables for the system fatigue reliability evaluation model. In a typical service mission, with a long flight time, the input power history of this kind of planetary system was collected, as shown in Figure 14a. The power peaks in the load history were extracted as load statistical characteristic and fitted with normal distribution, and the power peak distribution result is shown in Figure 14b. Based on the hierarchical finite element method, a mapping relationship of the same probability quantile from power distribution to stress distribution was established, and the power peak distribution was transformed into the tooth root stress peak distribution of each gear in the system. Final calculation results are shown in Figure 14d, where the PS represents the tooth root stress on the planet gear tooth (the tooth side meshed with the sun gear), and the other codes in the figure have similar meanings.

3. Tooth Probability Strength Fitting Based on Gear Fatigue Test

3.1. Gear Bending Fatigue Test

A gear bending fatigue accelerated life test was carried out by a power flow closed gear rotation testing machine, which provided strength input variables for the fatigue reliability prediction model of the large aviation planetary system. The overall layout of the test platform is shown in Figure 15a, and its working principle and verification standard are shown in the previous relevant research work [2]. For power flow input inside the testing machine, a conical friction surface-type mechanical loading device was adopted, as shown in Figure 15d, which ensured a reliable seal of the power flow under an ultra-high cycle. A vibration monitoring system (see Figure 15b,e) enabled the testing machine to have the function of real-time automatic shutdown, in case of sudden fatigue tooth breakage, which better ensured the consistency of the failure states for all the gear samples. The final shapes and sizes of each tooth root crack were basically the same, as shown in Figure 15h. At the same time, tooth root fatigue fracture surface morphology is shown in Figure 15i, and the macro morphological features, such as fatigue expansion zones and fatigue instantaneous fracture zones, can be clearly seen. A top-view of test gearbox inside is shown in Figure 15g. The gear pair adopted an assembly form of full tooth width contact, and the lubrication and cooling for the gear pair were realized by oil injection during the test. The structural details and specific parameters of the gear samples are shown in Figure 15f and Table 5, respectively, and the preparation process strictly followed the principle of “teeth fatigue strength equivalent”, that is, to make the sample parameters as equal, or similar, to the gear parameters in the planetary mechanism as possible, in terms of material properties, tooth geometry, machining, and heat treatment, etc. A reliable guarantee for this similarity level will help to improve fatigue reliability prediction accuracy for the planetary system. In addition, all the gear samples came from the same batch production process to minimize the dispersion of test data.

3.2. Tooth Probability Strength Fitting

The tooth root stress peak was taken as an evaluation index for the stress grade, and the bending fatigue performance of the tooth root was tested by the group method under four stress grades. The selected stress levels and the number of test points under various stress levels were 649 MPa (17 points), 618 MPa (22 points), 586 MPa (29 points), and 555 MPa (38 points), respectively. During the test, if any tooth on the gear sample failed first, the testing machine would automatically stop, and the direct data obtained from this was the gear life, rather than the tooth life. It represents the ability of the individual gear to maintain excellent transmission function under current stress level, so the gear life was also the “first broken tooth” life. From the perspective of probability, the more teeth on a gear, the more potential failure links. Therefore, under the same stress and revolution conditions, the failure risk of the gear will increase with the increase of the number of teeth. In the process of this test, a complete rotation of the sample was recorded as one gear life. In this conventional counting mode, the statistical characteristics of the number of teeth led to the difference in the probability life between the gear and tooth. In order to obtain direct strength input variables for the reliability prediction model, a probabilistic statistical transformation method was proposed to fit the tooth P-S-N curves, based on the gear life data.

The probability life relationship between the gear and tooth was established based on the concept of minimum order statistics. The fracture of any tooth on a gear will cause the gear to lose excellent transmission capacity. For this reason, it can be considered that the life of a gear depends on the minimum life of its teeth. Suppose ${\mathsf{X}}_1, {\mathsf{X}}_2, \dots , {\mathsf{X}}_n$ is a set of samples from a parent X, then ${\mathsf{X}}_{\min }=\min \left({\mathsf{X}}_1, {\mathsf{X}}_2, \dots , {\mathsf{X}}_n\right)$ is the minimum order statistics of the parent. This probability model will be applied to the life transformation calculation under the failure mode of tooth root bending fatigue, and then the gear probability life is equal to the minimum order statistics of tooth probability life.

Assuming that the cumulative distribution function of random variable X is $F\left(x\right)$ and its probability density function is $f\left(x\right)$, then the probability density function of the minimum order statistics of X can be expressed as

(8)$g_{\min }\left(x\right)=z\cdot {\left[1-F\left(x\right)\right]}^{z-1}\cdot f\left(x\right)$

If the two-parameter Weibull distribution is adopted to express tooth probability life, then the cumulative distribution function can be expressed as

(9)$F\left(x\right)=1-\exp \left[-{\left(x/\theta \right)}^{\beta }\right]$

and the probability density function is

(10)$f\left(x\right)=\left(\beta \cdot x^{\beta -1}/{\theta }^{\beta }\right)\exp \left[-{\left(x/\theta \right)}^{\beta }\right]$

where β and θ are, respectively, the shape parameter and scale parameter of the tooth life distribution.

Equations (9) and (10) are brought into Equation (8) to obtain the following equation

(11)$g_{\min }\left(x\right)=\left[\beta \cdot x^{\beta -1}/{\left(\theta /z^{1/\beta }\right)}^{\beta }\right]\exp \left[-{\left(x\cdot z^{1/\beta }/\theta \right)}^{\beta }\right]$

If the number of teeth on a gear is z, then $g_{\min }\left(x\right)$ directly represents the probability density function of the gear life distribution. From the expression of this function, it can be seen that the gear probability life also follows two-parameter Weibull distribution, and the shape parameter and scale parameter are as follows

(12)$\left\{\begin{array}{l}{\beta }_{\mathrm{Gear}}=\beta \\ {\theta }_{\mathrm{Gear}}=\theta /z^{1/\beta }\end{array}\right.$

In the statistical processing for the test data, the two-parameter Weibull distribution function was adopted to fit the probability distribution of gear life points under each stress level. The probability life transformation between the gear and tooth was then performed by model (12). Finally, a least square method was used to linearly fit the same probability quantiles of the tooth life distribution under each stress level in a single logarithmic coordinate system, and the results of tooth bending fatigue P-S-N curves obtained are shown in Figure 16. Under deterministic loading, the dispersion of fatigue life generally increases as stress level decreases; therefore, in a linear coordinate system, the P-S-N curve family will appear as an “umbrella” shape with a small upper opening and a large lower opening. However, due to the single logarithmic coordinate system, the curve family presents a corresponding inverted shape.

4. System Reliability Modeling Considering Planetary Transmission Characteristic

4.1. Conditional Probability Expectation Algorithm for Part Fatigue Reliability Calculation

The traditional “load and strength interference” analysis method was extended to establish a conditional probability expectation algorithm for calculating the fatigue reliability of the parts, based on the probability distribution of the stress level and the life distribution under the specified stress level. For the static strength failure of parts under one load, the reliability can be regarded as a function of stress, and the conditional reliability model under specified stress is established. That is, under the condition of stress $\sigma$, the probability calculation formula of static strength S greater than the stress is

(13)$\zeta \left(\sigma \right)={\displaystyle {\int }_{\sigma }^{\infty }f\left(S\right)\mathrm{d}S}$

where $f\left(S\right)$ is the probability density function of static strength S.

Furthermore, the static strength reliability model of parts that can reflect the effect of load uncertainty can be expressed as

(14)$R={\displaystyle {\int }_0^{\infty }h\left(\sigma \right)\cdot \zeta \left(\sigma \right)\mathrm{d}\sigma }$

where $h\left(\sigma \right)$ is the probability density function of stress $\sigma$.

If Equation (14) is the traditional “load and strength interference model”, then the basic meaning is extended from probability perspective. Based on the total probability calculation principle of continuous random variable, Equation (14) can be interpreted as the mathematical expectation of random function $\zeta \left(\sigma \right)$ in the definition domain of random variable $\sigma$.

For fatigue reliability, because it is difficult to obtain the fatigue strength distribution under specified life directly, a mathematical expression for calculating the fatigue reliability of parts directly, based on life distribution, can be constructed, according to above extended thinking. Assuming that $\zeta \left(n,\sigma \right)$ is the conditional fatigue reliability function under the specified stress level, that is, under the condition of stress level $\sigma$, the probability calculation formula of the life N greater than the number of load cycles n is

(15)$\zeta \left(n,\sigma \right)={\displaystyle {\int }_n^{\infty }f\left(N\left|\sigma \right.\right)\mathrm{d}N}$

where $f\left(N\left|\sigma \right.\right)$ is the life probability density function at stress level $\sigma$.

Correspondingly, under the action of random stress level $\sigma$, the fatigue reliability model of parts can be expressed as

(16)$R\left(n\right)={\displaystyle {\int }_0^{\infty }h\left(\sigma \right)\cdot \zeta \left(n,\sigma \right)\mathrm{d}\sigma }$

The fatigue reliability index of parts under random constant amplitude cyclic load can be directly calculated by Equation (16).

4.2. Fatigue Reliability Evaluation Model of Series System Considering Failure Dependence

The fracture of any tooth in a planetary system will affect the transmission capacity of the whole system. Therefore, if each tooth in the system is regarded as a potential failure unit, then the system is a typical series system. In the service process, the load borne by each tooth has significant dependence on the system input power, and the load dependence and general load randomness make the failure of each element not independent of each other, so the reliability of the series system cannot be simply considered as the reliability product of each unit.

The fatigue reliability evaluation model of a series system is established based on conditional probability expectation algorithm and considering the failure dependence among unit parts. First of all, only under the action of deterministic load, the failure of each part in the system is independent of each other, so the conditional reliability of the series system is equal to the conditional reliability product of each part. Then, considering the load uncertainty effect at the system-level, the probability that r parts in the system will not fail, that is, the fatigue reliability evaluation model of the series system is

(17)$R_{\mathrm{SYS}}\left(n\right)={\displaystyle {\int }_0^{\infty }{h}_0\left(\sigma \right){\displaystyle \prod _{i=1}^r\left[{\displaystyle {\int }_n^{\infty }{f}_i\left(N\left|\left(s_i\sigma +{\mu }_i\right)\right.\right)\mathrm{d}N}\right]}\mathrm{d}\sigma }$

where $h_0\left(\sigma \right)$ is the probability density function of standard normal distribution, and $f_i\left(N\left|\left(s_i\sigma +{\mu }_i\right)\right.\right)$ is the life probability density function of the ith part under stress level $s_i\sigma +{\mu }_i$.

The model assumes that the load follows normal distribution and realizes load normalization for each part through the mathematical relationship transformation between normal distribution and standard normal distribution, so as to consider the load uncertainty effect at the system-level.

4.3. Structural Optimization of Reliability Model Considering Sequence Characteristic

The kinematic equation of planetary transmission was deduced according to a periodic operation law in the planetary gear train; at the same time, the single tooth meshing times of various gears in the system within the same time interval were obtained, which matched the sequence characteristic attributes for the system fatigue reliability evaluation model. Assuming that the absolute angular velocities of the sun gear, planet gear, ring gear, and planet carrier are ${\omega }_{\mathrm{S}}$, ${\omega }_{\mathrm{P}}$, ${\omega }_{\mathrm{R}}$, and ${\omega }_{\mathrm{C}}$, respectively, then the mathematical relationship among them can be expressed as:

(18)$\left\{\begin{array}{l}{\omega }_{\mathrm{S}}=i_{\mathrm{SR}}^{\mathrm{C}}\cdot {\omega }_{\mathrm{R}}+i_{\mathrm{SC}}^{\mathrm{R}}\cdot {\omega }_{\mathrm{C}}\\ {\omega }_{\mathrm{R}}=i_{\mathrm{RS}}^{\mathrm{C}}\cdot {\omega }_{\mathrm{S}}+i_{\mathrm{RC}}^{\mathrm{S}}\cdot {\omega }_{\mathrm{C}}\\ {\omega }_{\mathrm{C}}=i_{\mathrm{CS}}^{\mathrm{R}}\cdot {\omega }_{\mathrm{S}}+i_{\mathrm{CR}}^{\mathrm{S}}\cdot {\omega }_{\mathrm{R}}\\ {\omega }_{\mathrm{P}}=i_{\mathrm{PC}}^{\mathrm{R}}\cdot {\omega }_{\mathrm{C}}+i_{\mathrm{PR}}^{\mathrm{C}}\cdot {\omega }_{\mathrm{R}}\end{array}\right.$

where $i_{\mathrm{ab}}^{\mathrm{c}}$ represents the ratio of the relative rotational speed of member a and member b, respectively, relative to member c, i.e., $i_{\mathrm{ab}}^{\mathrm{c}}=\left({\omega }_{\mathrm{a}}-{\omega }_{\mathrm{c}}\right)/\left({\omega }_{\mathrm{b}}-{\omega }_{\mathrm{c}}\right)$.

The kinematic equation of the planetary gear train can be derived from Equation (18)

(19)$\left\{\begin{array}{l}{\omega }_{\mathrm{S}}+p\cdot {\omega }_{\mathrm{R}}-\left(1+p\right)\cdot {\omega }_{\mathrm{C}}=0\\ {\omega }_{\mathrm{R}}+{\omega }_{\mathrm{S}}/p-\left(1+p\right)\cdot {\omega }_{\mathrm{C}}/p=0\\ {\omega }_{\mathrm{C}}-{\omega }_{\mathrm{S}}/\left(1+p\right)-p\cdot {\omega }_{\mathrm{R}}/\left(1+p\right)=0\\ {\omega }_{\mathrm{P}}-\left(1+p\right)\cdot {\omega }_{\mathrm{C}}/\left(1-p\right)+2p\cdot {\omega }_{\mathrm{R}}/\left(1-p\right)=0\end{array}\right.$

where p is the kinematic characteristic parameter of the planetary gear train, which is the ratio of the number of teeth between ring gear and sun gear, i.e., $p=z_{\mathrm{R}}/z_{\mathrm{S}}$.

Through Equation (19), the relative angular velocities of the sun gear, planet gear, and ring gear, relative to the planet carrier and their single tooth meshing times in the same time interval, can be obtained, and the kinematic parameters are shown in Table 6, where the positive and negative signs indicate that the rotation directions are opposite, and $k_{\mathrm{P}}$ is the number of planet gears in the system. System input speed ${\omega }_{\mathrm{S}}$ is a known condition, and $n_{\mathrm{PS}}\left(t\right)$ is the meshing times between the target single tooth of a planet gear and the sun gear within a time interval t, and the interpretation for other relevant parameters is similar.

The parameters of the single tooth meshing times in Table 6 were brought into the model (17), the calculation factors ${\zeta }_i$ of the tooth element conditional fatigue reliability of various gears were obtained, and the fatigue reliability evaluation model $R_{\mathrm{SYS}}\left(t\right)$ for the planetary system, considering failure dependence and meshing sequence, was as follows

(20)$\left\{\begin{array}{l}{\zeta }_{\mathrm{S}}={\displaystyle {\int }_{n_{\mathrm{S}}\left(t\right)}^{\infty }{f}_{\mathrm{S}}\left(N\left|\left(s_{\mathrm{S}}\sigma +{\mu }_{\mathrm{S}}\right)\right.\right)\mathrm{d}N}\\ {\zeta }_{\mathrm{PS}}={\displaystyle {\int }_{n_{\mathrm{PS}}\left(t\right)}^{\infty }{f}_{\mathrm{PS}}\left(N\left|\left(s_{\mathrm{PS}}\sigma +{\mu }_{\mathrm{PS}}\right)\right.\right)\mathrm{d}N}\\ {\zeta }_{\mathrm{PR}}={\displaystyle {\int }_{n_{\mathrm{PR}}\left(t\right)}^{\infty }{f}_{\mathrm{PR}}\left(N\left|\left(s_{\mathrm{PR}}\sigma +{\mu }_{\mathrm{PR}}\right)\right.\right)\mathrm{d}N}\\ {\zeta }_{\mathrm{R}}={\displaystyle {\int }_{n_{\mathrm{R}}\left(t\right)}^{\infty }{f}_{\mathrm{R}}\left(N\left|\left(s_{\mathrm{R}}\sigma +{\mu }_{\mathrm{R}}\right)\right.\right)\mathrm{d}N}\end{array}\right.$

(21)$R_{\mathrm{SYS}}\left(t\right)={\displaystyle {\int }_0^{\infty }{h}_0\left(\sigma \right)\cdot {\zeta }_{\mathrm{S}}^{z_{\mathrm{S}}}\cdot {\left({\zeta }_{\mathrm{PS}}\cdot {\zeta }_{\mathrm{PR}}\right)}^{k_{\mathrm{P}}\cdot z_{\mathrm{P}}}\cdot {\zeta }_{\mathrm{R}}^{z_{\mathrm{R}}}\mathrm{d}\sigma }$

It is assumed that the manufacturing and assembly process of such large aviation planetary mechanism products has well-met the accuracy and quality requirements in the product design. Based on the above model, the fatigue reliability evolution law of the planetary systems before the first overhaul was predicted, and the result is shown in Figure 17. The prediction result of the traditional model shows an excessively rapid downward trend, and the downward speed will raise with the increase of the number of units in the system, which is difficult to accept. The curve result, considering failure dependence (FD), was relatively accurate in overall shape, which could distinguish the evolution characteristics of product fatigue reliability during an accidental failure period and a depletion failure period. However, only the prediction result that considers both failure dependence and meshing sequence (MS) could achieve excellent agreement with the relevant historical information (based on the statistical processing results from the historical censored data [22]).

5. Reliability-Driven Optimization Design for Key Structural Elements

Based on the above reliability calculation method to redesign the key geometric features of the large aeronautical planetary mechanism, the reliability-sensitive structural elements of large thin-walled parts in the system were screened, and multi-objective dimensional optimization analysis was executed. Taking the rated working conditions of the system as the load boundary conditions of the simulation model, the uncertainty effect of the load was ignored in the reliability analysis processes, so the load distribution function in Equation (15) will be provided in the form of the determined tooth root stress peak. The simulation analysis shows that the ring gear rim and the planet carrier played a main supporting role for the planetary gear system, and their core structural parameters largely determined the meshing quality of the whole planetary gear system. Among them, insufficient rigidity of the ring gear rim led to excessive bending deformation for the gear teeth. Insufficient rigidity of the planet carrier baseplate will cause severe deflection for the planet shaft, and these design defects will increase the risk of fatigue tooth breakage in the planetary gear system, and the corresponding dimensional schematic is shown in Figure 18. In turn, both of these dimensional characteristics are core factors in the lightweight design of large aerospace planetary equipment, and their dimensional growth contributes significantly to the weight growth for the planetary mechanism. Therefore, taking the thickness of the rim and base plate as the design variables, based on the reliability analysis and calculation method proposed in this paper, under complex system elastic behavior coupling mechanism, the best stiffness matching result of the thickness of the rim and base plate that jointly meets the requirements of reliability and lightweight indicators is sought.

The elastic deformation results of the ring gear and the planet carrier, as well as the reliability results of the planetary gear system, are shown in Figure 19. A total of 30 data points are calculated in Figure—six points were taken in the range of 10~35 mm for the ring gear rim thickness, with 5 mm was taken for the point interval; 5 points were taken in the range of 12~36 mm for the planet carrier base plate thickness, with 6 mm for the point interval, and the spline curves were used to fit the data points. It can be found from Figure 19a that the maximum node resultant displacement almost stopped decreasing when the rim thickness reached 25.3 mm, indicating that the rigidity reserve of the ring gear raised by increasing the rim size can no longer be effectively utilized at this limit value. At the same time, the elastic deformation response of the ring gear is also influenced, to some extent, by the thickness dimension of the planet carrier base plate. Under the same rim thickness, the deformation of the ring gear slightly decreases with the increase of the base plate thickness, which indicates that the improvement of the planetary system stiffness conditions caused by the thickening of the planet carrier base plate also has a benign effect on the mechanical environment of the ring gear.

Figure 19b shows the variation pattern of the fatigue reliability index of the planetary gear system under the joint influence of rim and base plate sizes. As the thickness of the ring gear rim increases, the reliability of the system will continue to grow, and when the value of the rim thickness reaches 22.5 mm, the reliability growth is weak. At the same time, the reliability of the system also increases with the increase of the thickness of the planet carrier base plate, and the reliability growth is weak when the thickness reaches 30.5 mm. The results of the reliability curves in the figure reflect the mechanisms of the coupling effect of the rim and base plate dimensions on the fatigue reliability of the system. In the variation range of the rim thickness from 10 to 22.5 mm, the growth rate of each reliability curve increases with the increase of planet carrier base plate thickness, indicating that their scale growth in a certain range has a mutual promotion effect in optimizing the system reliability index. By setting the coordinate point with a reliability growth rate of less than 0.03% as the optimal design point, and considering the static strength requirements of the ring gear and the planet carrier, the best matching values of the two key structural dimensions of the large aero-planetary mechanism can be determined: 22.5 mm thickness for the ring gear rim and 30.5 mm thickness for the planet carrier base plate.

6. Concluding Remarks

This paper takes the large aviation planetary system as a research object, aims to accurately evaluate the reliability quality attributes formed in the system design process, and deeply excavates the inherent behavior characteristics for the planetary system in function realization. Based on this, a system fatigue reliability evaluation model was constructed by considering the failure dependence and meshing sequence among the functional logic units. At the same time, an advanced hierarchical finite element method was adopted to calculate the tooth dangerous load history under the influence of system global elastic behavior, and the tooth probability fatigue strength is statistically analyzed by the gear low-cycle fatigue test and life distribution transformation method. These will provide economic effective load and strength input variables for the reliability evaluation model, respectively, and the specific conclusions are as follows.

• (1). In the face of advanced simulation and analysis tasks for a large aviation planetary system, only considering the convenience of modeling and boundary condition setting, the computational efficiency of the hierarchical finite element method will be much higher than that of the general finite element method. Compared with the results of the tooth root stress obtained by the international standard method, the hierarchical finite element method had high calculation accuracy. Moreover, the stress results, without considering system elastic deformation, were 11.9%~17.3% higher than those considering this factor, which indicates that ignoring the flexible behavior characteristics of the large aviation planetary mechanism may directly lead to an over-conservative design scheme for the corresponding structural strength.

• (2). The gear bending fatigue accelerated life test was carried out by a power flow closed gear rotation testing machine; a total of 106 life data points were tested under four stress levels by a group method, and the probability life relationship between the gear and tooth was established based on the concept of minimum order statistics. The linear correlation of the P-S-N curves obtained by the statistical method in this paper was more than 96%, which ensures the effectiveness of strength input variables for the reliability model.

• (3). With the increase of the ring gear rim thickness of the large, thin-walled internal gear ring, its elastic deformation under the rated working condition gradually decreased, and the maximum node resultant displacement almost stopped decreasing when the rim thickness reached 25.3 mm, indicating that the rigidity reserve of the ring gear raised by increasing the rim size can no longer be effectively utilized at this limit value. At the same time, the elastic deformation response of the ring gear was also influenced, to some extent, by the thickness dimension of the planet carrier base plate. Under the same rim thickness, the deformation of the ring gear slightly decreased with the increase of the base plate thickness, which indicates that the improvement of the planetary system stiffness conditions caused by the thickening of the planet carrier base plate also has a benign effect on the mechanical environment of the ring gear.

• (4). Within a certain size range, the increase of the ring gear rim thickness and the planet carrier base plate thickness will improve the stiffness conditions of the planetary gear system and optimize the gear meshing performance, thus improving the fatigue reliability level of the planetary gear system. When the ring gear rim is thin, the increase in the thickness of the planet carrier base plate will accelerate the change in the reliability of the system caused by the thickening of the rim, while, when the ring gear rim is thick, the thickness of the planet carrier base plate will have the opposite effect. This “thin vs. thick” dimensional range will be determined by the other structural elements in the planetary mechanism, as well as the level of external loads. Based on the base plate thickness–rim thickness–reliability curve cluster, the best matching results for two key structural dimensions in a specified type of large aviation planetary system to meet the reliability and lightweight requirements were determined. The results of this study will provide important reference data for the structural optimization design of a large aviation planetary system.

Author Contributions

Conceptualization, M.L.; methodology, L.X.; software, M.L. and Y.L.; validation, M.L.; formal analysis, M.L.; investigation, M.L. and L.X.; resources, M.L.; data curation, M.L. and Y.L.; writing—original draft preparation, M.L.; writing—review and editing, M.L., Y.L. and L.X.; visualization, L.X.; supervision, M.L.; project administration, M.L.; funding acquisition, M.L. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank Liyang Xie for the reliability calculation model.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

$${\rho }_{fp}$$	Tooth root radius
CPU	Central processing unit
FE	Finite element
TRB	Tapered roller bearings
DRTRB	Double-row tapered roller bearing
RBB	Radial ball bearing
d1	Pitch circle diameter
PCD	Pitch circle diameter
TF	Tolerance factor
NRD	Node resultant displacement
PS	The tooth root stress on the planet gear tooth (the tooth side meshed with the sun gear)
P (kw)	The input power history of this kind of planetary system has been collected
Xi	A set of samples from a parent X
$F\left(\mathrm{x}\right)$	The cumulative distribution function of random variable X
$$f\left(x\right)$$	Probability density function
$$g_{\min }\left(x\right)$$	The probability density function of the minimum order statistics of X
β	The shape parameter of the tooth life distribution
θ	The Scale parameter of the tooth life distribution
z	Number of teeth on a gear
$$\sigma$$	Condition of stress
$$\zeta \left(\sigma \right)$$	Under the condition of stress $$\sigma$$, the probability calculation formula of static strength S greater than the stress
$$f\left(S\right)$$	The probability density function of static strength S
$$h\left(\sigma \right)$$	The probability density function of stress $$\sigma$$
$$R\left(n\right)$$	Under the action of random stress level $$\sigma$$, the fatigue reliability model of parts
$$h_0\left(\sigma \right)$$	The probability density function of standard normal distribution
$$\omega$$	Represents the angular velocity of the member
$$i_{\mathrm{ab}}^{\mathrm{c}}$$	Represents the ratio of the relative rotational speed of member a and member b respectively relative to member c
$$p$$	The kinematic characteristic parameter of the planetary gear train
$$k_{\mathrm{P}}$$	The number of planet gears in the system
$$n_{\mathrm{PS}}\left(t\right)$$	The meshing times between the target single tooth of a planet gear and the sun gear within time interval t
$${\zeta }_i$$	The calculation factors of tooth element conditional fatigue reliability of various gears
FD	Failure dependence
MS	Meshing sequence

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. System model of planetary mechanism.

Figure 2. Example of FE model and geometry of planet carrier.

Figure 3. Planet bearing element.

Figure 4. Geometrical parameter for internal gear modeling.

Figure 5. Bearing structure in planetary system.

Figure 6. Node search criteria: (a) shell search; (b) solid search. Where d1 is pitch circle diameter (PCD), and D and d2 indicate the search range. D = d1 × (1 + 2 × TF), d2 = d1 × (1 – 2 × TF).

Figure 7. Mechanical analysis result: (a) displacement nephogram of planet carrier; (b) roller force diagram of DRTRB.

Figure 8. Mechanical analysis result of planet gear: (a) displacement nephogram; (b) stress nephogram.

Figure 9. Elastic deformation result: (a) displacement nephogram of input shaft; (b) displacement nephogram of ring carrier.

Figure 10. Contour map of tooth surface modification: (a) planet gear tooth surface meshed with sun gear; (b) working tooth surface of sun gear; (c) planet gear tooth surface meshing with ring gear; (d) working tooth surface of ring gear.

Figure 11. Load line on tooth surface: (a) planet gear tooth surface meshed with sun gear; (b) working tooth surface of sun gear; (c) planet gear tooth surface meshing with ring gear; (d) working tooth surface of ring gear.

Figure 12. Evolution process of tooth root stress distribution on sun gear.

Figure 13. Maximum tooth root stress history: (a) sun gear; (b) planet gear.

Figure 14. Load input variable conversion: (a) input power history; (b) power peak distribution; (c) hierarchical finite element simulation; (d) stress peak distribution.

Figure 15. Gear bending fatigue test: (a) test platform; (b) vibration monitoring; (c) revolution counter; (d) loading device; (e) acceleration sensor; (f) gear sample; (g) gear installation; (h) root crack; (i) fracture morphology.

Figure 16. Statistical transformation of test data: (a) probability life transformation for tooth; (b) P-S-N curves fitting for tooth.

Figure 17. Prediction result of reliability model.

Figure 18. Sensitivity dimensions of system performance indexes: (a) rim thickness of the ring gear (b) base plate thickness of the planet carrier.

Figure 19. Elastic deformation and reliability results: (a) elastic deformation curves of the ring gear. (b) System reliability growth curves. Notes: BPT is the base plate thickness of the planet carrier.

Table 1

Geometric parameter of planetary gear train.

Parameters	Sun Gear	Planet Gear	Ring Gear
Module (mm)	5.012	5.012	5.012
Number of teeth	84	56	196
Number of gears	1	7	1
Pressure angle (°)	20	20	20
Helix angle (°)	0	0	0
Effective face width (mm)	120	120	120
Base circle diameter (mm)	394.671	263.114	920.899
Base circle pitch (mm)	14.761	14.761	14.761
Root fillet radius (mm)	2.757	2.870	2.657
Tooth surface hardness	60 HRC	60 HRC	60 HRC
Tooth core hardness	35 HRC	35 HRC	35 HRC
Elastic modulus (MPa)	2.07 × 105	2.07 × 105	2.07 × 105
Poisson ratio	0.3	0.3	0.3

Table 2

Bearing parameters in planetary system.

Parameters	TRB1	TRB2	DRTRB	RBB
External diameter (mm)	310	420	240	480
Internal diameter (mm)	200	300	160	360
Width (mm)	70	76	102	56
Number of rollers	31	40	72	24
Roller diameter (mm)	23	26	16	40
Roller length (mm)	50	55	38	-
Contact angle (°)	15.945	14.931	17.049	0

Table 3

Finite element model and mesh quality parameter.

Finite Element Models	Effective Size	Minimum Angle
[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]
Number of surface mesh node: 31,248	Number of surface mesh elements: 15,589
[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]
Number of surface mesh node: 9273	Number of surface mesh elements: 4606
[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]
Number of surface mesh node: 16,111	Number of surface mesh elements: 8021
[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]
Number of surface mesh node: 23,956	Number of surface mesh elements: 11,942

Table 4

Node connection between finite element model and solid model.

Connection Nodes	Connection Objects	Node Parameters
Input shaft	Meshing lines	Number of connection nodes: 1830Number of node groups: 7Node search criterion: solid searchNode selection method: spaceElement type: RBE2PCD = 10 mm, TF = 0.01
[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]
Planet gear	DRTRB	Number of connection nodes: 520Number of node groups: 2Node search criterion: solid searchNode selection method: surfaceElement type: RBE3PCD = 240 mm, TF = 0.01
[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]
Ring gear	Main reducer box	Number of connection nodes: 930Number of node groups: 34Node search criterion: shell searchNode selection method: surfaceElement type: RBE2PCD = 20 mm, TF = 0.01
[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]
Planet carrier	Planet shafts	Number of connection nodes: 1410Number of node groups: 7Node search criterion: solid searchNode selection method: surfaceElement type: RBE3PCD = 140 mm, TF = 0.01
[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]

Table 5

Parameter list of test gear.

Items	Parameters	Items	Parameters
Module (mm)	5	ISO quality grade	5
Number of teeth	25	Material brand	1Cr18Ni9Ti
Pressure angle (°)	20	Carburized depth (mm)	0.8 ± 0.13
Helix angle (°)	0	Tooth surface hardness	59–63 HRC
Face width (mm)	32	Tooth core hardness	35–48 HRC
Root fillet radius (mm)	2.7	Precision machining	Grinding

Table 6

Kinematic parameter of planetary gear train.

Components	Angular Velocity	Relative Angular Velocity	Single Tooth Meshing Times
Sun gear	$${\omega }_{\mathrm{S}}$$	$$p\cdot {\omega }_{\mathrm{S}}/\left(1+p\right)$$	$$n_{\mathrm{S}}\left(t\right)=p\cdot {\omega }_{\mathrm{S}}\cdot t\cdot k_{\mathrm{P}}/\left(1+p\right)$$
Planet gear	$${\omega }_{\mathrm{S}}/\left(1-p\right)$$	$$2p\cdot {\omega }_{\mathrm{S}}/\left(1-p^2\right)$$	$$n_{\mathrm{PS}}\left(t\right)=n_{\mathrm{PR}}\left(t\right)=2p\cdot {\omega }_{\mathrm{S}}\cdot t/\left(p^2-1\right)$$
Ring gear	0	$$-{\omega }_{\mathrm{S}}/\left(1+p\right)$$	$$n_{\mathrm{R}}\left(t\right)={\omega }_{\mathrm{S}}\cdot t\cdot k_{\mathrm{P}}/\left(1+p\right)$$
Planet carrier	$${\omega }_{\mathrm{S}}/\left(1+p\right)$$	0	-

References

1. Yin, Z.Y.; Fu, B.B.; Xue, T.B.; Wang, Y.H. Development of helicopter power transmission system technology. Appl. Mech. Mater.; 2011; 86, pp. 1-17. [DOI: https://dx.doi.org/10.4028/www.scientific.net/AMM.86.1]

2. Li, M.; Xie, L.; Ding, L. Load sharing analysis and reliability prediction for planetary gear train of helicopter. Mech. Mach. Theory; 2017; 115, pp. 97-113. [DOI: https://dx.doi.org/10.1016/j.mechmachtheory.2017.05.001]

3. Xue, H.L.; Liu, G.; Yang, X.H.; Han, B. Key technologies research of helicopter transmissions. Appl. Mech. Mater.; 2015; 743, pp. 55-60. [DOI: https://dx.doi.org/10.4028/www.scientific.net/AMM.743.55]

4. McFarland, J.M.; Riha, D.S. Uncertainty quantification methods for helicopter fatigue reliability analysis. Proceedings of the American Helicopter Society 65th Annual Forum; Grapevine, TX, USA, 27–29 May 2009; American Helicopter Society International: Grapevine, TX, USA, 2009.

5. Li, T.; Zhao, Z.; Sun, C.; Yan, R.; Chen, X. Adaptive channel weighted CNN with multisensor fusion for condition monitoring of helicopter transmission system. IEEE Sens. J.; 2020; 20, pp. 8364-8373. [DOI: https://dx.doi.org/10.1109/JSEN.2020.2980596]

6. Brown, M.A.; Chang, J.H. Analytical techniques for helicopter component reliability. Proceedings of the American Helicopter Society 64th Annual Forum; Montreal, QC, Canada, 29 April–1 May 2008; American Helicopter Society International: 2008.

7. Li, M.; Xie, L.Y.; Li, H.Y.; Ren, J.G. Life distribution transformation model of planetary gear system. Chin. J. Mech. Eng.; 2018; 31, pp. 1-8. [DOI: https://dx.doi.org/10.1186/s10033-018-0221-x]

8. Chen, J.; Li, W.; Sheng, L.; Jiang, S.; Li, M. Study on reliability of shearer permanent magnet semi-direct drive gear transmission system. Int. J. Fatigue; 2020; 132, 105387. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2019.105387]

9. Yan, Y. Load characteristic analysis and fatigue reliability prediction of wind turbine gear transmission system. Int. J. Fatigue; 2020; 130, 105259. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2019.105259]

10. Xie, L.; Wu, N.; Qian, W. Time domain series system definition and gear set reliability modeling. Reliab. Eng. Syst. Saf.; 2016; 155, pp. 97-104. [DOI: https://dx.doi.org/10.1016/j.ress.2016.06.009]

11. Wang, X.; Yang, Y.; Wang, W.; Chi, W. Simulating coupling behavior of spur gear meshing and fatigue crack propagation in tooth root. Int. J. Fatigue; 2020; 134, 105381. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2019.105381]

12. Zorko, D. Investigation on the high-cycle tooth bending fatigue and thermo-mechanical behavior of polymer gears with a progressive curved path of contact. Int. J. Fatigue; 2021; 151, 106394. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2021.106394]

13. Vučković, K.; Čular, I.; Mašović, R.; Galić, I.; Žeželj, D. Numerical model for bending fatigue life estimation of carburized spur gears with consideration of the adjacent tooth effect. Int. J. Fatigue; 2021; 153, 106515. [DOI: https://dx.doi.org/10.1016/j.ijfatigue.2021.106515]

14. Alshoaibi, A.M.; Fageehi, Y.A. Simulation of Quasi-Static Crack Propagation by Adaptive Finite Element Method. Metals; 2021; 11, 98. [DOI: https://dx.doi.org/10.3390/met11010098]

15. Mohsenzadeh, R.; Soudmand, B.H.; Shelesh-Nezhad, K. Failure analysis of POM ternary nanocomposites for gear applications: Experimental and finite element study. Eng. Fail. Anal.; 2022; 140, 106606. [DOI: https://dx.doi.org/10.1016/j.engfailanal.2022.106606]

16. Li, Q.; Xie, L. Analysis and optimization of tooth surface contact stress of gears with tooth profile deviations, meshing errors and lead crowning modifications based on finite element method and Taguchi method. Metals; 2020; 10, 1370. [DOI: https://dx.doi.org/10.3390/met10101370]

17. Chen, K.; Ma, H.; Che, L.; Li, Z.; Wen, B. Comparison of meshing characteristics of helical gears with spalling fault using analytical and finite-element methods. Mech. Syst. Signal Process.; 2019; 121, pp. 279-298. [DOI: https://dx.doi.org/10.1016/j.ymssp.2018.11.023]

18. Concli, F.; Maccioni, L.; Fraccaroli, L.; Bonaiti, L. Early crack propagation in single tooth bending fatigue: Combination of finite element analysis and critical-planes fatigue criteria. Metals; 2021; 11, 1871. [DOI: https://dx.doi.org/10.3390/met11111871]

19. Craig, R.R.; Bampton, M.C.C. Coupling of Substructures for Dynamic Analysis. AIAA J.; 1968; 6, pp. 1313-1319. [DOI: https://dx.doi.org/10.2514/3.4741]

20. Ericson, T.M.; Parker, R.G. Experimental measurement and finite element simulation of elastic-body vibration in planetary gears. Mech. Mach. Theory; 2021; 160, 104264. [DOI: https://dx.doi.org/10.1016/j.mechmachtheory.2021.104264]

21. Wang, C.; Dong, B.; Parker, R.G. Impact of planet mesh phasing on the vibration of three-dimensional planetary/epicyclic gears. Mech. Mach. Theory; 2021; 164, 104422. [DOI: https://dx.doi.org/10.1016/j.mechmachtheory.2021.104422]

22. Henry, Z.S. Bell Helicopter Advanced Rotorcraft Transmission (ART) Program; Textron Bell Helicopter: Fort Worth, TX, USA, 2013.