Analytical Investigation on Load Sharing Performance of Planetary Gear Transmission under Loop Maneuver

: In order to explore the influence of airframe space motion on the load sharing performance (LSP) of planetary gear transmission (PGT) inside an aircraft, the combined influence of the rotation of the carrier (internal non-inertial system, INIS) and space motion of the fuselage (external non-inertial system, ENIS) is considered. Firstly, the motion equations of PGT members under the air-frame loop maneuver are derived, and then the time-varying meshing stiffness, backlash, and tooth profile error are considered; a system-level dynamics model of PGT with casing is established and solved by numerical integration. The results show that the loop maneuver intensifies the fluctuation of the dynamic meshing force in the time domain and does not change its spectral distribution, but it will affect the peak value of the low frequency part. After considering the loop maneuver, there are differences in the load sharing coefficient (LSC) of each branch, and the fluctuation range of LSC increases. The interactive influence of working conditions, bearing stiffness, manufacturing eccentricity errors and loop motion parameters on the LSP is investigated and compared in depth. This research can be utilized to guide the design of aircraft’s drivetrain.


Introduction
Compared with ordinary gear boxes, planetary gear boxes have the characteristics of a large transmission ratio, compact structure, strong bearing capacity and stable operation, so they are widely used in the aerospace, shipbuilding and wind power generation fields [1].The uniformity of load distribution between each planet gear is an important prerequisite to give full play to these advantages.However, due to the influence of installation and manufacturing errors, component elastic deformation, temperature and other factors, making the load distribution on each transmission branch of PGT uniform is difficult, which greatly weakens its transmission advantages.Therefore, many scholars have studied the influence of various parameters on the LSP of PGT.
Bodas and Kahraman [2] researched the effects of pin hole position error, gear tooth thickness deviation and transmission error on the static load sharing of PGT.Singh et al. [3] explained the load-sharing mechanism of PGT under floating and non-floating conditions, and gave the generalized calculation form of LSC, but only considered the case of a certain planet gear with position errors.In this regard, subsequently, a uniform epicyclic load sharing map [4] was presented to express the LSP of the system at any position error and torque level.A non-linear bending-torsional coupled model for double-row PGT was proposed by Sheng and Zhu et al. [5]; the influence of eccentricities errors, bearing stiffness, torsional stiffness of the first stage carrier and input rotation rate on the behaviors of the dynamic load sharing characteristics (LSCs) were investigated.Cooley and Parker [6] investigated the stability of high-speed planetary gears considering the gyroscopic effects.Chu et al. [7] studied the influence of gravity, installation angle and ring support stiffness on the LSC of the horizontal axis wind turbine.Park et al. [8] analyzed the effect of helix modification on load distribution over the gear tooth surface and LSCs by establishing the overall system model of the wind power gearbox.Wei and Zhang et al. [9] constructed a dynamic model of PGT supported by journal bearings, and studied the influence of journal bearing characteristics, working conditions and flexible structure on LSCs.From this knowledge, many studies have discussed the LSCs of PGT and considered many factors, like errors, support stiffness, working conditions, gravity, modification, structure, etc.And the method of analyzing the LSCs can be divided into three aspects: numerical theory, virtual prototype simulation and experimental method [10].
Gu and Velex [11,12] established a lumped mass model to explore the influence of errors on LSCs of PGT under different supporting modes.Ren and Qin [13] et al. proposed a new generalized herringbone planetary gear train dynamics model; the lumped parameter method was used to research the influence of errors and floating on LSCs.A new nonlinear dynamic model was established by Han et al. [14] to study the dynamic response of the planetary gear system under nonlinear parameter excitation.The equations of motion are derived by the Lagrangian method and solved by the numerical integration method.Singh et al. [15][16][17] explored the influence of errors, number of planets and support stiffness on LSP by using the finite-contact model.The load-sharing behavior experiment for the PGT with flexible ring gear was formulated by Shen et al. [18], and the test was performed under different working conditions.Li et al. [19] established the calculation model of the NGW planetary gear train and obtained several methods to improve the load sharing, which were verified by experiments.From the foregoing, research on dynamic modeling and the behaviors of the PGT is extensive and in-depth.However, all these studies assume that PGT is supported on a static basis; only the following few academics have considered the motion of the frame.
In 2018, Chu et al. [20] proposed a rotational-translational-axial dynamic model of the planetary gear under pitching base motion.And the influence rules of parameters on response spectra and load sharing conditions were derived, including the amplitude and frequency of the pitching base motion, and the rotating speed of the carrier.In 2019, a dynamic model of the wind turbine drivetrain with platform motions was developed by Zhu et al. [21] who found that the platform motions not only introduce additional excitation frequencies and increased system vibration but also increased the risk of the resonance.In 2020, based on the hovering motion of the airframe, Wei et al. [22] investigated the dynamic behavior of the PGT in a compound non-inertial system, and found that the non-inertial effect could dramatically affect the bearing force, vibration and load sharing performance.The above results indicate that it is necessary to consider the influence of basic motion when studying the dynamics characteristics of PGT.Although the above studies consider the PGT supported on a motion basis, the interactive influence of various factors and the airframe's space movements on LSCs is not clear.
Subject to the capacity of the human body, the acceleration of manned aircraft will not be too large, and modern fighter aircraft can generally withstand an overload of 9 G.In the future, with the further development of science and technology, the application of helicopters, UAVs and other aircrafts will be more extensive and continue to develop towards high-maneuvering flight (as shown in Figure 1); by then, the aircraft may be able to withstand accelerations of 20 G or more.So it is necessary to explore how the space motion of the airframe will affect the dynamic characteristics of the driveline inside the aircraft.For example, how does the space motion affect the meshing of the gear pair, and how does it affect the work of the bearing and impact on the lubrication system, so as to guide the design of the aircraft transmission system.This paper takes the PGT as the research object.Based on the previous research [22], considering the combined effects of the INIS caused by the planet carrier and ENIS caused by the space motion of the airframe, the motion equation of the PGT components under loop motion is deduced, and the time-varying meshing stiffness, backlash, and tooth profile error are considered.Considering the lightweight design of the casing and its rich geometric features, the casing is coupled in the dynamic model, a system-level dynamic model of the PGT with casing is established.Finally, the changes in dynamic meshing force and LSC under different loop motion parameters are analyzed, and the interactive effects of working conditions, bearing stiffness, gear manufacturing eccentricity errors and loop motion parameters on the LSP are investigated.

