计入结构柔性的风电齿轮箱行星轮系动力学特性研究

谭建军，李 浩，冉 峯，朱才朝，宋朝省，李祖锋; TAN Jian?jun; LI Hao; RAN Feng; ZHU Cai?chao; SONG Chao?sheng; LI Zu?feng

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

计入结构柔性的风电齿轮箱行星轮系动力学特性研究 PDF

- ORCID：
谭建军 ¹
✉
- ORCID：
李浩 ¹
- ORCID：
冉峯 ¹
- ORCID：
朱才朝 ¹
- ORCID：
宋朝省 ¹
- ORCID：
李祖锋 ²

1. 重庆大学机械传动国家重点实验室,重庆 400044； 2. 重庆工商大学智能制造服务国际科技合作基地,重庆 400067

中图分类号： TH113.1； TH132.422

最近更新：2023-11-09

DOI：10.16385/j.cnki.issn.1004-4523.2023.05.002

摘要

针对常规梁/壳单元无法考虑复杂构件几何特征、计算精度低和大规模有限元计算量大、系统级建模与动态分析困难的问题，提出一种计入结构柔性的行星轮系变速动力学建模方法。以某型5 MW级风电齿轮箱低速级行星轮系为研究对象，根据内齿圈结构及边界特征，采用有限元缩聚理论建立内齿圈轮齿与弹性支撑的耦合关系，通过引入啮合副变速表征变量和内齿圈虚拟振动线位移，将驱动轮转角与行星轮系多轮齿啮合状态进行关联，并计入行星架和太阳轮轴柔性，利用界面位移协调条件将各构件耦合，建立行星轮系变速动力学模型。研究结果表明，当5个行星轮受力平衡时，内齿圈轮齿节点瞬态变形整体呈“五角星”形状，并随着行星架旋转；行星轮系轮齿动载荷先增大后减小，变化趋势与单对齿啮合刚度相似；稳态工况下的内齿圈轮齿节点振动位移呈现大幅低频波动和高频振动的叠加特征；输入扭矩突变会破坏多个行星轮之间的动载荷平衡，而柔性内齿圈可在一定程度上吸收部分因冲击载荷引起的构件振动，提高了均载性能。

关键词

机械动力学; 风电齿轮箱; 行星轮系; 结构柔性; 内齿圈

引言

在“十四五”期间，中国风电补贴国家政策将全面取消，对风电机组度电成本提出了严苛要求，而发展10 MW及以上超大功率风电机组是降低度电成本的有效措施之一^［

1］。具有多行星轮功率分流的风电齿轮箱被广泛应用于风电机组中传递兆瓦级功率，是极为重要的传动装置。为了满足更大功率的传递需求，风电齿轮箱行星轮系内齿圈和行星架等构件尺寸将会设计得更大，行星轮数量也将增多，容易在随机气动转矩作用下产生过大的系统振动噪声和动载荷，增大疲劳失效风险。因此，开展计入结构柔性的风电齿轮箱行星轮系动力学特性研究具有重要意义。

国内外学者围绕行星轮系动力学建模、固有特性和动态响应等开展了深入研究，取得了诸多有益的研究成果。Lin等^［

2］利用行星轮系动力学模型分析了轴承支撑刚度、齿轮副啮合刚度等设计参数对系统固有频率与模态动能的影响，发现了模态跃迁现象。Eritenel等^{［参考文献 3

百度学术}3］分析了行星轮系自由振动模态，将行星轮系振动模态归纳为三类振动模式，即中心构件扭转、中心构件平移和行星轮振动模式。Guo等^{［参考文献 4‑5}4‑5］考虑时变啮合刚度、非线性齿面接触和轴承游隙等时变参数激励，分析了行星轮系非线性振动特性。Zhai等^{［参考文献 6

百度学术}6］考虑时变行星架装配误差，建立了多级行星轮系动力学模型，分析了行星轮系动态啮合力。Öztürk等^{［参考文献 7

百度学术}7］以弹性力学为基础建立了行星轮系纯扭转动力学模型，研究了轮齿修形对行星轮系振动特性的影响。为了提高行星轮系动力学模型计算精度，Parker等^{［参考文献 8

百度学术}8］将内齿圈视为弹性连续体，建立行星轮系弹性‑集中参数混合动力学模型，分析了行星轮非均布对系统固有特性的影响，在常见的三类振动模式^{［参考文献 3

百度学术}3］基础之上，发现了内齿圈的弯曲振动模式。Fan等^{［参考文献 9

百度学术}9］、Guan等^{［参考文献 10

百度学术}10］采用壳理论和Timoshenko梁理论计算内齿圈结构柔性，建立行星轮系壳/梁单元‑集中参数混合动力学模型，分析了内齿圈动态变形和系统固有特性。Kahraman等^{［参考文献 11

百度学术}11］采用有限元法建立了行星轮系准静态分析模型，分析了内齿圈柔性对构件动应力和变形的影响。张俊等^{［参考文献 12

百度学术}12］将连续体柔性内齿圈离散为多段等效虚拟弹簧连接的刚性内齿圈段，研究了内齿圈柔性对系统固有特性的影响。许华超等^{［参考文献 13

百度学术}13］采用平面梁单元建立弹性边界柔性直齿内齿圈的振动分析模型，分析了组合激励下系统固有特性。魏静等^{［参考文献 14

百度学术}14］采用轴系单元将内齿圈和行星架进行离散化建模，研究了内齿圈和行星架结构柔性对多级行星轮系动态啮合力的影响。为了考虑复杂结构几何特征，Abousleiman等^{［参考文献 15

百度学术}15］、Portron等^{［参考文献 16

百度学术}16］考虑内齿圈柔性，采用有限元法建立了行星轮系有限元‑集中参数混合动力学模型，分析了系统振动特性。Bettaïeb等^{［参考文献 17

百度学术}17］、Guilbert等^{［参考文献 18

百度学术}18］将模态综合法引入到定轴齿轮动力学建模中，有效提高了建模精度，并降低了系统自由度。随后，许华超等^{［参考文献 19

百度学术}19］采用有限元‑集中参数混合方法建立计入传动轴和机匣结构柔性的直升机主减速器混合动力学模型，研究了机匣结构对系统振动特性的影响。

目前，行星轮系动力学建模方法可以大致划分为集中参数模型^［

2‑7，12］、混合模型^{［参考文献 8‑10}8‑10，13‑19］和有限元模型^{［参考文献 11

百度学术}11］。集中参数模型侧重于对系统固有特性、激励机理及动态载荷分配等进行初步分析，在建模时常将齿轮、轴和轴承简化为一体，以简单径向、弯曲和扭转刚度代替弹性轴的复杂受载状况。此类模型未计入内齿圈、行星架和箱体等构件的弹性变形，整体计算精度不高。混合模型包括梁/壳单元‑集中参数混合模型和有限元‑集中参数混合模型，由于此类建模方法在集中参数模型的基础之上计入了内齿圈、行星架和传动轴等构件的弹性变形，整体计算精度得到提高。相比于梁/壳单元‑集中参数混合模型，有限元‑集中参数混合模型可以考虑复杂的结构几何特征，适用性更好，但多用于恒定转速下的振动分析，少有文献关注计入内齿圈结构柔性的行星轮系变速变载动力学建模。此外，集中参数模型和混合模型常将齿轮副多轮齿啮合简化为等效单点啮合，以简单的综合啮合刚度代替复杂的多轮齿啮合过程，造成在求解轮齿动载荷时精度不高。有限元模型可以考虑复杂结构几何特征，并能较好地模拟轮齿接触状态，全面反映行星轮系各构件受载状况，但此类模型建模过程复杂、计算量庞大、系统级建模分析困难，一般不适用于动态设计场合。

