Dynamic characteristics of aircraft gear transmission systems under overload level-flight conditions

1. Introduction

Gear transmission systems are widely used in core power systems of aircraft due to their high-power density, strong load-bearing capacity, and precise motion transmission characteristics [1]-[5]. Due to the complexity of internal excitations such as time-varying meshing stiffness, backlash, and transmission errors, coupled with long-term service under extreme conditions such as high-speed and heavy loads, gear systems exhibit significant nonlinear vibration characteristics and cause prominent dynamic problems. In severe cases, this can lead to gear failure, system malfunction, and ultimately affect the reliability of the entire system [6-8]. Notably, international aviation giants such as Sikorsky and Boeing have continuously increased their research and development investment in transmission systems, focusing on three core technological areas: high-speed operational stability, lightweight design, and vibration and noise control [9]. This reflects increasing industry requirements for dynamic stability, structural reliability, and acoustic performance of transmission systems. Therefore, investigating the dynamic characteristics of aviation gear transmission systems is essential for enhancing the technical performance and safe operation level of aviation equipment.

Regarding dynamic modeling methodologies for gear transmission systems, Wang et al. [10] established a fully coupled bending-torsional-lateral-pendular model and a bending-torsion coupled model for two-stage helical gear transmission systems, considering shaft flexibility. By comparing the effects of these two models on the system’s vibrational characteristics, their research findings demonstrate that lateral and pendular motions cannot be ignored in the precise vibration analysis of helical gear pairs. Yuan et al. [11] established a full-degree-of-freedom model for helical gear transmission systems using Timoshenko beam elements. Dynamic response results indicate that: Manufacturing errors induce sidebands in the vibration response spectrum; under light-load conditions, sidebands dominate the spectrum, while increasing torque enhances the prominence of meshing frequency and its harmonics relative to the side-bands. Li et al. [12] proposed a nonlinear dynamics model for spur gear pairs incorporating time-varying backlash, multi-state meshing and contact ratio. This model established nine safety domains to quantify system stability levels, enabling safety-basin-based analysis of bifurcation effects on operational security.

Beyond dynamic modeling methodologies, internal excitation parameters stand out as key factors influencing the dynamic response of gear systems and have received considerable scholarly attention. Core internal excitations in gear systems primarily encompass time-varying mesh stiffness, meshing damping, transmission error, and backlash. Current mathematical representations show relative consensus for meshing damping and backlash: meshing damping is predominantly modeled via Coulomb damping theory [13], while backlash is universally characterized by piecewise nonlinear displacement functions [14]. In contrast, precise modeling of mesh stiffness and transmission error remains an active research frontier. Cooley et al. [15] compared the averaged slope method and local deflection slope method for calculating mesh stiffness, investigating their respective applicability in static and dynamic contexts. Wei et al. [16] developed an analytical mesh stiffness model for helical gears incorporating axial deflection, enhancing computational efficiency through sliced tooth contact analysis. Duan et al. [17] investigated time-varying transmission error under varying load, rotational speed, and housing mass conditions, considering the effects of flexible shafts and housing, with experimental validation. For refined contact behavior characterization, Abruzzo et al. [18] demonstrated viable integration of FEA-derived stiffness/error data into lumped-parameter models, achieving strong correlation between nonlinear response predictions and experimental validation; the study concurrently elucidated torque and speed dependencies on dynamic transmission error through experimental investigation. Benaicha et al. [19] created a system-level flexible multibody gear contact model incorporating Lagrange formulation, face-to-face contact detection, tooth compliance, and profile deviations. Through combined consideration of global and local deformations, this model achieves smoother peak-to-peak transitions behavior in predicted stiffness and transmission error.

Current dynamic analyses of gear transmission systems predominantly rely on a fixed coordinate system, assuming the gearbox is rigidly mounted to the aircraft fuselage while neglecting additional effects from base motion. However, modern aircraft face increasingly stringent performance demands requiring frequent overload flight maneuvers such as steep turns, zoom climb, and roll escape. These maneuvers subject transmission components to significant inertial excitation loads, triggering strong nonlinear dynamic coupling between aircraft spatial motions and elastic deformations. Under maneuver overloads, generalized inertia forces acting on critical transmission members can induce flexible-body deformations, thereby modifying gear system operational characteristics. Considering its significance, research on non-inertial frame gear dynamics has witnessed a notable surge in recent years.

Yi et al. [20] established a gear-rotor system dynamics model for roll, pitch, and yaw maneuver environments using Lagrange energy method, investigating transient time-frequency responses under maneuvering parameters and mass imbalance excitation. Qiu et al. [21] developed a planetary gear dynamics model during base pitch motion, employing Runge-Kutta methods to solve system responses and identify spectral characteristics. Wei et al. [22] investigated the load-sharing characteristics, contact dynamic stress, and bending dynamic stress of planetary gear transmission systems under aircraft maneuvers including level flight, rolling, and looping. The results indicate that maneuvering flight conditions exacerbate load imbalance in planetary gear trains and elevate dynamic stress levels. Yue et al. [23] constructed a non-inertial frame finite element dynamics model incorporating gear-rotor-bearing assemblies with Timoshenko beam elements, identifying natural frequencies and critical speeds via Campbell diagrams while systematically analyzing complex influences of translational/rotational carrier motion on inherent characteristics and dynamic responses. Ye et al. [24] validated their theoretical model through experimental tests and further investigated the dynamic behavior of gear transmission systems during aircraft hover maneuvers. Specifically, they analyzed the effects of additional inertial forces and gyroscopic moments induced by airframe maneuvering motion on the dynamic responses of critical system components including casings, bearings and shafts under non-inertial frames.