System Coordinate System Settings
The object of this paper is the PGT system of the main reducer of a certain type of helicopter, as shown in Figure 2. The ring gear is fixed, the 4 planet gears are evenly distributed along the circumference, and the power is input by the sun gear and transmitted to the rotor through the planet carrier shaft.The coordinate system settings are shown in Figure 3.The degrees of freedom (DOF) of the sun gear, ring gear and translational DOF of the planet carrier are defined in the moving coordinate system oc-xcyczc rotating with the planet carrier.The DOF of the planet gears and pins are defined in the coordinate system op-ξηzp; the torsional DOF θc of the planet carrier and DOF of the casing are defined in the coordinate system o0-x0y0z0.Set the coordinate system o1-x1y1z1 to describe the aircraft's rotation around itself, and set the Earth coordinate system O-XYZ to describe the complex attitude of the airframe in space.Among them, the axes of the coordinate system o1-x1y1z1 are, respectively, parallel to the axes of the coordinate system o0-x0y0z0, where o1 is the center of gravity of the aircraft, and The coordinate system of the central member node (taking the sun gear as an example) is shown in Figure 4a; the coordinate system of the pin shaft node is shown in Figure 4b, which is similar to that of the planet gear node.The coordinate system of the casing node is shown in Figure 4c.ω represents the angular velocity of the coordinate system oc-xcyczc around point oc, ω = ωc • kc, ωc denotes the speed of the planet carrier.Ω is the rotational angular velocity of the coordinate system o1-x1y1z1.The radius vector of point o1 in the Earth coordinate system O-XYZ is r0, the radius vector of point oc in o1-x1y1z1 is r1.rs0 = zs0 • kc denotes the initial radial vector caused by the distribution of the nodes of the sun gear shaft along its own axis, zs0 represents the coordinates of node Ms relative to the origin oc, and the radius vector of point Ms in O-XYZ is rMs.The pin shaft node Mpp is located on the equivalent base circle of the planet carrier, and there is an initial vector radius rbc.Pin shaft nodes are distributed in the direction of their own axis opzp, so the moving point Mpp has an initial radial rpp0 similar to rs0, and the radius vector of point Mpp in O-XYZ is rMpp.
Different from the central member, planet gear and pin, when inertial effects under maneuvering flight conditions are not considered, the DOF of the casing nodes are all x 0 y 0 o 0

Airframe coordinate system
Earth coordinate system The geometric relationship of the vectors is shown in Equation ( 1), taking the radius vector of the sun gear node Ms as an example.where i, j and k are the three unit vectors of the coordinate axis O-XYZ, i1, j1 and k1 are the three unit vectors of the coordinate axis o1-x1y1z1, and the three unit vectors of the coordinate axis oc-xcyczc are ic, jc and kc.