综上所述，为了提高变速变载工况下风电齿轮箱行星轮系动力学性能的预测精度，针对常规梁/壳单元无法考虑复杂构件几何特征、计算精度低和大规模有限元计算量大、系统级建模分析困难的问题，基于有限元‑集中参数混合模型，提出一种计入结构柔性的行星轮系变速动力学建模方法。以某型5 MW级风电齿轮箱低速级行星轮系为研究对象，采用有限元缩聚理论对内齿圈和行星架进行建模，并将齿轮啮合副常规等效单点啮合细化至多对轮齿啮合，建立能够反映复杂构件几何特征、轮齿动载荷且适用于变速运行的行星轮系动力学模型，为风电齿轮箱行星轮系动态设计提供理论基础。

1 行星轮系动力学建模

1.1 子结构缩聚单元

由于实际工程结构往往较为复杂，整体结构有限元模型节点自由度数量庞大，造成计算效率低、收敛慢等问题，因此常采用子结构缩聚方法消除“整体单元”中除需要保留节点以外的所有节点自由度，将子结构系统矩阵规模缩减至可接受的维度。

建立如图1所示的行星架（c）和内齿圈（r）有限元模型，根据其结构特征，在其连接位置、轴承支撑和内齿圈轮齿处设置缩聚点，包括行星架的驱动力矩（缩聚点1）、轴承支撑（缩聚点2和3）和销轴（缩聚点4~8）；内齿圈的轮齿（缩聚点1~93）、箱体连接（缩聚点94和95）和螺栓连接（缩聚点96~120）。定义如图2所示的 O‑XYZ为原点固定在内齿圈中心的绝对坐标系；o_r_j‑x_r_jy_r_jz_r_j为原点固定在内齿圈轮齿缩聚点和螺栓连接缩聚点的绝对坐标系，且与O‑XYZ平行；o_c‑x_cy_cz_c为原点固定在行星架中心并随其转动的随动坐标系；o_cp_i‑x_cp_iy_cp_iz_cp_i为原点固定在销轴缩聚点并随行星架转动的随动坐标系，且与o_c‑x_cy_cz_c平行。此外，内齿圈缩聚点94和95参照O‑XYZ；行星架缩聚点1~3参照o_c‑x_cy_cz_c。R_c为太阳轮到行星轮的中心距离，ϕ_p_i为行星轮i的位置角。

图1 有限元缩聚模型

Fig.1 Finite element condensation model

图2 缩聚点坐标系定义

Fig.2 Definition of coordinate systems of condensation points

内齿圈轮齿缩聚点设置在每个轮齿节圆且位于齿宽中心位置处，则根据角度可以计算出任意内齿圈轮齿节点在O‑XYZ中的坐标如下式所示：

X_{r j} = \{\begin{matrix} x_{r j} \\ y_{r j} \\ z_{r j} \end{matrix}\} = \{\begin{matrix} R_{m r} c o s ϕ_{r j} \\ R_{m r} s i n ϕ_{r j} \\ 0 \end{matrix}\}

（1）

式中 $R_{m r}$ 为内齿圈节圆半径； $ϕ_{r j}$ 为内齿圈轮齿节点 $j$ 的位置角。

同理可得螺栓连接缩聚点坐标，不再详述。

通过柔性多点约束（MPC）^［

20］将缩聚点与对应的界面节点进行连接，采用固定界面模态综合法对行星架和内齿圈进行子结构缩聚，可得缩聚后构件

χ

（

χ

=r， c）在自身参考坐标系下的位移向量为：

u^{χ} = \{\begin{matrix} u_{q}^{χ} \\ u_{w}^{χ} \end{matrix}\} = [\begin{matrix} Φ_{q q}^{χ} & Φ_{q w}^{χ} \\ 0 & I \end{matrix}] \{\begin{matrix} q_{q}^{χ} \\ q_{w}^{χ} \end{matrix}\} = Φ^{χ} q_{q w}^{χ}

（2）

式中 $u_{q}^{χ}$ 和 $u_{w}^{χ}$ 分别为内部节点位移和缩聚点位移的物理坐标； $Φ_{q q}^{χ}$ 和 $Φ_{q w}^{χ}$ 分别为内部节点模态矩阵和约束模态矩阵； $q_{q}^{χ}$ 和 $q_{w}^{χ}$ 分别为内部节点和约束模态的模态坐标； $Φ^{χ}$ 为缩聚后的模态矩阵； $q_{q w}^{χ}$ 为对应的模态坐标；I为单位矩阵。

可得缩聚后行星架和内齿圈的自由振动方程为：

\tilde{M} \ddot{q} + \tilde{C} \dot{q} + \tilde{K} q = 0

（3）

式中 $\tilde{M} = {(Φ^{χ})}^{T} M Φ^{χ}$ ， $\tilde{K} = {(Φ^{χ})}^{T} K Φ^{χ}$ ，其中 $M$ 和 $K$ 分别为全局质量矩阵和全局刚度矩阵； $\tilde{C}$ 为阻尼矩阵，采用Rayleigh阻尼计算^［

21］；q=

\{q_{q}^{r}

q_{w 1}^{r}

…

q_{w 120}^{r}

q_{q}^{c}

q_{w 1}^{c}

…

{q_{w 8}^{c}\}}^{T}

，其中

q_{w k}^{χ} = {\{\underset{u χ k}{\underset{︸}{x_{χ k}, y_{χ k}, z_{χ k}}} \underset{θ χ k}{\underset{︸}{θ_{χ k}^{x}, θ_{χ k}^{y}, θ_{χ k}^{z}}}\}}^{T}

为缩聚点

k

在自身参考坐标系下的广义位移向量。

1.2 啮合与支撑单元

采用如图3所示的广义坐标系对太阳轮和行星轮进行建模，包括太阳轮随动坐标系o_s‑x_sy_sz_s和行星轮随动坐标系o_p_i‑x_p_iy_p_iz_p_i，x_p_i和y_p_i分别为行星轮径向和切向方向。行星轮中心与对应的行星架销轴缩聚点重合，行星架缩聚点转动自由度 $θ_{c υ}^{z}$ （ $υ$ 取值参考图1（b））参照O‑XYZ，其余自由度参照o_c‑x_cy_cz_c和o_cp_i‑x_cp_iy_cp_iz_cp_i（见1.1节）。R_bs，R_bp和R_br分别为太阳轮、行星轮和内齿圈的基圆半径。图4所示为行星轮系啮合副中轮齿动态啮合点 $M_{c}$ 在啮合平面上的位置。

图3 行星轮系啮合副