Existing studies predominantly employ lumped-mass or rigid-flexible coupling methodologies, which fail to model all critical transmission components as flexible bodies, neglect real-time dynamic coupling between substructures, and inadequately address the effects of flexible deformations and component interactions. To overcome these limitations, this work establishes a fully-coupled multibody dynamics model encompassing flexible shaft assemblies, deformable gear pairs, deformable bearings, and flexible casing based on multi-flexible-body dynamics modeling. Furthermore, we derive analytical expressions for additional inertial excitations induced by overload level flight and integrate them into the model, systematically investigating dynamic responses and critical component load characteristics under varying overload values. This research provides crucial theoretical underpinnings for high-accuracy dynamic design and performance assessment of gear transmission systems under extreme overload maneuvering conditions.

2. Dynamic modeling of gear transmission systems under overload level flight maneuvers

2.1. Substructure condensation model

To establish a fully flexible multibody dynamics model for the gear transmission system, the system is first divided into three substructures: the driving gear shaft, the driven gear shaft, and the casing. Each substructure is discretized and a component-level finite element model is established, as illustrated in Fig. 1.

Fig. 1Substructure division

Such dense meshing inevitably leads to an excessive number of degrees of freedom (DOFs), resulting in inefficient calculation and incompatibility with system dynamics analysis requirements. Therefore, Guyan condensation [25-27] is adopted to reduce the original finite element model, which greatly lowers the calculation cost while ensuring the necessary accuracy.

The governing equations of motion for each substructure are expressed as Eq. (1):

1

$M\ddot{X}+C\dot{X}+KX=F,$

where $X$ denotes the node displacement vector, $F$ represents the node load vector, and $M$, $C$, $K$ stand for the mass matrix, damping matrix, and stiffness matrix respectively. Based on substructure characteristics, DOFs are partitioned into master DOFs $X_m$ and slave DOFs $X_s$. The master DOFs $X_m$ correspond to external nodes, while the slave DOFs $X_s$ correspond to internal nodes, and the dimension of $X_m$ is much smaller than that of $X_s$.

Derived from static equilibrium conditions, the partitioned stiffness matrix equation is expressed as Eq. (2):

2

$KX=\left[\begin{array}{cc}{K}_{mm}& K_{ms}\\ K_{sm}& K_{ss}\end{array}\right]\left[\begin{array}{c}{X}_m\\ X_s\end{array}\right]=F=\left[\begin{array}{c}{F}_m\\ 0\end{array}\right],$

where, $K_{mm}$ and $K_{ss}$ are formulated as the stiffness matrices corresponding to the DOFs of master and slave nodes, respectively. $K_{ms}={\left(K_{sm}\right)}^T$ is designated as the coupling stiffness matrix. The external load applied to the master node is denoted by $F_m$. Since no external load is applied to the slave nodes, the load vector is set to a zero vector.

Derived from the second equation of Eq. (2), the displacement of slave nodes can be solved as:

3

$K_{sm}{X}_m+K_{ss}{X}_s=0\Rightarrow X_s=-{\left(K_{ss}\right)}^{-1}{K}_{sm}{X}_m.$

The displacement vector can be expressed in terms of the master node displacement vector as:

4

$X=\left[\begin{array}{c}{X}_m\\ X_s\end{array}\right]=\left[\begin{array}{c}I\\ -{\left(K_{ss}\right)}^{-1}{K}_{sm}\end{array}\right]X_m.$

Under the static equilibrium assumption neglecting the inertia of slave nodes, Eq. (4) is substituted into Eq. (1) and premultiplied by $T^T$:

5

$\underset{M^{*}}{\underbrace{T^{T}MT}}{\ddot{X}}_m+\underset{C^{*}}{\underbrace{T^{T}CT}}{\dot{X}}_m+\underset{K^{*}}{\underbrace{T^{T}KT}}{X}_m=\underset{F^{*}}{\underbrace{T^{T}F}}.$

Thus, from Eq. (5), the condensed substructure dynamic equations are obtained as follows:

6

$M^{*}{\ddot{X}}_m+C^{*}{\dot{X}}_m+K^{*}{X}_m=F^{*}.$

To achieve further reduction in DOFs, condensation nodes are established at substructure coupling interfaces. Kinematic constraints between these condensation nodes and the master nodes are enforced using the Multi-point Constraint (MPC) method [28]:

7

$X_c={\sum }_{i=1}^N\left[{\omega }_i{\left(X_m\right)}_i\right],$

where $X_c$ denotes the DOFs of the condensation node, ${\omega }_i$ represents the weight coefficient of the master node $i$, and $N$ indicates the number of coupled master nodes. Functionally, condensation nodes are categorized into five types, as illustrated in Fig. 2:

1) Gear Mesh Condensation Node, located at the center of the gear bore;

2) Power Input Condensation Node, located at the center of input positions in the driving gear shaft core;

3) Power Output Condensation Node, located at the center of output positions in the driven gear shaft core;

4) Bearing Coupling Condensation Node, located at the center of the corresponding shaft hole and the casing bearing bore;

5) Geometric Feature Condensation Node, located at the shaft segments with high slenderness ratios or complex local geometry.

Through Eq. (7), coupling relations between condensation nodes and master nodes are established. This allows the model parameters to be represented by condensation nodes in dynamic modeling, significantly reducing the DOFs while maintaining sufficient precision, thereby enhancing computational efficiency in multibody dynamics.

To account for the effects of gear body flexibility on transmission performance, a flexible gear model is constructed using the tooth surface condensation method to accurately simulate tooth bending deformation, shear deformation, and other effects. Three condensation nodes are equidistantly arranged along the tooth width direction at the intersection line between the tooth width mid-plane and the pitch cylinder surface for each tooth of the helical gear. As shown in Fig. 3, a local Cartesian coordinate system is established on the helical gear end face, with its $y$-axis coinciding with the tooth symmetry axis and its $z$-axis aligned with the shaft bore centroidal axis. The angular displacement $\mathrm{\Delta }\theta$ between adjacent condensation nodes within a single tooth is thus defined as:

