Nonlinear analysis of a geared rotor system supported by fluid film journal bearings

doi:10.1016/j.jsv.2020.115269

Journal of Sound and Vibration
声音与振动》杂志

Volume 475, 9 June 2020, 115269
第 475 卷，2020 年 6 月 9 日，115269

https://doi.org/10.1016/j.jsv.2020.115269 Get rights and content 获取权利和内容

Abstract 摘要

This paper presents a novel approach for modeling and analyzing a geared rotor-bearing system including nonlinear forces in the gear set and the supporting fluid film journal bearings. The rotordynamics system model has five degrees of freedom that define the transverse displacements of the shaft-gear centerlines and the relative displacement of the gear tooth contact point. The journal bearing nonlinear forces are obtained via a solution of Reynolds equation for lubricant film pressure utilizing the finite element method. Co-existing, steady-state, autonomous and non-autonomous responses are obtained in an accurate and computationally efficient manner utilizing the multiple shooting and continuation algorithms. This yields the full manifolds of the multiple bifurcation system. Chaos is identified with maximum Lyapunov exponents, frequency spectra, Poincaré attractors, etc. The results reveal a dependence of the gear set contact conditions and system nonlinear response characteristics, i.e. jump, co-existing responses, subharmonic resonances and chaos on the choice of journal bearing parameters. The results also show that Hopf bifurcations, which occur along with oil whirl in a journal bearing system, can be attenuated by increasing the gear torque.
本文提出了一种新方法，用于模拟和分析齿轮转子轴承系统，包括齿轮组和支撑流体薄膜轴颈轴承中的非线性力。转子动力学系统模型有五个自由度，分别定义了轴-齿轮中心线的横向位移和齿轮齿接触点的相对位移。轴颈轴承的非线性力是通过利用有限元法求解润滑油膜压力的雷诺方程得到的。利用多重拍摄和延续算法，以精确且计算效率高的方式获得共存、稳态、自主和非自主响应。这产生了多重分岔系统的完整流形。混沌可通过最大 Lyapunov 指数、频率谱、Poincaré 吸引子等进行识别。结果表明，齿轮组接触条件和系统非线性响应特征，即跳跃、共存响应、亚谐共振和混沌，取决于轴颈轴承参数的选择。结果还表明，在轴颈轴承系统中伴随油旋发生的霍普夫分岔可以通过增加齿轮扭矩来减弱。

Keywords 关键词

Nonlinear gear dynamics

Multiple shooting method

Jump phenomena

Multiple co-existing response

Chaos

Journal bearing

非线性齿轮动力学多重射击法跳跃现象多重共存响应混沌期刊轴承

1. Introduction 1.导言

Gearing system speeds and operating torques continue to increase in high-performance machinery. This amplifies the effects of nonlinearities in the gears including tooth backlash and time-varying mesh stiffness. Backlash describes the intentional clearance provided between mating teeth to prevent binding and to include a thin lubricant film between the teeth for heat removal and reduced wear. Backlash causes intermittent loss of contact between the teeth creating a nonlinear force and torque. The mesh stiffness varies periodically with time due to the variation of the number of tooth pairs in contact, and the variation of the point of contact along with the tooth profiles. The time-varying stiffness of the meshing teeth may lead to parametric resonances, which are principal sources of internal excitations and vibrations in gear transmission systems. The backlash forces and time-varying stiffness interact yielding a complex nonlinear, parametrically excited system with both torsional and lateral vibration. Accurate and computationally efficient gear dynamic models, including nonlinear forces and parametric excitations, are required for the effective design of gear sets and the machinery in which they form a critical component.
在高性能机械中，齿轮系统的速度和工作扭矩不断提高。这放大了齿轮中的非线性效应，包括齿隙和随时间变化的啮合刚度。齿隙指的是啮合齿之间故意留出的间隙，以防止咬合并在齿间形成一层薄薄的润滑油膜，从而达到散热和减少磨损的目的。齿隙会导致齿间间歇性失去接触，从而产生非线性力和扭矩。由于接触齿对数量的变化以及接触点随齿廓的变化，啮合刚度随时间周期性变化。啮合齿的时变刚度可能导致参数共振，而参数共振是齿轮传动系统内部激励和振动的主要来源。反向间隙力和时变刚度相互作用，产生了一个复杂的非线性参数激励系统，并伴有扭转和横向振动。要有效设计齿轮组及其关键部件机械，就必须建立精确且计算效率高的齿轮动态模型，包括非线性力和参数激励。

Significant prior research has been performed on the nonlinear dynamic response of geared systems. Kahraman and Singh [1] analyzed the effect of backlash on a single-degree-of-freedom gear model employing both analytical and numerical simulations. They validated their model by comparison with experimental results and found that the nonlinear characteristics caused chaotic and subsynchronous resonance responses. Kahraman and Singh [2] examined interactions between gear backlash nonlinearity and the bearing clearances and identified chaotic and subharmonic responses. Kahraman and Singh [3] included time-varying stiffness and clearance nonlinearity in their numerical model of geared systems and identified strong coupling effects between these characteristics. Blankenship and Kahraman [4] presented an experimental – analytical correlation study of a geared system including backlash nonlinearity and parametric excitation. Their predictions of co-existing solutions with the harmonic balance method were confirmed experimentally. Kahraman and Blankenship [5] observed subharmonic resonances in a geared system experiment, which were demonstrated to be strongly dependent on damping ratio and stiffness variation of the gear mesh. Kahraman and Blankenship [6] experimentally observed chaotic vibration, jump phenomena, and subharmonic response due to parametric and backlash excitations. Ranghothama and Narayanan [7] employed an incremental harmonic balance method, arc-length continuation, Floquet theory, and Lyapunov exponents to examine the bifurcation characteristics of a three-degree-of-freedom geared rotor-bearing model. Theodossiades and Natsiavas [8] introduced a new analytical method for a gear system with time-varying stiffness and backlash using perturbations techniques. Al-shyyab and Kahraman [9,10] investigated the nonlinear response of a multi-mesh gear system using a multi-term harmonic balance method. The effects of gear parameters on the nonlinear behavior were studied for both period-one and sub-harmonic motions. Liu et al. [11] analyzed the effect of gear mesh damping and backlash amplitude on the states of gear meshing and nonlinear behaviors of a gear pair. Yang et al. [12,13] predicted the nonlinear vibration of a gear system subjected to multi-frequency excitations utilizing a multiple time scales method. They confirmed the interaction between different harmonic excitations and the complex nonlinear behaviors caused by the multi-frequency excitations. Yang et al. [14] performed parametric studies to investigate the influence of the contact ratio, spacing error, transmitted load and mesh damping of a gear using a fifth-order Runge-Kutta method. Wang et al. [15] analyzed the effect of modulation internal excitation on the gear system and verified the accuracy of the prediction by comparing its results with the experimental measurements.
之前已经对齿轮系统的非线性动态响应进行了大量研究。Kahraman 和 Singh [1] 通过分析和数值模拟，分析了反向间隙对单自由度齿轮模型的影响。他们通过与实验结果的对比验证了自己的模型，并发现非线性特性会导致混乱和次同步共振响应。Kahraman 和 Singh [2] 研究了齿轮反向间隙非线性与轴承间隙之间的相互作用，确定了混乱和次谐波响应。Kahraman 和 Singh [3] 在他们的齿轮系统数值模型中加入了时变刚度和间隙非线性，并确定了这些特性之间的强耦合效应。Blankenship 和 Kahraman [4] 对包括反向间隙非线性和参数激励的齿轮系统进行了实验-分析关联研究。他们用谐波平衡法预测的共存解决方案在实验中得到了证实。Kahraman 和 Blankenship [5] 在齿轮系统实验中观察到了次谐波共振，并证明其与齿轮啮合的阻尼比和刚度变化密切相关。Kahraman 和 Blankenship [6] 通过实验观察到了参数和反向间隙激励引起的混乱振动、跳跃现象和次谐波响应。Ranghothama 和 Narayanan [7] 采用增量谐波平衡法、弧长延续、Floquet 理论和 Lyapunov 指数研究了三自由度齿轮转子轴承模型的分岔特性。Theodossiades 和 Natsiavas [8] 利用扰动技术为具有时变刚度和反向间隙的齿轮系统引入了一种新的分析方法。Al-shyyab 和 Kahraman [9,10] 采用多期谐波平衡法研究了多啮合齿轮系统的非线性响应。研究了齿轮参数对周期一运动和次谐波运动的非线性行为的影响。Liu 等人[11] 分析了齿轮啮合阻尼和齿隙振幅对齿轮啮合状态和齿轮副非线性行为的影响。Yang 等人[12,13]利用多时间尺度方法预测了齿轮系统在多频激励下的非线性振动。他们证实了不同谐波激励之间的相互作用以及多频激励引起的复杂非线性行为。Yang 等人[14]采用五阶 Runge-Kutta 方法对齿轮的接触比、间距误差、传递载荷和啮合阻尼的影响进行了参数研究。Wang 等人[15 [15] 分析了调制内激振对齿轮系统的影响，并将其结果与实验测量结果进行比较，验证了预测的准确性。

Nonlinear vibration in different types of gears has been investigated. Motahar et al. [16] performed a numerical, nonlinear dynamics study of a bevel gear system. Tip and root modifications were introduced to study their influence on gear vibration. Yang and Lim [17] developed a hypoid gear model considering time varying mesh stiffness, backlash nonlinearity and time-varying bearing stiffness. They showed that the backlash nonlinearity could suppress parametric instability induced by the time-varying bearing stiffness, under certain operating conditions. Wang and Lim [18] studied the effect of gear mesh stiffness asymmetry for the drive and coast sides of the hypoid gear system, and confirmed that the mesh stiffness at the drive side has more significant effect on the nonlinear dynamics. Ambarisha and Parker [19] investigated nonlinear dynamics of a planetary gear system. They applied the profile of the time-varying mesh stiffness obtained from a finite element analysis to improve accuracy. Zhao and Ji [20] performed numerical simulations of a wind turbine gearbox having two planetary gear trains. Complex nonlinear responses of the gearbox were shown to result from a time-varying mesh stiffness, backlash nonlinearity and static transmission error. Xinghui et al. [21] analyzed parametric resonance of a planetary gear subjected to speed fluctuations. The gear model considers time-varying mesh stiffness, and the instability boundaries for the fundamental and combinations resonances were derived based on a perturbation analysis.
对不同类型齿轮的非线性振动进行了研究。Motahar 等人[16] 对锥齿轮系统进行了非线性动力学数值研究。他们引入了齿尖和齿根修正，研究它们对齿轮振动的影响。Yang 和 Lim [17] 建立了一个准双曲面齿轮模型，考虑了时变啮合刚度、反隙非线性和时变轴承刚度。他们的研究表明，在一定的工作条件下，反齿隙非线性可以抑制时变轴承刚度引起的参数不稳定性。Wang 和 Lim [18] 研究了准双曲面齿轮系统驱动侧和沿岸侧齿轮啮合刚度不对称的影响，证实驱动侧的啮合刚度对非线性动力学的影响更大。Ambarisha 和 Parker [19] 研究了行星齿轮系统的非线性动力学。他们应用从有限元分析中获得的时变啮合刚度轮廓来提高精度。Zhao 和 Ji [20] 对具有两个行星齿轮系的风力涡轮机齿轮箱进行了数值模拟。结果表明，齿轮箱的复杂非线性响应是由时变啮合刚度、反向间隙非线性和静态传动误差造成的。Xinghui 等人[21] 分析了速度波动下行星齿轮的参数共振。该齿轮模型考虑了时变啮合刚度，并根据扰动分析得出了基本共振和组合共振的不稳定边界。

Some researchers have explored approaches to suppress vibrations induced by gear nonlinearities. Cheon's [22] simulation study investigated the effect of a one-way clutch to reduce the dynamic transmission error of a geared system. Cheon [23] employed a phasing approach to reduce time-varying mesh stiffness and the resulting vibration, especially at the fundamental resonance.
一些研究人员探索了抑制齿轮非线性引起的振动的方法。Cheon [22] 的模拟研究调查了单向离合器对减少齿轮系统动态传输误差的影响。Cheon[23]采用了一种相位方法来减少时变啮合刚度和由此产生的振动，尤其是在基本共振时。

Stochastic methods have been applied to study the effects of uncertainty in gear parameters. Bonori and Pellicano [24] utilized a stochastic model to analyze the effect of manufacturing error on nonlinear gear dynamics and showed that this could induce chaotic vibrations in the gear system. Wei et al. [25] included modeling uncertainties of a gear system, such as mesh stiffness and damping, and determined the resulting response levels using an interval harmonic balance method.
随机方法已被用于研究齿轮参数不确定性的影响。Bonori 和 Pellicano [24] 利用随机模型分析了制造误差对非线性齿轮动力学的影响，结果表明这会引起齿轮系统的混沌振动。Wei 等人[25] 将啮合刚度和阻尼等齿轮系统的建模不确定性包括在内，并使用区间谐波平衡法确定了由此产生的响应水平。

Various analytical and modeling methods have been applied to gear dynamics simulation. Kim et al. [26] investigated the effect of smoothing functions on clearance nonlinearity of an oscillator and showed how the adjustment of a regulating factor associated with the smoothing functions yielded more reliable predictions. Farshidianfar and Saghafi [27] applied a Melnikov type analysis to investigate homoclinic bifurcations and chaotic responses in a geared system. Gou et al. [28] employed a cell mapping theory to analyze the multi-parameter coupling characteristics of gear parameters. Li et al. [29] used an incremental harmonic balance method to analyze gear systems with internal and external periodic excitations.
各种分析和建模方法已被应用于齿轮动力学模拟。Kim 等人[26]研究了平滑函数对振荡器间隙非线性的影响，并展示了调整与平滑函数相关的调节因子如何产生更可靠的预测结果。Farshidianfar 和 Saghafi [27] 应用梅尔尼科夫分析法研究了齿轮系统中的同轴分岔和混沌响应。Gou 等人[28] 采用单元映射理论分析了齿轮参数的多参数耦合特性。Li 等人[29] 采用增量谐波平衡法分析了具有内部和外部周期性激励的齿轮系统。

Hydrodynamic journal bearings are widely employed in geared systems with high speed and load requirements due to their relatively high stiffness and damping. Theodossiades and Natsiavas [30] investigated the effect of gear and journal bearing parameters on bifurcation, chaos and oil whirl. They represented the journal bearing force with a finite-length impedance method. Baguet and Jacquenot [31] developed a finite element shaft model to study the interactions between a helical gear and a finite-length bearing and showed that a linearized bearing coefficient model does not provide accurate predictions of gear vibrations, especially at high speed and load conditions. Fargère and Velex [32] investigated the effects of the bearing oil inlet location and thermal response on the gear system dynamics. These effects change the journal static equilibrium position, which in turn alters the dynamic response of the system. Liu et al. [33] studied the interactions between tooth wedging effect and journal bearing clearance using the approximate short journal bearing theory. Simulation results showed that varying the operating speed or applied torque may cause the occurrence of oil whirl response of the rotordynamic systems. The effect of tooth wedging on the vibration level of the geared-rotor system is also presented.
流体动力轴颈轴承因其相对较高的刚度和阻尼，被广泛应用于具有高速和高负载要求的齿轮系统中。Theodossiades 和 Natsiavas [30] 研究了齿轮和轴颈轴承参数对分叉、混沌和油旋的影响。他们用有限长度阻抗法表示了轴颈轴承力。Baguet 和 Jacquenot [31] 建立了一个有限元轴模型来研究斜齿轮和有限长度轴承之间的相互作用，结果表明线性化轴承系数模型不能准确预测齿轮振动，尤其是在高速和负载条件下。Fargère 和 Velex [32] 研究了轴承进油位置和热反应对齿轮系统动力学的影响。这些影响改变了轴颈的静态平衡位置，进而改变了系统的动态响应。Liu 等人[33]利用近似短轴颈轴承理论研究了齿楔效应和轴颈轴承间隙之间的相互作用。仿真结果表明，改变运行速度或施加扭矩可能会导致旋转动力系统出现油旋响应。此外，还介绍了齿楔对齿轮-转子系统振动水平的影响。

Kim and Palazzolo [34,35] employed shooting with deflation to study the nonlinear response of a Jeffcott rotor supported by floating ring bearings. The effects of changing parameters such as bearing length-to-diameter ( $L/D$ ) ratio and including the thermal effect of the lubricant were presented. Kim and Palazzolo [36] studied the bifurcation of a heavily loaded rotor with five-pad tilting pad bearings. A shooting/arc-length continuation approach was utilized to obtain quasi-periodic and chaotic motions, the latter being confirmed by maximum Lyapunov exponents.
Kim 和 Palazzolo [34,35]采用放气射击法研究了由浮动环轴承支撑的 Jeffcott 转子的非线性响应。研究介绍了改变轴承长径比 ( $L/D$ ) 等参数以及润滑剂热效应的影响。Kim 和 Palazzolo [36] 研究了装有五片倾斜垫轴承的重载转子的分叉问题。利用射击/弧长延续方法获得了准周期和混沌运动，后者由最大 Lyapunov 指数证实。

Prior models for coupled gearset-bearing vibration generally utilized lower fidelity or steady-state bearing models, and presented results in less rigorous nonlinear dynamics formats. This may have been motivated by the high computational expense of employing higher fidelity bearing models and presenting results in advanced nonlinear dynamics formats. Bearing forces were typically represented with linear spring and damping constants, or were obtained using short bearing theory, with highly simplified oil film cavitation models. The simplified approaches may lead to significant prediction error especially for steady-state responses with orbits that are relatively large (>15%) with respect to the bearing clearance. The present approach provides a highly accurate, finite element-based solution of the finite-length, Reynold's equation accounting for cavitation at each time step in the numerical integration. Additionally, results are presented in advanced nonlinear dynamics formats including bifurcation diagrams, maximum Lyapunov exponent plots, and Poincaré attractor plots. Computation time is held within practical limits utilizing multiple shooting and continuation algorithms, and with the use of embedded C++ components in the MATLAB code, and parallel processing.
之前的齿轮组-轴承耦合振动模型一般采用较低保真度或稳态轴承模型，并以不太严格的非线性动力学格式呈现结果。这可能是由于采用保真度较高的轴承模型并以高级非线性动力学格式呈现结果的计算成本较高。轴承力通常使用线性弹簧和阻尼常数表示，或使用短轴承理论和高度简化的油膜气蚀模型获得。简化方法可能会导致显著的预测误差，尤其是对于相对于轴承游隙而言轨道相对较大（>15%）的稳态响应。本方法提供了基于有限元的高精度有限长度雷诺方程解决方案，在数值积分的每个时间步长上都考虑了气蚀。此外，计算结果以先进的非线性动力学形式呈现，包括分岔图、最大 Lyapunov 指数图和 Poincaré 吸引子图。利用多重射击和延续算法，以及在 MATLAB 代码中使用嵌入式 C++ 组件和并行处理，将计算时间控制在实际范围内。

The highlight and original contribution of the work is to provide a computationally efficient, high fidelity and rigorously presented modeling approach for the dynamics of the five-degree-of-freedom dual shaft-gear pair system supported on fluid film bearings. This approach involves finite-length bearing models, advanced multiple shooting and continuation methods, gear flexibility and transmission error effects, bifurcation and Poincaré attractor diagrams, and maximum Lyapunov exponents for identifying chaotic behavior. Finally, this approach is applied to parametric studies with varying journal bearing and gear mesh stiffness parameters.
这项工作的亮点和原创性贡献在于为流体薄膜轴承支撑的五自由度双轴齿轮对系统的动力学提供了一种计算效率高、保真度高且严谨的建模方法。该方法涉及有限长度轴承模型、先进的多重射击和延续方法、齿轮柔性和传动误差效应、分岔和泊恩卡莱吸引子图，以及用于识别混沌行为的最大 Lyapunov 指数。最后，该方法被应用于不同轴颈轴承和齿轮啮合刚度参数的参数研究。

2. Modeling of a geared rotor system supported by fluid film journal bearings
2.由流体薄膜轴颈轴承支撑的齿轮转子系统建模

2.1. Five-degree-of-freedom gear-bearing-rotor model
2.1.五自由度齿轮-轴承-转子模型

Fig. 1 shows a centered gear pair attached to parallel rotors that are each supported by fluid film journal bearings. The model is composed of two rigid rotors having mass elements $m_{i}$ , radii $R_{i}$ and polar moments of inertia $J_{i}$ . The subscript $i$ denotes the driving ( $i = 1$ ) and driven ( $i = 2$ ) geared rotors.
图 1 显示了连接到平行转子上的对中齿轮，每个转子都由流体薄膜轴颈轴承支撑。该模型由两个刚性转子组成，其质量元素为 $m_{i}$ ，半径为 $R_{i}$ ，极惯性矩为 $J_{i}$ 。下标 $i$ 表示驱动转子 ( $i = 1$ ) 和从动转子 ( $i = 2$ ) 。

An external torque $T_{1}$ is applied to the driving gear. A nonlinear mesh coupling consisting of tooth backlash and time-varying stiffness is modeled to transmit torque between driving and driven gears. The motion coordinates for the model include ( $θ_{1}$ , $θ_{2}$ , $x_{1}$ , $x_{2}$ , $y_{1}$ , $y_{2}$ ) as shown in Fig. 2.
外部扭矩 $T_{1}$ 作用于驱动齿轮。由齿隙和时变刚度组成的非线性啮合耦合被建模用于在驱动齿轮和从动齿轮之间传递扭矩。模型的运动坐标包括 ( $θ_{1}$ , $θ_{2}$ , $x_{1}$ , $x_{2}$ , $y_{1}$ , $y_{2}$ ) ，如图 2 所示。

The dynamic transmission error (DTE), $δ (t)$ is given by(1) $δ (t) = R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)$ where $e_{r} (t)$ represents the static transmission error. The analytical description of the time-varying mesh stiffness and the static transmission error can be expressed in the form of Fourier series as [3](2) $\begin{array}{l} k_{m} (t) = k_{0} + \sum_{i = 1}^{\infty} s_{i} k_{0} \cos (i ω_{g} t - φ_{i}) \\ e_{r} (t) = e_{0} + \sum_{j = 1}^{\infty} p_{j} e_{0} \cos (j ω_{g} t - ψ_{j}) \end{array}$ where $k_{0}$ is a mean mesh stiffness, $e_{0}$ is a mean static transmission error, and $s_{i}$ and $p_{j}$ are the amplitude of the Fourier series components. The phase angles of the Fourier series are represented by $φ_{i}$ and $ψ_{j}$ , respectively.
$δ (t)$ 的动态传输误差 (DTE) 由 (1) $δ (t) = R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)$ 得出，其中 $e_{r} (t)$ 代表静态传输误差。对时变网格刚度和静态传输误差的分析描述可以用傅里叶级数的形式表示 [3] (2) $\begin{array}{l} k_{m} (t) = k_{0} + \sum_{i = 1}^{\infty} s_{i} k_{0} \cos (i ω_{g} t - φ_{i}) \\ e_{r} (t) = e_{0} + \sum_{j = 1}^{\infty} p_{j} e_{0} \cos (j ω_{g} t - ψ_{j}) \end{array}$ ，其中 $k_{0}$ 是平均网格刚度， $e_{0}$ 是平均静态传输误差， $s_{i}$ 和 $p_{j}$ 是傅里叶级数分量的振幅。傅里叶级数的相位角分别用 $φ_{i}$ 和 $ψ_{j}$ 表示。