Fig.3 Meshing gear pairs in the planetary gear train

图4 行星轮系啮合平面

Fig.4 Meshing plane of the planetary gear train

1.2.1 啮合单元

定义太阳轮和行星轮i在自身参考坐标系下的广义位移向量分别为 $X_{s}$ = ${\{u_{s}^{T} θ_{s}^{T}\}}^{T}$ = $\{x_{s}$ $y_{s}$ $z_{s}$ $θ_{s}^{x}$ $θ_{s}^{y}$ ${θ_{s}^{z}\}}^{T}$ 和 $X_{p i} = {\{u_{p i}^{T} θ_{p i}^{T}\}}^{T} = {\{x_{p i} y_{p i} z_{p i} θ_{p i}^{x} θ_{p i}^{y} θ_{p i}^{z}\}}^{T}$ 。根据矢量叠加原理，可得太阳轮‑行星轮i啮合副在啮合点 $M_{c}$ 处的相对位移为 $δ_{s p i} (M_{c})$ ^［

22］。然而，由于行星架缩聚点转动属于刚体转动而非扭角，为了使内齿圈轮齿节点扭振位移可以反馈至内齿圈‑行星轮啮合线上，在文献［22］基础之上，在内齿圈轮齿节点j‑行星轮i啮合点

M_{c}

处的相对位移中引入内齿圈虚拟振动线位移

θ_{c 1}^{z} R_{b r}

，如下式所示：

\begin{array}{l} δ_{r p i}^{j} (M_{c}) = [u_{p i} + θ_{p i} \times (\bar{o_{p i} T_{r p ‑ p}} + \bar{T_{r p ‑ p} M_{c}})] {\bar{I}}_{r p i}^{j} - \\ [u_{r j} + θ_{r j} \times (\bar{O T_{r p ‑ r}} + \bar{T_{r p ‑ r} M_{c}} - X_{r j})] {\bar{I}}_{r p i}^{j} + \\ θ_{c 1}^{z} R_{b r} - δ_{r p i}^{j e} (M_{c}) \end{array}

（4）

式中 $δ_{r p i}^{j e} (M_{c})$ 为啮合误差； ${\bar{I}}_{r p i}^{j}$ 为内齿圈轮齿节点j‑行星轮i啮合副的啮合线方向矢量； $\bar{T_{r p ‑ p} M_{c}}$ 和 $\bar{T_{r p ‑ r} M_{c}}$ 见文献［

22］。

基于 $δ_{s p i} (M_{c})$ 和 $δ_{r p i}^{j} (M_{c})$ ，可得太阳轮‑行星轮i啮合副、内齿圈轮齿节点j ‑行星轮i啮合副在啮合点 $M_{c}$ 处时，各构件6个自由度广义位移向啮合线方向转化的投影矢量分别如下式所示^［

22］：

V_{s p i}^{l} (M_{c}) = \{\begin{matrix} V_{s}^{l} |_{6 \times 1} \\ V_{p i}^{l} |_{6 \times 1} \end{matrix}\}

（5）

V_{r p i}^{j} (M_{c}) = \{\begin{matrix} V_{p i}^{l} |_{6 \times 1} \\ V_{r}^{j} |_{6 \times 1} \end{matrix}\}

（6）

式中 $V_{s p i}^{l}$ 和 $V_{r p i}^{j}$ 为啮合向量； $V_{s}^{l}$ 和 $V_{p i}^{l}$ 分别为在太阳轮‑行星轮i啮合副中对应第 $l$ 对啮合轮齿的太阳轮和行星轮广义位移向量； $V_{p i}^{l}$ 和 $V_{r}^{j}$ 分别为在内齿圈轮齿节点j ‑行星轮i啮合副中对应第 $l$ 对啮合轮齿的行星轮和内齿圈广义位移向量。

为了简化表述，将由内齿圈虚拟振动线位移 $θ_{c 1}^{z} R_{b r}$ （式（4））产生的时变参数激励处理为激振力。

然而，在文献［

2‑14，17‑19］中太阳轮‑行星轮啮合副和内齿圈‑行星轮啮合副常采用等效综合啮合刚度将啮合副中多轮齿啮合过程简化为等效单点啮合，其啮合刚度矩阵

K_{s p} (M_{c})

和

K_{r p} (M_{c})

分别如下式所示：

K_{s p} (M_{c}) = \sum_{i = 1}^{N_{p}} Λ_{s p i}^{r o w} K_{s p i} (M_{c}) V_{s p i} V_{s p i}^{T} Λ_{s p i}^{c o l}

（7）

K_{r p} (M_{c}) = \sum_{i = 1}^{N_{p}} Λ_{r p i}^{r o w} K_{r p i} (M_{c}) V_{r p i} V_{r p i}^{T} Λ_{r p i}^{c o l}

（8）

式中 $N_{p}$ 为行星轮数量； $Λ$ 为定位矩阵，即将啮合刚度矩阵扩展至整个系统矩阵； $K_{s p i}$ 和 $K_{r p i}$ 分别为太阳轮‑行星轮i啮合副、内齿圈‑行星轮i啮合副的综合啮合刚度； $V_{s p i}$ 和 $V_{r p i}$ 为对应的啮合向量，并未细化至每对啮合轮齿。

啮合阻尼矩阵 $C_{s p} (M_{c})$ 和 $C_{r p} (M_{c})$ 的形式分别与式（7）和式（8）相同，其中啮合阻尼采用经验公式计算，见文献［

23］。

基于式（7）和（8），按照“化整为零”的思想，将齿轮副综合啮合刚度按照其重合度等效分解到每对啮合轮齿上，同时为了模拟每对啮合轮齿从啮入到啮出的状态变化，引入轮齿啮合判断系数 $Θ_{s p i}^{l}$ 和 $Θ_{r p i}^{j \to l}$ ，可得：

K_{s p} (M_{c}) =

\sum_{i = 1}^{N_{p}} \{\sum_{l = 1}^{c e i l (ξ_{s p i})} Λ_{s p i}^{l r o w} Θ_{s p i}^{l} k_{s p i}^{l} (M_{c}) V_{s p i}^{l} {(V_{s p i}^{l})}^{T} Λ_{s p i}^{l c o l}\}

（9）

\begin{array}{l} K_{r p} (M_{c}) = \\ \sum_{i = 1}^{N_{p}} \sum_{j = 1}^{N_{r}} \{\begin{matrix} \sum_{l = 1}^{c e i l (ξ_{r p i})} Λ_{r p i}^{j r o w} Θ_{r p i}^{j \to l} k_{r p i}^{j \to l} (M_{c}) \end{matrix} V_{r p i}^{j} {(V_{r p i}^{j})}^{T} Λ_{r p i}^{j c o l}\} \end{array}

（10）

式中 $N_{r}$ 为内齿圈齿数； $k_{s p i}^{l}$ 为太阳轮‑行星轮i啮合副中参与啮合的第 $l$ 对齿啮合刚度； $k_{r p i}^{j \to l}$ 为第j个内齿圈轮齿节点在内齿圈‑行星轮i啮合副中对应参与啮合的第 $l$ 对齿啮合刚度； $c e i l (\cdot)$ 为朝正无穷大方向取整； $ξ$ 为重合度； $V_{s p i}^{l}$ 为太阳轮‑行星轮i啮合副中参与啮合的第 $l$ 对轮齿啮合向量； $V_{r p i}^{j}$ 为内齿圈轮齿节点j‑行星轮i啮合副中参与啮合的轮齿啮合向量。