8

$\mathrm{\Delta }\theta =\frac{2\pi b}{3p_z},$

where: $b$ is the face width, $p_z$ is the lead. For a right-hand helical gear, the angle between any tooth’s condensation node and the $y$-axis is given by Eq. (9):

9

${\phi }_{i,j}=\frac{2\pi }{z}\left(i-1\right)+\mathrm{\Delta }\theta \left(j-\frac{1}{2}\right),$

where $i$ denotes the tooth index ($i=$1, 2…, $n$), $j$ is the index of the condensation node within the same tooth ($j=$1, 2, 3), and $n$ represents the number of teeth. The Cartesian coordinates of any condensation node are calculated from the angular position as:

10

$x_{i,j}=r_p\cos {\phi }_{i,j}, y_{i,j}=r_p\sin {\phi }_{i,j}, z_{i,j}=\frac{b}{3}\left(j-1\right)+\frac{b}{6},$

where $r_p$ represents the pitch circle radius. The coordinates of gear tooth condensation nodes are then used to generate the condensation node set through computational programming, enabling full-flexibility modeling of the gear system.

Fig. 2Condensation nodes configuration

Fig. 3Condensation nodes of helical gear

2.2. Substructure coupling relationships

The gear transmission system dynamically couples the shaft substructures through gear mesh elements connected to gear mesh condensation nodes. Specifically, this interaction is characterized by equivalent spring-damper elements along the line of action (LoA), where the kinematic constraints are enforced. To quantify the displacement projection mechanism, the time-varying mesh vector $V\left(t\right)$ projects the multi-degree-of-freedom displacement of the gear pair onto the LoA direction. This projection relationship is explicitly defined by the following vector expression:

11

$V{\left(t\right)}^T=\left[\begin{array}{ccc}\sin \left(\alpha -\gamma \right)\cos \beta & \cos \left(\alpha -\gamma \right)\cos \beta & \sin \beta \\ -r_{b1}\sin \left(\alpha -\gamma \right)\sin \beta & -r_{b1}\cos \left(\alpha -\gamma \right)\sin \beta & r_{b1}\cos \beta \\ -\sin \alpha \cos \beta & -\cos \alpha \cos \beta & -\sin \beta \\ -r_{b2}\sin \alpha \sin \beta & -r_{b2}\cos \alpha \sin \beta & -r_{b2}\cos \beta \end{array}\right],$

where: $r_{b1}$ and $r_{b2}$ are respectively defined as the base circle radii of the driving and driven gears, $\alpha$ represents the tooth pressure angle, $\gamma$ denotes the position angle of the driving gear, and $\beta$ is the base circle helix angle.

Taking transmission error $e\left(t\right)$ and backlash $b\left(t\right)$ into account, the relative deformation ${\mathrm{\Delta }}_m$ along the normal direction of the meshing line is corrected to:

12