Kinematic Analysis
The loop maneuver covers maneuvering actions such as zooming, inverted flight, diving, and Kulbit, so it is an important indicator to measure the maneuverability of an aircraft.As shown in Figure 5, when the aircraft is doing loop motion, it flips around the OX axis of the Earth-fixed coordinate system with the angular velocity Ωd and the turning radius R. The motion parameters are shown in Equation (2).The coordinates of o0 in the coordinate system o1-x1y1z1 are (l1 = 1 m, l2 = 1 m, l3 = 1 m).
Each component of the planetary gear system is defined under its own local coordinate system.In order to obtain the additional force vector excitation of each component derived from the absolute motion of the body, it is necessary to establish a transformation relationship between the coordinate systems, and the base vector relationship of multiple coordinate systems as shown in Equation (3).
The absolute acceleration of each node in dual non-inertial frame can be obtained by taking the derivative of rMs, rMpp, rMx.The author has made a more detailed derivation of the motion equation of the component when the airframe performs any movement in space [22], and will not repeat them here.Then, the additional term aa (remove terms when only INIS condition is considered) derived from the basic motion can be obtained by combining Equations ( 2) and ( 3).Under the loop motion, the additional term of the absolute acceleration of the sun gear shaft node in the coordinate system oc-xcyczc is shown in Equation (4), which is the same as that of the carrier shaft node.The additional term of the absolute acceleration of the ring gear node is also similar to Equation ( 4), but does not contain the term resulting from the axial distribution of the nodes.where aa represents the additional term of absolute acceleration caused by only the ENIS, and the subscript s represents the sun gear.
2 sin sin 0.5sin( 2) The additional terms of absolute acceleration of pin nodes in the coordinate system op-ξηzp are shown in Equation (5), which is parallel to that of the planet gear nodes, except for the term resulting from the axial distribution of the nodes.
where the subscript pp represents the pin, θcp = φN + θc, and φN denotes the position angle of the Nth planet gear.
Since the effects of the casing is not previously considered, the casing nodes are now detailed.When only considering the inertial action in the INIS, the absolute acceleration anx of node Mx is expressed as Equation ( 6), which is obtained by two derivations of the radius vector rnx of Mx versus time.According to Equation ( 6), the initial radius vector rx0 of the casing node has no actual effect.where an denotes the absolute acceleration caused by only the INIS, the subscript x indicates the casing.
Affected by the space motion of the airframe, the radius vector of the casing node is rMx = r 0 + r 1 + r x0 + r x , as shown in Figure 4c.The derivative term of absolute acceleration is obtained by two derivations of the radius vector rMx versus time, as shown in Equation (7), which can continue to derive Equation ( 8) by substituting Equation ( 2) and Equation (3). ) where the v denotes the velocity vector.
Combined with Equation ( 6), the additional terms of the absolute acceleration of the casing node are all generated by ENIS.
Additionally, the gravity component of each part changes under the condition of loop motion.The time-varying gravity component of each part is shown in the following Equation (9).
where the subscript j1 = s, c, r represents the sun gear, planet carrier and ring gear, respectively, j2 = p, pp, indicates planet gear and pin.

Substructure Modeling
The node finite element models [23] of the shafting substructures considering the flexibility are established by Timoshenko beam elements [24].Among them, the sun gear shafting is divided into 5 nodes, the pin shafting is divided into 4 nodes, and the planet carrier shafting is divided into 11 nodes, while planet gear and ring gear are modeled as rigid bodies.According to the structure and size of the shaft segment, the beam elements can be divided into equal section beam elements and variable section beam elements (this article is mainly a tapered beam), as shown in Figure 6.The casing is the basic bearing part of the helicopter's main reducer.The designer's strict pursuit of a high power-to-weight ratio and light weight make the assumption of large stiffness of the casing no longer applicable.In addition, to meet the needs of load bearing, installation, vibration reduction, lubrication, and light weight, the casing is arranged with rich geometric features such as bearing holes, mounting ears, restraint holes, and reinforcing ribs, which are typically complex special-shaped components.
Due to the complex structure and rich geometric features of the casing, the finite element method is usually used to model the casing.However, the original degrees of freedom of the finite element model is large, which leads to the low efficiency of the dynamics calculation.Therefore, using polycondensation technology [25], a high-precision model of the casing substructure is established, as shown in Figure 7.The upper and lower casings are set with 38 and 16 nodes, respectively, and the frequency error of the first 12 orders of the casing super-unit model is within 5.2%.The basic parameters of casing materials are shown in Table 1.