将式（9）和（10）中 $k_{s p i}^{l}, k_{r p i}^{j \to l}$ ， $Θ_{s p i}^{l}$ 和 $Θ_{r p i}^{j \to l}$ 称为啮合副变速表征变量。为了使行星轮系变速运行，需要建立啮合副变速表征变量与行星轮（驱动轮）转角之间的函数关系。考虑太阳轮‑行星轮1与太阳轮‑行星轮i啮合副的相位差 $γ_{s i}$ 、内齿圈‑行星轮1与内齿圈‑行星轮i啮合副的相位差 $γ_{r i}$ 、太阳轮‑行星轮1与内齿圈‑行星轮1啮合副的相位差 $γ_{s r}$ ^［

24］，以及相邻轮齿啮合顺序^{［参考文献 25

百度学术}25］，可得在太阳轮‑行星轮i啮合副、内齿圈‑行星轮i啮合副中参与啮合的第

l

对轮齿啮合点M_c在啮合平面上的位置分别为：

u_{s p i}^{l} = R_{b p} m o d [θ_{p i}^{z} - γ_{s i} + \frac{2 π (l_{s p i} - 1)}{z_{p}}, \frac{c e i l (ξ_{s p i}) 2 π}{z_{p}}]

（11）

u_{r p i}^{l} = R_{b p} \cdot m o d [θ_{p i}^{z} - γ_{r i} - γ_{s r} + \frac{2 π (l_{r p i} - 1)}{z_{p}}, \frac{c e i l (ξ_{r p i}) 2 π}{z_{p}}]

（12）

式中 $l_{s p i} \in [1, c e i l (ξ_{s p i})]$ ， $l_{r p i} \in [1, c e i l (ξ_{r p i})]$ 。

在啮合区内（见图4），啮合点从啮入到啮出刚好经历一个单对齿啮合刚度变化周期，可得：