${\mathrm{\Delta }}_m=\left\{\begin{array}{ll}V{\left(t\right)}^T{X}_{mesh}-e\left(t\right)-b\left(t\right),& V{\left(t\right)}^T{X}_{mesh}-e\left(t\right)>b\left(t\right),\\ 0,& -b\left(t\right)\le V{\left(t\right)}^T{X}_{mesh}-e\left(t\right)\le b\left(t\right),\\ V{\left(t\right)}^T{X}_{mesh}-e\left(t\right)+b\left(t\right),& V{\left(t\right)}^T{X}_{mesh}-e\left(t\right)<-b\left(t\right),\end{array}\right.$

where, $X_{mesh}$ denotes the gear mesh condensation node vector, the formulations of $e\left(t\right)$ and $b\left(t\right)$ are defined per reference [29].

Consequently, the normal dynamic meshing force $F_n$ is governed by:

13

$F_n=k_m\cdot {\mathrm{\Delta }}_m+c_m\cdot {\dot{\mathrm{\Delta }}}_m,$

where, $k_m$ and $c_m$ are defined as the time-varying meshing stiffness and meshing damping coefficient of the gear pair, respectively.

For helical gears, $F_n$ is decomposed into three-dimensional force components, which include the tangential force $F_T$, radial force $F_r$, and axial force $F_{\alpha }$, as shown in Eq. (14):

14

$\left\{\begin{array}{l}{F}_T=F_n\cos \alpha \cos \beta ,\\ F_{\gamma }=F_n\sin \alpha ,\\ F_{\alpha }=F_n\cos \alpha \sin \beta .\end{array}\right.$

The coupling relationship is depicted in Fig. 4.

Fig. 4The gear meshing element

The bearing support element is implemented to establish the coupling relationship between the shafting system and the casing, as illustrated in Fig. 5.

Fig. 5The bearing support element

To characterize its dynamic behavior, the six-degree-of-freedom bearing support stiffness matrix $K_{bi}$ and damping matrix $C_{bi}$ must be established. The reaction force equation for bearing $i$ is expressed as:

15

$F_{bi}=K_{bi}{q}_{bi}+C_{bi}{\dot{q}}_{bi},$

where, the force vector $F_{bi}$, the velocity vector ${\dot{q}}_{bi}$, and displacement vector $q_{bi}$ of bearing $i$ are defined as follows:

16

$\left\{\begin{array}{l}{F}_{bi}={\left[\begin{array}{cccccc}{F}_{b{i}_x}& F_{b{i}_y}& F_{b{i}_z}& F_{b{i}_{\theta x}}& F_{b{i}_{\theta y}}& F_{b{i}_{\theta z}}\end{array}\right]}^T,\\ {\dot{q}}_{bi}={\left[\begin{array}{cccccc}{\dot{x}}_{bi}& {\dot{y}}_{bi}& {\dot{z}}_{bi}& {\dot{\theta }}_{b{i}_x}& {\dot{\theta }}_{b{i}_y}& {\dot{\theta }}_{b{i}_z}\end{array}\right]}^T,\\ q_{bi}={\left[\begin{array}{cccccc}{x}_{bi}& y_{bi}& z_{bi}& {\theta }_{b{i}_x}& {\theta }_{b{i}_y}& {\theta }_{b{i}_z}\end{array}\right]}^T.\end{array}\right.$

2.3. Additional excitation during overload level flight

To characterize the kinematic behavior of aircraft gear transmission systems under maneuvering overload conditions, a hierarchical coordinate framework is established as depicted in Fig. 6. The Earth inertial coordinate system $OXYZ$ serves as the global reference basis. The gearbox-body coordinate system $O_c{X}_c{Y}_c{Z}_c$ is affixed to the aircraft airframe, maintaining continuous coincidence with the aircraft body axes in both position and orientation throughout all maneuvering states. Local coordinate systems $o_1{x}_1{y}_1{z}_1$ and $o_2{x}_2{y}_2{z}_2$ define the instantaneous spatial poses of the driving and driven gear shafts respectively.

Based on the theory of flexible multibody system dynamics, the position vector $r_p$ of an arbitrary point $p$ on flexible substructures in the inertial coordinate system $OXYZ$ is expressed by Eq. (17):

17

$r_p=r_o+A\left(u_o+u_r\right),$

where: $r_0$ denotes the position vector of the origin $o_c$ of the gearbox-body coordinate system $O_c{X}_c{Y}_c{Z}_c$ with respect to the inertial coordinate system $OXYZ$; $A$ represents the rotation transformation matrix from $O_c{X}_c{Y}_c{Z}_c$ to $OXYZ$; $u_o$ designates the position vector of point p in the undeformed state within $O_c{X}_c{Y}_c{Z}_c$, coinciding with point $p_o$; $u_r$ indicates the additional displacement vector from $p_o$ to $p$ due to elastic deformation.

Fig. 6Kinematic analysis

a) Node position vector

b) Hierarchical coordinate framework

By differentiating Eq. (17), the absolute velocity vector $v_p$ of point $p$ in the inertial coordinate system $OXYZ$ is obtained:

18

$v_p={\dot{r}}_o+\mathrm{\Omega }\times \left[A\left(u_o+u_r\right)\right]+A{\dot{u}}_r,$

where: ${\dot{r}}_o$ denotes the translational velocity of the origin $o_c$ in the gearbox-body coordinate system $O_c{X}_c{Y}_c{Z}_c$; $\mathrm{\Omega }$ represents the angular velocity vector of $O_c{X}_c{Y}_c{Z}_c$; $\mathrm{\Omega }\times Au$ corresponds to the rotational convective velocity, where $u=u_o+u_r$ is the position vector of $p$ in $O_c{X}_c{Y}_c{Z}_c$; $A{\dot{u}}_r$ defines the relative velocity due to elastic deformation.

Differentiating Eq. (18) gives acceleration vector $a_p$ of point $p$ on the flexible body in the inertial coordinate system $OXYZ$:

19

$a_p={\ddot{r}}_o+\dot{\mathrm{\Omega }}\times Au+\mathrm{\Omega }\times \mathrm{\Omega }\times Au+2\mathrm{\Omega }\times A{\dot{u}}_r+A{\ddot{u}}_r$

where: ${\ddot{r}}_o$ is defined as the translational acceleration of the origin $o_c$ in the gearbox-body coordinate system $O_c{X}_c{Y}_c{Z}_c$; $\dot{\mathrm{\Omega }}\times Au$ denotes the tangential acceleration due to angular acceleration; $\mathrm{\Omega }\times \mathrm{\Omega }\times Au$ represents the centripetal acceleration; $2\mathrm{\Omega }\times A{\dot{u}}_r$ is the Coriolis acceleration; $A{\ddot{u}}_r$ defines the relative acceleration induced by elastic deformation.

According to d’Alembert’s principle, the expression for the additional inertial force $F_{ft}$ at any node $i$ in the gear transmission system is given by:

20

$F_{fi}=-m_i\cdot a_{fi},$

where $m_i$ is the mass of node $i$, and $a_{ft}$ is the absolute acceleration vector of node $i$ in the inertial system $OXYZ$, which determined by Eq. (18). Three-dimensional spatial maneuvers can be decomposed into a combination of level flight, circular motion, and rolling motion [30]. Under variable-speed level flight condition ($\mathrm{\Omega }=$ 0), the additional inertial acceleration comprises solely translational acceleration ${\ddot{r}}_o$ and the relative acceleration of elastic deformation $A{\ddot{u}}_r$. Considering gravitational force $G_i$, the resultant force components at node $i$ are expressed as Eq. (21):

21

$\left\{\begin{array}{l}{F}_{fi}=-m_i\cdot \left({\ddot{r}}_o+A{\ddot{u}}_r\right)=-m_i\cdot a_{ei},\\ G_{fi}=m_i\cdot g\cdot k,\end{array}\right.$

where, $g$ is defined as the gravitational acceleration, and $k$ represents the vertically upward unit vector.

2.4. System coupling dynamic model

Considering gear transmission systems operating under overload level flight conditions, the complete coupled dynamic equations are established in the inertial coordinate system $OXYZ$. This model integrates substructure condensation matrices, time-varying gear mesh characteristics, bearing support stiffness and maneuver flight effects, expressed as:

22

$\left[\begin{array}{cc}{M}_s& 0\\ 0& M_j\end{array}\right]\left[\begin{array}{c}{\ddot{X}}_s\\ {\ddot{X}}_j\end{array}\right]+\left[\begin{array}{cc}{C}_s+C_m\left(t\right)+C_{b11}& -C_{b12}\\ -C_{b21}& C_j+C_{b22}\end{array}\right]\left[\begin{array}{c}{\dot{\tilde{X}}}_s\\ {\dot{X}}_j\end{array}\right]+\left[\begin{array}{cc}{K}_s+K_m\left(t\right)+K_{b11}& -K_{b12}\\ -K_{b21}& K_j+K_{b22}\end{array}\right]\left[\begin{array}{c}{\tilde{X}}_s\\ X_j\end{array}\right]=T\left(t\right)+G\left(t\right)+F_f\left(t\right).$

In the left-side dynamic terms, $M_s$, $C_s$, and $K_s$ are denoted as the condensed mass, damping, and stiffness matrices of the gear shaft substructure; $M_j$, $C_j$, and $K_j$ are represented as the condensed matrices of the casing substructure. $C_m\left(t\right)$ and $K_m\left(t\right)$ are defined as the time-varying meshing damping matrix and time-varying meshing stiffness matrix of the gear pair, respectively. $K_{b11}$, $K_{b22}$, $-K_{b12}$ and $-K_{b21}$ are denoted as the corresponding block matrices of the bearing stiffness, respectively; ${\dot{\tilde{X}}}_s$ and ${\tilde{X}}_s$ are defined as the modified velocity and displacement vectors with transmission error coupling.

In the right-side excitation terms, $T\left(t\right)=\left[\dots 0{ T}_{in}\left(t\right)\dots 0 T_{out}\left(t\right)\right]$ is denoted as the driving/driven torque vector, in which the non-zero terms correspond to the input and output nodes. $G\left(t\right)$ is defined as the overall system gravity vector, and $F_f\left(t\right)$ represents the overall system additional inertial force, whose nodal component values are shown in Eq. (21).

The gear shaft systems are coupled through gear meshing elements, and the shaft systems are coupled with the casing through bearing elements, forming a coupled dynamic model of a fully flexible gear-shaft-bearing-casing multibody system. The coupling topological relationships are illustrated in Fig. 7.

Fig. 7Coupled topology of fully flexible gear multibody dynamic model

3. Internal excitation dynamics in aircraft transmissions under overload level flight maneuvers

Focused on a specific transmission stage within an aeroengine accessory gearbox (AGB), this study investigates its dynamic characteristics under different maneuvering overload conditions during level flight, with the system operating at rated conditions of 1500 r/min input speed and 300 kW input power. As shown in Fig. 8, without loss of generality, the maneuvering overload acceleration $a_e$ is applied tangentially to the radial gear mesh force direction, and relevant parameters are provided in Table 1.

Fig. 8Overload level flight maneuver schematic diagram

The vibration response of the power input node at the driving shaft under overload level flight conditions (–50 g to +50 g) is presented to exhibit typical dynamic characteristics. It should be noted that the ±50 g range represents theoretical boundaries designed to magnify the system's sensitivity to inertial loads for parametric analysis, rather than practical aircraft operational envelopes. Extending the acceleration limits theoretically allows for a clearer elucidation of the dynamic response trends.

Time-domain analysis (Fig. 9(a)) shows that the mean value of vibration displacement is symmetric about the static condition (0 g), and a quasi-linear offset is exhibited with the overload acceleration. The maximum dynamic offset is found to reach 70.59 μm, with a 128 % increase achieved compared to the static reference value, while the fluctuation range of vibration displacement is maintained stably within ±3.83 μm. Frequency-domain analysis (Fig. 9(b)) further reveals that the dominant vibration energy is concentrated on the second order harmonic ($2f_m=$1150 Hz) and third order harmonic ($3f_m=$1725 Hz) of the meshing frequency, while the fundamental meshing frequency $f_m$ (575 Hz) is suppressed.

The distribution characteristics of the time-domain displacement of the driven shaft (Fig. 10(a)) are shown to be identical to those of the driving shaft, yet the vibration amplitude is significantly increased. The mean value of the maximum offset under 50 g is found to reach 81.13 μm, exceeding that of the driving shaft by 14.9 %. The frequency-domain spectrum (Fig. 10(b)) is shown to reveal the coexistence of the shaft rotational frequency $f_s$ (25 Hz), the fundamental meshing frequency $f_m$ (575 Hz), and their second-/third order harmonics (1150 Hz/1725 Hz), exhibiting broader spectral energy distribution than the driving shaft. Comparative analysis demonstrates that maneuvering overload substantially modifies the mean displacement characteristics of the shafts. The larger mass of the driven shaft generates stronger inertial forces, resulting in more severe deformation offsets. Spectral analysis demonstrates that although the overload acceleration does not alter the spectral distribution, higher-order harmonics become the predominant frequencies in the vibrational displacement response.

Fig. 9Vibration displacement at power input location of driving gear shaft

a) Time domain

b) Frequency domain