Substructure Coupling
Internal coupling of shafting: Figure 8 shows that the sun gear, planet gear and ring gear are coupled through the meshing relationship; the planet gear and pin are coupled through the bearing; the side plate and pin are coupled through the spatial relationship, as shown in Figure 9.  where N denotes the installation phase angle of the Nth planet gear.αsp and αrp represent the mesh angle of the sun-planet and planet-ring, respectively.rbs, rbp, and rbr are the base circle radius of the sun gear, planet gear, and ring gear, respectively.And rbc denotes the equivalent radius of the planet carrier.
Coupling of shafting and casing: the sun gear shaft and lower casing; the carrier and upper casing are respectively coupled through bearings; the ring gear is coupled with the upper casing and lower casing through bolts.The bolt and bearing connection equation are as follows.2; Fi and Fj denote the force vector of the casing and shafting, respectively.Finally, the system-level dynamic model of PGT with the 472-DOF (37 × 4 + 54 × 6) shown in Figure 10 is formed, and the dynamic meshing force of each gear pair is obtained by utilizing the numerical integration [23].The system dynamics equation is as follows: where M, C and K are the mass matrices, damping matrices and stiffness matrices of the system as a whole; A(t) represents the absolute acceleration; () X t denotes the displace- ment vector considering nonlinear factors; T(t) indicates the external excitation; G(t) represents the gravity vector.When there is only internal non-inertial action, A(t) and G(t) are shown in Equation (13).The gravity influence is constant under such conditions.

Internal Dynamic Excitations
The basic parameters of the gear are shown in Table 3.According to the principle of potential energy [26], the comprehensive time-varying mesh stiffness of external meshing and internal meshing can be obtained, which includes bending stiffness, shear stiffness, axial compression stiffness, equivalent stiffness of gear body deformation and Hertz contact stiffness of the gear pair.The total cumulative pitch error and single tooth tangential deviation are fitted to the superposition form of simple harmonic functions to simulate the tooth profile error excitation in the meshing process.The minimum and maximum backlash and tooth profile error are shown in Table 4.

Analysis of LSP under Airframe Loop Maneuver
Based on the system dynamics model and kinematic model of the airframe loop motion, the meshing forces Fspi(t) and Frpi(t) of the external and internal meshing can be obtained by solving the system dynamics equation with mathematical software Matlab2019.Let LSCspik1 and LSCrpik2 denote the LSC of the external and internal meshing in each meshing frequency cycle, respectively; then, the LSC of each meshing frequency cycle is as follows [5].
The LSC of the i-th external and internal meshing gear pairs over a period of time are respectively: The LSC of the PGT is shown in Equation (17).The larger the LSC, the more uneven the load distribution of the transmission system.

( ) ( )
It should be noted that the calculated data are processed in the software Matlab2019, Origin2021 and Microsoft Visio2019 to obtain the following figures.