\{\begin{matrix} k_{s p i}^{l} = k_{s p i}^{l} (\frac{u_{s p i}^{l}}{ξ_{s p i} p_{b t}} T_{s p i}) \\ k_{r p i}^{j \to l} = k_{r p i}^{j \to l} (\frac{u_{r p i}^{l}}{ξ_{r p i} p_{b t}} T_{r p i}) \end{matrix}

（13）

式中 $T_{s p i}$ 和 $T_{r p i}$ 为单对齿啮合刚度变化周期； $p_{b t}$ 为基圆节距。

同时，结合式（11）和（12），可得啮合副变速表征变量 $Θ_{s p i}^{l}$ 和 $Θ_{r p i}^{j \to l}$ 为：

\{\begin{matrix} Θ_{s p i}^{l} = \{\begin{array}{l} 1, u_{s p i}^{l} \in [0, ξ_{s p i} p_{b t}] \\ 0, 其他 \end{array} \\ Θ_{r p i}^{j \to l} = \{\begin{array}{l} 1, u_{r p i}^{l} \in [0, ξ_{r p i} p_{b t}] \\ 0, 其他 \end{array} \end{matrix}

（14）

1.2.2 支撑单元

行星架销轴缩聚点与行星轮通过轴承支撑单元耦合，但销轴缩聚点的参考坐标系与行星轮节点不一致，需要采用旋转矩阵 $T_{c p}$ ^［

26］将销轴缩聚点的广义位移向量

q_{j k}^{c}

转换至坐标系o_p_i‑x_p_iy_p_iz_p_i：

(q_{j k}^{c})' = T_{c p} q_{j k}^{c}

（15）

基于式（15），可得销轴缩聚点与对应行星轮节点的耦合刚度矩阵 $K_{c p}$ 为：

K_{c p} = [\begin{matrix} K_{b}^{c i} & K_{b}^{c i ‑ p i} \\ K_{b}^{p i ‑ c i} & K_{b}^{p i} \end{matrix}]

（16）

式中 $K_{b}$ 为 $K_{c p} 的$ 子矩阵，仅保留对角线项，其中， $K_{b}^{c i} = d i a g ([k_{x}^{p} 0 k_{z}^{p} k_{θ x}^{p} k_{θ y}^{p} k_{θ z}^{p}])$ ， $K_{b}^{p i ‑ c i} = - K_{b}^{c i}$ ， $K_{b}^{p i} = d i a g ([k_{x}^{p} k_{y}^{p} k_{z}^{p} k_{θ x}^{p} k_{θ y}^{p} k_{θ z}^{p}])$ ， $K_{b}^{c i ‑ p i} = - K_{b}^{p i}$ ，对应的阻尼矩阵 $C_{c p}$ 形式与 $K_{c p}$ 相同，采用Rayleigh阻尼计算^［

21］。太阳轮、内齿圈和行星架的轴承支撑单元形式与

K_{b}^{p i}

一致，不再详述。

1.3 轴段单元

结合如图5所示的太阳轮轴结构特征，采用计入剪切变形影响的修正Euler‑Bernouli 梁单元^［

27］，建立轴段梁单元动力学方程。其中，第

J

个梁单元包括2个节点，每个节点考虑6个自由度，定义在自身参考坐标下的广义位移向量为：

图5 太阳轮轴段单元建模

Fig.5 Modeling of sun shaft segment

X_{s J} = {\{x_{J} y_{J} z_{J} θ_{J}^{x} θ_{J}^{y} θ_{J}^{z} x_{J + 1} y_{J + 1} z_{J + 1} θ_{J + 1}^{x} θ_{J + 1}^{y} θ_{J + 1}^{z}\}}^{T}

（17）

可得第 $J$ 个梁单元自由振动方程为：

M_{J} {\ddot{X}}_{s J} + C_{J} {\dot{X}}_{s J} + K_{J} X_{s J} = 0

（18）

式中 $M_{J}$ 为12×12质量矩阵，采用集中质量形式，即M_J= $\frac{ρ A l_{e n}}{2}$ diag（1，1，1， $\frac{I}{2 A}$ ， $\frac{I}{2 A}$ ， $\frac{I}{A}$ ，1，1，1， $\frac{I}{2 A}$ ， $\frac{I}{2 A}$ ， $\frac{I}{A}$ ），其中， $ρ$ 为材料密度， $A$ 为单元截面面积， $l_{e n}$ 为单元长度， $I$ 为极惯性矩； $K_{J}$ 为12×12刚度矩阵；第 $J$ 个梁单元阻尼矩阵 $C_{J}$ 采用Rayleigh阻尼计算^［

21］。

1.4 激振力单元

齿形误差是齿轮制造加工中产生的典型高频误差，在齿轮啮合过程中表现为周期性位移激励，可以通过三坐标测量机等仪器实测获取。本文仅考虑单齿切向偏差，其与齿轮啮频相关，将其视作简谐函数^［

23］，可得齿轮副啮合误差表达式为：

\{\begin{matrix} δ_{s p i}^{e} (M_{c}) = 0.5 (f_{i_s}^{'} + f_{i_p i}^{'}) s i n (w_{m_s p i} t + φ_{m_s p i}) \\ δ_{r p i}^{j e} (M_{c}) = 0.5 (f_{i_r}^{'} + f_{i_p i}^{'}) s i n (w_{m_r p i} t + φ_{m_r p i}) \end{matrix}

（19）

式中 $f_{i}^{'}$ 为单齿切向偏差^［

28］；

w_{m}

和

φ_{m}

分别为啮频及其初相位；

t

为时间。

此外，结合式（9）和（10），可得由内齿圈虚拟振动线位移 $θ_{c 1}^{z} R_{b r}$ 产生的激振力矩阵为：

\begin{array}{l} F_{r p}^{θ} (M_{c}) = - θ_{c 1}^{z} R_{b r} \cdot \\ \sum_{i = 1}^{N_{p}} \sum_{j = 1}^{N_{r}} \{\begin{matrix} \sum_{l = 1}^{c e i l (ξ_{r p i})} Λ_{r p i}^{j r o w} Θ_{r p i}^{j \to l} k_{r p i}^{j \to l} (M_{c}) V_{r p i}^{j} Λ_{r p i}^{j c o l} \end{matrix}\} \end{array}

（20）

同理，可得由齿轮啮合误差产生的误差激振力矩阵，不再详述。

1.5 系统耦合模型

根据如图6所示的太阳轮、内齿圈、行星架和行星轮节点自由度及其耦合关系，定义系统节点在自身参考坐标下的广义位移向量为：

X_{s y s} = {\{X_{c}^{T} X_{r}^{T} X_{p 1}^{T} \dots X_{p N_{p}}^{T} X_{s}^{T} q_{q}^{c} q_{q}^{r}\}}^{T}

（21）

式中 $X_{c}$ = $\{X_{c 1}^{T}$ $X_{c 2}^{T}$ … ${X_{c N_{c}}^{T}\}}^{T}$ ， $X_{r}$ = $\{X_{r 1}^{T}$ $X_{r 2}^{T}$ … ${X_{r N_{r}}^{T}\}}^{T}$ ， $X_{s} = {\{X_{s 1}^{T} X_{s 2}^{T} \dots X_{s N_{s}}^{T}\}}^{T}$ ，其中下标“ $N_{c}$ ”和“ $N_{s}$ ”分别为行星架和太阳轮轴节点数量； $q_{q}^{c}$ 和 $q_{q}^{r}$ 分别表示行星架和内齿圈的内部节点模态坐标，分别保留前50阶和70阶模态（截止模态频率选为额定工况下3倍高速级啮频）。

图6 行星轮系耦合模型

Fig.6 Coupling model of the planetary gear train

将太阳轮、内齿圈、行星架和行星轮的质量矩阵、刚度矩阵与阻尼矩阵按照式（21）所示顺序进行组装，建立行星轮系动力学模型：

M_{s y s} {\ddot{X}}_{s y s} + C_{s y s} {\dot{X}}_{s y s} + K_{s y s} X_{s y s} = F_{s y s}

（22）

式中 $M_{s y s}$ ， $K_{s y s}$ 和 $C_{s y s}$ 分别为系统质量矩阵、刚度矩阵和阻尼矩阵； $F_{s y s}$ 为系统激振力矩阵，包括输入扭矩（ $T_{i n}$ ）、误差激励、内齿圈虚拟振动线位移激励和负载（ $T_{o u t}$ ），其中负载根据太阳轮转速对转速‑转矩运行特性曲线进行插值得到^［

29］。

行星轮系动力学模型，共计138个节点+120阶内部节点保留模态，948个自由度。

由于模型自由度数量较大，造成行星轮系动力学模型求解困难，本文采用精细积分法（PIM）^［

30］求解式（22），计算各节点振动位移和振动速度，步长取为1×10^‒5。

2 结果分析与讨论

某型5 MW级海上风电齿轮箱低速级行星轮系齿轮参数和轴承支撑刚度分别如表1和2所示，时变啮合刚度和啮合误差如图7所示。其中，啮合刚度和轴承支撑刚度采用ROMAX软件在额定工况下计算得到，齿轮选用5级精度。

表1 行星轮系齿轮参数

Tab.1 Gear parameters in the planetary gear train

部件	齿数	法向模数/mm	螺旋角/(°)	压力角/(°)	齿宽/mm	变位系数
太阳轮	32	24	5	20	530	0.8591
行星轮	29	24	5	20	530	0.8360
内齿圈	93	24	5	20	530	0.7973

图7 时变啮合刚度和啮合误差

Fig.7 Time‑varying meshing stiffness and meshing error

表2 轴承支撑刚度

Tab.2 Bearing supporting stiffness

部件	$k_{x} /$ (N⋅m^-1)	$k_{y}$ /(N⋅m^-1)	$k_{z} /$ (N⋅m^-1)	$k_{θ x} /$ (N⋅m⋅rad^-1)	$k_{θ y} /$ (N⋅m⋅rad^-1)	$k_{θ z}$ /(N⋅m⋅rad^-1)
太阳轮	5×10⁹	5×10⁹	5×10⁹	2×10⁹	2×10⁹	—
行星架(上)	2.2×10⁹	2.2×10⁹	2.4×10⁹	2.5×10⁹	2.5×10⁹	—
行星架(下)	9.7×10⁹	9.7×10⁹	2.8×10⁹	6.1×10⁹	6.1×10⁹	—
行星轮	1×10¹⁰	1×10¹⁰	2.59×10⁹	3.4×10⁸	3.4×10⁸	—
内齿圈	1×10⁹	1×10⁹	1×10⁹	1×10⁹	1×10⁹	1×10⁹

2.1 模型验证

在额定工况下（T_in=4×10⁶ N⋅m），当行星架分别动态旋转至2π/10 rad和2π/5 rad时，行星轮、内齿圈和太阳轮的瞬态受力平衡位置如图8所示，其中红色线表示变形后的状态，蓝色线表示原始状态。图9所示为在行星架3D有限元模型中，当A‑A处缩聚点固定时，在行星架销轴缩聚点施加顺时针（正视B‑B平面）切向力后行星架静态变形，模拟当行星架逆时针旋转时销轴缩聚点变形后的状态。

图8 行星轮系各构件瞬态受力平衡位置

Fig.8 Transient forced equilibrium position of each component in the planetary gear train

注：在X‑Y平面中，内齿圈径向变形放大2000倍，行星架径向变形放大500倍，太阳轮径向变形放大10000倍；在Y‑Z平面中，内齿圈轴向变形放大50000倍，行星架轴向变形放大50000倍，太阳轮轴向变形放大5000倍。

图9 行星架静态变形

Fig.9 Static deformation of the carrier

从图8中可知，行星轮切向位置变化与行星架旋转方向（逆时针旋转）相反，这与图9所示的有限元模型销轴切向变形趋势较为吻合，其主要原因是由于行星轮主要传递周向力矩，造成行星轮沿切向的受力较大。此外，内齿圈轮齿节点的受力平衡位置整体呈“五角星”形状，且会随着行星架逆时针旋转。这是由于内齿圈同时与5个行星轮啮合，导致内齿圈受挤压变形。太阳轮径向的受力平衡位置变化很小，其主要原因是由于行星轮径向对称布置，同时由于斜齿轮轴向分力相对于径向分力较小，内齿圈、行星架和太阳轮的轴向受力平衡位置变化较小。

2.2 啮合副变速表征变量

在恒定输入扭矩（T_in=4×10⁶ N⋅m）和扭矩突变工况（T_in=4×10⁶→2×10⁶→4×10⁶ N⋅m，T_in=4×10⁶→0.2×10⁶→4×10⁶ N⋅m）下，行星架转速变化如图10所示，对应的恒定输入扭矩工况下内齿圈‑行星轮啮合副变速表征变量（ $k_{r p i}^{j \to l}$ 和 $Θ_{r p i}^{j \to l}$ ）如图11所示。从图10中可知，当输入扭矩恒定时，初始阶段行星架转速逐渐增加，随后达到额定转速（输入转速11.27 r/min）并趋于稳定波动；当输入扭矩骤减时，行星架转速随之降低，而当扭矩恢复时，其逐渐增速至额定转速，同时还可以观察到扭矩跌落幅值越大，行星架转速波动也越大。从图11中可知，随着行星架转角增大，内齿圈齿数序号逐渐增大，表明内齿圈轮齿依序参与啮合；并且由于多行星轮啮合存在相位差异，造成啮合副变速表征变量的初始位置沿内齿圈均匀分布。此外，啮合副变速表征变量曲线的两侧边缘出现了明显的周期性“毛刺”，这是由于内齿圈‑行星轮啮合副重合度为非整数时造成的同时参与啮合的齿数周期性变化。

图10 行星架转速变化

Fig.10 Changes of carrier rotation speed

图11 内齿圈‑行星轮啮合副变速表征变量

Fig.11 Characteristic variables for describing the variable speed process of meshing gear pairs of ring gear‑planetary gear

注： No.i表示内齿圈‑行星轮i啮合副

为了进一步分析内齿圈轮齿啮合状态的变化规律，分别提取了内齿圈‑行星轮1啮合副载荷、内齿圈轮齿沿x向的振动位移，分别如图12和13所示。从图12中可知，随着行星架转角增加大，内齿圈‑行星轮1同时参与啮合的齿数在2和3之间周期性波动，并且各轮齿动载荷先增大后减小，整体变化趋势与单对齿啮合刚度相似（如图7所示）。从图13中可知，随着行星架转角增加，内齿圈轮齿振动位移出现明显的循环波动，间隔相位为2π/5 rad，这是由于内齿圈同一轮齿先后与不同行星轮啮合；当行星架转角一定时，内齿圈所有轮齿振动位移整体呈原点对称分布，其主要原因是内齿圈轮齿动载荷方向（沿啮合线方向）与对应的局部坐标系之间的夹角会随着行星架转动而周期性变化。

图12 内齿圈‑行星轮1啮合副载荷

Fig.12 Meshing gear pair loads of ring gear‑planetary gear 1

图13 内齿圈轮齿节点沿x向的振动位移

Fig.13 Vibration displacements of ring gear tooth nodes in the x‑direction

2.3 稳态工况

2.3.1 构件振动

当T_in=4×10⁶ N⋅m时，太阳轮、行星轮和内齿圈振动位移分别如图14~16所示。其中，全柔模型表示本文模型（式（22）），刚性模型则是将式（22）中太阳轮轴、内齿圈等构件视作刚体，但仍保留行星轮系多轮齿动态啮合关系。从图14~16可知，考虑构件结构柔性后，太阳轮、行星轮和内齿圈振动位移幅值出现了不同程度增大。其中，内齿圈振动位移变化最为明显，呈现大幅低频波动（相邻行星轮相位与图13相似）和高频振动（齿轮啮频）的叠加特征；内齿圈轮齿节点1的最大振动位移峰值出现在y向，其次依次为x向和z向，分别为118.4，36.8和12.5 μm，并且观察到x向振动位移峰值位置超前于y向和z向。其主要原因是随着行星架逆时针旋转，当行星轮与内齿圈轮齿节点1刚好啮合时，沿着啮合线方向的轮齿动载荷在内齿圈轮齿节点1处x向的投影量达最大，使得x向振动位移最早出现峰值；随着行星架转角继续增加，轮齿动载荷在y向的投影量逐渐增至最大值，随后出现振动位移峰值，并且在整个啮合过程中，轮齿动载荷沿y向的投影分量大于x向的投影分量。

图14 太阳轮振动位移

Fig.14 Vibration displacements of sun gear

图15 行星轮1振动位移

Fig.15 Vibration displacements of planetary gear 1

图16 内齿圈轮齿节点1振动位移

Fig.16 Vibration displacements of ring gear tooth node 1

2.3.2 轮齿载荷

图17和18分别为不同建模方式对太阳轮‑行星轮1和内齿圈‑行星轮1啮合副载荷的影响。其中，刚度分解表示根据啮合刚度幅值（如图7所示），将全柔模型计算得到的动态啮合力分解到各轮齿上，进而得到轮齿动载荷 ${F_{s p i (r p i)}^{l}|}_{d e}$ ，如下式所示：

图17 太阳轮‑行星轮1啮合副载荷

Fig.17 Meshing gear pair loads of sun gear‑planetary gear 1

图18 内齿圈‑行星轮1啮合副载荷

Fig.18 Meshing gear pair loads of ring gear‑planetary gear 1

{F_{s p i (r p i)}^{l}|}_{d e} = \frac{{[\sum_{l = 1}^{c e i l [ξ_{s p i (r p i)}]} F_{s p i (r p i)}^{l}]|}_{i}}{{[\sum_{l = 1}^{c e i l [ξ_{s p i (r p i)}]} k_{s p i (r p i)}^{l}]|}_{i}} {[k_{s p i (r p i)}^{l}]|}_{i}

（23）

式中 $F_{s p i (r p i)}^{l}$ 为内齿圈（太阳轮）‑行星轮i啮合副中参与啮合的第 $l$ 对轮齿动载荷。

从图17，18中可知，相较于刚度分解模式下的轮齿动载荷，考虑多轮齿动态啮合状态后，其轮齿动载荷峰值更大，且更偏向于左侧。因此若将动态啮合力直接按照啮合刚度幅值简单地分解到各轮齿上，容易低估轮齿动载荷幅值。在同一对内齿圈‑行星轮i啮合副中，多轮齿啮合点沿啮合线方向等间隔分布，此时啮合点越靠近从动轮（内齿圈轮齿节点），内齿圈轮齿节点弹性变形在该处的投影分量将会越小（即式（4）中 $θ_{r j}$ 叉乘系数更小），使得对应的啮合副相对位移更大，加之时变啮合刚度作用，造成相较于刚性模型，考虑构件柔性后，内齿圈‑行星轮啮合副的轮齿动载荷峰值会朝左侧偏移（啮入起始点侧）。因此，考虑内齿圈结构柔性有助于获取更为准确的内齿圈‑行星轮的齿轮动载荷分布。

2.3.3 均载系数

为了描述整个行星轮系因外载荷变动、构件变形、时变啮合刚度与传动误差等动态参数激励引起的均载性能变化，行星轮系均载系数的计算式为^［

31］：

K_{γ} (θ_{c}) = m a x \{\frac{N_{p} \sum_{l = 1}^{c e i l (ξ_{r p i})} F_{r p i}^{l}}{\sum_{i = 1}^{N_{p}} \sum_{l = 1}^{c e i l (ξ_{r p i})} F_{r p i}^{l}}\}

（24）

计算结果如图19所示。从图19中可知，相较于刚性模型，考虑构件柔性后，行星轮均载系数增大。这是由于相较于全柔模型，刚性模型的内齿圈被简化为集中质量点，其振动位移幅值偏小（如图16所示），且无法计算因轮齿动载荷造成的内齿圈局部挤压变形（如图8（a）所示），造成计算得到的内齿圈与5个行星轮之间的动态啮合力较为理想。分析结果表明，若采用内齿圈‑行星轮啮合副的动态啮合力计算行星轮均载系数（式（24）），建议在模型中考虑内齿圈结构柔性。

图19 行星轮均载系数

Fig.19 Load‑sharing coefficient of planetary gear

2.4 变速变载工况

由于风速的随机性，风电齿轮箱需要在变速变载工况下长期运行。为了验证本文模型适用于变速变载工况，通过预设输入扭矩产生突变，获取构件振动位移、轮齿载荷以及均载系数。输入扭矩预设过程：① t ≤1 s时，T_in=4×10⁶ N‧m；② 1 s<t ≤1.2 s时，T_in=4×10⁵ N⋅m；③ 1.2 s<t时，T_in=4×10⁶ N⋅m。

2.4.1 构件振动

图20~22分别为扭矩突变工况下太阳轮、行星轮和内齿圈的振动位移。从图20~22中可知，当输入扭矩骤跌时，各构件的瞬态振动位移幅值明显降低；全柔模型和刚性模型均可以较好地计算太阳轮和行星轮的瞬态振动响应；相较于刚性模型，全柔模型可以计算得到更为详细的内齿圈轮齿节点瞬态振动响应，有助于指导内齿圈结构优化设计。

图20 太阳轮振动位移

Fig.20 Vibration displacements of sun gear

图21 行星轮1振动位移

Fig.21 Vibration displacements of planetary gear 1

图22 内齿圈轮齿节点1振动位移

Fig.22 Vibration displacements of ring gear tooth node 1

2.4.2 轮齿载荷

图23和24分别为扭矩突变工况下太阳轮‑行星轮1和内齿圈‑行星轮1啮合副载荷的变化情况。从图23，24可知，在全柔模型中，输入扭矩跌落会降低轮齿动载荷，且不同轮齿动载荷降幅存在差异，太阳轮‑行星轮1的轮齿动载荷分别降至-177.3 和-188 kN⋅m，内齿圈‑行星轮1的轮齿动载荷分别降至-202.5和-110.8 kN⋅m，这意味对应的轮齿已发生脱齿，且容易产生反向接触。刚性模型计算得到的太阳轮‑行星轮啮合副瞬态载荷变化与全柔模型计算结果差异较小；由于内齿圈结构柔性会吸收部分因冲击载荷引起的轮齿节点振动，造成刚性模型计算得到的内齿圈‑行星轮啮合副中部分轮齿瞬态载荷略大于全柔模型计算结果。此外，相比于全柔模型计算结果，刚度分解模式下各轮齿动载荷均大于零，表明该方法不适用于计算变速变载工况下的轮齿动载荷。

图23 太阳轮‑行星轮1啮合副载荷

Fig.23 Meshing gear pair loads of sun gear‑planetary gear 1

图24 内齿圈‑行星轮1啮合副载荷

Fig.24 Meshing gear pair loads of ring gear‑planetary gear 1

2.4.3 均载系数

图25所示为扭矩突变工况下行星轮均载系数的变化情况。从图25可知，扭矩突变会显著增大行星轮均载系数，表明此时内齿圈与5个行星轮之间的动态啮合力出现了明显差异。此外，由于柔性内齿圈的吸振作用，造成在扭矩突变期间内，柔性模型计算得到的行星轮均载系数明显小于刚性模型计算结果。因此，在计算变速变载工况下的行星轮均载系数时，建议考虑内齿圈结构柔性，可以提高计算结果精度。

图25 行星轮均载系数

Fig.25 Load‑sharing coefficient of planetary gear

3 结论

本文将行星轮系啮合副常规等效单点啮合细化至多对轮齿啮合，通过引入啮合副变速表征变量和内齿圈虚拟振动线位移，构建了驱动轮转角与行星轮系多轮齿啮合状态的映射关系，并基于模态缩减理论，将内齿圈轮齿与弹性支撑进行耦合，建立计入结构柔性的风电齿轮箱行星轮系变速动力学模型，分析了工况变化对啮合副变速表征变量、构件振动、轮齿载荷和均载系数的影响，主要结论如下：

（1）当5个行星轮受力平衡时，行星轮切向变形量方向与行星架旋转方向相反；内齿圈轮齿节点瞬态变形整体呈“五角星”形状，且会随着行星架旋转；行星轮系各构件径向受力平衡位置变化最大，轴向受力平衡位置变化较小。

（2）在稳态工况下，内齿圈轮齿节点振动位移呈现大幅低频波动和高频振动的叠加特征；行星轮系轮齿动载荷先增大后减小，整体变化趋势与单对齿啮合刚度相似；轮齿动载荷会使柔性内齿圈产生局部挤压变形，影响均载性能。

（3）当输入扭矩跌落时，行星轮系各构件瞬态振动位移幅值将骤减；扭矩突变会破坏多个行星轮之间的动载荷平衡，恶化均载性能；柔性内齿圈可在一定程度上吸收部分因冲击载荷引起的构件振动，提高行星轮均载性能。

本文方法可以直接考虑传动构件复杂结构几何特征与多轮齿啮合状态，避免因简化处理带来的精度不足，为实现重载风电齿轮箱行星轮系变速动力学建模、提高动态特性分析精度提供了一种有效技术手段。

参考文献

朱才朝，周少华，张亚宾，等.滑动轴承在风电齿轮箱中的应用现状与发展趋势［J］.风能，2021（9）：38‑42. [百度学术]

ZHU Caichao， ZHOU Shaohua， ZHANG Yabin， et al. Application status and development trend of sliding bearings in wind turbine gearboxes［J］. Wind Energy， 2021（9）：38‑42. [百度学术]

Lin J， Parker R G. Sensitivity of planetary gear natural frequencies and vibration modes to model parameters［J］. Journal of Sound and Vibration， 1999， 228（1）： 109‑128. [百度学术]

Eritenel T， Parker R G. Modal properties of three‑dimensional helical planetary gears［J］. Journal of Sound and Vibration， 2009， 325（1‑2）： 397‑420. [百度学术]

Guo Y， Parker R G. Dynamic analysis of planetary gears with bearing clearance［J］. Journal of Computational and Nonlinear Dynamics， 2012， 7（4）：041002. [百度学术]

Guo Y， Keller J， Parker R G. Nonlinear dynamics and stability of wind turbine planetary gear sets under gravity effects［J］. European Journal of Mechanics‑A/Solids， 2014， 47： 45‑57. [百度学术]

Zhai H， Zhu C， Song C， et al. Influences of carrier assembly errors on the dynamic characteristics for wind turbine gearbox［J］. Mechanism and Machine Theory， 2016， 103： 138‑147. [百度学术]

Öztürk V Y， Cigeroglu E， Özgüven H N. Ideal tooth profile modifications for improving nonlinear dynamic response of planetary gear trains［J］. Journal of Sound and Vibration， 2021， 500： 116007. [百度学术]

Parker R G， Wu X. Vibration modes of planetary gears with unequally spaced planets and an elastic ring gear［J］. Journal of Sound and Vibration， 2010， 329（11）： 2265‑2275. [百度学术]

Fan Z， Zhu C， Song C. Dynamic analysis of planetary gear transmission system considering the flexibility of internal ring gear［J］. Iranian Journal of Science and Technology， Transactions of Mechanical Engineering， 2020， 44（3）： 695‑706. [百度学术]

Guan X， Tang J， Hu Z， et al. Dynamic analysis of spur gear pair established by flexible ring and time‑varying mesh model［J］. Journal of the Brazilian Society of Mechanical Sciences and Engineering， 2022， 44（4）： 1‑21. [百度学术]

Kahraman A， Kharazi A A， Umrani M. A deformable body dynamic analysis of planetary gears with thin rims［J］. Journal of Sound and Vibration， 2003， 262（3）： 752‑768. [百度学术]

张俊，宋轶民，王建军.计入齿圈柔性的直齿行星传动动力学建模［J］.机械工程学报，2009，45（12）：29‑36. [百度学术]

ZHANG Jun， Song Yimin， WANG Jianjun. Dynamic modeling for spur planetary gear transmission with flexible ring gear［J］. Journal of Mechanical Engineering， 2009，45（12）：29‑36. [百度学术]

许华超，秦大同.内激励下弹性边界柔性直齿内齿圈振动响应研究［J］.机械工程学报，2018，54（9）：161‑167. [百度学术]

XU Huachao， QIN Datong. Vibration response of flexible spur ring gear with elastic foundation under internal excitation［J］. Journal of Mechanical Engineering， 2018，54（9）：161‑167. [百度学术]

魏静，张爱强，秦大同，等.考虑结构柔性的行星轮系耦合振动特性研究［J］.机械工程学报，2017，53（1）：1‑12. [百度学术]

WEI Jing， ZHANG Aiqiang， QIN Datong， et al. Coupling vibration analysis for planetary gear system considering flexible structure［J］. Journal of Mechanical Engineering， 2017，53（1）：1‑12. [百度学术]

Abousleiman V， Velex P， Becquerelle S. Modeling of spur and helical gear planetary drives with flexible ring gears and planet carriers［J］. Journal of Mechanical Design， 2007，129（1）： 95‑106. [百度学术]

Portron S， Velex P， Abousleiman V. A hybrid model to study the effect of tooth lead modifications on the dynamic behavior of double helical planetary gears［J］. Proceedings of the Institution of Mechanical Engineers， Part C： Journal of Mechanical Engineering Science， 2019， 233（21‑22）： 7224‑7235. [百度学术]

Bettaïeb M N， Velex P， Ajmi M. A static and dynamic model of geared transmissions by combining substructures and elastic foundations—applications to thin‑rimmed gears［J］. Journal of Mechanical Design， 2007，129（2）：184‑194. [百度学术]

Guilbert B， Velex P， Dureisseix D， et al. A mortar‑based mesh interface for hybrid finite‑element/lumped‑parameter gear dynamic models‑applications to thin‑rimmed geared systems［J］. Journal of Mechanical Design， 2016， 138（12）：123301. [百度学术]

许华超，秦大同，刘长钊，等.计入结构柔性的直升机主减速器振动特性分析［J］.航空动力学报，2019，34（5）：1020‑1028. [百度学术]

XU Huachao， QIN Datong， LIU Changzhao， et al. Vibration characteristics for helicopter main gearbox considering structural flexibility［J］. Journal of Aerospace Power，2019，34（5）：1020‑1028. [百度学术]

Helsen J， Vanhollebeke F， Marrant B， et al. Multibody modelling of varying complexity for modal behaviour analysis of wind turbine gearboxes［J］. Renewable Energy， 2011， 36（11）： 3098‑3113. [百度学术]

Wang J. Rayleigh coefficients for series infrastructure systems with multiple damping properties［J］. Journal of Vibration and Control， 2015， 21（6）： 1234‑1248. [百度学术]

Abousleiman V， Velex P. A hybrid 3D finite element/lumped parameter model for quasi‑static and dynamic analyses of planetary/epicyclic gear sets［J］. Mechanism and Machine Theory，2006，41（6）： 725‑748. [百度学术]

Wei J， Zhang A， Qin D， et al. A coupling dynamics analysis method for a multistage planetary gear system［J］. Mechanism and Machine Theory， 2017， 110： 27‑49. [百度学术]

Parker R G， Lin J. Mesh phasing relationships in planetary and epicyclic gears［J］. Journal of Mechanical Design， 2004， 126（2）： 365‑370. [百度学术]

Liu C， Qin D， Liao Y. Dynamic model of variable speed process for herringbone gears including friction calculated by variable friction coefficient［J］. Journal of Mechanical Design， 2014， 136（4）： 041006. [百度学术]

Tan J， Zhu C， Song C， et al. Dynamic modeling and analysis of wind turbine drivetrain considering platform motion［J］. Mechanism and Machine Theory， 2019， 140： 781‑808. [百度学术]

常乐浩，刘更，吴立言.齿轮综合啮合误差计算方法及对系统振动的影响［J］.机械工程学报，2015，51（1）：123‑130. [百度学术]

CHANG Lehao， LIU Geng， WU Liyan. Determination of composite meshing errors and its influence on the vibration of gear system［J］. Journal of Mechanical Engineering， 2015，51（1）：123‑130. [百度学术]

成大先. 机械设计手册［M］. 六版. 北京：化学工业出版社， 2017. [百度学术]

CHENG Daxian. Mechanical Design Manual［M］. 6th ed. Beijing： Chemical Industry Press， 2017. [百度学术]

谭建军. 双馈风电机组传动链系统耦合动力学研究［D］. 重庆：重庆大学， 2019. [百度学术]

TAN Jianjun. Coupling dynamics research of a doubly‑fed wind turbine drivetrain［D］. Chongqing： Chongqing University，2019. [百度学术]

Wan X， Zhi Q. Precise integration method for LQG optimal measurement feedback control problem［J］. Applied Mathematics and Mechanics， 2000， 21（12）： 1417‑1422. [百度学术]

Qiu X， Han Q， Chu F. Load‑sharing characteristics of planetary gear transmission in horizontal axis wind turbines［J］. Mechanism and Machine Theory， 2015， 92： 391‑406. [百度学术]