Fig. 10Vibration displacement at power output location of driven gear shaft

a) Time domain

b) Frequency domain

It is revealed by these findings that the fluctuation of the center distance is the key mechanism through which the dynamic response of the system is affected by the overload inertial force. Quantification of this mechanism requires characterizing the center-distance offset evolution against overload acceleration profiles.

The root-mean-square (RMS) center distance offset is shown in Fig. 11. Under static condition (0 g), meshing separation forces induce bidirectional elastic deformation of the gear shafts, yielding a baseline offset of 23.63 μm. With maneuvering overload application, two response regimes emerge: From –50 g to 0 g, the large gear (driven gear) experiences co-directional inertial centrifugal force and radial meshing force ($-x$ direction), while force opposition occurs at the small gear (driving gear). The asymmetric shaft deformation induced by inertial force disparity between gears drives a peak center distance offset of 35.47 μm at –50 g. From 0 g to 50 g, force vector reversal creates counter-directional alignment at the large gear but co-directional at the small gear (+$x$ direction). This inertial-mesh force cancellation, amplified by mass asymmetry, generates a linear reduction in offset to 11.81 μm at 50 g. The 66.7 % reduction in offset amplitude from 35.47 μm to 11.81 μm demonstrates that inertial forces dominate center-distance fluctuations.