Influence of Loop Maneuver on LSC
The input speed of the sun gear is constant at 4000 r/min, the power is 285 kW, and the angular acceleration of the loop motion is 0 = d Ω .When discussing the effect of loop angular velocity (LAV) on LSC, the keep ?? loop radius (LR) R = 10 m; when discussing the effect of LR on LSC, the ?? keep Ωd = 3 rad/s (similarly hereinafter), and "Static" in legend indicates that the aircraft is immobile without space movement.
Figure 11 presents the variation of the dynamic meshing forces of sun-planet1 with gear mesh time.When the aircraft is at rest, the dynamic meshing force fluctuation is small, and its peak-to-peak value is 3554.When the body is in loop motion (3 rad/s,10 m), the peak-to-peak value of the meshing force becomes 3952, with a change rate of 11.12%, and the fluctuation is intensified.Comparing Figure 11a,b, the LAV has a greater influence on the volatility of the dynamic meshing force than the LR.When Ωd = 6 rad/s and R = 10 m, the centripetal acceleration is 360 m/s 2 , which is equal to the situation of Ωd = 3 rad/s and R = 40 m.The peakto-peak value of 3958 in the case of Ωd = 3 rad/s and R = 10 m was taken as the reference value, and the peak-to-peak value of the dynamic meshing force in the former case was 4446, which increased by 12.33% compared with the reference value.In the latter case, the peak-to-peak value is 3994, an increase of 0.91% over the reference value.
Figure 12 shows the maximum and minimum values of the dynamic meshing force within 0.3s under different loop parameters.As the LAV increases, the maximum and minimum values of the meshing force tend to increase and decrease, respectively.However, the increase in the LR has little effect on the maximum and minimum values of the meshing force.This difference can be seen in Equation ( 4), when the angular acceleration d Ω is 0, the R only plays a role in the z term, while Ωd exists in all three terms, x, y and z.The meshing force is applied in the x-y plane.Figure 13 shows the frequency spectra for the dynamic mesh force of the sun-planet1, which is obtained by using the fast Fourier transform.Table 5 shows the various excitation frequencies at an input speed of 4000 r/min.It can be observed that some main frequencies 16.67, 1716, 3432 Hz, and so on, emerge in the frequency spectra of the sun-planet1.At the meshing frequency (1716 Hz), the peak value of the frequency spectra is the highest, which indicates that the gear meshing frequency is the direct current component of the dynamic mesh force.Compared with when the airframe is stationary, after considering the body's loop movement, the frequency distribution does not significantly change, and no new frequency components appear.However, the peak value of the frequency spectrum increases significantly in the low frequency part.When the aircraft is immobile, the peak value at frequency 16.67 Hz is 59.99, and when doing loop motion (Ωd = 3 rad/s, R = 10 m), the peak value at this frequency becomes 253.8, an increase of 323.07%.At the same centripetal acceleration, the LAV has a greater effect on the peak at low frequencies than the LR.
Figure 13c shows the frequency spectrum of the low-frequency part at different sun gear rotating speeds under the same loop motion parameters (Ωd = 3 rad/s, R = 10 m).The low frequency part mainly has two frequencies; one is the sun gear shaft frequency, whike the other is close to the carrier shaft frequency.The peak value corresponding to the shaft frequency of the sun gear does not change with the change in the loop motion parameters.
At a frequency of 16.67Hz, the peak value increases with the LAV as shown in Figure 14a, showing a parabolic growth trend, and the peak value increases almost linearly with LR, as shown in Figure 14b.Using Equation ( 15), the LSC of each mesh pair in each meshing frequency cycle is calculated, and the LSC changes with the number of meshing frequency cycles as shown in Figures 15 and 16.The example of external meshing is analyzed.When the airplane is stationary, the LSC shows periodic changes, and the LSC curve of each mesh pair is basically the same, but there is a phase difference.After considering the loop maneuver, it is obvious that there are differences in the LSC of each branch.The maximum value of the LSC of sun-planet3 is larger than that of the other routes, while the maximum value of the LSC of sun-planet4 is the smallest.
As the LAV increases, the fluctuation of each curve becomes more severe.The range of LSC in the number of meshing frequency cycles shown in the figure changes from 0.9564~1.0410when the aircraft is at rest to 0.9202~1.0765(Ωd = 3 rad/s), 0.8763~1.1204(Ωd = 6rad/s), 0.8105~1.1888(Ωd = 9rad/s).?? Keeping Ωd = 3 rad/s.The range of LSC changes as follows with the increase of LR: 0.9171~1.0807(R = 40 m), 0.9126~1.0849(R = 70 m), 0.9081~1.0892(R = 100 m).Increasing the LAV has an obvious effect on the load sharing coefficient, while changing the LR has a small effect on the load sharing coefficient.For ease of analysis and simplicity, the load distribution balance of PGT is generally expressed in terms of the LSC over a period of time, see Equation (17).
Figure 17a shows the effect of LAV on the root-mean-square value of the dynamic meshing force and LSC.When LAV is increased from 0 rad/s to 12 rad/s with an interval of 0.5, the meshing force of sun-planet2 and sun-planet3 increases in the of a parabola; while meshing force of sun-planet4 decreases nonlinearly.The four curves tend to deviate from each other, which inevitably destroys the load sharing performance.The LSCsp increase nonlinearly with the increase of LAV.As the LR increases with an interval of 5, the meshing force of sun-planet1 and sun-planet4 decreases nonlinearly; on the contrary, the meshing force of sun-planet2 and sun-planet3 increases nonlinearly.The difference in the meshing force of each gear pair gradually increases, which leads to an almost linear LSCsp increase versus the increase of LR, as shown in Figure 17b.
Under the aircraft's loop maneuver, each component is subjected to multiple effects of INIS, ENIS and gravity G.In order to further clarify the influence of each factor on LSC, each factor is separately applied to the system on the basis of only considered INIS.As shown in Figure 18, gravity has little influence on the system LSC, while ENIS has a significant influence on the LSP.Observing the enlarged view in the figure, when INIS and gravity are considered simultaneously, LSC changes with LAV, but it remains constant with LR.This is because at the same running time and different LAV, the θd value is different, resulting in different components of gravity on each coordinate axis at the same time.