The term $ω_{g}$ is the gear mesh frequency represented by(3) $ω_{g} = N_{i} ω_{i}$ where $N_{i}$ is the number of gear teeth, $ω_{i}$ is rotor operating frequency and $i = 1,2$ . For this study $ω_{1} = ω_{2} = ω$ is used, which follows the convention in the related literature [33]. The pressure angle $α$ is assumed to remain constant during operation. Plain journal bearings support both rigid shafts, and their nonlinear fluid film force models are explained in section 2.2.
$ω_{g}$ 是齿轮啮合频率，由 (3) $ω_{g} = N_{i} ω_{i}$ 表示，其中 $N_{i}$ 是齿轮齿数， $ω_{i}$ 是转子工作频率， $i = 1,2$ 。本研究采用 $ω_{1} = ω_{2} = ω$ ，这与相关文献[33]中的惯例一致。压力角 $α$ 假设在运行过程中保持恒定。滑动轴承支撑两个刚性轴，其非线性流体膜力模型将在 2.2 节中解释。

As noted in Ref. [26], a tooth backlash model defined with a piecewise linear function in the governing nonlinear differential equations may result in convergence difficulties when employing a Newton-Raphson method. Therefore, the following smoothening function presented in the same reference is also used in the present study.(4) $ρ (t) = \frac{1}{2} {(δ (t) - b) [1 + tanh (σ (δ (t) - b))]} + \frac{1}{2} {(δ (t) + b) [1 + tanh (- σ (δ (t) + b))]}$ where $ρ (t)$ represents relative gear mesh displacement considering backlash, $b$ is the half-length of the tooth backlash amplitude ( $b = \frac{b_{0}}{2}$ ) and $σ$ is a modulating factor which affects the accuracy of the backlash representation and convergence [26]. The value $σ = 100$ is selected for this study.
正如参考文献[26]所指出的，在使用牛顿-拉夫逊方法时，如果在非线性微分方程中使用片断线性函数来定义齿隙模型，可能会导致收敛困难。[26] 所述，在非线性微分方程中使用片断线性函数定义的齿隙模型在使用牛顿-拉夫逊方法时可能会导致收敛困难。因此，本研究也采用了同一参考文献中的以下平滑函数。 (4) $ρ (t) = \frac{1}{2} {(δ (t) - b) [1 + tanh (σ (δ (t) - b))]} + \frac{1}{2} {(δ (t) + b) [1 + tanh (- σ (δ (t) + b))]}$ ，其中 $ρ (t)$ 表示考虑了齿隙的齿轮啮合相对位移， $b$ 是齿隙振幅的半长 ( $b = \frac{b_{0}}{2}$ )， $σ$ 是影响齿隙表示精度和收敛性的调节因子 [26]。本研究选择的值为 $σ = 100$ 。

The coupling force between the driving and driven gear mesh is given by(5) $F_{m0} = k_{m} (t) δ (t) + c_{m} \dot{δ} (t)$ where $c_{m}$ represents mesh damping, and it is assumed to be constant in this study. $δ (t)$ represents the dynamic transmission error in Eq. (1), and $\dot{δ} (t)$ is its derivative.
驱动和从动齿轮啮合之间的耦合力由 (5) $F_{m0} = k_{m} (t) δ (t) + c_{m} \dot{δ} (t)$ 给出，其中 $c_{m}$ 代表啮合阻尼，在本研究中假定为常数。 $δ (t)$ 代表公式 (1) 中的动态传输误差， $\dot{δ} (t)$ 是其导数。

The equations of motion for the six-degree-of-freedom gear-bearing rotor system are(6) $\begin{array}{l} J_{1} {\ddot{θ}}_{1} + k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) R_{1} + c_{m} \dot{δ} (t) R_{1} = T_{1} \\ J_{2} {\ddot{θ}}_{2} - k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) R_{2} - c_{m} \dot{δ} (t) R_{2} = - T_{2} \\ m_{1} {\ddot{x}}_{1} + k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) \sin (α) + c_{m} \dot{δ} (t) \sin (α) = F_{b 1 x} \\ m_{1} {\ddot{y}}_{1} + k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) \cos (α) + c_{m} \dot{δ} (t) \cos (α) = F_{b 1 y} - m_{1} g \\ m_{2} {\ddot{x}}_{2} - k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) \sin (α) - c_{m} \dot{δ} (t) \sin (α) = F_{b 2 x} \\ m_{2} {\ddot{y}}_{2} - k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) \cos (α) - c_{m} \dot{δ} (t) \cos (α) = F_{b 2 y} - m_{2} g \end{array}$ By replacing the term $k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r}) + c_{m} \dot{δ} (t)$ with $F_{m0}$ , the equations become(7) $\begin{array}{l} J_{1} {\ddot{θ}}_{1} + R_{1} F_{m0} = T_{1} \\ J_{2} {\ddot{θ}}_{2} - R_{2} F_{m0} = - T_{2} \\ m_{1} {\ddot{x}}_{1} + F_{m0} \sin (α) = F_{b 1 x} \\ m_{1} {\ddot{y}}_{1} + F_{m0} \cos (α) = F_{b 1 y} - m_{1} g \\ m_{2} {\ddot{x}}_{2} - F_{m0} \sin (α) = F_{b 2 x} \\ m_{2} {\ddot{y}}_{2} - F_{m0} \cos (α) = F_{b 2 y} - m_{2} g \end{array}$ where $F_{bix}$ and $F_{biy}$ represents the $i$ th bearing forces in the $x$ and $y$ directions, and $m_{1} g$ and $m_{2} g$ terms represent gravity forces. By multiplying each of equations with $R_{1}$ and $R_{2}$ , the first two become(8) $\begin{array}{l} J_{1} R_{1} {\ddot{θ}}_{1} + {R_{1}}^{2} F_{m 0} = R_{1} T_{1} \\ J_{2} R_{1} {\ddot{θ}}_{2} - {R_{2}}^{2} F_{m 0} = - R_{1} T_{2} \end{array}$
六自由度齿轮轴承转子系统的运动方程为 (6) $\begin{array}{l} J_{1} {\ddot{θ}}_{1} + k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) R_{1} + c_{m} \dot{δ} (t) R_{1} = T_{1} \\ J_{2} {\ddot{θ}}_{2} - k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) R_{2} - c_{m} \dot{δ} (t) R_{2} = - T_{2} \\ m_{1} {\ddot{x}}_{1} + k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) \sin (α) + c_{m} \dot{δ} (t) \sin (α) = F_{b 1 x} \\ m_{1} {\ddot{y}}_{1} + k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) \cos (α) + c_{m} \dot{δ} (t) \cos (α) = F_{b 1 y} - m_{1} g \\ m_{2} {\ddot{x}}_{2} - k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) \sin (α) - c_{m} \dot{δ} (t) \sin (α) = F_{b 2 x} \\ m_{2} {\ddot{y}}_{2} - k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r} (t)) \cos (α) - c_{m} \dot{δ} (t) \cos (α) = F_{b 2 y} - m_{2} g \end{array}$ 将 $k_{m} (t) (R_{1} θ_{1} - R_{2} θ_{2} + x_{1} \sin (α) - x_{2} \sin (α) + y_{1} \cos (α) - y_{2} \cos (α) - e_{r}) + c_{m} \dot{δ} (t)$ 项替换为 $F_{m0}$ 后，方程变为 (7) $\begin{array}{l} J_{1} {\ddot{θ}}_{1} + R_{1} F_{m0} = T_{1} \\ J_{2} {\ddot{θ}}_{2} - R_{2} F_{m0} = - T_{2} \\ m_{1} {\ddot{x}}_{1} + F_{m0} \sin (α) = F_{b 1 x} \\ m_{1} {\ddot{y}}_{1} + F_{m0} \cos (α) = F_{b 1 y} - m_{1} g \\ m_{2} {\ddot{x}}_{2} - F_{m0} \sin (α) = F_{b 2 x} \\ m_{2} {\ddot{y}}_{2} - F_{m0} \cos (α) = F_{b 2 y} - m_{2} g \end{array}$ ，其中 $F_{bix}$ 和 $F_{biy}$ 代表 $x$ 和 $y$ 方向上的 $i$ th 轴承力， $m_{1} g$ 和 $m_{2} g$ 项代表重力。将每个方程与 $R_{1}$ 和 $R_{2}$ 相乘，前两个方程变为 (8) $\begin{array}{l} J_{1} R_{1} {\ddot{θ}}_{1} + {R_{1}}^{2} F_{m 0} = R_{1} T_{1} \\ J_{2} R_{1} {\ddot{θ}}_{2} - {R_{2}}^{2} F_{m 0} = - R_{1} T_{2} \end{array}$

Divide the two equation with $J_{1}$ and $J_{2}$ respectively, and then subtracting the second equation from the first one, to obtain(9) $R_{1} {\ddot{θ}}_{1} - R_{1} {\ddot{θ}}_{2} + (\frac{{R_{1}}^{2}}{J_{1}} + \frac{{R_{2}}^{2}}{J_{2}}) F_{m0} = \frac{R_{1}}{J_{1}} T_{1} + \frac{R_{2}}{J_{2}} T_{2}$ Substituting $p = R_{1} θ_{1} - R_{2} θ_{2}$ and manipulating the equation yields(10) $\ddot{p} + (\frac{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}}{J_{1} J_{2}}) F_{m0} = \frac{R_{1}}{J_{1}} T_{1} + \frac{R_{2}}{J_{2}} T_{2}$ Dividing through by $(\frac{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}}{J_{1} J_{2}})$ yields(11) $(\frac{J_{1} J_{2}}{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}}) \ddot{p} + F_{m0} = (\frac{J_{1} J_{2}}{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}}) (\frac{R_{1}}{J_{1}} T_{1} + \frac{R_{2}}{J_{2}} T_{2})$ Substitute $J_{e}$ for $(\frac{J_{1} J_{2}}{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}})$ to obtain(12) $J_{e} \ddot{p} + F_{m0} = T_{e}$ where $T_{e} = J_{e} (\frac{R_{1}}{J_{1}} T_{1} + \frac{R_{2}}{J_{2}} T_{2})$ is an equivalent input torque term.
将两个等式分别与 $J_{1}$ 和 $J_{2}$ 相除，然后将第二个等式减去第一个等式，得出 (9) $R_{1} {\ddot{θ}}_{1} - R_{1} {\ddot{θ}}_{2} + (\frac{{R_{1}}^{2}}{J_{1}} + \frac{{R_{2}}^{2}}{J_{2}}) F_{m0} = \frac{R_{1}}{J_{1}} T_{1} + \frac{R_{2}}{J_{2}} T_{2}$ 将 $p = R_{1} θ_{1} - R_{2} θ_{2}$ 代入并运算等式，得出 (10) $\ddot{p} + (\frac{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}}{J_{1} J_{2}}) F_{m0} = \frac{R_{1}}{J_{1}} T_{1} + \frac{R_{2}}{J_{2}} T_{2}$ 将 $(\frac{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}}{J_{1} J_{2}})$ 除以 (11) $(\frac{J_{1} J_{2}}{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}}) \ddot{p} + F_{m0} = (\frac{J_{1} J_{2}}{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}}) (\frac{R_{1}}{J_{1}} T_{1} + \frac{R_{2}}{J_{2}} T_{2})$ 将 $J_{e}$ 代入 $(\frac{J_{1} J_{2}}{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}})$ ，得出 (12) $J_{e} \ddot{p} + F_{m0} = T_{e}$ $T_{e} = J_{e} (\frac{R_{1}}{J_{1}} T_{1} + \frac{R_{2}}{J_{2}} T_{2})$ ，其中是等效输入扭矩项。