Fig. 11Centre distance offset of gear transmission system

The above analysis confirms that maneuvering acceleration loads alter the operating center distance of the gear pair, directly modifying its meshing stiffness characteristics. As acceleration shifts from –50 g to 50 g, the time varying center distance progressively decreases, extending the effective contact length. This geometric change increases the contact ratio, enabling more tooth pairs to participate in load sharing. Consequently, time domain analysis (Fig. 12(a)) demonstrates that the mean time varying meshing stiffness increases from 3.66×10⁹ N/m to 3.68×10⁹ N/m. Frequency domain reveals (Fig. 12(b)) that the fundamental mesh frequency $f_m$ (575 Hz) identified as the dominant stiffness excitation component, which exhibits amplitude attenuation under rising acceleration loads. Conversely, the second harmonic 2$f_m$ (1150 Hz) and third harmonic 3$f_m$ (1725 Hz) are progressively increased with the change of maneuvering overload (Fig. 12(c)), indicating that maneuvering-induced inertial forces suppresses the fundamental meshing frequency, driving spectral energy redistribution from the main frequency to higher-order harmonics.

Under level flight dynamic accelerations, dynamic center distance fluctuations within the gear pair induce kinematic deviation of the actual contact point from the theoretical mesh position, constituting the dominant excitation mechanism for transmission error (TE). As quantified in Fig. 13(a), at –50 g maneuvering overload acceleration, TE is found to increase by an average of 5.13 μm compared to the aircraft stationary condition. As the maneuvering acceleration changes from –50 g to 50 g, the absolute value of transmission error is reduced. At 50 g, TE registers a 5.15 μm decrease versus stationary condition. Spectral decomposition of TE time-series data (Fig. 13(b)) reveals three dominant components: the fundamental meshing frequency $f_m$(575 Hz), second harmonic 2$f_m$ (1150 Hz), and third harmonic 3$f_m$ (1725 Hz). While $f_m$ maintains persistent spectral dominance, its amplitude undergoes progressive attenuation across the –50 g→ 50 g acceleration envelope. Concurrently, harmonic energies intensify, showcasing vibrational energy migration from the fundamental meshing frequency to higher-order harmonics (Fig. 13(c)).

Analysis of Fig. 14(a) confirms that progressive reduction in center distance offset intensifies gear flank engagement, resulting in backlash reduction from 34.3 μm at –50 g to 11.3 μm at 50 g. The resultant magnitude reduction reaches 66.7 %, validating overload-induced suppression of tooth clearance. Spectral decomposition in Fig. 14(b) reveals three dominant components: shaft rotational frequency $f_s$ (25 Hz), fundamental meshing frequency $f_m$ (575 Hz) and associated harmonics. Crucially, $f_m$maintains persistent spectral dominance while preserving identical amplitude and distribution characteristics across all overload conditions. Spectral analysis confirms that maneuvering acceleration negligibly affects the spectral characteristics and amplitude magnitudes in backlash measurements.

Fig. 12Time varying gear mesh stiffness

a) Time domain

b) Frequency domain

c) Amplitude trend

Fig. 13Time varying gear transmission error

a) Time domain

b) Frequency domain

c)

Fig. 14Time varying backlash

a) Time domain

b) Frequency domain

4. Effects of differential overload maneuvers on load characteristics

Analysis establishes that aircraft base motion-induced inertial loads represent external operating conditions for the transmission system, while input/output torque profiles define internal loading constraints. This dual-condition loading model enables the investigation of overload level-flight maneuvers on the dynamic load transfer properties of shaft-bearing systems. By solving the fully coupled dynamic equations of the gear transmission system accounting for the additional effects of maneuvering flight, the dynamic bearing loads at each bearing node under different level-flight acceleration $a_e$ are obtained. All load vectors are transformed into the inertial coordinate system $OXYZ$, with bearing positions illustrated in Fig. 15.

Fig. 15Bearing mounting locations and designations

Fig. 16(a-d) quantify bearing loads dynamics at Bearing 1 (power-input end) of the driving shaft, while Fig. 16(e-h) characterize Bearing 2 (non-input end). Results indicate that the bearing loads of the driving shaft exhibit nonlinear growth with varying accelerations from –50 g to 50 g. The $x$-direction component aligned with the maneuver direction exhibits dominant response: for Bearing 1, the $x$-direction load increases from 2.28 kN at –50 g to 5.17 kN at 50 g, marking 127 % growth (Fig. 16(a)); for Bearing 2, the $x$-direction load linearly increases from 5.58 kN to 8.27 kN, increasing by 48 % (Fig. 16(e)).

$y$-direction bearing loads show negligible sensitivity to acceleration changes, as confirmed by Fig. 16(b) and 16(f). Conversely, $z$-direction load exhibit acceleration-dependent increases (Fig. 16(c) and 16(g)), albeit with significantly lower growth rates versus $x$-direction.

Fig. 16Load characteristics of driving shaft bearings

a) Bearing 1 $x$-direction load

b) Bearing 1 $y$-direction load

c) Bearing 1 $z$-direction load

d) Bearing 1 radial load

e) Bearing 2 $x$-direction load

f) Bearing 2 $y$-direction load

g) Bearing 2 $z$-direction load

h) Bearing 2 radial load

Resultant radial load analysis reveals that the radial load on Bearing 1 increases from 11.47 kN (–50 g) to 12.58 kN (50 g), an increase of 9.7 %, with the load orientation angle rotating clockwise from 78.64° to 65.07°, while confined to Quadrant I (Fig. 16(d)). For Bearing 2, the radial load increases from 12.61 kN (0 g) to 13.82 kN (50 g), an increase of 9.6 %, with the load orientation angle shifting from 63.01° to 52.93° in Quadrant I (Fig. 16(h)). The radial resultant load on Bearing 2 is greater compared to Bearing 1, but Bearing 1 exhibits more pronounced changes in the $x$-direction response.