Influence of Working Conditions on LSC
The change trend of LSC with input speed during loop is shown in Figure 19, where the load torque is 2800 Nm when the input speed is changed, and the input speed is 4000 r/min when the load torque is changed.The variation of LSC with input speed is very complex, and there are many peaks and troughs.In particular, there is a mutation in the LSC around 4500 r/min, which may produce resonance.The change trend of LSC with the raise of input speed is affected by loop maneuver.When the aircraft is stationary, LSC reduces with the increase in rotational speed after the input rotational speed reaches 5500 rpm; while when Ωd = 3 rad/s, LSC increases with the enhancement of rotational speed after the rotational speed reaches 5700 rpm.?? Keep Ωd = 3 rad/s, change R, the trend of each curve is shown in Figure 19b.The LSC reduces nonlinearly with the raise in load torque, as shown in Figure 20.Additionally, the difference between the LSC corresponding to different LAV and LR decreases with the enhancement of the load torque.When the torque equals 1.441 kNm, the difference in LSC between "static" and "6 rad/s" is 0.223; while when the torque increases to 7.824 kNm, the difference decreases to 0.034.So the rise in the load torque suppresses the influence of the additional effect caused by loop maneuver on the LSC.Therefore, the input torque of the reducer can be designed higher.

Influence of Bearing Stiffness on LSC
The PGT contains a sun bearing, four planet bearings and two carrier bearings, and the positions of each bearing are shown in Figure 2. By setting the bearing stiffness kbii multiplied by the stiffness ratio Kkb (0 < Kkb ≤ 1), the effect of the bearing stiffness variation on the LSC under loop maneuver is studied, as shown in Equation ( 18), and kbii is shown in Table 2. Since changing the LR alone has little effect on the LSC, the subsequent study will not analyze the change in the LR.
where ii = s, p, c1, c2.kbii and kbiir represent the initial bearing stiffness and reduced bearing stiffness, respectively.Explaining the legend in Figure 21, "Sun bearing Kkb = 0.7" refers to reducing only the sun bearing stiffness kbs; "Planet bearing Kkb = 0.7" means simultaneously reducing the stiffness of the four planet bearings (kbp1, kbp2, kbp3, and kbp4); "Sun + Planet bearing Kkb = 0.7" denotes reducing both the sun bearing stiffness kbs and planet bearing stiffness kbp; "Carrier bearing Kkb = 0.7" refers to reducing only the carrier bearing stiffness kbc (kbc1 and kbc2); "All bearings Kkb = 0.7" means all bearing stiffness multiplied by the same stiffness ratio.According to the enlarged figure in Figure 21, when the LAV is at 0-3 rad/s, reducing the stiffness of different bearings to 0.7 times the original can reduce the LSC to a certain extent.The effect of reducing the stiffness of the sun gear bearing is not very obvious, and the effect of reducing the stiffness of all bearings is the best.But with the increase in the LAV, in addition to reducing the stiffness of the sun gear bearing, other ways will make the LSC greater than the original.It shows that reducing the sun gear bearing stiffness has a better application range for the optimization of load sharing performance under different loop motion parameters.
The variation trend in LSC versus bearing stiffness under different LAV is shown in Figure 22. Figure 22a shows that the LSC increases with the enhancement of the stiffness of the sun gear bearing, and the change rate of the curve slightly increases with the increase in the LAV.
1.035  The change trend of LSC with kbp and "kbs + kbp" under loop motion is significantly changed in contrast with that when the aircraft is stationary according to Figure 22b,c.LSC undergoes a change of decrease, increase, decrease and increase as kbp and "kbs + kbp" increase when the airframe is static.While Ωd = 3 rad/s, R = 10 m, the LSC increases slightly at first and then slowly decreases as the stiffness of the planet gear bearing increases; when Ωd = 5 rad/s, 6 rad/s, the LSC hardly changes between Kkb = 0.5~0.6,which is very different from that when the aircraft is stationary.
When Ωd = 3 rad/s, 4 rad/s, the LSC increases first with the increase of "kbs + kbp", then decreases, and finally slowly increases; when Ωd = 5 rad/s, 6 rad/s.LSC increases with the raise of "kbs + kbp", but when Kkb = 0.7~0.95, the curve change rate is very small.