The term $δ (t)$ was inserted into Eq. (4) to include the backlash nonlinearity effect. Then, from Eqs. (4), (5), the gear meshing force including the backlash nonlinearity effect becomes(13) $F_{m} = k_{m} ρ (t) + c_{m} \dot{δ} (t)$
在公式 (4) 中插入了 $δ (t)$ 一词，以包含反向间隙非线性效应。然后，根据公式 (4)、(5)，包括反向间隙非线性效应的齿轮啮合力变为 (13) $F_{m} = k_{m} ρ (t) + c_{m} \dot{δ} (t)$

Finally, the equations including the backlash nonlinearity, time-varying mesh stiffness and the static transmission error become(14) $\begin{array}{l} J_{e} \ddot{p} + F_{m} = T_{e} \\ m_{1} {\ddot{x}}_{1} + F_{m} \sin (α) = F_{b 1 x} \\ m_{1} {\ddot{y}}_{1} + F_{m} \cos (α) = F_{b 1 y} - m_{1} g \\ m_{2} {\ddot{x}}_{2} - F_{m} \sin (α) = F_{b 2 x} \\ m_{2} {\ddot{y}}_{2} - F_{m} \cos (α) = F_{b 2 y} - m_{2} g \end{array}$ where J_e is the equivalent inertia of two gears, i.e. $(\frac{J_{1} J_{2}}{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}})$ . The torsional natural frequency $ω_{n}$ of the system is defined as $\sqrt{\frac{k_{0}}{J_{e}}}$ .
最后，包括齿隙非线性、时变啮合刚度和静态传动误差的方程变为 (14) $\begin{array}{l} J_{e} \ddot{p} + F_{m} = T_{e} \\ m_{1} {\ddot{x}}_{1} + F_{m} \sin (α) = F_{b 1 x} \\ m_{1} {\ddot{y}}_{1} + F_{m} \cos (α) = F_{b 1 y} - m_{1} g \\ m_{2} {\ddot{x}}_{2} - F_{m} \sin (α) = F_{b 2 x} \\ m_{2} {\ddot{y}}_{2} - F_{m} \cos (α) = F_{b 2 y} - m_{2} g \end{array}$ ，其中 J _e 是两个齿轮的等效惯性，即 $(\frac{J_{1} J_{2}}{J_{2} {R_{1}}^{2} + J_{1} {R_{2}}^{2}})$ 。系统的扭转固有频率 $ω_{n}$ 定义为 $\sqrt{\frac{k_{0}}{J_{e}}}$ 。

For validation purposes, the simulation result are compared with experimental measurements [4] in Fig. 3. The experiment was conducted using relatively stiff ball bearings so the $x$ and $y$ journal motions are assumed fixed in the simulation. The gear parameters for backlash, Fourier coefficients of time varying mesh stiffness and amplitude of static transmission error from Ref. [4] are employed in the simulation. The root mean square (RMS) value of the dynamic transmission error is plotted with respect to operating speed, showing good agreement between prediction and test results.
出于验证目的，模拟结果与实验测量结果[4]进行了比较，如图 3 所示。实验是使用相对较硬的球轴承进行的，因此在模拟中假定 $x$ 和 $y$ 轴颈运动是固定的。参考文献[4]中的反向间隙、啮合刚度时变傅里叶系数和静态传输误差振幅等齿轮参数被用于仿真。[4] 中的反向间隙傅里叶系数和静态传输误差振幅。动态传动误差的均方根值随运行速度而变化，显示出预测与测试结果之间的良好一致性。

Table 1 provides a second validation case through comparison of the five-degree-of-freedom gear-bearing dynamic model's predicted natural frequencies with those provided in Ref. [15]. Since natural frequencies are characteristics of a linear model, the backlash and time varying stiffness were omitted and the bearing forces were represented by the stiffness and damping provided in the reference. The correlation is shown in the table and confirms excellent agreement.
表 1 提供了第二个验证案例，将五自由度齿轮轴承动态模型的预测自然频率与参考文献[15]中提供的自然频率进行了比较。[15].由于自然频率是线性模型的特征，因此省略了反向间隙和时变刚度，轴承力由参考文献中提供的刚度和阻尼表示。表中显示了相关性，并证实两者非常吻合。

Table 1. Comparison of calculated natural frequencies with [15].
表 1.计算的自然频率与 [15] 的比较。

Empty Cell	Calculated natural frequencies in Ref. [15]	Natural frequencies based on current model
1st mode 第一模式	0 Hz 0 赫兹	0 Hz 0 赫兹
2nd mode 第二模式	1149 Hz 1149 赫兹	1149.3 Hz 1149.3 赫兹
3rd mode 第三模式	1293 Hz 1293 赫兹	1293.6 Hz 1293.6 赫兹
4th mode 第四模式	1604 Hz 1604 赫兹	1604.3 Hz 1604.3 赫兹
5th mode 第 5 模式	1799 Hz 1799 赫兹	1799.1 Hz 1799.1 赫兹
6th mode 第六模式	5043 Hz 5043 赫兹	5043.7 Hz 5043.7 赫兹

2.2. Finite element model of plain journal bearing
2.2.滑动轴承的有限元模型

The Reynolds equation [34] for an incompressible lubricant combines the fluid continuity and momentum equations into a partial differential equation for film pressure, and is given by(15) $\frac{\partial}{\partial θ} (\frac{h^{3}}{12 μ} \frac{\partial p}{\partial θ}) + \frac{\partial}{\partial z} (\frac{h^{3}}{12 μ} \frac{\partial p}{\partial z}) = \frac{R_{J} ω_{J}}{2} \frac{\partial h}{\partial θ} + \frac{\partial h}{\partial t}$ where $ω_{J}$ is the rotating speed of the journal, and $R_{J}$ and $μ$ represent the radius of the journal and the viscosity of the lubricant, respectively. The centers of the bearing and the journal are $O_{B}$ and $O_{J}$ in Fig. 4, respectively.
不可压缩润滑油的雷诺方程 [34] 将流体连续性方程和动量方程合并为油膜压力偏微分方程，其式为 (15) $\frac{\partial}{\partial θ} (\frac{h^{3}}{12 μ} \frac{\partial p}{\partial θ}) + \frac{\partial}{\partial z} (\frac{h^{3}}{12 μ} \frac{\partial p}{\partial z}) = \frac{R_{J} ω_{J}}{2} \frac{\partial h}{\partial θ} + \frac{\partial h}{\partial t}$ ，其中 $ω_{J}$ 为轴颈的旋转速度， $R_{J}$ 和 $μ$ 分别代表轴颈半径和润滑油粘度。轴承和轴颈的中心分别为图 4 中的 $O_{B}$ 和 $O_{J}$ 。

The displacements of the journal center relative to the bearing center in the $x$ and $y$ directions are $x_{J}$ and $y_{J}$ , respectively, and $p$ is the pressure in the lubricant film. Expressions for fluid film thickness h and its derivative $\frac{\partial h (θ)}{\partial t}$ at $θ$ are given by(16) $h (θ) = C_{B} - x_{J} \cos θ - y_{J} \sin θ \frac{\partial h (θ)}{\partial t} = - {\dot{x}}_{J} \cos θ - {\dot{y}}_{J} \sin θ$ where $C_{B}$ represents the bearing radial clearance.
轴颈中心相对于轴承中心在 $x$ 和 $y$ 方向上的位移分别为 $x_{J}$ 和 $y_{J}$ ， $p$ 为润滑油膜中的压力。液膜厚度 h 及其在 $θ$ 处的导数 $\frac{\partial h (θ)}{\partial t}$ 的表达式为 (16) $h (θ) = C_{B} - x_{J} \cos θ - y_{J} \sin θ \frac{\partial h (θ)}{\partial t} = - {\dot{x}}_{J} \cos θ - {\dot{y}}_{J} \sin θ$ ，其中 $C_{B}$ 代表轴承径向游隙。

The mathematical model assumes rigid shafts and rigid attachments between the bearings and ground. Therefore, the journal motions $x_{J}$ and $y_{J}$ are identical to their respective gear centerline motions. Thus $x_{1}$ , $y_{1}$ are identical to $x_{J 1}$ and $y_{J 1}$ , and $x_{2}$ , $y_{2}$ are identical to $x_{J 2}$ and $y_{J 2}$ .
数学模型假设轴是刚性的，轴承和地面之间的连接也是刚性的。因此， $x_{J}$ 和 $y_{J}$ 的轴颈运动与其各自的齿轮中心线运动相同。因此， $x_{1}$ , $y_{1}$ 与 $x_{J 1}$ 和 $y_{J 1}$ 相同， $x_{2}$ , $y_{2}$ 与 $x_{J 2}$ 和 $y_{J 2}$ 相同。

The finite element mesh of a fluid film is illustrated in Fig. 5. The coordinate $θ$ corresponds to the circumferential direction of the film and the direction of rotation is from the left ( $θ_{B}$ ) to the right ( $θ_{E}$ ). The axial coordinate is represented with $z$ and only a half-length ( $\frac{L}{2}$ ) of the film is modeled due to its symmetry. The pressure on the bottom ( $z = 0$ ) side of the mesh are set to ambient pressure $P_{ambient}$ . Continuous pressure and flow condition are imposed on the left and right sides of the mesh. The zero-flow condition at the symmetric side ( $z = \frac{L}{2}$ ) and the continuous pressure and pressure gradient conditions at left/right sides are applied as follows(17) $v_{z = L/ 2} = 0, P_{θ = π} = P_{θ = - π}, \frac{\partial p}{\partial θ_{θ = π}} = \frac{\partial p}{\partial θ_{θ = - π}}$
流体薄膜的有限元网格如图 5 所示。坐标 $θ$ 对应薄膜的圆周方向，旋转方向是从左 ( $θ_{B}$ ) 到右 ( $θ_{E}$ )。轴向坐标用 $z$ 表示，由于薄膜的对称性，只对半长 ( $\frac{L}{2}$ ) 的薄膜进行建模。网格底部 ( $z = 0$ ) 的压力设置为环境压力 $P_{ambient}$ 。在网格的左右两侧施加了连续的压力和流动条件。对称侧 ( $z = \frac{L}{2}$ ) 的零流动条件和左右两侧的连续压力和压力梯度条件如下所示 (17) $v_{z = L/ 2} = 0, P_{θ = π} = P_{θ = - π}, \frac{\partial p}{\partial θ_{θ = π}} = \frac{\partial p}{\partial θ_{θ = - π}}$

The Reynolds equation is solved with a mesh of triangular simplex finite elements, which interpolate the two-dimensional pressure distribution in the film domain. The instantaneous reaction force on a journal is obtained by integrating the pressure distribution. Considering that symmetry condition the journal reaction force becomes(18) $F_{bi} = {\begin{cases} F_{bix} \\ F_{biy} \end{cases}} = 2 \int_{0}^{L / 2} \int_{- π}^{π} p {\begin{cases} \cos θ \\ \sin θ \end{cases}} d θ d z$
雷诺方程是通过三角简支有限元网格求解的，该网格对膜域中的二维压力分布进行了插值。通过对压力分布进行积分，可获得轴颈上的瞬时反作用力。考虑到对称条件，轴颈反作用力变为 (18) $F_{bi} = {\begin{cases} F_{bix} \\ F_{biy} \end{cases}} = 2 \int_{0}^{L / 2} \int_{- π}^{π} p {\begin{cases} \cos θ \\ \sin θ \end{cases}} d θ d z$

3. Nonlinear steady-state solution methods
3.非线性稳态求解方法

3.1. Multiple shooting method
3.1.多重拍摄法

The shooting method (SM) is a numerical procedure that utilizes an iterative algorithm and numerical integration of the nonlinear differential equations to locate co-existing, periodic equilibrium states. The SM provides a guided iterative search to locate the state vectors that repeat after a specified, or unknown in the case of autonomous systems, period. The single shooting method (SSM) is widely used in nonlinear dynamics research because of its simplicity. However, SSM may experience convergence problems, especially at saddle-node points. Multiple shooting methods (MSM) improves the numerical stability of the SSM by dividing time intervals into smaller ones. Compared to SSM, the MSM shows more robust convergence to periodic states and is less sensitive to the selection of initial state guesses. In addition, MSM is more suitable for parallel computing, thus making it desirable for systems with a large number of degrees of freedom.
射击法（SM）是一种数值程序，它利用非线性微分方程的迭代算法和数值积分来定位共存的周期性平衡状态。单射法提供了一种引导性迭代搜索，以找到在指定周期后重复出现的状态矢量，或在自主系统中的未知周期后重复出现的状态矢量。单次射击法（SSM）因其简便性被广泛应用于非线性动力学研究。然而，SSM 可能会遇到收敛问题，尤其是在鞍节点处。多重射击法（MSM）通过将时间间隔划分为更小的时间间隔，提高了 SSM 的数值稳定性。与 SSM 相比，MSM 能更稳健地收敛到周期状态，对初始状态猜测的选择也不那么敏感。此外，MSM 更适合并行计算，因此适用于具有大量自由度的系统。