Analysis of the driven shaft bearing loads (Fig. 17) confirms that acceleration variations from –50 g to 50 g induce significant load redistribution. At Bearing 3 (non-output end), absolute magnitude of the $x$-direction bearing load decreases from 16.82 kN to 2.24 kN, a reduction of 86.68 %. Meanwhile, the $x$-direction load at Bearing 4 (power-output end) changes direction, from –7.79 kN to 6.03 kN. Similar to the driving shaft, the $y$-direction load remains unaffected by acceleration variations.

In the $z$-direction, the absolute load at Bearing 3 decreases from 3.63 kN to 2.98 kN, a reduction of 17.91 %, and at Bearing 4, the absolute $z$-direction load reduces from 3.12 kN to 2.47 kN, a reduction of 20.83 %. The evolution of radial resultant load shows a 41.87 % reduction at Bearing 3 from 20.23 kN to 11.76 kN, with the load orientation angle in Quadrant I, increasing from 33.13° to 78.11°. At Bearing 4, the radial load initially decreases from 13.82 kN (–50 g) to 11.23 kN (10 g), then rises again to 12.63 kN (50 g). Correspondingly, the load angle first rotates from 54.01° (–50 g) to 88.02° (10 g), then further shifts to 61.50° (50 g), consistently remaining in Quadrant I.

In conclusion, bearing loads of both driving and driven shafts undergo significant changes under maneuvering overload level flight conditions. Therefore, the influence of additional inertial loads should not be neglected in bearing selection and strength verification.

Maneuvering overload accelerations from –50 g to 50 g induce only minimal quantitative alterations in gear meshing force characteristics. As illustrated in Fig. 18, the RMS of the radial component exhibits a linear micro-reduction from 10.55 kN to 10.54 kN, representing a marginal decay of 0.015 kN and a 0.14 % decay rate. Radial peak-to-peak fluctuations remain constrained within 0.4-0.8 kN throughout the acceleration envelope.

Fig. 17Load characteristics of driven shaft bearings

a) Bearing 3 $x$-direction load

b) Bearing 3 $y$-direction load

c) Bearing 3 $z$-direction load

d) Bearing 3 radial load

e) Bearing 4 $x$-direction load

f) Bearing 4 $y$-direction load

g) Bearing 4 z-direction load

h) Bearing 4 radial load

Fig. 18Gear mesh force

a) Radial force

b) Axial force

Contrastingly, the axial component RMS maintains near-invariant stability, varying between 5.850 kN and 5.852 kN with maximum deviation below 0.002 kN. Axial peak-to-peak oscillations are confined to a 0.1-0.3 kN magnitude range, showing that the axial component is barely affected by acceleration.

This confirms that maneuvering accelerations induce only micro-scale radial mean load attenuation without substantial alteration to dynamic impact characteristics (peak-to-peak signatures) or axial load-bearing stability.

As shown in Fig. 19, quantitative analysis of the time varying gear Hertz contact stress under maneuvering accelerations (–50 g to 50 g) confirms a 0.02 % reduction in RMS values, decreasing from 715.72 MPa to 715.55 MPa, while the mean contact stress shows a 0.022 % decline , reducing from 715.49 MPa to 715.33 MPa. These minimal variations demonstrate that inertial loads negligibly alter the stress characteristics in aerospace gear contacts. Frequency-domain analysis reveals that the dominant excitation frequencies of Hertzian contact stress are the gear meshing frequency $f_m$ (575 Hz) and its harmonics. The fundamental meshing frequency amplitude $f_m$ (575 Hz) exhibits a linear reduction with variations in level-flight acceleration, decreasing from 17.26 MPa at –50 g to 17.15 MPa at 50 g.

Fig. 19Gear Hertz contact stress

a) Time domain

b) Frequency domain

The results indicate that the contact behavior at the mesh interface is strictly governed by the torque and speed. Under the rated operating conditions, the meshing force remains the primary excitation. Consequently, the incremental dynamic effects induced by maneuvering overloads are essentially overwhelmed by the main transmission load, leading to the observed stability in gear mesh forces and Hertzian contact stresses.

5. Experimental verification and discussion

While experimental validation under sustained translational maneuvering acceleration is the ideal approach, to perform an effective steady-state spectral analysis on high-frequency gear mesh vibrations, the transmission system must be maintained in a strictly constant overload state for a sufficient steady-state duration. Quantitatively, maintaining a continuous translational acceleration even for a few seconds would necessitate a long high-precision linear track and result in an extremely high terminal velocity, which exceeds the physical and safety constraints of conventional laboratory linear motion platforms. Furthermore, conventional high-g shock equipment, such as drop towers, can only generate transient accelerations lasting a few milliseconds, which is incapable of capturing the steady-state dynamic evolution required for this study.

To validate the proposed modeling methodology despite these physical constraints, a kinematic equivalence verification is conducted based on a generalized non-inertial framework. The fully flexible coupled multibody dynamics model is established upon a generalized arbitrary non-inertial reference frame. Mathematically and physically, pure translational overload represents a specific case where the base's rotational angular velocity is zero. The introduction of base rotation transitions the system into a rotating non-inertial environment. Consequently, validating the model under complex rotational conditions inherently substantiates the theoretical formulations.