Influence of Manufacturing Errors on LSC
Manufacturing error refers to the deviation in the geometric center, center of mass and the actual center of rotation of the component generated during the processing of the gear.The offset direction of the actual center of rotation and center of mass relative to the theoretical geometric center of the component rotates with the rotation of the gear.The deviation between the geometrical center and actual center of rotation will produce displacement excitation, and the deviation between the center of mass and actual center of rotation will produce centrifugal force, resulting in an unbalanced excitation.
Taking the meshing pair of the sun-planet as an example, as shown in Figure 23.Os is the geometric center of the sun gear, Os1 is the actual center of rotation, and Os2 is the center of mass.ems denotes the error of the actual rotation center of the sun gear and γs is its phase angle; eds represents the manufacturing error for the center of mass of the sun gear, and ∅s is its phase angle.Both are defined in the moving coordinate system Osf-xsf ysf zsf rotating with the sun gear.The projection of the manufacturing error of the actual rotation center of the sun gear, planet gear and ring gear on the external and internal meshing line is shown in Equation ( 19) [27].The planet gear manufacturing errors include the carrier manufacturing eccentricity errors.
where emp and emr represent the manufacturing error values of the planet gear and ring gear, respectively; γp and γr denote the initial phase angle of the manufacturing error of each gear, respectively.ωsrc = (ωs − ωc), ωprc = (ωp − ωc) and ωrrc = (ωr − ωc) indicate the rotational angular velocity of each gear relative to the carrier, respectively.
When only INIS is considered, the manufacturing error for the center of mass of the member will lead to an unbalanced force Fe1, as shown in Equation (20).where mw1 represents the mass of the member, and w1 = s, p, r, and ae1 denotes the absolute acceleration, as shown in Equation (23).Substituting Equation (22) into Equation ( 21) can be further expressed as Equation (23).
( 2 ) When ENIS is considered, the manufacturing error will add an unbalanced force Fe2, as shown in Equation (24).It is assumed that the actual center of rotation coincides with the geometric center, and only the deviation between the center of mass and actual center of rotation is considered.When studying the influence of the manufacturing error of a certain gear on the LSC, the manufacturing error value of other gears is kept at 5 μm and the initial phase angle is 20°.When studying the influence of the planet gear manufacturing error on LSC, the four planet gears are set to the same.
Figure 24 shows the time domain diagram of dynamic meshing force when the sun gear manufacturing eccentricity error and loop motion are considered or not.When the airplane is stationary and manufacturing eccentricity error is 0 μm, the peak-to-peak value of the dynamic meshing force is 3605.When the manufacturing eccentricity error is set to 120 μm, the peak-to-peak value of the dynamic meshing force is 4201, an increase of 16.53%.When the influence of the loop movement is reconsidered, the peak-to-peak value of the dynamic meshing force is 4927, which increases by 36.67%.It can obviously be found that the mesh force fluctuations with manufacturing errors are significantly larger than those without manufacturing errors.The time-domain data in Figure 24 were subjected to fast Fourier transform to obtain the spectrum diagram in Figure 25, which adds a data case, that is, considering the loop motion, the manufacturing eccentricity error of the sun gear is 11 μm.This is because in the previous analysis, the manufacturing eccentricity error of the sun gear is set to this value.Considering the eccentricity error of the sun gear, a new frequency is introduced, which is close to the relative shaft frequency frs1 of the sun gear.When the manufacturing eccentricity error is very small, the corresponding peak value of this frequency is tiny.As shown in Figure 26, when the aircraft is immobile, LSC increases with the raise of the manufacturing error of the sun gear and planet gear, while it first decreases and then increases with the enhancement of the manufacturing error of the ring gear.The increase in LAV changes the trend of LSC versus the gear manufacturing error.Under the condition of loop motion, LSC decreases first and then increases with the enhancement of the manufacturing error of the sun gear, and the extreme point will move to the right as the LAV increases, as shown in Figure 26a.Therefore, after considering the loop movement, the eccentricity error value of the sun gear is not as small as possible.As shown in Figure 26b, although LSC increases as ep increases at different LAV on the whole, when the manufacturing eccentricity error is in a small range (ep = 0-15 μm), LSC basically remains unchanged.At each LAV, LSC increases almost linearly with the raise of ring gear manufacturing error, as shown in Figure 26c.Therefore, the eccentricity error of planet gear and ring gear should be in a small range.