The MSM algorithm is explained in this section. The non-autonomous nonlinear equations of motion can be represented by the first order form as(19) $\frac{d}{d τ} x = x^{'} = h (x, τ, p) (n \times 1)$ where $x$ is a state vector, $τ$ is an explicit time variable in the forcing term and $p$ represents the physical parameters of a system. The period of the steady-state harmonic response defined by a user is represented by the minimum period $P_{F}$ and a rational number $R$ in the non-autonomous case(20) $P_{R} = R P_{F}$
本节将解释 MSM 算法。非自主非线性运动方程可以用一阶形式表示为 (19) $\frac{d}{d τ} x = x^{'} = h (x, τ, p) (n \times 1)$ ，其中 $x$ 是状态向量， $τ$ 是强迫项中的显式时间变量， $p$ 表示系统的物理参数。用户定义的稳态谐波响应周期由最小周期 $P_{F}$ 和非自治情况下的有理数 $R$ 表示 (20) $P_{R} = R P_{F}$

The solution at the end of the period $P_{R}$ is represented(21) $x_{T_{R}} = x (τ = P_{R}, x (0) = x_{0})$ where $x_{0}$ is the initial condition state vector. If $x_{0}$ is a solution on an orbital equilibrium state of the period $P_{R}$ , it will result that(22) $x_{T_{R}} = x (τ = P_{R}, x (0) = x_{0}) = x_{0} x_{0} = g + e$ where $g$ is a user-defined guess of initial conditions, and $e$ is an error term.
$P_{R}$ 周期结束时的解表示为 (21) $x_{T_{R}} = x (τ = P_{R}, x (0) = x_{0})$ ，其中 $x_{0}$ 是初始条件状态向量。如果 $x_{0}$ 是周期为 $P_{R}$ 的轨道平衡状态的解，则结果为 (22) $x_{T_{R}} = x (τ = P_{R}, x (0) = x_{0}) = x_{0} x_{0} = g + e$ ，其中 $g$ 是用户自定义的初始条件猜测， $e$ 是误差项。

Unlike SSM which requires the single end-point constraint ( $x_{T_{R}} = x_{0}$ ) as explained above, MSM divides $P_{R}$ into smaller intervals and generates multiple constraints as follows:(23) $\begin{array}{l} x_{P_{1}} = x (τ = P_{1}, x (0) = x_{0}) = x_{1} \\ x_{P_{2}} = x (τ = P_{2}, x (0) = x_{1}) = x_{2} \\ ⋮ \\ x_{P_{m - 1}} = x (τ = P_{m - 1}, x (0) = x_{m - 2}) = x_{m - 1} \\ x_{P_{m}} = x (τ = P_{m}, x (0) = x_{m - 1}) = x_{m} = x_{0} \end{array}$ where divided time intervals are $0 = P_{0} < P_{1} < \dots < P_{m} = P_{R}$ . Note that $m$ is the number of time intervals defined by a user. A multi-dimensional Newton-Raphson method is applied to update $x_{0}$ .
如上所述，SSM 需要单个端点约束 ( $x_{T_{R}} = x_{0}$ )，而 MSM 则不同，它将 $P_{R}$ 划分为更小的时间间隔，并生成多个约束，如下所示： (23) $\begin{array}{l} x_{P_{1}} = x (τ = P_{1}, x (0) = x_{0}) = x_{1} \\ x_{P_{2}} = x (τ = P_{2}, x (0) = x_{1}) = x_{2} \\ ⋮ \\ x_{P_{m - 1}} = x (τ = P_{m - 1}, x (0) = x_{m - 2}) = x_{m - 1} \\ x_{P_{m}} = x (τ = P_{m}, x (0) = x_{m - 1}) = x_{m} = x_{0} \end{array}$ ，其中划分的时间间隔为 $0 = P_{0} < P_{1} < \dots < P_{m} = P_{R}$ 。请注意， $m$ 是用户定义的时间间隔数。采用多维牛顿-拉夫逊方法更新 $x_{0}$ 。

For an autonomous system, the additional phase condition should be defined since the point of the periodic solution at a specific time is not unique. In this study, a phase condition that sets the DTE from Eq. (1) as zero is used [42].(24) $H (x_{0}, P_{A}) = g (x_{0}, P_{A}) - x_{0} = 0 c (x_{0}) = 0$ where $P_{A}$ is a period of an autonomous system orbit to be identified along with an initial condition $x_{0}$ .
对于自主系统，由于特定时间的周期解点并非唯一，因此应定义附加的相位条件。本研究采用的相位条件是将式 (1) 中的 DTE 设为零[42]。 (24) $H (x_{0}, P_{A}) = g (x_{0}, P_{A}) - x_{0} = 0 c (x_{0}) = 0$ ，其中 $P_{A}$ 是待识别的自主系统轨道周期，以及初始条件 $x_{0}$ 。

3.2. Arc-length continuation
3.2.弧长延续

The shooting method may take considerable computation time, especially when plotting co-existing solution loci versus system parameters. “Continuation algorithms” have been developed to generate the loci (branch plots) with significantly increased efficiency relative to conducting independent SM searches for each parameter value. The Arc-length Continuation (AC) method [42] is applied in this research. AC provides robust solution searches even in high curvature regions by utilizing the trajectory of the solution curves along an arc-length as shown in Fig. 6.
射击法可能会耗费大量计算时间，尤其是在绘制共存解位置与系统参数的关系图时。相较于对每个参数值进行独立的 SM 搜索，"连续算法 "的开发大大提高了生成位置（分支图）的效率。本研究采用了弧长延续（AC）方法[42]。如图 6 所示，AC 通过利用求解曲线沿弧长的轨迹，即使在高曲率区域也能提供稳健的求解搜索。

An additional unknown $ω$ is involved in the iteration search. The next solution is(25) ${\begin{cases} {\tilde{g}}_{n}^{i + 1} \\ ω_{n}^{i + 1} \end{cases}} = {\begin{cases} {\tilde{g}}_{n}^{i} \\ ω_{n}^{i} \end{cases}} + [\begin{matrix} {\tilde{J}}_{g}^{i} & J_{ω}^{i} \\ \frac{\partial k_{n}^{i^{T}}}{\partial x} & \frac{\partial k_{n}^{i}}{\partial ω} \end{matrix}] {\begin{cases} - f ({\tilde{g}}_{n}^{i}, ω_{n}^{i}) \\ - q ({\tilde{g}}_{n}^{i}, ω_{n}^{i}, s) \end{cases}}$ where $ω$ and $s$ are operating parameters and arc-length of the solution curve, respectively, $n$ is a current step number, $J_{ω}$ is Jacobian matrix with regard to $ω$ , and $k$ is the constraint imposed on the solution procedures as(26) $k (x, ω, s) = ν {‖ {\tilde{g}}_{n + 1}^{i} - {\tilde{g}}_{n} ‖}_{2}^{2} + {(ω_{n + 1}^{i} - ω_{n})}^{2} - {(Δ s)}^{2}$ where $ν$ is a relaxation factor, and $Δ s$ is an arc-length.
迭代搜索中还涉及到一个未知数 $ω$ 。下一个解决方案是 (25) ${\begin{cases} {\tilde{g}}_{n}^{i + 1} \\ ω_{n}^{i + 1} \end{cases}} = {\begin{cases} {\tilde{g}}_{n}^{i} \\ ω_{n}^{i} \end{cases}} + [\begin{matrix} {\tilde{J}}_{g}^{i} & J_{ω}^{i} \\ \frac{\partial k_{n}^{i^{T}}}{\partial x} & \frac{\partial k_{n}^{i}}{\partial ω} \end{matrix}] {\begin{cases} - f ({\tilde{g}}_{n}^{i}, ω_{n}^{i}) \\ - q ({\tilde{g}}_{n}^{i}, ω_{n}^{i}, s) \end{cases}}$ ，其中 $ω$ 和 $s$ 分别是运行参数和求解曲线的弧长， $n$ 是当前步数， $J_{ω}$ 是与 $ω$ 有关的雅各布矩阵， $k$ 是对 (26) $k (x, ω, s) = ν {‖ {\tilde{g}}_{n + 1}^{i} - {\tilde{g}}_{n} ‖}_{2}^{2} + {(ω_{n + 1}^{i} - ω_{n})}^{2} - {(Δ s)}^{2}$ 的求解程序施加的约束，其中 $ν$ 是松弛因子， $Δ s$ 是弧长。

Then Newton-Raphson iteration is performed until the convergence criteria are satisfied. A continuation of periodic solution searches is carried out to plot a frequency response in the excitation frequency range of interest, including bifurcation points. The phase condition in Eq. (24) is added to Eq. (25) for the autonomous system continuation algorithm.
然后进行牛顿-拉夫逊迭代，直到满足收敛标准。进行周期性解的延续搜索，以绘制相关激励频率范围内的频率响应，包括分岔点。式 (24) 中的相位条件被添加到式 (25) 中，用于自主系统延续算法。

3.3. Stability identification based on the shooting method
3.3.基于射击法的稳定性识别

The Jacobian matrix of the shooting method is calculated to determine the local stability of periodic solutions. More specifically, the eigenvalues of the Jacobian matrix at a steady-state solution identifies the solution's stability and its bifurcation type. Perturbed solutions are computed to generate the Jacobian matrix entries. The system is considered unstable if the maximum magnitude of the eigenvalues is larger than unity.
计算射影法的雅各布矩阵可确定周期解的局部稳定性。更具体地说，雅各布矩阵在稳态解处的特征值可确定解的稳定性及其分岔类型。通过计算扰动解来生成雅各布矩阵项。如果特征值的最大值大于 1，则系统被认为是不稳定的。

3.4. Lyapunov exponents for identifying chaos
3.4.识别混沌的李亚普诺夫指数

Various approaches are used to identify the presence of chaos in the response of a nonlinear dynamical system. The most widely used approach is to calculate Maximum Lyapunov exponent (MLE, $μ_{\max}$ ). Lyapunov exponents indicate the rate of separation of two infinitesimally close trajectories in the local phase space [42]. A total of $n$ initial separation vectors with different directions are used for a system with $n$ states, to obtain a spectrum of $n$ Lyapunov exponents (LE, $μ_{i} (i = 1,2, ..., n)$ ) for calculating the rates of separation. The simultaneous numerical integrations of nonlinear differential equations and linearized form of them are required for the MLE calculation.(27) $Nonlinear differential equations : \dot{x} = g (x)$ (28) $Linearized form of Eq. (27) : \dot{β} = L (β)$
有多种方法可用于识别非线性动力系统响应中是否存在混沌。最广泛使用的方法是计算最大李亚普诺夫指数（MLE， $μ_{\max}$ ）。Lyapunov 指数表示两个无限接近的轨迹在局部相空间中的分离率[42]。对于具有 $n$ 状态的系统，总共使用了 $n$ 不同方向的初始分离矢量，从而得到 $n$ 李亚普诺夫指数（LE, $μ_{i} (i = 1,2, ..., n)$ ）的频谱，用于计算分离率。计算 MLE 时需要同时对非线性微分方程及其线性化形式进行数值积分。 (27) $Nonlinear differential equations : \dot{x} = g (x)$ (28) $Linearized form of Eq. (27) : \dot{β} = L (β)$

The actual and linearized trajectories with perturbed initial conditions are calculated from Eq. (27) and Eq. (28), respectively, at various times along the nonlinear system trajectory. Deviation distances are obtained from the difference between the nonlinear trajectory and linearized trajectories.(29) $Δ (t) = \sqrt{δ x_{1}^{2} + δ x_{1}^{2} + \dots + δ x_{n}^{2}}$
根据公式 (27) 和公式 (28) 分别计算出非线性系统轨迹上不同时间的扰动初始条件下的实际轨迹和线性化轨迹。根据非线性轨迹和线性化轨迹之间的差值得出偏差距离。 (29) $Δ (t) = \sqrt{δ x_{1}^{2} + δ x_{1}^{2} + \dots + δ x_{n}^{2}}$

An appropriate time interval $t_{f}$ is selected for the numerical integrations to avoid a numerical error. Sets of orthonormalized perturbed vectors are obtained from a Gram-Schmidt procedure as(30) ${\tilde{β}}_{1} = \frac{β_{1} (t_{f})}{‖ β_{1} (t_{f}) ‖}, {\tilde{β}}_{2} = \frac{β_{2} (t_{f}) - (β_{2} (t_{f}) \cdot {\tilde{β}}_{1}) {\tilde{β}}_{1}}{‖ β_{2} (t_{f}) - (β_{2} (t_{f}) \cdot {\tilde{β}}_{1}) {\tilde{β}}_{1} ‖}, and {\tilde{β}}_{m} = \frac{β_{m} (t_{f}) - \sum_{i = 1}^{m - 1} (β_{m} (t_{f}) \cdot {\tilde{β}}_{i}) {\tilde{β}}_{i}}{‖ β_{m} (t_{f}) - \sum_{i = 1}^{m - 1} (β_{m} (t_{f}) \cdot {\tilde{β}}_{i}) {\tilde{β}}_{i} ‖}$
为避免数值误差，数值积分选择了适当的时间间隔 $t_{f}$ 。正交化扰动向量集通过格拉姆-施密特程序得到，即 (30) ${\tilde{β}}_{1} = \frac{β_{1} (t_{f})}{‖ β_{1} (t_{f}) ‖}, {\tilde{β}}_{2} = \frac{β_{2} (t_{f}) - (β_{2} (t_{f}) \cdot {\tilde{β}}_{1}) {\tilde{β}}_{1}}{‖ β_{2} (t_{f}) - (β_{2} (t_{f}) \cdot {\tilde{β}}_{1}) {\tilde{β}}_{1} ‖}, and {\tilde{β}}_{m} = \frac{β_{m} (t_{f}) - \sum_{i = 1}^{m - 1} (β_{m} (t_{f}) \cdot {\tilde{β}}_{i}) {\tilde{β}}_{i}}{‖ β_{m} (t_{f}) - \sum_{i = 1}^{m - 1} (β_{m} (t_{f}) \cdot {\tilde{β}}_{i}) {\tilde{β}}_{i} ‖}$

After conducting the integrations of Eqs. (28), (29), (30) for $r$ times, the Lyapunov exponents are obtained as(31) $μ_{i} = \frac{1}{r t_{f}} \sum_{k = 1}^{r} \ln (Δ_{i}^{k} (t_{k}))$ where $Δ_{i}^{k} (t_{k})$ is the denominator of the orthonormal vector $β_{i}^{k}$ , $k$ denotes the $k_{t h}$ time step and $i$ represents $i$ th vector element.
对 $r$ 次进行公式 (28)、(29)、(30) 的积分后，Lyapunov 指数为 (31) $μ_{i} = \frac{1}{r t_{f}} \sum_{k = 1}^{r} \ln (Δ_{i}^{k} (t_{k}))$ ，其中 $Δ_{i}^{k} (t_{k})$ 是正交向量 $β_{i}^{k}$ 的分母， $k$ 表示 $k_{t h}$ 的时间步长， $i$ 表示 $i$ 的第 1 个向量元素。

The MLE is used as a quantitative measure to determine the chaotic response of a nonlinear dynamical system as follows.(32) $\begin{array}{l} • μ_{max} > 0 : System is chaotic (necessary but not a sufficient condition for chaos) \\ • μ_{\max} < 0 : System attracts to a fixed point or a stable periodic orbit (asymptotically stable) \\ • μ_{\max} = 0 : The orbit is quasi - periodic \end{array}$
MLE 用作确定非线性动力系统混沌响应的定量指标，其计算方法如下。 (32) $\begin{array}{l} • μ_{max} > 0 : System is chaotic (necessary but not a sufficient condition for chaos) \\ • μ_{\max} < 0 : System attracts to a fixed point or a stable periodic orbit (asymptotically stable) \\ • μ_{\max} = 0 : The orbit is quasi - periodic \end{array}$

4. Numerical example 4.数字示例

The five-degree-of-freedom geared rotor-bearing model of Fig. 2 and Eq. (14) is utilized to demonstrate various nonlinear phenomena induced by gear and journal bearing nonlinearities. The multiple shooting/continuation method and Lyapunov exponents discussed in section 3 are utilized along with direct numerical integration. MATLAB ODE15s was used with a relative tolerance of $10^{- 9}$ for computing the Jacobian matrix in the shooting/continuation procedure. Embedded C++ coding and parallel processing are utilized in the MATLAB program to accelerate the execution. The results are divided into 3 sections (1) parametric resonances/jump phenomena and the effect of journal bearing parameter variations on those phenomena, and (2) chaotic responses due to gear nonlinearity and the effect of journal bearing parameters on the responses, and (3) the effect of gear mesh stiffness on the hydrodynamic stability of the gear system supported by journal bearings. Solid and dashed lines indicate stable and unstable responses, respectively in all figures. Table 2 summarizes the parameters of the spur gear pair and journal bearings used in this study.
图 2 中的五自由度齿轮转子轴承模型和公式 (14) 用于演示齿轮和轴颈轴承非线性引起的各种非线性现象。在直接数值积分的同时，还使用了第 3 节中讨论的多重射击/延续法和 Lyapunov 指数。使用 MATLAB ODE15s 计算射击/延续过程中的雅各布矩阵，相对容差为 $10^{- 9}$ 。在 MATLAB 程序中使用了嵌入式 C++ 编码和并行处理，以加快执行速度。结果分为三部分：（1）参数共振/跳动现象以及轴颈轴承参数变化对这些现象的影响；（2）齿轮非线性引起的混沌响应以及轴颈轴承参数对响应的影响；（3）齿轮啮合刚度对轴颈轴承支撑的齿轮系统流体动力学稳定性的影响。在所有图表中，实线和虚线分别表示稳定和不稳定响应。表 2 概述了本研究中使用的正齿轮副和轴颈轴承的参数。

Table 2. Parameters of the spur gear pair and journal bearings.
表 2.正齿轮副和轴颈轴承的参数。

Gear parameters 齿轮参数
Mass 质量	$m_{1} = m_{2} = 9.3276$ kg $m_{1} = m_{2} = 9.3276$ 公斤
Moment of inertia 惯性矩	$J_{1} = J_{2} = 0.03187$ $kg \cdot m^{2}$
Radius of gears 齿轮半径	$R_{1} = R_{2} = 0.0982$ m
Pressure angle 压力角	$α = 20 °$
Backlash amplitude 反向间隙振幅	$b_{0} = 100 μm$
Mean mesh (tooth) stiffness 平均网格（齿）刚度	$k_{0} = 1 e 8 N/m$
Mesh damping ratio 网格阻尼比	$ζ_{m} = 0.01 - 0.025$
Number of gear teeth 齿轮齿数	$N_{1} = N_{2} = 28$
Applied torque ( $T_{1}$ ) 应用扭矩 ( $T_{1}$ )	100–3000 N m 100-3000 牛米
Journal bearing parameters 滚动轴承参数
Bearing diameter 轴承直径	$D_{B} = 0.092$ $m$
Bearing clearance 轴承间隙	$C_{B} = 74 - 184$ μm
Bearing $L / D$ ratio 轴承 $L / D$ 比率	$0.3 - 2$
Lubricant viscosity 润滑油粘度	$μ = 10 - 90$ mPa s

4.1. Parametric resonances and jump phenomenon
4.1.参数共振和跳跃现象

This section presents results for parametric instability, including fundamental and subharmonic resonances induced by time varying mesh stiffness. Both the fundamental and subharmonic resonances are parametric resonances since they are removed when the time varying component of the gear mesh stiffness is removed. Parametric resonance and jump phenomena of a spur gear system with backlash nonlinearity and time-varying mesh stiffness were treated in Refs. [[1], [2], [3], [4], [5], [6], [7], [8]] and more recently in Refs. [[9], [10], [11], [12], [13], [14], [15]]. However, the effects of journal bearing parameters on the nonlinear response of the geared system needs further investigation. Parametric instability of the fundamental and double-period subharmonic frequencies is presented below including journal bearing effects. Applied torque input is considered as an excitation source, and time-varying mesh stiffness with the frequency corresponding to gear mesh frequency $ω_{g}$ in Eq. (3) is included. No imbalance excitation and static transmission error $e_{r} (t)$ in Eq. (2) are applied in this section. The non-autonomous, multiple shooting/continuation methods in section 3 are applied to the system equations to analyze the influence of journal bearing parameters on nonlinear responses. Only the first Fourier coefficient of the time-varying mesh stiffness $s_{1}$ in Eq. (2) was included, and it was set equal to 0.2, since values from 0.1 to 0.3 are employed in the literatures [[1], [2], [3], [4], [5], [6], [7], [8], [9], [10], [11], [12], [13], [14], [15]]. The bearing $L/D$ ratio is 1, the radial clearance is 105 μm, the lubricant viscosity is 30 mPa s, and an applied torque of 1250 N m is used for the nominal values of the simulation. Each parameter is varied to investigate its effect on the softening effect due to gear nonlinearities.
本节介绍了参数不稳定性的结果，包括由啮合刚度时变引起的基谐和次谐波共振。基谐波和次谐波共振都是参数共振，因为在消除齿轮啮合刚度的时变分量时，它们都会消失。文献[[1]、[2]、[3]、[4]、[5]、[6]、[7]、[8]]和最近的文献[[9]、[10]、[11]、[12]、[13]、[14]、[15]]讨论了具有反齿隙非线性和时变啮合刚度的直齿轮系统的参数共振和跳跃现象。然而，轴颈轴承参数对齿轮系统非线性响应的影响还需要进一步研究。下面介绍基频和双周期次谐波频率的参数不稳定性，包括轴颈轴承的影响。外加扭矩输入被视为激励源，并包含与公式（3）中齿轮啮合频率 $ω_{g}$ 相对应频率的时变啮合刚度。本节不使用不平衡激励和公式 (2) 中的静态传输误差 $e_{r} (t)$ 。第 3 节中的非自主、多重射击/延续方法应用于系统方程，以分析轴颈轴承参数对非线性响应的影响。公式 (2) 中只包含时变啮合刚度 $s_{1}$ 的第一个傅立叶系数，并将其设置为 0.2，因为文献[[1]、[2]、[3]、[4]、[5]、[6]、[7]、[8]、[9]、[10]、[11]、[12]、[13]、[14]、[15]]中采用的值为 0.1 至 0.3。轴承 $L/D$ 比率为 1，径向游隙为 105 μm，润滑剂粘度为 30 mPa s，模拟的名义值为 1250 N m 的外加扭矩。通过改变每个参数来研究其对齿轮非线性软化效应的影响。

Direct numerical integrations with initial conditions determined from the shooting/continuation procedure are performed to obtain peak-to-peak displacement amplitude vs. rotor operating speed plots. At the same time, the stability of the periodic solutions is presented and obtained from the shooting method's Jacobian matrix eigenvalues.
根据射击/延续程序确定的初始条件进行直接数值积分，以获得峰-峰位移振幅与转子运行速度的关系图。同时，还介绍了周期解的稳定性，并从射击法的雅各布矩阵特征值中获得了稳定性。

4.1.1. Fundamental resonance
4.1.1.基本共振

Fig. 7 shows the response amplitude regarding operating speed for the non-dimensional peak-to-peak displacements of DTE ( $\frac{δ}{b_{0}}$ ) where $δ$ is defined from Eq. (1). The applied torque is 1100 N m in this case. The linear system resonance will occur when the gear mesh frequency $ω_{g} = N ω = ω_{n} = \sqrt{\frac{k_{0}}{J_{e}}}$ . This occurs at around 2500 rpm in this case.
图 7 显示了 DTE ( $\frac{δ}{b_{0}}$ ) 的非尺寸峰峰位移在运行速度方面的响应振幅，其中 $δ$ 由公式 (1) 定义。本例中的应用扭矩为 1100 N m。当齿轮啮合频率为 $ω_{g} = N ω = ω_{n} = \sqrt{\frac{k_{0}}{J_{e}}}$ 时，线性系统将发生共振。在本例中，共振频率约为 2500 rpm。

Including backlash causes a softening effect related left-leaning backbone curve with an unstable branch occurring between 1800 and 2300 rpm in Fig. 7(a). Multiple co-existing solutions are seen to occur when the rotor speed is in the vicinity of $\frac{ω_{n}}{N}$ where $N$ represents the number of gear teeth. The presence of co-existing solutions is clearly seen to result from including backlash, and the peak DTE severity is seen to nearly double when backlash is included. Note that the time-varying stiffness effect is included in both the with and without backlash cases. The increased DTE caused by the gear nonlinearity is also confirmed by the time response of both cases at 2000 rpm in Fig. 7(b).
在图 7(a)中，包括反向间隙会导致与左倾主干曲线相关的软化效应，不稳定分支出现在 1800 至 2300 rpm 之间。当转子转速在 $\frac{ω_{n}}{N}$ 附近时，会出现多个并存解，其中 $N$ 代表齿轮齿数。可以清楚地看到，并存解的出现是由于加入了反向间隙，当加入反向间隙时，DTE 的峰值严重程度几乎增加了一倍。请注意，有齿隙和无齿隙情况下都包含了时变刚度效应。图 7(b) 中两种情况在 2000 rpm 时的时间响应也证实了齿轮非线性导致的 DTE 增加。

The results in Fig. 8 are obtained for three different applied torques $T_{1}$ , i.e., 1,000, 1250 and 1500 N m ( $T_{2} = 0$ ), where $T_{1}$ is for driving gear and $T_{2}$ is for driven gear. All three cases exhibit a fundamental parametric resonance caused by the time-varying mesh stiffness, with a softening effect starting from near 2300 rpm. The amplitudes of the forced steady-state harmonic response increase with increased torque inputs. In all cases of Fig. 8(a), the unstable responses emerge through saddle-node bifurcations around 2300 rpm. Three or five multiple, co-existing steady-state responses occur between 1700 and 2300 rpm depending on the applied torque values. The critical rpm (2300 rpm) for jump-up behavior is nearly independent of applied torque variation, but in contrast, the critical rpm for jump-down behavior is more sensitive to applied torque variation. The critical rpm of jump-down behavior tends to increase with increasing drive torque. The separation between the jump-up and jump-down rpm is reduced with increasing applied torque. The emergence of the right-leaning portion of the response curves, which observed at the highest applied torque 1500 N m may be explained in terms of the backlash forces. Three different meshing states can exist depending on the maximum DTE; no impact, single-sided impact, and double-sided impact. In Fig. 8(a), the system shows only a single-sided contact with the 1000 N m applied torque, which is the source of the primary softening effect. Increasing the applied input torque gives rise to a hardening effect, which introduces additional co-existing responses and jump-up/down frequencies. The peak-peak response at 1500 N m shows a clear hardening effect along with the softening effect. The number of multiple co-existing responses increased from three to five with the applied torque increased from 1000 to 1500 N m. The figure also shows a tendency of the response curve to move rightwards towards the zero backlash response as the DTE increases. This is consistent with a greater engagement of teeth as the amplitude increases. These results are consistent with the experimental and numerical results in Ref. [4], which used a single-degree-of-freedom gear model considering only torsional motion. Prior research utilized a simple rigid or linear stiffness and damping bearing model, or analytical short bearing theory, or finite-length impedance method, which precluded the accurate investigation of nonlinear bearing parameter effects on the response. In this study, utilizing the finite element method for the nonlinear fluid film force of the finite-length bearing, the bearing $L/D$ ratio is varied from 0.3 to 2 in Fig. 8(b). The peak-peak DTE displacements in the frequency range away from the resonance region are not significantly affected by the bearing $L/D$ ratio variation. However, similar to the applied torque input case, the jump-down frequency is influenced by the $L/D$ ratio variation, since the jump-down event occurs at a relatively lower frequency range with higher $L/D$ ratio. The jump down speed is lowered by about 100 rpm for $L/D$ = 2 compared with $L/D$ = 0.3. Fig. 9 presents the bifurcation diagrams corresponding to the input torques 1000 and 1500 N m in Fig. 8(a), in a manner that highlights the jump phenomena and multiple co-existing solutions.
图 8 中的结果是在三种不同的应用扭矩 $T_{1}$ 下得到的，即 1,000、1250 和 1500 N m ( $T_{2} = 0$ )，其中 $T_{1}$ 代表驱动齿轮， $T_{2}$ 代表从动齿轮。所有三种情况都表现出由时变啮合刚度引起的基本参数共振，从接近 2300 rpm 开始出现软化效应。强迫稳态谐波响应的振幅随着扭矩输入的增加而增大。在图 8(a) 的所有情况下，不稳定响应都是在 2300 rpm 附近通过鞍节点分岔出现的。根据所应用的扭矩值，在 1700 至 2300 rpm 之间会出现三个或五个多重并存的稳态响应。跃升行为的临界转速（2300 rpm）几乎不受外加扭矩变化的影响，但相比之下，跃降行为的临界转速对外加扭矩变化更为敏感。随着驱动扭矩的增加，跳降行为的临界转速也会增加。随着应用扭矩的增加，上跳和下跳转速之间的间隔也会缩小。在最高应用扭矩 1500 N m 时出现的响应曲线右倾部分可以用反冲力来解释。根据最大 DTE，可存在三种不同的啮合状态：无冲击、单侧冲击和双侧冲击。在图 8(a)中，系统在施加 1000 N m 扭矩时仅显示单面接触，这是主要软化效应的来源。增加所施加的输入扭矩会产生硬化效应，从而带来额外的共存响应和上跳/下跳频率。1500 N m 时的峰-峰响应显示出明显的硬化效应和软化效应。随着施加扭矩从 1000 牛米增加到 1500 牛米，多重共存响应的数量从三个增加到五个。这与振幅增大时齿的啮合增大是一致的。这些结果与参考文献 [4] 中的实验和数值结果一致。[4] 中的实验和数值结果一致，后者使用的是单自由度齿轮模型，只考虑了扭转运动。之前的研究使用了简单的刚性或线性刚度和阻尼轴承模型，或分析性短轴承理论，或有限长度阻抗法，这些方法无法精确研究非线性轴承参数对响应的影响。在本研究中，利用有限元法计算有限长度轴承的非线性流体膜力，在图 8(b) 中，轴承 $L/D$ 的比率在 0.3 到 2 之间变化。在远离共振区的频率范围内，DTE 的峰-峰位移受轴承 $L/D$ 比率变化的影响不大。然而，与外加扭矩输入情况类似，跳降频率受 $L/D$ 轴承比变化的影响，因为当 $L/D$ 轴承比越高时，跳降事件发生在相对较低的频率范围内。与 $L/D$ = 0.3 相比，当 $L/D$ = 2 时，跳降速度降低了约 100 rpm。图 9 显示了与图 8（a）中输入扭矩 1000 和 1500 N m 相对应的分岔图，突出显示了跳跃现象和多个并存解。

Fig. 10 illustrates the effect of varying bearing lubricant viscosity on steady-state, nonlinear harmonic response. Lubricant viscosities of 10 and 90 mPa s, in addition to the nominal value of 30 mPa s, are simulated for two $L/D$ ratio ( $L/D$ = 0.3 and 1). Increasing the viscosity decreases the jump-down frequencies and increases the peak resonant amplitude in both $L/D$ cases.
图 10 显示了轴承润滑油粘度变化对稳态非线性谐波响应的影响。除了 30 mPa s 的标称值外，还模拟了 10 和 90 mPa s 的润滑油粘度对两个 $L/D$ ratio ( $L/D$ = 0.3 和 1) 的影响。在两个 $L/D$ 的情况下，增加粘度都会降低下跃频率并增加共振振幅峰值。

Fig. 11 shows the frequency-amplitude diagram with radial bearing clearances of $C_{B}$ = 184 μm and $C_{B}$ = 74 μm, along with nominal $C_{B}$ = 105 μm. The results with the smallest clearance, i.e., $C_{B}$ = 74 μm, show a broader range of co-existing solutions region compared to other values in both the $L/D$ = 0.3 and the $L/D$ = 1 cases. The $C_{B}$ = 74 μm and $L/D$ = 1 case displays a double-sided impact and resulting in hardening effect in Fig. 11(b). The $L/D$ = 1 case shows more reduction in jump-down frequency as compared to the $L/D$ = 0.3 case. This reduction may be attributed to the fact that the $L/D$ = 1 geometry has more fluid film area than for $L/D$ = 0.3, and hence greater force to affect the parametric resonance.
图 11 显示了径向轴承间隙为 $C_{B}$ = 184 μm 和 $C_{B}$ = 74 μm 以及标称值 $C_{B}$ = 105 μm 时的频率-振幅图。在 $L/D$ = 0.3 和 $L/D$ = 1 的情况下，与其他值相比，最小间隙（即 $C_{B}$ = 74 μm）的结果显示出更大范围的共存解区域。图 11(b) 中， $C_{B}$ = 74 μm 和 $L/D$ = 1 的情况显示出双面影响，并导致硬化效应。与 $L/D$ = 0.3 的情况相比， $L/D$ = 1 的情况显示出更大的下跳频率减小。这种降低可能是由于 $L/D$ = 1 的几何形状比 $L/D$ = 0.3 的几何形状有更大的流体膜面积，因此影响参数共振的力更大。

Fig. 12 shows five co-existing steady-state responses in the frequency and phase plane domain for $C_{B}$ = 74 μm case in Fig. 11 (b). For the phase portrait, the rotor spin speed equals 1850 rpm, and it is obtained via the non-autonomous MSM developed in section 3. The excitations include the time-varying mesh stiffness and the applied torque. Periodic responses 1, 3 and 5 are stable and 2 and 4 are unstable as predicted by the eigenvalues of the Jacobian matrix of the MSM. Note that the unstable forced harmonic responses cannot be obtained by direct numerical integration, but only with a directed search approach such as the MSM.
图 12 显示了图 11 (b) 中 $C_{B}$ = 74 μm 情况下频率和相位平面域中的五个共存稳态响应。在相位图中，转子旋转速度等于 1850 rpm，它是通过第 3 节中开发的非自主 MSM 获得的。激励包括时变网格刚度和外加扭矩。根据 MSM 的雅各布矩阵特征值预测，周期响应 1、3 和 5 是稳定的，2 和 4 是不稳定的。请注意，不稳定的受迫谐波响应无法通过直接数值积分获得，只能通过 MSM 等定向搜索方法获得。

The mesh deformation $ρ (t)$ in Eq. (4) corresponding to the multiples responses in Fig. 12 are shown in Fig. 13. The mesh deformation in Fig. 13(a) has static and dynamic components and never loses contact. The response in Fig. 13(b) has zero mesh deformation for most of the period and demonstrates the effect of single-sided contact induced by the backlash nonlinearity. Fig. 13(c) shows double-sided contact cycling between positive, zero and negative mesh deformation states. The result of this contact behavior is a net hardening effect as evidenced by the right-leaning secondary bend in the response 5 curve of Fig. 12.
图 13 显示了公式 (4) 中与图 12 中的多重响应相对应的网格变形 $ρ (t)$ 。图 13(a)中的网格变形包含静态和动态部分，并且从未失去接触。图 13(b)中的响应在大部分时间内网格变形为零，显示了反向间隙非线性引起的单边接触效应。图 13(c) 显示了在正、零和负网格变形状态之间循环的双面接触。这种接触行为的结果是净硬化效应，图 12 的响应曲线 5 中右倾的二次弯曲就是证明。

Fig. 14 illustrates the repelling/attracting motions among stable and unstable forced harmonic response orbits at 1850 rpm. The unstable Response 2 is repelled towards the stable attractors, i.e. Response 1 (Fig. 14(a)) or Response 3 (Fig. 14(b)). Thus, the unstable manifold (the dotted line) initiated by a saddle-node bifurcation in Fig. 12 acts as a border manifold, which provides information about the convergence route in phase space.
图 14 显示了 1850 rpm 转速下稳定和不稳定强迫谐波响应轨道之间的排斥/吸引运动。不稳定的响应 2 被排斥向稳定的吸引子，即响应 1（图 14(a)）或响应 3（图 14(b)）。因此，图 12 中由鞍节点分岔引发的不稳定流形（虚线）作为边界流形，提供了有关相空间收敛路线的信息。

4.1.2. Subharmonic resonances
4.1.2.次谐波共振

Prior work [3,38,39] revealed that subharmonic resonance of geared systems is highly dependent on the amplitude of time-varying mesh stiffness, mesh damping ratio and static applied torque. The focus here is to demonstrate the effect of journal bearing parameters on the subharmonic vibrations occurring for a rotor speed in the vicinity of $\frac{2 ω_{n}}{N}$ . The maximum rpm is extended from 3100 rpm for the fundamental resonance case to 5500 rpm to observe the subharmonic resonances. Fig. 15(a) shows that increasing torque increases the peak-peak DTE amplitude, the backlash nonlinearity generally causes a softening (left-leaning resonance) effect, but a double-sided contact is evidenced by the right-leaning, hardening peak at 750 N m torque, similar with the fundamental resonance in Fig. 8.
先前的研究[3,38,39]表明，齿轮系统的次谐振与时变啮合刚度、啮合阻尼比和静态应用扭矩的振幅有很大关系。这里的重点是演示转子转速在 $\frac{2 ω_{n}}{N}$ 附近时，轴颈轴承参数对次谐振的影响。最大转速从基频共振情况下的 3100 rpm 扩展到 5500 rpm，以观察次谐波共振。图 15(a) 显示，增加扭矩会增加 DTE 的峰值振幅，反向间隙非线性通常会导致软化（左倾共振）效应，但在 750 N m 扭矩时出现的右倾硬化峰值证明了双面接触，这与图 8 中的基本共振类似。

Fig. 16 provides waterfall diagrams corresponding to the 750 N m torque case in Fig. 15(a) and the bearing $L/D$ = 0.6 case in Fig. 15(b). Fig. 16(a) and (b) correspond to run-up and run-down in speed for the 750 N m torque case with $L/D$ = 1. Fig. 16(c) and (d) correspond to run-up and run-down in speed for the bearing $L/D$ = 0.6 case with 625 N m torque input. The figures show jump-up and jump-down bifurcations near 2500 rpm, and the presence of a sub-harmonic resonance and a 0.5x gear mesh frequency response.
图 16 提供了与图 15(a) 中 750 N m 扭矩情况和图 15(b) 中轴承 $L/D$ = 0.6 情况相对应的瀑布图。图 16(a)和(b)对应于 $L/D$ = 1 时 750 N m 扭矩情况下的速度上升和下降。图 16(c)和(d)对应于输入 625 N m 扭矩时轴承 $L/D$ = 0.6 的速度上升和下降。图中显示了 2500 rpm 附近的上跳和下跳分岔，以及亚谐波共振和 0.5 倍齿轮啮合频率响应。

In Fig. 15, unstable branches appear as a period-doubling bifurcation emerges near 5000 rpm, which corresponds to twice the fundamental resonance frequency rotor speed ( $\frac{2 ω_{n}}{N}$ ). A saddle-node bifurcation occurs as the period-doubled unstable branches reach their peak amplitude, yielding a stable branch. As a result, one stable solution with the gear mesh frequency ( $ω_{g}$ ) and two unstable/stable solutions with half gear mesh frequency ( $0.5 ω_{g}$ ) coexist in the operating speed range between 3000 and 5000 rpm.
在图 15 中，当周期加倍分岔出现在 5000 rpm 附近时，不稳定分支就会出现，这相当于基频共振频率转子转速的两倍 ( $\frac{2 ω_{n}}{N}$ )。当周期加倍的不稳定分支达到峰值振幅时，会出现鞍节点分岔，产生一个稳定分支。因此，在 3000 至 5000 rpm 的工作转速范围内，一个齿轮啮合频率的稳定解 ( $ω_{g}$ ) 和两个齿轮啮合频率为一半的不稳定/稳定解 ( $0.5 ω_{g}$ ) 同时存在。

Variation of the journal bearing parameters significantly affects the subharmonic resonance as shown in Fig. 15, Fig. 17. The sub-harmonic resonance is seen to vanish when considering a drop in $L/D$ from 0.8 to 0.6. Physically this may result from the more significant force of the larger $L/D$ bearing inducing the subharmonic vibrations. Fig. 17(a) shows that increasing the lubricant viscosity extends the frequency overlap range by decreasing the jump-down frequency. Similar trends are also observed in the case of bearing clearance in Fig. 17(b). The subharmonic is seen to disappear for the low viscosity (10 mPa s) case and the low clearance (105 μm) case. Fig. 18 shows multiple co-existing DTE versus time plots for the 30 mPa s case in Fig. 17(a): 3575 rpm, $L/D$ = 1, an applied torque of 625 N m, a lubricant viscosity of 30 mPa s and a bearing clearance of 36.8 μm. The $t ω_{n}$ in Fig. 18 is the non-dimensional time variable where the time $t$ is multiplied by the torsional natural frequency $ω_{n}$ . Note that the Response 1 in Fig. 18 (a) has the frequency equal to $ω_{g}$ , while the Response 2 and 3 in Fig. 18(b) and (c) have the half frequency of $ω_{g}$ .
如图 15 和图 17 所示，轴颈轴承参数的变化对次谐波共振有很大影响。当 $L/D$ 从 0.8 降到 0.6 时，亚谐波共振消失。从物理角度看，这可能是由于较大的 $L/D$ 轴承产生了更大的力，从而诱发了次谐波振动。图 17(a)显示，润滑油粘度的增加降低了下跳频率，从而扩大了频率重叠范围。图 17（b）中轴承游隙的情况也有类似趋势。在低粘度（10 mPa s）和低游隙（105 μm）情况下，次谐波消失了。图 18 显示了图 17（a）中 30 mPa s 情况下多个并存的 DTE 与时间关系图：转速为 3575 rpm， $L/D$ = 1，应用扭矩为 625 N m，润滑油粘度为 30 mPa s，轴承游隙为 36.8 μm。图 18 中的 $t ω_{n}$ 是非量纲时间变量，其中时间 $t$ 乘以扭转固有频率 $ω_{n}$ 。请注意，图 18 (a) 中响应 1 的频率等于 $ω_{g}$ ，而图 18 (b) 和 (c) 中响应 2 和 3 的半频率为 $ω_{g}$ 。

4.2. Chaotic response 4.2.混沌响应

The effect of linearized bearing model stiffness and damping coefficients, static transmission error, tooth friction, etc. on the chaotic response of geared systems was investigated in Refs. [2,6,7,24,40,41], and chaotic behavior was observed experimentally in Ref. [6]. The effects of varying journal bearing parameters using a nonlinear bearing model, on the geared system's chaotic response are discussed below. The model parameter values include a low torque load of 200 N m, a light mesh damping ratio of 0.01, $L/D$ = 0.3, and the mean static transmission error $e_{0}$ is set to $b_{0}$ , the coefficient of static transmission error $p_{j}$ in Eq. (2) is set to $0.15$ . Unbalance force is not included in this simulation. The nonlinear equations are solved using MATLAB's ODE 15s with a relative tolerance of $10^{- 5}$ . The time response corresponding to the first 1500 gear mesh periods of the system was discarded from the sampled data to ensure that the responses reached steady-state conditions. The steady-state responses during the last 500 gear mesh periods were employed for the plots shown below. The MLEs of the spur gear system were plotted to identify the onset of chaotic motion. The MLEs converged after 600-time intervals with 0.25 gear mesh periods per interval.
文献[2,6,7,24,40,41]研究了线性化轴承模型的刚度和阻尼系数、静态传动误差、轮齿摩擦等对齿轮系统混沌响应的影响，并在实验中观察到了混沌行为。[6].下文讨论了利用非线性轴承模型改变轴颈轴承参数对齿轮系统混沌响应的影响。模型参数值包括 200 N m 的低扭矩载荷、0.01 的轻啮合阻尼比、 $L/D$ = 0.3、平均静态传输误差 $e_{0}$ 设为 $b_{0}$ ，公式 (2) 中的静态传输误差系数 $p_{j}$ 设为 $0.15$ 。本模拟不包括不平衡力。非线性方程使用 MATLAB 的 ODE 15s 求解，相对容差为 $10^{- 5}$ 。从采样数据中剔除了与系统前 1500 个齿轮啮合周期相对应的时间响应，以确保响应达到稳态条件。下图采用了最后 500 个齿轮啮合周期的稳态响应。绘制正齿轮系统的 MLE 以确定混沌运动的开始。MLE 在 600 个时间间隔（每个间隔 0.25 个齿轮啮合周期）后收敛。

Fig. 19 shows DTE bifurcation diagrams versus rotor rpm while varying torque, bearing lubricant viscosity and radial clearance. The viscosity and clearance are held fixed at 40 mPa s and 105 μm, respectively, while varying torque in Fig. 19(a). The torque and clearance are held fixed at 1100 N m and 105 μm respectively while varying viscosity in Fig. 19(b). The torque and viscosity are held fixed at 1100 N m and 40 mPa s respectively while varying radial clearance in Fig. 19(c). The Poincaré dots are sampled at the gear mesh frequency period ( $Δ t = \frac{2 π}{ω_{g}}$ ).
图 19 显示了在改变扭矩、轴承润滑油粘度和径向游隙时，DTE 分岔图与转子转速的关系。在图 19(a)中，在改变扭矩时，粘度和间隙分别固定为 40 mPa s 和 105 μm。在图 19(b)中，扭矩和间隙分别固定为 1100 牛米和 105 微米，同时改变粘度。在图 19(c)中，扭矩和粘度分别固定为 1100 牛米和 40 毫帕秒，同时改变径向间隙。Poincaré 点取样于齿轮啮合频率周期 ( $Δ t = \frac{2 π}{ω_{g}}$ )。

With a relatively low applied torque (100 N m), the system performs 1× synchronous motions at the gear mesh frequency until the operating speed reaches about 1800 rpm. Then chaotic responses emerge and are maintained at higher speeds. Note that with increasing applied torque from 1000 to 2000 N m, the operating speeds where jump occurs are slightly increased from 2050 to 2150 rpm, which implies that the variation of applied torque has an impact on the natural frequencies of the geared system supported by journal bearings, via a torsional-lateral coupling. Comparison of three cases reveals that the increasing torques tend to mitigate the chaotic responses.
在施加相对较低的扭矩（100 牛米）时，系统以齿轮啮合频率执行 1 倍同步运动，直到运行速度达到约 1800 转/分钟。然后出现混沌响应，并在更高转速下保持不变。需要注意的是，随着施加扭矩从 1000 牛米增加到 2000 牛米，发生跳跃的运行速度从 2050 转/分钟略微增加到 2150 转/分钟，这意味着施加扭矩的变化通过扭转侧向耦合对轴颈轴承支撑的齿轮系统的固有频率产生了影响。对三种情况进行比较后发现，扭矩的增加往往会减轻混乱响应。

Fig. 19(b) shows that the chaos in DTE response starts at lower speed and ranges over a larger set of values as the bearing lubricant viscosity increase from 10 to 70 mPa s. In all three diagrams of Fig. 19(b), with higher operating speeds over 2500 rpm, the chaotic responses turn into period four and two motions, and eventually into period one motion. Fig. 19(c) shows how the chaotic behavior in DTE occurs with the radial bearing clearance values from 74 to 184 μm. Period-doubling bifurcations occur with decreasing speeds from 4500 rpm to lower speeds in the 74 and 105 μm cases, as exemplified by the period one motion turning into period two motion around 3800 rpm, followed by additional period doubling into chaos. In Fig. 19(b) and (c), the lower viscosity and high bearing radial clearance tend to suppress the chaotic behaviors and show relatively stable motions.
图 19(b)显示，随着轴承润滑剂粘度从 10 mPa s 增加到 70 mPa s，DTE 响应从较低转速开始出现混沌，并在较大的数值范围内变化。在图 19(b)的所有三个图中，随着工作转速超过 2500 rpm，混沌响应变成周期四和周期二运动，最终变成周期一运动。图 19(c) 显示了当径向轴承游隙值从 74 μm 到 184 μm 时，DTE 的混沌行为是如何发生的。在 74 μm 和 105 μm 的情况下，随着转速从 4500 rpm 下降到更低的转速，会出现周期加倍分岔，例如，周期一运动在 3800 rpm 左右转变为周期二运动，随后周期加倍进入混沌状态。在图 19(b)和(c)中，较低的粘度和较高的轴承径向游隙倾向于抑制混沌行为，并显示出相对稳定的运动。

Fig. 20 is included to more clearly illustrate the occurrence of chaotic and period-doubling bifurcations corresponding to the 40 mPa s case in Fig. 19(b). Attractors are presented for the four different operating speeds (2300 rpm, 2900 rpm, 3500 rpm and 4100 rpm). The strange attractor and corresponding positive MLE value of 0.056 both confirm chaotic behavior at 2300 rpm. The number of attractors (4 → 2→1) confirms the occurrence of period-doubling bifurcations as the speed decreases from 4100 to 3500 to 2900 rpm.
图 20 更清楚地说明了与图 19（b）中 40 mPa s 情况相对应的混沌和周期加倍分岔的发生。图中显示了四种不同运行速度（2300 rpm、2900 rpm、3500 rpm 和 4100 rpm）下的吸引子。奇异吸引子和相应的正 MLE 值 0.056 都证实了 2300 rpm 时的混沌行为。吸引子的数量（4 → 2→1）证实，当转速从 4100 转/分下降到 3500 转/分再下降到 2900 转/分时，出现了周期加倍分岔。

Fig. 21 shows DTE bifurcation and MLE diagrams versus applied torque, journal bearing radial clearance and viscosity, at the fixed operating speed 1540 rpm. The viscosity is held fixed at 30 mPa s, and the clearance is held fixed at 100 μm while varying torque in Fig. 21(a). The Fourier coefficient of the static transmission error in Eq. (2) is set to $0.15$ and applied to the system. The sampling frequency of Poincaré dots is the gear mesh frequency $ω_{g}$ .
图 21 显示了在固定运行速度 1540 rpm 下，DTE 分岔图和 MLE 图与应用扭矩、轴颈轴承径向游隙和粘度的关系。粘度固定为 30 mPa s，游隙固定为 100 μm，同时在图 21（a）中改变扭矩。将式（2）中静态传输误差的傅立叶系数设为 $0.15$ 并应用于系统。波恩卡莱点的采样频率为齿轮啮合频率 $ω_{g}$ 。

The MLEs quantitatively confirm the existence of chaotic behavior. Fig. 21(a) presents the bifurcation diagram and MLE plot with applied torque ranging from 0 to 2400 N m. Increasing the applied torque is seen to eliminate chaotic motion, and synchronous 1× motion appears at applied torques above 400 N m. The MLEs show positive values, indicating chaos, with low applied torque (<400 N m). For instance, at 200 N m, the corresponding MLE has the positive value of $+$ 0.07.

The radial bearing clearance is varied from 60 to 185 μm in Fig. 21(b) with the rotor speed at 1540 rpm and the lubricant viscosity of 30 mPa s. Chaotic motion occurs in three separate clearance ranges (i.e., 60–70 μm, 95–105 μm, 175–185 μm), and period-doubling routes to chaos with decreasing and increasing parameter values are observed in the ranges 60–80 μm and 130–180 μm, respectively. In Fig. 21(c), the applied torque is held fixed at 100 N m, and the clearance is held fixed at 100 μm while varying viscosity. The figure shows chaos appearing in two ranges of lubricant viscosity: over the ranges of 0–7 mPa s and 55–85 mPa s. All bifurcation diagrams in Fig. 21 display the period-doubling route to chaos, which indicates that the system goes into chaotic states by doubling its period with increasing or decreasing control parameters.

The periodic, period-doubling and chaotic DTE behaviors implied in Fig. 21(b) for bearing clearance variation, are further confirmed in the frequency spectrum, phase portrait, and Poincaré attractor plots in Fig. 22. Three radial bearing clearance values, i.e., 80, 74 and 61 μm, are examined. As can be seen by comparing Fig. 21, Fig. 22, the geared system with MLE = $-$ 0.005 (80 μm case) and $-$ 0.002 (74 μm case) has period two and period four DTE responses, respectively, confirmed by the Poincaré attractors. The 61 μm bearing clearance case in Fig. 22(c) was selected to illustrate chaotic DTE response based on an MLE ( $+$ 0.06), from Fig. 21(b). The phase portrait orbit has a clear aperiodic response, the frequency spectrum has a broadband character, and the corresponding Poincaré dots form a strange attractor. These results show that the system experiences period-doubling bifurcation with decreasing bearing clearance as the system transitions into chaotic responses.

4.3. Effect of gear mesh stiffness on oil whirl

Oil whirl is a rotor dynamics term describing a self-excited, shaft vibration anomaly common to rotating machinery supported by oil film, fixed pad, journal bearings. The dominant symptom of oil whirl is a sub-synchronous limit cycle vibration that is sustained by the journal bearing forces, as opposed to the external excitation such as rotor mass imbalance. In this section, the effect of gear mesh stiffness on the oil whirl instability is investigated. The bearing supported geared rotor pair is modeled as an autonomous (unforced) system with perfectly balanced rotors, and as such is analyzed utilizing the autonomous shooting method. The multiple shooting/continuation method for autonomous cases developed in section 3 is used for identifying periodic solutions. This approach for drawing the bifurcation diagrams was incapable of finding the solutions near saddle-node bifurcations without manual assistance. The continuation algorithm was restarted near the saddle-node bifurcation points and the step size was reduced to help the solutions converge. This typically requires two or three restarts to plot one complete bifurcation diagram.

Fig. 23 shows a bifurcation diagram of the non-dimensional maximum and minimum vertical position of the rotor center in the load direction versus operating speed. Results are shown for two mesh stiffness values, i.e., k₀ = 1e7 N/m and 1e9 N/m. The input torque applied to the gear pair is 1000 N m. Run-up/run-down simulations are conducted separately, using direct numerical integration, and the two results are combined to generate the figure. In Fig. 23(a), as the operating speed increases, the jump-up phenomenon is observed at 37,500 rpm, which corresponds to the onset speed of oil whirl. Jump-down frequency is also predicted at 34,500 rpm with run-down simulation, and consequently, the amplitude of vertical rotor response has been abruptly reduced at the same speed. The result with high mesh stiffness value, i.e., 1e9 N/m in Fig. 23(b), shows a delayed onset speed of oil whirl compared with 1e7 N/m case, such that the jump-up and jump-down frequencies are predicted at 39,500 and 37,000 rpm, respectively. From this result, it is verified that mesh stiffness variation not only affects the high frequency vibration in a geared system but also may change the characteristics of hydrodynamic stability of a geared system supported by journal bearings.

It should be noted that the result with direct numerical integration is incapable of identifying co-existing solutions and the stability of responses. The bifurcation diagram is obtained with the continuation algorithm and the result is presented in Fig. 24. In the case of mesh stiffness 1e7 N/m, in Fig. 24(a), the stable equilibrium position EP which is maintained over low operating speeds switches to unstable motion after crossing Hopf bifurcation point at 37,210 rpm. The transition from stable to unstable response is verified with the continuation algorithm, which provides the eigenvalues of the Jacobian matrix moving out of the unit circle in the complex plane. The point ① corresponds to a Hopf bifurcation, where an unstable subsynchronous response PS 1 emerges and approaches the saddle bifurcation point ②. The amplitudes of the maximum and minimum responses of PS 1 quickly grow as the speed decreases until the saddle point is reached. After the saddle bifurcation, the subsynchronous branch becomes the stable response PS 2 and its maximum and minimum amplitudes slowly approach the bearing clearance limit. The results using the continuation algorithm are drawn in the same figure with the direct numerical integration results in Fig. 24(b). Two results agree well throughout the all speed ranges, except the unstable subsynchronous response PS 1 is only predicted with the continuation algorithm. The vibration frequencies of the branches are shown in Fig. 24(c) and its zoomed version in (d). The frequency range where subsynchronous vibration was observed is located between 17,000 and 23,000 rpm, which corresponds to the 45–50% of operating speed.

Bifurcation diagrams with three different mesh stiffness (1e7, 1e8 and 1e9 N/m) at 1000 N m torque are drawn using the continuation algorithm in Fig. 25. The onset speed of oil whirl is increased from 37,210 to 39,310 rpm as the mesh stiffness increases from 1e7 to 1e8 N/m. In contrast, the onset speed shift is relatively insignificant (110 rpm) between higher mesh stiffness cases, i.e., 1e8 and 1e9 N/m. This result shows that the gear and the journal bearing parameters are coupled and the variation of gear parameter may affect the hydrodynamic stability characteristics of fluid film bearings.

The effect of applied torque on oil whirl onset speed was presented in Refs. [30,33], and the results showed that higher applied torque on the geared system delayed the oil whirl onset speed. This aspect of applied torque effect on oil whirl onset speed may be affected by the amplitude of gear mesh stiffness. The interaction between the applied torque and gear mesh stiffness is investigated with various torque and mesh stiffness values in Table 3. Six applied torque amplitudes from 50 to 1000 N m were applied as mesh stiffness varies from 1e7 to 1e9 N/m. The onset speed of oil whirl is identified by employing the continuation algorithm to find the speed where the jump-up phenomenon occurs. As observed in the reference [30,33], the onset speed generally increases with higher applied torques for all three stiffness cases. From the table, it is also evident that the effect of mesh stiffness magnitude on the onset speed is closely related to the amplitude of torque values. When the lowest applied torque (50 N m) is applied, the onset speed delay caused by the mesh stiffness variation from 1e7 to 1e9 N/m is only 40 rpm. Meanwhile, the speed delay increases up to 2210 rpm with the application of the highest torque value (1000 N m). This result confirms that the transition of oil whirl onset speed with different mesh stiffness is less significant at lower applied torques.

Table 3. Oil whirl onset speed with different torque and mesh stiffness (Identified utilizing Continuation algorithm).

Applied torque	1e7 $N / m$	1e8 $N / m$	1e9 $N / m$	Onset speed difference
50 N m	10,450 rpm	10,480 rpm	10,490 rpm	40 rpm
200 N m	18,760 rpm	18,890 rpm	18,910 rpm	150 rpm
400 N m	25,360 rpm	25,780 rpm	25,800 rpm	440 rpm
600 N m	30,750 rpm	31,390 rpm	31,410 rpm	660 rpm
800 N m	34,780 rpm	36,070 rpm	36,130 rpm	1350 rpm
1000 N m	37,210 rpm	39,310 rpm	39,420 rpm	2210 rpm

5. Conclusions

The nonlinear behavior and bifurcation of a geared rotor system supported by fluid film journal bearings were investigated employing a multiple shooting/continuation algorithm. Nonlinear effects included in the model are nonlinear fluid film force in journal bearing, gear backlash, and time-varying mesh stiffness. The present study confirms that the nonlinearities in a gear pair may induce nonlinear behaviors such as the jump phenomenon, co-existing responses, subharmonic resonances and chaotic responses in the five-degree-of-freedom, gear-journal bearing system model.

The effect of the gear applied torque and journal bearing paramaters on the nonlinear phenomena were investigated. The simulation with varying gear input torque showed that the separation between jump-up and jump-down speed is reduced with high input torques. The high input torque also induced a hardening effect which is not observed in low torque values. It was also confirmed that as bearing $L/D$ ratio and bearing lubricant viscosity are increased, or bearing clearances are decreased, the frequency where the gear nonlinearity-induced jump phenomenon occurs is lowered, and the number of multiple responses is increased, along with the double-sided contact of meshes. In addition, the influence of the input torque and journal bearing parameters on the subharmonic responses were investigated. The simulation results revealed that the high input torque gives rise to the hardening effect as well as the softening effect in the subharmonic resonance region. It is also shown that small bearing $L/D$ ratio, lubricant viscosity and bearing clearance suppress the subharmonic resonances.

The impact of journal bearing parameters on the chaotic response were investigated via direct numerical integration, bifurcation diagrams, spectrums, Poincaré attractors and maximum Lyapunov exponents. As compared to the former studies, which used the linearized bearing stiffness and damping coefficients, the present study utilized the nonlinear journal bearing modeled with the finite element method. The results showed that the high-applied input torques to the gear suppress the chaotic response in the system. Chaotic motions and period-doubling bifurcations were observed at constant operating speed, as the value of lubricant viscosity and bearing clearance varied.

The effect of gear mesh stiffness on the oil whirl phenomenon of the journal bearing was also studied. Using the continuation algorithm, it was verified that the increased gear mesh stiffness delays the onset speed of oil whirl. In addition, the mesh stiffness effect on the oil whirl phenomenon was sensitive to the magnitude of the gear input torque, as the amount of the onset speed delay was found to be more significant with high input torques.

Future investigations for bifurcation and nonlinear dynamics of a gear supported by journal bearings will include thermal effect in the bearing lubricant, other types of hydrodynamic journal bearings such as pressure dam or tilting pad bearings. A more detailed gear-rotor model including a finite element shaft and a lubricant between gear meshes will also be developed. Experimental verification will be conducted to obtain validation results for the theoretical models.

CRediT authorship contribution statement

Dongil Shin: Investigation, Visualization, Writing - original draft, Data curation, Software, Methodology, Conceptualization. Alan Palazzolo: Methodology, Supervision, Validation, Writing - review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgment

The authors gratefully acknowledge the financial support on this research from the Texas A&M Turbomachinery Research Consortium (TRC).

References

[1]
A. Kahraman, R. Singh
Non-linear dynamics of a spur gear pair
J. Sound Vib. (1990), 10.1016/0022-460X(90)90582-K
Google Scholar
[2]
A. Kahraman, R. Singh
Non-linear dynamics of a geared rotor-bearing system with multiple clearances
J. Sound Vib. (1991), pp. 469-506
View PDF View article View in Scopus Google Scholar
[3]
A. Kahraman, R. Singh
Interactions between time-varying mesh stiffness and clearance non-linearities in a geared system
J. Sound Vib. (1991), 10.1016/0022-460X(91)90527-Q
Google Scholar
[4]
G.W. Blankenship, A. Kahraman
Steady state forced response of a mechanical oscillator with combined parametric excitation and clearance type nonlinearity
J. Sound Vib., 185 (1995), pp. 743-765
View PDF View article View in Scopus Google Scholar
[5]
A. Kahraman, G.W. Blankenship
Interactions between commensurate parametric and forcing excitations in a system with clearance
J. Sound Vib. (1996), 10.1006/jsvi.1996.0361
Google Scholar
[6]
A. Kahraman, G.W. Blankenship
Experiments on nonlinear dynamic behavior of an oscillator with clearance and periodically time-varying parameters
J. Appl. Mech., 64 (1997), p. 217, 10.1115/1.2787276
View in Scopus Google Scholar
[7]
A. Raghothama, S. Narayanan
Bifurcation and chaos in geared rotor bearing system by incremental harmonic balance method
J. Sound Vib. (1999), 10.1006/jsvi.1999.2264
Google Scholar
[8]
S. Natsiavas Theodossiades
Non-linear dynamics of gear-pair systems with periodic stiffness and backlash
J. Sound Vib. (2000), 10.1006/jsvi.1999.2490
Google Scholar
[9]
A. Al-Shyyab, A. Kahraman
Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: period-one motions
J. Sound Vib., 284 (1–2) (2005), pp. 151-172
View PDF View article View in Scopus Google Scholar
[10]
A. Al-Shyyab, A. Kahraman
Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: sub-harmonic motions
J. Sound Vib., 279 (1–2) (2005), pp. 417-451
View PDF View article View in Scopus Google Scholar
[11]
S. Zhou, J. Liu, C. Li, B. Wen
Nonlinear behavior of a spur gear pair transmission system with backlash
J. Vibroeng., 16 (2014) (1465), pp. 3850-3861
Google Scholar
[12]
Y. Yang, H. Li, Y. Dai
Nonlinear vibration characteristics of spur gear system subjected to multiple harmonic excitations
Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci., 233 (2019), pp. 6026-6050, 10.1177/0954406219858171
View in Scopus Google Scholar
[13]
Y. Yang, M. Xu, Y. Du, P. Zhao, Y. Dai
Dynamic analysis of nonlinear time-varying spur gear system subjected to multi-frequency excitation
JVC/J. Vib. Contr., 25 (2019), pp. 1210-1226, 10.1177/1077546318814951
View in Scopus Google Scholar
[14]
Y. Yang, L. Cao, H. Li, Y. Dai
Nonlinear dynamic response of a spur gear pair based on the modeling of periodic mesh stiffness and static transmission error
Appl. Math. Model., 72 (2019), pp. 444-469, 10.1016/j.apm.2019.03.026
View PDF View article View in Scopus Google Scholar
[15]
Z. Wang, L. Zhang, Y.Q. Luo, C.Z. Chen
Dynamic model of spur gear pair with modulation internal excitation
Int. J. Rotating Mach. (2017), 10.1155/2017/1264904
Google Scholar
[16]
H. Motahar, F.S. Samani, M. Molaie
Nonlinear vibration of the bevel gear with teeth profile modification
Nonlinear Dynam., 83 (2016), pp. 1875-1884, 10.1007/s11071-015-2452-z
View in Scopus Google Scholar
[17]
J. Yang, T. Lim
Dynamics of coupled nonlinear hypoid gear mesh and time-varying bearing stiffness systems
SAE Int. J. Passeng. Cars - Mech. Syst., 4 (2011), pp. 1039-1049, 10.4271/2011-01-1548
View in Scopus Google Scholar
[18]
J. Wang, T.C. Lim
Effect of tooth mesh stiffness asymmetric nonlinearity for drive and coast sides on hypoid gear dynamics
J. Sound Vib., 319 (2009), pp. 885-903, 10.1016/j.jsv.2008.06.021
View PDF View article View in Scopus Google Scholar
[19]
V.K. Ambarisha, R.G. Parker
Nonlinear dynamics of planetary gears using analytical and finite element models
J. Sound Vib. (2007), 10.1016/j.jsv.2006.11.028
Google Scholar
[20]
M. Zhao, J.C. Ji
Nonlinear torsional vibrations of a wind turbine gearbox
Appl. Math. Model., 39 (2015), pp. 4928-4950, 10.1016/j.apm.2015.03.026
View PDF View article View in Scopus Google Scholar
[21]
X. Qiu, Q. Han, F. Chu
Investigation of parametric instability of the planetary gear under speed fluctuations
Shock Vib., 2017 (2017), 10.1155/2017/6851903
Google Scholar
[22]
C. Gill-Jeong
Nonlinear behavior analysis of spur gear pairs with a one-way clutch
J. Sound Vib. (2007), 10.1016/j.jsv.2007.01.044
Google Scholar
[23]
C. Gill-Jeong
Numerical study on reducing the vibration of spur gear pairs with phasing
J. Sound Vib. (2010), 10.1016/j.jsv.2010.04.005
Google Scholar
[24]
G. Bonori, F. Pellicano
Non-smooth dynamics of spur gears with manufacturing errors
J. Sound Vib. (2007), 10.1016/j.jsv.2007.05.013
Google Scholar
[25]
S. Wei, Q.K. Han, X.J. Dong, Z.K. Peng, F.L. Chu
Dynamic response of a single-mesh gear system with periodic mesh stiffness and backlash nonlinearity under uncertainty
Nonlinear Dynam. (2017), 10.1007/s11071-017-3435-z
Google Scholar
[26]
T.C. Kim, T.E. Rook, R. Singh
Effect of smoothening functions on the frequency response of an oscillator with clearance non-linearity
J. Sound Vib. (2003), 10.1016/S0022-460X(02)01469-4
Google Scholar
[27]
A. Farshidianfar, A. Saghafi
Global bifurcation and chaos analysis in nonlinear vibration of spur gear systems
Nonlinear Dyn. Int. J. Nonlinear Dyn. Chaos Eng. Syst. ISSN Nonlinear Dyn. (2014), pp. 783-806, 10.1007/s11071-013-1104-4
0924
View in Scopus Google Scholar
[28]
X.F. Gou, L.Y. Zhu, D.L. Chen
Bifurcation and chaos analysis of spur gear pair in two-parameter plane
Nonlinear Dynam., 79 (2014), pp. 2225-2235, 10.1007/s11071-014-1807-1
Google Scholar
[29]
Y. Li, T. Chen, X. Wang
Non-linear dynamics of gear pair with dynamic backlash subjected to combined internal and external periodic excitations
JVC/J. Vib. Contr., 22 (2016), pp. 1693-1703, 10.1177/1077546314544350
View in Scopus Google Scholar
[30]
S. Theodossiades, S. Natsiavas
On geared rotordynamic systems with oil journal bearings
J. Sound Vib., 243 (2001), pp. 721-745, 10.1006/jsvi.2000.3430
View PDF View article View in Scopus Google Scholar
[31]
S. Baguet, G. Jacquenot
Nonlinear couplings in a gear-shaft-bearing system
Mech. Mach. Theor., 45 (2010), pp. 1777-1796, 10.1016/j.mechmachtheory.2010.08.009
View PDF View article View in Scopus Google Scholar
[32]
R. Fargère, P. Velex
Influence of clearances and thermal effects on the dynamic behavior of gear-hydrodynamic journal bearing systems
J. Vib. Acoust. Trans. ASME, 135 (2013), 10.1115/1.4025018
Google Scholar
[33]
Z. Liu, Z. Liu, J. Zhao, G. Zhang
Study on interactions between tooth backlash and journal bearing clearance nonlinearity in spur gear pair system
Mech. Mach. Theor., 107 (2017), pp. 229-245, 10.1016/j.mechmachtheory.2016.09.024
View PDF View article Google Scholar
[34]
S. Kim, A.B. Palazzolo
Shooting with deflation algorithm-based nonlinear response and Neimark-Sacker bifurcation and chaos in floating ring bearing systems
J. Comput. Nonlinear Dynam., 12 (2017), pp. 1-15, 10.1115/1.4034733
Google Scholar
[35]
S. Kim, A.B. Palazzolo
Effects of thermo hydrodynamic (THD) floating ring bearing model on rotordynamic bifurcation
Int. J. Non Lin. Mech., 95 (2017), pp. 30-41, 10.1016/j.ijnonlinmec.2017.05.003
View PDF View article View in Scopus Google Scholar
[36]
S. Kim, A.B. Palazzolo
Bifurcation analysis of a rotor supported by five-pad tilting pad journal bearings using numerical continuation
J. Tribol., 140 (2018), pp. 1-12, 10.1115/1.4037699
Google Scholar
[37]
D. Shin, A. Palazzolo
Multiple Steay State Response Solver for Nonlinear Gear, Bearing and Mafunction Simulation
Texas A&M University, College Station, TX, USA (2017)
Techincal report no. TRC-400124-00073
Google Scholar
[38]
Q. Ma, A. Kahraman
Subharmonic resonances of a mechanical oscillator with periodically time-varying, piecewise-nonlinear stiffness
J. Sound Vib. (2006), 10.1016/j.jsv.2005.11.026
Google Scholar
[39]
J. Wang, T.C. Lim, M. Li
Dynamics of a hypoid gear pair considering the effects of time-varying mesh parameters and backlash nonlinearity
J. Sound Vib., 308 (2007), pp. 302-329, 10.1016/j.jsv.2007.07.042
View PDF View article View in Scopus Google Scholar
[40]
L. Xiang, Y. Jia, A. Hu
Bifurcation and chaos analysis for multi-freedom gear-bearing system with time-varying stiffness
Appl. Math. Model., 40 (2016), pp. 10506-10520, 10.1016/j.apm.2016.07.016
View PDF View article View in Scopus Google Scholar
[41]
S. Zhou, G. Song, M. Sun, Z. Ren
Nonlinear dynamic response analysis on gear-rotor-bearing transmission system
JVC/J. Vib. Contr., 24 (2018), pp. 1632-1651, 10.1177/1077546316667178
View in Scopus Google Scholar
[42]
Ali H. Nayfeh, Balakumar Balachandran
Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods
John Wiley & Sons (2008)
Google Scholar

Cited by (32)

Nonlinear parametric vibration and stability analysis of worm drive system with time-varying meshing stiffness and backlash
2024, Journal of Sound and Vibration
The worm drive system with time-varying meshing stiffness may produce parametric resonance, resulting in structural damage. However, few studies focus on the parametric vibration of the drive system. Under certain working conditions, parametric excitation from time-varying meshing stiffness could cause instability and severe vibration. Therefore, determining the working conditions that cause parameter instability is the key to the design of the drive system. A parametric vibration model of the drive system with time-varying meshing stiffness and backlash is developed to simulate actual working conditions. The instantaneous contact line between meshing teeth is obtained based on the geometric coordinate relationship. With the help of the tooth slice model, the static meshing stiffness is derived, which changes with the rotation angle of the worm shaft. The time-varying meshing stiffness model is established by Fourier series expansion. The influence of backlash on the system is integrated into the dynamic driving error, and the analytical function is obtained by polynomial fitting. The proposed meshing stiffness model was verified by static loading experiments, and experimental errors from different sources, such as base deflection and shaft key deformation, were explored. A dynamic response test rig with shaker excitation was developed to verify the proposed vibration model. The Floquet multiplier is applied to judge the stability. The dynamic snap-through phenomenon in the meshing drive and the influence of the damping ratio on the stability boundary are studied. The results show that it is valuable to use dynamic analysis to find measures to avoid parametric resonance.
Vibration stability and bifurcation analysis of spur gear systems based on time-varying backlash and random load
2024, International Journal of Non-Linear Mechanics
Gear backlash and random load affect the stability of gear transmission, which are the main factors causing meshing impact. This study presents the establishment of a six-degree-of-freedom model for the spur gear system, employing the lumped parameter theory and D'Alembert's principle. The model comprehensively incorporates the effects of time-varying meshing stiffness (TVMS), static transmission error (STE), time-varying gear backlash, and random load. The system's equations of motion are solved using the fourth-order Runge-Kutta method. Subsequently, the dynamic response of the constant backlash constant load model is compared with that of the time-varying backlash constant load model and the constant backlash random load model. This research quantitatively dissects the impact of the nonlinear characteristics of backlash and load on the dynamic response of the gear system. This analysis is conducted utilizing a variety of methods, including bifurcation diagrams, Floquet characteristic multipliers, Lyapunov exponents, phase diagrams, and Poincaré sections. The study demonstrates that fluctuations in backlash and load during meshing compromise the stability of the periodic solution, thereby escalating the complexity of predicting system motion. Moreover, observations include the jump phenomena, quasi-periodic and periodic motions, and the trajectory of the system transitioning into a chaotic state via bifurcation. This investigation serves as a valuable reference for future explorations into spur gear systems.
Vibration attenuation characteristics of squeeze film dampers in spiral bevel gear systems
2023, Applied Mathematical Modelling
A central drive bevel gear system of aero-engines extracts power from a high-pressure rotor. The high speed of the high-pressure rotor aggravates the gear system's vibration. This work adopts squeeze film dampers (SFD) as damping devices to suppress excessive vibration of the spiral bevel gear rotor system. This paper establishes the general Reynolds equation of SFD and calculates the nonlinear oil-film force of SFD in real-time based on the finite element method. Besides, the finite element method and Timoshenko beam element are employed to model the flexible gear shaft. Then, the dynamic model of the spiral bevel gear system supported by SFDs is established, including the nonlinear oil-film force of SFD and time-varying meshing stiffness of the gear system. The dynamic response of the spiral bevel gear system is analyzed theoretically and experimentally to study the vibration attenuation characteristics of SFD. The findings indicate that SFD can effectively suppress the vibration amplitude of the system at high speeds. The theoretical analysis and experimental amplitude are within an order of magnitude, and the changing trend is consistent.
Nonlinear dynamics of non-orthogonal offset face gear-bearing transmission system
2023, International Journal of Non-Linear Mechanics
A nonlinear dynamic model of a non-orthogonal offset face gear-bearing transmission system is established, in which nonlinear factors such as input load fluctuation, gear clearance, bearing clearance, and transmission error are considered. The dynamic response of the system is analyzed by solving the differential equation with the numerical method. The subharmonic resonance and primary resonance characteristics of the torsional vibration displacement are analyzed by the multi-scale method. The results show that the change of meshing damping will make the system present rich nonlinear characteristics, and the equivalent static load and meshing damping ratio will affect the resonance characteristics of the system. Choosing the appropriate system parameters can make the system in a stable state and suppress the amplitude of primary resonance and subharmonic resonance.
Computational prediction of structure-borne noise in a two-stage double-helical gearbox using harmonic balance method
2023, Mechanical Systems and Signal Processing
The acoustics of gearboxes contribute as a quality feature to the success of a gear design by increasing consumer satisfaction. Particularly in special machine engineering, ensuring quiet operation is a challenging task because the production and testing of a prototype is usually not a feasible option. A suitable calculation methodology is thus crucial for finding the optimal design. Unfortunately, an increase in accuracy most often results in high simulation times and thus complicates a continuous design support. A promising approach regarding computational efficiency for nonlinear periodic systems is the harmonic balance method (HBM). However, its application to gearboxes is still quite rare or even non-existent in terms of this acoustic motivation. Convergence behavior and memory requirements are ambitious aspects for the HBM. Consequently, tracking the vibration propagation from the gear mesh to the housing in one simulation run requires special treatment. This article examines appropriate implementation recommendations and further investigates the applicability of the HBM for gear acoustics from a practical point of view. Experimental data serves as the basis for validation purposes.
Multi-meshing-state and disengaging-proportion analyses of a gear-bearing system considering deterministic-random excitation based on nonlinear dynamics
2023, Journal of Sound and Vibration
Citation Excerpt :
The gear drivetrain is a non-smooth nonlinear system due to the backlash, time-varying meshing stiffness, tooth surface friction, internal and external excitations and so on [12–14]. In the past decades, many scholars have modeled and analyzed the nonlinear dynamics of gear systems with the consideration of various influencing factors and revealed many nonlinear phenomena, such as chaos, bifurcation, multi-stability, attraction domain, etc [15–17]. Cirelli et al. [18] proposed a refined methodology to simulate the nonlinear dynamic response of spur gears using a multibody contact model with flexible teeth.
Backlash induces gear-teeth separation or back-side tooth mesh (BTM), resulting in gear disengagement (GD), which heavily aggravates the dynamic performance and transmission quality of gear drivetrains. This study focuses on the disengaging proportion (DP) characteristics of a gear- bearing systems based on a multi-meshing-state (MMS) model. A new nonlinear dynamic model of gear-bearing system with MMS and deterministic-random external excitation is established considering the synergistic effect of backlash and non-integer contact ratio. A modified time-varying mesh stiffness with temperature stiffness is included in the model. Afterwards, the effect of the gear pair meshing vibration on bearing radial vibration is studied. The bifurcation and evolution of the MMS behavior are investigated changing the system parameters. An algorithm of DP is proposed in the time domain based on the established MMS model. The characteristics of DP are systematically discussed with the effects of dynamic response type, phase trajectory topology, bifurcation and chaos, and correlation between parameters. The results demonstrate that grazing bifurcation reduces the number of GD and chaos is prone to BTM. The disengaging proportion is highly dependent on the phase trajectory topology and generally increases from period-1 to quasi-periodic and then to chaotic responses. Bifurcation changes the variation trend of DP, thus there is a peak near period-doubling bifurcation and a trough near grazing bifurcation. The chaotic response leads to a jump in DP. The correlation between parameters considerably dominates the DP. This study provides a reference for revealing the mechanism of gear system disengagement and improving the dynamic performance.

View all citing articles on Scopus

View Abstract

[1] [1]
A. Kahraman, R. Singh
Non-linear dynamics of a spur gear pair
J. Sound Vib. (1990), 10.1016/0022-460X(90)90582-K
Google Scholar

[2] [2]
A. Kahraman, R. Singh
Non-linear dynamics of a geared rotor-bearing system with multiple clearances
J. Sound Vib. (1991), pp. 469-506
View PDF View article View in Scopus Google Scholar

[3] [3]
A. Kahraman, R. Singh
Interactions between time-varying mesh stiffness and clearance non-linearities in a geared system
J. Sound Vib. (1991), 10.1016/0022-460X(91)90527-Q
Google Scholar

[4] [4]
G.W. Blankenship, A. Kahraman
Steady state forced response of a mechanical oscillator with combined parametric excitation and clearance type nonlinearity
J. Sound Vib., 185 (1995), pp. 743-765
View PDF View article View in Scopus Google Scholar

[5] [5]
A. Kahraman, G.W. Blankenship
Interactions between commensurate parametric and forcing excitations in a system with clearance
J. Sound Vib. (1996), 10.1006/jsvi.1996.0361
Google Scholar

[6] [6]
A. Kahraman, G.W. Blankenship
Experiments on nonlinear dynamic behavior of an oscillator with clearance and periodically time-varying parameters
J. Appl. Mech., 64 (1997), p. 217, 10.1115/1.2787276
View in Scopus Google Scholar

[7] [7]
A. Raghothama, S. Narayanan
Bifurcation and chaos in geared rotor bearing system by incremental harmonic balance method
J. Sound Vib. (1999), 10.1006/jsvi.1999.2264
Google Scholar

[8] [8]
S. Natsiavas Theodossiades
Non-linear dynamics of gear-pair systems with periodic stiffness and backlash
J. Sound Vib. (2000), 10.1006/jsvi.1999.2490
Google Scholar

[9] [9]
A. Al-Shyyab, A. Kahraman
Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: period-one motions
J. Sound Vib., 284 (1–2) (2005), pp. 151-172
View PDF View article View in Scopus Google Scholar

[10] [10]
A. Al-Shyyab, A. Kahraman
Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: sub-harmonic motions
J. Sound Vib., 279 (1–2) (2005), pp. 417-451
View PDF View article View in Scopus Google Scholar

[11] [11]
S. Zhou, J. Liu, C. Li, B. Wen
Nonlinear behavior of a spur gear pair transmission system with backlash
J. Vibroeng., 16 (2014) (1465), pp. 3850-3861
Google Scholar

[12] [12]
Y. Yang, H. Li, Y. Dai
Nonlinear vibration characteristics of spur gear system subjected to multiple harmonic excitations
Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci., 233 (2019), pp. 6026-6050, 10.1177/0954406219858171
View in Scopus Google Scholar

[13] [13]
Y. Yang, M. Xu, Y. Du, P. Zhao, Y. Dai
Dynamic analysis of nonlinear time-varying spur gear system subjected to multi-frequency excitation
JVC/J. Vib. Contr., 25 (2019), pp. 1210-1226, 10.1177/1077546318814951
View in Scopus Google Scholar

[14] [14]
Y. Yang, L. Cao, H. Li, Y. Dai
Nonlinear dynamic response of a spur gear pair based on the modeling of periodic mesh stiffness and static transmission error
Appl. Math. Model., 72 (2019), pp. 444-469, 10.1016/j.apm.2019.03.026
View PDF View article View in Scopus Google Scholar

[15] [15]
Z. Wang, L. Zhang, Y.Q. Luo, C.Z. Chen
Dynamic model of spur gear pair with modulation internal excitation
Int. J. Rotating Mach. (2017), 10.1155/2017/1264904
Google Scholar

[16] [16]
H. Motahar, F.S. Samani, M. Molaie
Nonlinear vibration of the bevel gear with teeth profile modification
Nonlinear Dynam., 83 (2016), pp. 1875-1884, 10.1007/s11071-015-2452-z
View in Scopus Google Scholar

[17] [17]
J. Yang, T. Lim
Dynamics of coupled nonlinear hypoid gear mesh and time-varying bearing stiffness systems
SAE Int. J. Passeng. Cars - Mech. Syst., 4 (2011), pp. 1039-1049, 10.4271/2011-01-1548
View in Scopus Google Scholar

[18] [18]
J. Wang, T.C. Lim
Effect of tooth mesh stiffness asymmetric nonlinearity for drive and coast sides on hypoid gear dynamics
J. Sound Vib., 319 (2009), pp. 885-903, 10.1016/j.jsv.2008.06.021
View PDF View article View in Scopus Google Scholar

[19] [19]
V.K. Ambarisha, R.G. Parker
Nonlinear dynamics of planetary gears using analytical and finite element models
J. Sound Vib. (2007), 10.1016/j.jsv.2006.11.028
Google Scholar

[20] [20]
M. Zhao, J.C. Ji
Nonlinear torsional vibrations of a wind turbine gearbox
Appl. Math. Model., 39 (2015), pp. 4928-4950, 10.1016/j.apm.2015.03.026
View PDF View article View in Scopus Google Scholar

[21] [21]
X. Qiu, Q. Han, F. Chu
Investigation of parametric instability of the planetary gear under speed fluctuations
Shock Vib., 2017 (2017), 10.1155/2017/6851903
Google Scholar

[22] [22]
C. Gill-Jeong
Nonlinear behavior analysis of spur gear pairs with a one-way clutch
J. Sound Vib. (2007), 10.1016/j.jsv.2007.01.044
Google Scholar

[23] [23]
C. Gill-Jeong
Numerical study on reducing the vibration of spur gear pairs with phasing
J. Sound Vib. (2010), 10.1016/j.jsv.2010.04.005
Google Scholar

[24] [24]
G. Bonori, F. Pellicano
Non-smooth dynamics of spur gears with manufacturing errors
J. Sound Vib. (2007), 10.1016/j.jsv.2007.05.013
Google Scholar

[25] [25]
S. Wei, Q.K. Han, X.J. Dong, Z.K. Peng, F.L. Chu
Dynamic response of a single-mesh gear system with periodic mesh stiffness and backlash nonlinearity under uncertainty
Nonlinear Dynam. (2017), 10.1007/s11071-017-3435-z
Google Scholar

[26] [26]
T.C. Kim, T.E. Rook, R. Singh
Effect of smoothening functions on the frequency response of an oscillator with clearance non-linearity
J. Sound Vib. (2003), 10.1016/S0022-460X(02)01469-4
Google Scholar

[27] [27]
A. Farshidianfar, A. Saghafi
Global bifurcation and chaos analysis in nonlinear vibration of spur gear systems
Nonlinear Dyn. Int. J. Nonlinear Dyn. Chaos Eng. Syst. ISSN Nonlinear Dyn. (2014), pp. 783-806, 10.1007/s11071-013-1104-4
0924
View in Scopus Google Scholar

[28] [28]
X.F. Gou, L.Y. Zhu, D.L. Chen
Bifurcation and chaos analysis of spur gear pair in two-parameter plane
Nonlinear Dynam., 79 (2014), pp. 2225-2235, 10.1007/s11071-014-1807-1
Google Scholar

[29] [29]
Y. Li, T. Chen, X. Wang
Non-linear dynamics of gear pair with dynamic backlash subjected to combined internal and external periodic excitations
JVC/J. Vib. Contr., 22 (2016), pp. 1693-1703, 10.1177/1077546314544350
View in Scopus Google Scholar

[30] [30]
S. Theodossiades, S. Natsiavas
On geared rotordynamic systems with oil journal bearings
J. Sound Vib., 243 (2001), pp. 721-745, 10.1006/jsvi.2000.3430
View PDF View article View in Scopus Google Scholar

[31] [31]
S. Baguet, G. Jacquenot
Nonlinear couplings in a gear-shaft-bearing system
Mech. Mach. Theor., 45 (2010), pp. 1777-1796, 10.1016/j.mechmachtheory.2010.08.009
View PDF View article View in Scopus Google Scholar

[32] [32]
R. Fargère, P. Velex
Influence of clearances and thermal effects on the dynamic behavior of gear-hydrodynamic journal bearing systems
J. Vib. Acoust. Trans. ASME, 135 (2013), 10.1115/1.4025018
Google Scholar

[33] [33]
Z. Liu, Z. Liu, J. Zhao, G. Zhang
Study on interactions between tooth backlash and journal bearing clearance nonlinearity in spur gear pair system
Mech. Mach. Theor., 107 (2017), pp. 229-245, 10.1016/j.mechmachtheory.2016.09.024
View PDF View article Google Scholar

[34] [34]
S. Kim, A.B. Palazzolo
Shooting with deflation algorithm-based nonlinear response and Neimark-Sacker bifurcation and chaos in floating ring bearing systems
J. Comput. Nonlinear Dynam., 12 (2017), pp. 1-15, 10.1115/1.4034733
Google Scholar

[35] [35]
S. Kim, A.B. Palazzolo
Effects of thermo hydrodynamic (THD) floating ring bearing model on rotordynamic bifurcation
Int. J. Non Lin. Mech., 95 (2017), pp. 30-41, 10.1016/j.ijnonlinmec.2017.05.003
View PDF View article View in Scopus Google Scholar

[36] [36]
S. Kim, A.B. Palazzolo
Bifurcation analysis of a rotor supported by five-pad tilting pad journal bearings using numerical continuation
J. Tribol., 140 (2018), pp. 1-12, 10.1115/1.4037699
Google Scholar

[37] [37]
D. Shin, A. Palazzolo
Multiple Steay State Response Solver for Nonlinear Gear, Bearing and Mafunction Simulation
Texas A&M University, College Station, TX, USA (2017)
Techincal report no. TRC-400124-00073
Google Scholar

[38] [38]
Q. Ma, A. Kahraman
Subharmonic resonances of a mechanical oscillator with periodically time-varying, piecewise-nonlinear stiffness
J. Sound Vib. (2006), 10.1016/j.jsv.2005.11.026
Google Scholar

[39] [39]
J. Wang, T.C. Lim, M. Li
Dynamics of a hypoid gear pair considering the effects of time-varying mesh parameters and backlash nonlinearity
J. Sound Vib., 308 (2007), pp. 302-329, 10.1016/j.jsv.2007.07.042
View PDF View article View in Scopus Google Scholar

[40] [40]
L. Xiang, Y. Jia, A. Hu
Bifurcation and chaos analysis for multi-freedom gear-bearing system with time-varying stiffness
Appl. Math. Model., 40 (2016), pp. 10506-10520, 10.1016/j.apm.2016.07.016
View PDF View article View in Scopus Google Scholar

[41] [41]
S. Zhou, G. Song, M. Sun, Z. Ren
Nonlinear dynamic response analysis on gear-rotor-bearing transmission system
JVC/J. Vib. Contr., 24 (2018), pp. 1632-1651, 10.1177/1077546316667178
View in Scopus Google Scholar

[42] [42]
Ali H. Nayfeh, Balakumar Balachandran
Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods
John Wiley & Sons (2008)
Google Scholar

Outline 概要

Cited by (32) 引用 (32)

Figures (25) 数字 (25)