As illustrated in Fig. 20(a), a large-scale constant acceleration rotating test bench is utilized to simulate maneuvering overloads. The test bench has a single-arm length of 7.8 m and is driven by a YB2 series flameproof three-phase asynchronous motor. By adjusting the arm’s rotational speed to 33.8 r/min, 47.9 r/min, and 58.6 r/min, maneuvering accelerations of 10 g, 20 g, and 30 g are generated in the $x$-direction of the gear transmission system under test. The tested gear transmission system primarily consists of a driving motor, a gearbox, and a magnetic powder brake, as shown in Fig. 20(b). The experimental system utilizes a spur gear pair, characterized by a modulus of 2 mm and a pressure angle of 20°. The number of teeth for the driving and driven gears are 17 and 42, with corresponding tooth widths of 12 mm and 10 mm, respectively. During the test, the input speed is maintained at 1500 r/min, and the load torque is adjusted to 1.85 N·m by regulating the excitation current of the magnetic powder brake. Notably, although the specific physical parameters of the experimental setup differ from those in the numerical simulation, the setup is structurally analogous to the theoretical model, sharing an identical fundamental configuration and underlying dynamic mechanisms.

Fig. 20Test equipment

a) Acceleration rotating test equipment

b) Test gearbox

The vibration acceleration of the system was measured using a PCB-352C03 type integrated circuit piezoelectric (ICP) accelerometer, with the sampling frequency configured at 10240 Hz. The time-domain signals and amplitude spectra under various maneuvering acceleration conditions were compared with the results derived from the theoretical model.

For different magnitudes of maneuvering overloads, the vibration acceleration responses and spectral characteristics in the $x$-direction at measurement point 1 (located on the bearing housing at the input end) are illustrated in Figs. 21-24. The time-domain signals exhibit distinct periodic oscillation patterns, where the characteristic period of 0.04s corresponds to the 25 Hz shaft rotation frequency. Spectral analysis reveals that the dominant frequency of the acceleration response predicted by the theoretical model is the 4th harmonic of the mesh frequency (1700 Hz), followed by the sum of the 3rd mesh harmonic and the shaft rotation frequency, and the difference between the 6th mesh harmonic and twice the shaft rotation frequency. The distribution of characteristic frequencies in the experimental signals remains consistent with the theoretical analysis. The phenomenon where the amplitudes of the 8th and 9th harmonics in the experimental results exceed the theoretical predictions is primarily attributed to the simplifications in the model, such as bearing support stiffness and system damping. Nevertheless, the evolutionary trends of the overall spectral characteristics demonstrate consistency.

Comparative analysis of the amplitude spectra under various maneuvering acceleration conditions reveals that the amplitudes of the acceleration frequency components predicted by the theoretical model increase progressively with the maneuvering acceleration. Simultaneously, the amplitude of the dominant frequency component captured in the experimental signals also exhibits a near-linear increase as the maneuvering acceleration rises. These results demonstrate that the proposed full-flexible coupled dynamics modeling method can effectively capture the dynamic evolution trends of the gear system under maneuvering overloads. Since rotational maneuvers and translational accelerations can be described within the framework of the generalized non-inertial reference frame, the validation of the dynamic characteristics under rotational conditions further substantiates the reliability of the modeling methodology for predicting the system's dynamic responses under aircraft overload level-flight conditions.

Fig. 21Vibration acceleration under 0 g maneuvering acceleration

a) Theoretical time domain

b) Theoretical frequency domain

c) Experimental time domain

d) Experimental frequency domain

Fig. 22Vibration acceleration under 10 g maneuvering acceleration

a) Theoretical time domain

b) Theoretical frequency domain

c) Experimental time domain

d) Experimental frequency domain

Fig. 23Vibration acceleration under 20 g maneuvering acceleration

a) Theoretical time domain

b) Theoretical frequency domain

c) Experimental time domain

d) Experimental frequency domain

Fig. 24Vibration acceleration under 30 g maneuvering acceleration

a) Theoretical time domain

b) Theoretical frequency domain

c) Experimental time domain

d) Experimental frequency domain

6. Conclusions

This study focuses on clarifying the dynamic response mechanism of aircraft gear transmission systems under overload level-flight conditions, based on a fully flexible coupled multibody dynamics model that integrates flexible shafts, deformable gear pairs, flexible bearings, and deformable casings. By incorporating additional inertial excitations induced by overload maneuvers and conducting systematic analysis, the regulatory effects of maneuvering acceleration on the system’s internal excitation and load characteristics are elucidated, with key conclusions as follows:

1) Maneuvering overload induces a linear variation in vibration response by causing shafts deformation. This shaft deformation further intensifies the dynamic fluctuation of center distance offset, which in turn leads to changes in time-varying mesh stiffness, transmission error, and tooth backlash. This cascade effect clarifies that the fluctuation of center distance offset acts as the immediate driver for the evolution of the gear system’s core internal excitation parameters, while shaft deformation serves as the critical precursor linking maneuvering overload to center distance offset variation.

2) Maneuvering overload alters the dynamic force characteristics of critical components, primarily manifesting in changes to bearing support forces and force angles. For the shaft-bearing system, the redistribution of inertial loads driven by maneuvering acceleration directly modifies the magnitude and direction of bearing support forces, which is consistent with the variation law of shafts deformation and center distance offset.

3) Within the –50 g to +50 g maneuvering acceleration range, the root mean square (RMS) and mean values of gear meshing force and Hertzian contact stress exhibit slight linear attenuation under maneuvering overload, while their peak-to-peak values show no significant fluctuation. This indicates that maneuvering overload has a negligible impact on the dynamic impact characteristics of gear meshing contact, and the stability of gear meshing behavior is maintained despite changes in internal excitation and bearing force.

4) Experimental validation on a rotating test bench substantiates the model’s accuracy via kinematic equivalence. The amplitude of the 4th mesh harmonic in the vibration acceleration response increases progressively with maneuvering load, aligning with the theoretical trend and confirming the methodology’s reliability for predicting gear dynamics under overload level-flight conditions.

This study provides a theoretical framework and quantitative benchmark for the dynamic reliability design of transmission systems in high-maneuverability aircraft, with findings that can guide the structural optimization and performance assessment of aircraft gear transmissions under extreme maneuvering conditions.