Conclusions
In this study, a dynamic model of a PGT system with a casing under the non-inertial system of loop maneuver is established, considering several excitations such as time-varying mesh stiffness, backlash and tooth profile error.The influence of the loop maneuver on the LSP is studied, the interactive effects of working conditions, bearing stiffness, gear manufacturing errors and loop motion parameters on the LSP are investigated, and the following conclusions are drawn: (1) The loop movement intensifies the fluctuation of the dynamic meshing force in the time domain, and the LAV has a greater influence on the volatility of the dynamic meshing force than the LR.Loop motion does not change the spectral distribution of the dynamic meshing force but will significantly affect the peak value of the low frequency part.After considering the loop maneuver, it is obvious that there are differences in the LSC of each branch, and the fluctuation range of LSC increases.The LSC increases nonlinearly with the rise in the LAV, while it increases almost linearly with the enhancement of the LR.
(2) The LSC decreases with the rise in load torque, and the increase in the load torque will reduce the difference between the LSC under each loop motion parameter.Therefore, the input torque of the reducer can be designed higher.In general, the LSC increases with the raise in the input speed.But the variation trend in a certain input speed range is changed by the loop maneuver.
(3) Only when the LAV is in a relatively small range can reducing the bearing stiffness improve load sharing.It is worth noting that reducing the sun gear bearing stiffness has a better application range for the optimization of load sharing performance under different loop motion parameters.The change trend of LSC with kbp and "kbs + kbp" under loop maneuver is changed in contrast with that when the aircraft was static.
(4) Unbalanced force is added to the gear manufacturing errors under the INIS and ENIS.The manufacturing eccentricity errors aggravate the fluctuation of dynamic meshing force and introduce a new excitation frequency.The variation law of the LSC versus the manufacturing errors is affected by the loop maneuver.After considering the loop movement, the eccentricity error value of the sun gear is not as small as possible, but the eccentricity error of the planet gear and ring gear should be in a small range.

Conflicts of Interest:
The authors declare that they have no conflict of interest regarding the publication of this paper.

Figure 3 .
Figure 3. Multi-coordinate system.(a) The internal coordinate system of the PGT; (b) External coordinate system.

Figure 4 .
Figure 4.The coordinate system of the component nodes: (a) Sun gear node; (b) Pin node; (c) Casing node.

Figure 9 .
Figure 9. Spatial coupling relationship between pin and side plate.

Figure 10 .
Figure 10.System-level coupled dynamics model of PGT.

Figure 11 .
Figure 11.The time domain responses of the dynamic meshing force of sun-planet1.(a) At different LAV; (b) At different LR; (c) Partial enlargement.

Figure 12 .
Figure 12.The maximum and minimum values of dynamic meshing force within 0.3 s.(a) At different LAV; (b) At different LR.

Figure 13 .
Figure 13.Frequency spectra of the dynamic meshing force of sun-planet1.(a) At different LAV; (b) At different LR;.(c) At different input speeds of sun gear.

Figure 14 .
Figure 14.Spectra amplitude of the dynamic meshing force of sun-planet1 at 16.67 Hz.(a) At different LAV; (b) At different LR.

Figure 17 .
Figure 17.LSC and dynamic meshing force mean versus different loop motion parameters.(a) At different LAV; (b) At different LR.

Figure 18 .
Figure 18.The influence of different terms on LSC.(a) At different LAV; (b) At different LR.

Figure 19 .
Figure 19.Influence of input speed on LSC under loop maneuver.(a) At different (b) different LR.

Figure 20 .
Figure 20.Influence of load torque on LSC under loop maneuver.(a) At different LAV; (b) At different LR.

Figure 21 .
Figure 21.The effect of reducing the stiffness of each bearing on the LSC at different LAV.

Figure 22 .
Figure 22.Effect of reducing bearing stiffness on LSC at different LAV.(a) Sun bearing; (b) Planet bearing; (c) Sun + Planet bearing.

Figure 24 .
Figure 24.The time domain responses of the dynamic meshing force of sun-planet1 with manufacturing errors.

Figure 25 .
Figure 25.Frequency spectra of the dynamic mesh force of sun-planet1 with manufacturing errors.

Figure 26 .
Figure 26.Influence of gear manufacturing eccentricity error on LSC at different LAV.(a) Sun gear; (b) Planet gear; (c) Ring gear.

i0, j0 and k0 are
the inertial coordinate system o0-x0y0z0.The radius vector of the casing node Mx: rnx = rx0 + rx.rx0 = xx0 • i0 + yx0 • j0 + zx0 • k0 represents the initial constant radius vector generated by the spatial distribution of nodes Mx, and xx0, yx0, zx0 indicate the initial coordinates of the node Mx; the unit vectors of the coordinate axis.rx = xx • i0 + yx • j0 + zx • k0, xx, yx and zx denote the real-time dynamic displacement response value of node Mx, and the radius vector of point Mx in O-XYZ is rMx. )

Table 1 .
Material properties of the casing.

Table 3 .
Basic parameters of PGT.