CN106354975A - 一种获取行星齿轮错位量的有限元方法 - Google Patents

Abstract

本发明涉及一种获取行星齿轮错位量的有限元方法，包括以下步骤：1)建立轴部件的有限元模型；2)建立行星架的有限元模型；3)建立滚动轴承的有限元模型；4)建立行星齿轮有限元模型；5)建立行星齿轮传动系统静力学模型：根据轴部件、滚动轴承、行星架、齿轮之间的连接关系，采用有限元方法组集各部件的刚度矩阵，建立完整的行星齿轮传动系统静力学模型；6)传动系统静力学求解：采用牛顿‑拉弗森方法求解系统非线性静力学方程；7)计算行星齿轮错位量：即太阳轮和行星轮的错位量以及齿圈和行星轮的错位量。

Description

一种获取行星齿轮错位量的有限元方法

技术领域

本发明涉及一种获取行星齿轮错位量的有限元方法，属于机械传动技术领域。

背景技术

行星齿轮广泛应用于驱动桥、变速箱等机械传动系统，为了确保行星齿轮设计满足性能要求，国际上普遍采用ISO6336对行星齿轮的承载能力进行计算校核。行星齿轮传动系统受载变形，会引起齿轮错位，导致齿轮受载不均，错位量越大，齿轮的承载能力越差。在ISO6336中，齿轮错位量f_sh对齿轮性能的影响体现为齿间载荷分配系数K_Fα、K_Hα和齿向载荷分布系数K_Fβ、K_Hβ，其中，齿间载荷分配系数体现了齿轮错位量对轮齿之间载荷分配不均匀程度的影响，齿向载荷分布系数体现了齿轮错位量对齿宽方向载荷分布不均匀程度的影响。

行星齿轮错位量f_sh由传动系统的结构形式和载荷工况共同决定，目前有关行星齿轮错位量计算方法的研究鲜有报道，式(1)为ISO6336(参考文献：ISO 6336-1-2007Calculation of load capacity of spur and helical gears—Part 1:Basicprinciples,introduction and general influence factors)给出的单对齿轮副错位量f_sh的近似计算公式：

$f_{sh}=\frac{F_m}{b}0.023\[\left|B^{*}+K^{\′}\frac{ls}{d_1^2}{\left(\frac{d_1}{d_{sh}}\right)}^4-0.3\right|+0.3\]{\left(\frac{b}{d_1}\right)}^2---\left(1\right)$

上式中，F_m为分度圆上的平均端面切向力；b为齿宽；B^*为功率系数；K′为主动轮结构系数；l为轴承跨距；s为主动轮齿宽中点至轴承跨距中点的距离；d₁为主动轮分度圆直径；d_sh为主动轮弯曲变形当量直径。

ISO6336齿轮错位量计算方法主要存在以下不足：1)公式过于简化，仅适用于简单结构形式和简单受力状态下单对齿轮副错位量的近似计算，无法准确体现行星齿轮传动系统的结构特征，尤其没有考虑滚动轴承的非线性刚度特性对行星齿轮传动系统受载变形的影响。2)无法准确体现行星齿轮传动系统在承受复杂载荷工况时各行星轮分支齿轮错位量的差异。3)无法准确体现行星轮在不同方位时的齿轮错位量变化。

发明内容

针对上述问题，本发明的目的是提供一种获取行星齿轮错位量的有限元方法，该方法能够克服上述技术问题。

为实现上述目的，本发明采取以下技术方案：一种获取行星齿轮错位量的有限元方法，包括以下步骤：1)建立轴部件的有限元模型：采用欧拉-伯努利空间梁单元建立轴部件的有限元模型，获取轴部件的刚度矩阵；2)建立行星架的有限元模型：采用刚性梁单元对行星架单元进行模拟，获得行星架的刚度矩阵；3)建立滚动轴承的有限元模型：采用具有非线性耦合刚度特性的轴承单元对滚动轴承进行模拟，获得滚动轴承的刚度矩阵；4)建立行星齿轮有限元模型：分别建立太阳轮与行星轮的等效啮合模型和齿圈与行星轮的等效啮合模型，获得太阳轮与行星轮的等效啮合刚度矩阵和齿圈与行星轮的等效啮合刚度矩阵；5)建立行星齿轮传动系统静力学模型：根据轴部件、滚动轴承、行星架、齿轮之间的连接关系，采用有限元方法组集各部件的刚度矩阵，建立完整的行星齿轮传动系统静力学模型；6)传动系统静力学求解：采用牛顿-拉弗森方法求解系统非线性静力学方程；7)计算行星齿轮错位量：即太阳轮和行星轮的错位量以及齿圈和行星轮的错位量。

所述步骤1)中的轴部件包括输入轴、行星架轴、齿圈轴、行星轮销轴和行星轮轴，其中，所述行星轮销轴和行星轮轴的数量均等于行星轮的个数。

所述步骤4)中建立太阳轮与行星轮的等效啮合模型并获得太阳轮与行星轮的等效啮合刚度矩阵的过程如下：

太阳轮与行星轮为外啮合齿轮副，建立系统全局坐标系OXYZ和行星轮局部坐标系o₂x₂y₂z₂，其中，系统全局坐标系OXYZ的Z轴正方向与X、Y轴满足右手定则，坐标原点O定义为太阳轮中心；定义行星轮方位角θ为x₂与X的夹角，太阳轮中心梁单元节点为o₁，行星轮中心梁单元节点为o₂，太阳轮等效啮合节点为p₁，行星轮在太阳轮一端的等效啮合节点为p₂₁，p₁和p₂₁均有6个自由度；o₁与p₁、o₂与p₂₁之间采用刚性梁单元连接，刚度矩阵表示为K_g21，p₁和p₂₁之间采用沿齿轮等效啮合力作用线方向的空间弹簧单元连接，坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₁表示为式(3)：

n₀₂₁＝[n_x21,n_y21,n_z21]^T (3)

上式中，n_x21、n_y21、n_z21分别为n₀₂₁在各坐标轴方向上的分量系数；

等效啮合节点p₁和p₂₁的坐标位置相同，根据式(4)求得：

$\left\{\begin{array}{c}\stackrel{\‾}{o_1{p}_1}=\frac{N_1}{N_1+N_2}L\\ \stackrel{\‾}{o_2{p}_{21}}=\frac{N_2}{N_1+N_2}L\end{array}\right.---\left(4\right)$

上式中，N₁为太阳轮齿数；N₂为行星轮齿数；L为太阳轮与行星轮的中心距；

太阳轮与行星轮齿轮副的名义啮合分力如式(5)所示：

$\left\{\begin{array}{c}{F}_{t021}=2T_1/d_1\\ F_{r021}=F_{t021}{\mathrm{tan\α}}_n/cos\mathrm{\β}\\ F_{a021}=F_{t021}tan\mathrm{\β}\end{array}\right.---\left(5\right)$

上式中，F_t021为名义切向力；F_r021为名义径向力；F_a021为名义轴向力；T₁为行星轮分支输入转矩大小，设太阳轮总输入转矩T平均分配给各行星轮分支，则T₁＝T/N_p；d₁为太阳轮分度圆直径；α_n为法向压力角；β为螺旋角；

如果齿轮存在变位，齿轮副的真实啮合分力如式(6)所示：

$\left\{\begin{array}{c}{F}_{tw21}=T/\stackrel{\‾}{o_1{p}_1}\\ F_{nv21}=\sqrt{F_{t021}^2+F_{r021}^2-F_{tw21}^2}\\ F_{aw21}=F_{a021}\end{array}\right.---\left(6\right)$

上式中，F_tw21为真实切向力；F_rw21为真实径向力；F_aw21为真实轴向力；

太阳轮与行星轮齿轮副的法向啮合力如式(7)所示：

$F_{n21}=\sqrt{F_{tw21}^2+F_{rw21}^2+F_{aw21}^2}---\left(7\right)$

太阳轮和行星轮在坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₁的各分量系数如式(8)所示：

$\left\{\begin{array}{c}{n}_{x21}=-F_{rw21}/F_{n21}\\ n_{y21}=k_L{F}_{tw21}/F_{n21}\\ n_{z21}=-k_L{k}_R{F}_{aw21}/F_{n21}\end{array}\right.---\left(8\right)$

上式中，k_L为工况系数，当输入转矩T的方向为+Z时，k_L＝1，当输入转矩T的方向为-Z时，k_L＝-1；k_R为齿轮旋向系数，当太阳轮为右旋时，k_R＝1，当太阳轮为左旋时，k_R＝-1；

全局坐标系下太阳轮与行星轮齿轮副等效啮合力作用线方向矢量如式(9)所示：

$n_{21}^T=n_{021}^{T}H\left(\mathrm{\θ}\right)---\left(9\right)$

上式中，H(θ)为行星轮方位角θ对应的坐标变换矩阵，其由式(10)得到：

$H\left(\mathrm{\θ}\right)=\left[\begin{array}{ccc}cos\mathrm{\θ}& -\sin \mathrm{\θ}& 0\\ sin\mathrm{\θ}& \cos \mathrm{\θ}& 0\\ 0& 0& 1\end{array}\right]---\left(10\right)$

太阳轮和行星轮等效啮合节点平动自由度之间的刚度耦合如式(11)所示：

$K_{m21}=k_{m21}{\left[\begin{array}{cc}{n}_{21}^T{n}_{21}& -n_{21}^T{n}_{21}\\ -n_{21}^T{n}_{21}& n_{21}^T{n}_{21}\end{array}\right]}_{6\×6}---\left(11\right)$

上式中，k_m21为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得；

太阳轮和行星轮等效啮合节点转动自由度之间的刚度耦合如式(12)所示：

$K_{m\mathrm{\θ}21}=k_{m\mathrm{\θ}21}{\left[\begin{array}{cc}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& -n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\\ -n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\end{array}\right]}_{6\×6}---\left(12\right)$

上式中，n_θ21为转动自由度耦合矢量，如式(13)所示：

n_θ21＝[-n_y21,n_x21,0]H(θ)； (13)

k_mθ21为太阳轮和行星轮啮合弯曲刚度系数，如式(14)所示：

$\begin{array}{c}{k}_{m\mathrm{\θ}21}=\frac{{\mathrm{\ΔM}}_{21}}{{\mathrm{\Δ\γ}}_{21}}=\frac{1}{{\mathrm{\Δ\γ}}_{21}}\[\underset{b/2}{\overset{b}{\∫}}{\mathrm{\Δ\γ}}_{21}z\frac{k_{m21}}{b}\left(z-\frac{b}{2}\right)dz-\underset{0}{\overset{b/2}{\∫}}{\mathrm{\Δ\γ}}_{21}z\frac{k_{m21}}{b}\left(\frac{b}{2}-z\right)dz\]\\ =\frac{b^2}{12}{k}_{m21}\end{array}---\left(14\right)$

上式中，ΔM₂₁表示太阳轮和行星轮接触线产生单位转角Δγ₂₁时，对齿宽中点的弯矩；b为有效齿宽，取太阳轮和行星轮齿宽的较小值；z为齿宽方向的局部坐标；dz为齿宽方向的微变量；

式(15)为太阳轮和行星轮等效啮合节点自由度对应的完整的齿轮等效啮合刚度矩阵：

$K_{m21}={\left[\begin{array}{cccc}{k}_{m21}{n}_{21}^T{n}_{21}& 0& -k_{m21}{n}_{21}^T{n}_{21}& 0\\ 0& k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& 0& -k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\\ -k_{m21}{n}_{21}^T{n}_{21}& 0& k_{m21}{n}_{21}^T{n}_{21}& 0\\ 0& -k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& 0& k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\end{array}\right]}_{12\×12}---\left(15\right)$

上式中，n₂₁为全局坐标系下的等效啮合力作用线方向矢量，如式(9)所示；n_θ21为转动自由度耦合矢量，如式(13)所示；k_m21为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得；k_mθ21为齿轮啮合弯曲刚度系数，如式(14)所示。

所述步骤4)中建立齿圈与行星轮的等效啮合模型并获得齿圈与行星轮的等效啮合刚度矩阵的过程如下：

齿圈与行星轮为内啮合齿轮副，定义齿圈中心梁单元节点为o₃，其坐标位置与全局坐标系原点O、太阳轮中心梁单元节点o₁位置相同，行星轮中心梁单元节点为o₂，齿圈等效啮合节点为p₃，行星轮在齿圈一端的等效啮合节点为p₂₃，p₃和p₂₃均有6个自由度；o₃与p₃、o₂与p₂₃之间采用刚性梁单元连接，刚度矩阵表示为K_g23，p₃和p₂₃之间采用沿齿轮等效啮合力作用线方向的空间弹簧单元连接，坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₃表示为式(16)：

n₀₂₃＝[n_x23,n_y23,n_z23]^T (16)

上式中，n_x23、n_y23、n_z23分别为n₀₂₃在各坐标轴方向上的分量系数；

等效啮合节点p₃和p₂₃的坐标位置相同，根据式(17)求得：

$\left\{\begin{array}{c}\stackrel{\‾}{o_3{p}_3}-\stackrel{\‾}{o_2{p}_{23}}=L\\ \frac{\stackrel{\‾}{o_2{p}_{23}}}{\stackrel{\‾}{o_3{p}_3}}=\frac{N_2}{N_3}\end{array}\right.---\left(17\right)$

上式中，N₂为行星轮齿数；N₃为齿圈齿数；L为齿圈与行星轮的中心距。

齿圈与行星轮齿轮副的名义啮合分力如式(18)所示：

$\left\{\begin{array}{c}{F}_{t023}=2T_3/d_3\\ F_{r023}=F_{t023}{\mathrm{tan\α}}_n/cos\mathrm{\β}\\ F_{a023}=F_{t023}tan\mathrm{\β}\end{array}\right.---\left(18\right)$

上式中，F_t023为名义切向力；F_r023为名义径向力；F_a023为名义轴向力；T₃为齿圈承受的总转矩在每个行星轮分支上的转矩分量；d₃为齿圈分度圆直径；α_n为法向压力角；β为螺旋角；

如果齿轮存在变位，齿轮副的真实啮合分力如式(19)所示：

$\left\{\begin{array}{c}{F}_{tw23}=T_3/\stackrel{\‾}{o_3{p}_3}\\ F_{w23}=\sqrt{F_{t023}^2+F_{r023}^2-F_{N23}^2}\\ F_{aw23}=F_{a023}\end{array}\right.---\left(19\right)$

上式中，F_tw23为真实切向力；F_rw23为真实径向力；F_aw23为真实轴向力；

齿圈与行星轮齿轮副法向啮合力如式(20)所示：

$F_{n23}=\sqrt{F_{tw23}^2+F_{rw23}^2+F_{aw23}^2}---\left(20\right)$

齿圈和行星轮在坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₃的各分量系数如式(21)所示：

$\left\{\begin{array}{c}{n}_{x23}=-F_{rw23}/F_n\\ n_{y23}=k_L{F}_{tw23}/F_n\\ n_{z23}=-k_L{k}_R{F}_{aw23}/F_n\end{array}\right.---\left(21\right)$

上式中，工况系数k_L和齿轮旋向系数k_R的取值与式(8)相同；

全局坐标系下齿圈与行星轮齿轮副等效啮合力作用线方向矢量如式(22)所示：

$n_{23}^T=n_{023}^{T}H\left(\mathrm{\θ}\right)---\left(22\right)$

上式中，H(θ)为行星轮方位角θ对应的坐标变换矩阵，由式(10)得到。

齿圈和行星轮等效啮合节点平动自由度之间的刚度耦合如式(23)所示：

$K_{m23}=k_{m23}{\left[\begin{array}{cc}{n}_{23}^T{n}_{23}& -n_{23}^T{n}_{23}\\ -n_{23}^T{n}_{23}& n_{23}^T{n}_{23}\end{array}\right]}_{6\×6}---\left(23\right)$

上式中，k_m23为齿圈和行星轮啮合刚度系数，采用ISO6336给出的计算方法求得；

齿圈和行星轮等效啮合节点转动自由度之间的刚度耦合如式(24)所示：

$K_{m\mathrm{\θ}23}=k_{m\mathrm{\θ}23}{\left[\begin{array}{cc}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& -n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\\ -n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\end{array}\right]}_{6\×6}---\left(24\right)$

上式中，n_θ23为转动自由度耦合矢量，如式(25)所示：

n_θ23＝[-n_y23,n_x23,0]H(θ) (25)

k_mθ23为齿轮啮合弯曲刚度系数，如式(26)所示：

$\begin{array}{c}{k}_{m\mathrm{\θ}23}=\frac{{\mathrm{\ΔM}}_{23}}{{\mathrm{\Δ\γ}}_{23}}=\frac{1}{{\mathrm{\Δ\γ}}_{23}}\[\underset{b/2}{\overset{b}{\∫}}{\mathrm{\Δ\γ}}_{23}z\frac{k_{m23}}{b}\left(z-\frac{b}{2}\right)dz-\underset{0}{\overset{b/2}{\∫}}{\mathrm{\Δ\γ}}_{23}z\frac{k_{m23}}{b}\left(\frac{b}{2}-z\right)dz\]\\ =\frac{b^2}{12}{k}_{m23}\end{array}---\left(26\right)$

上式中，ΔM₂₃表示齿圈和行星轮接触线产生单位转角Δγ₂₃时，对齿宽中点的弯矩；b为有效齿宽，取齿圈和行星轮齿宽的较小值；z为齿宽方向的局部坐标；dz为齿宽方向的微变量；

式(27)为齿圈和行星轮等效啮合节点自由度对应的完整的齿轮等效啮合刚度矩阵：

$K_{m23}={\left[\begin{array}{cccc}{k}_{m23}{n}_{23}^T{n}_{23}& 0& -k_{m23}{n}_{23}^T{n}_{23}& 0\\ 0& k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& 0& -k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\\ -k_{m23}{n}_{23}^T{n}_{23}& 0& k_{m23}{n}_{23}^T{n}_{23}& 0\\ 0& -k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& 0& k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\end{array}\right]}_{12\×12}---\left(27\right)$

上式中，n₂₃为全局坐标系下的等效啮合力作用线方向矢量，如式(22)所示；n_θ23为转动自由度耦合矢量，如式(25)所示；k_m23为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得；k_mθ23为齿轮啮合弯曲刚度系数，如式(26)所示。

所述步骤5)中，行星齿轮传动系统静力学方程为：

K(δ)δ＝F (28)

上式中，K(δ)为行星齿轮传动系统的非线性刚度矩阵，由轴模型刚度矩阵K_s、行星架单元刚度矩阵K_r、滚动轴承非线性刚度矩阵K_b、齿轮刚性梁单元刚度矩阵K_g21和K_g23、齿轮等效啮合单元刚度矩阵K_m21和K_m23组集而成；δ为系统模型节点自由度的静变形矢量；F为系统外载荷矢量。

所述步骤7)中太阳轮和行星轮的错位量的求解过程如下：

设静力学计算求得的太阳轮、行星轮和齿圈的等效啮合节点绕全局坐标系X轴和Y轴的转角变形分别为和根据行星轮方位角θ，将全局坐标系下的转角变形按照式(31)变换到行星轮局部坐标系o₂x₂y₂z₂：

上式中，H(-θ)为行星轮方位角-θ对应的坐标变换矩阵；

太阳轮和行星轮各自在端面啮合力作用线方向上的错位量如式(32)所示：

$\left\{\begin{array}{c}{f}_{sh1}={\mathrm{\Φ}}_{x1}b{\mathrm{cos\α}}_{tw21}+{\mathrm{\Φ}}_{y1}b{\mathrm{sin\α}}_{tw21}\\ f_{sh21}={\mathrm{\Φ}}_{x2}b{\mathrm{cos\α}}_{tw21}+{\mathrm{\Φ}}_{y2}b{\mathrm{sin\α}}_{tw21}\end{array}\right.---\left(32\right)$

上式中，b为有效齿宽，取太阳轮和行星轮齿宽的较小值；α_tw21为考虑齿轮变位的真实端面压力角，如式(33)所示：

${\mathrm{\α}}_{tw21}=acos\left(\frac{d_1-d_2}{2L}{\mathrm{cos\α}}_{t21}\right)---\left(33\right)$

上式中，d₁为太阳轮分度圆直径；d₂为行星轮分度圆直径；L为齿轮副中心距；α_t21为分度圆端面压力角。

太阳轮和行星轮的总错位量如式(34)所示：

f_sh12＝|f_sh1-f_sh21|。 (34)

所述步骤7)中齿圈和行星轮的错位量的求解过程如下：

齿圈和行星轮各自在端面啮合力作用线方向上的错位量如式(35)所示：

$\left\{\begin{array}{c}{f}_{sh3}={\mathrm{\Φ}}_{x3}b{\mathrm{cos\α}}_{tw23}+{\mathrm{\Φ}}_{y3}b{\mathrm{sin\α}}_{tw23}\\ f_{sh23}={\mathrm{\Φ}}_{x2}b{\mathrm{cos\α}}_{tw23}+{\mathrm{\Φ}}_{y2}b{\mathrm{sin\α}}_{tw23}\end{array}\right.---\left(35\right)$

上式中，b为有效齿宽，取齿圈和行星轮齿宽的较小值；α_tw23为考虑齿轮变位的真实端面压力角，如式(36)所示：

${\mathrm{\α}}_{tw23}=acos\left(\frac{d_3-d_2}{2L}{\mathrm{cos\α}}_{t23}\right)---\left(36\right)$

上式中，d₂为行星轮分度圆直径；d₃为齿圈分度圆直径；L为齿轮副中心距；α_t23为分度圆端面压力角。

齿圈和行星轮的总错位量如式(37)所示：

f_sh32＝|f_sh3-f_sh23|。 (37)

本发明由于采取以上技术方案，其具有以下优点：1、本发明方法能够全面准确地体现行星齿轮传动系统的结构特征，采用的考虑轴承刚度耦合性和非线性的轴承单元能够准确体现滚动轴承刚度非线性特性对行星齿轮错位量的影响。2、本发明能够准确计算各行星轮分支齿轮副在任意工况和任意行星轮方位下的齿轮错位量。3、本发明基于经典非线性轴承理论和有限元方法，具有可靠的理论基础，采用空间梁单元模型和齿轮等效啮合模型实现行星齿轮传动系统建模，易于在各类常用的编程语言环境下编程实现，采用牛顿-拉弗森方法实现快速求解，具有较高的计算效率。本发明可广泛应用于行星齿轮传动系统和同类齿轮传动系统的设计建模和计算分析。

附图说明

图1是本发明方法的流程示意图；

图2是典型行星齿轮传动系统的机构简图；

图3是行星架模型的空间结构示意图；

图4是行星齿轮啮合模型平面示意图；

图5是齿轮错位量示意图；

图6是行星齿轮传动系统平面示意图；

图7是行星齿轮初始方位平面示意图；

图8是输入轴平面示意图；

图9是行星架轴平面示意图；

图10是齿圈轴平面示意图；

图11是行星轮销轴平面示意图；

图12是行星轮轴平面示意图；

图13是各行星轮分支齿轮错位量随行星轮方位角变化曲线；

图14是齿轮错位量随输入转矩大小的变化曲线。

其中，图2和图5中：1、输入轴；2、行星架轴；3、齿圈轴；4、行星轮销轴；5、行星轮轴；6、行星架单元；7、滚子轴承；8、太阳轮轮齿；9、行星轮轮齿；10、齿圈轮齿；11、输入端；12、输出端。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

如图1所示，本发明提供的获取行星齿轮错位量的有限元方法包括以下步骤：

1)建立轴部件的有限元模型：采用欧拉-伯努利空间梁单元建立轴部件的有限元模型，如图2所示，行星齿轮传动系统中的轴部件包括：输入轴1、行星架轴2、齿圈轴3、行星轮销轴4、行星轮轴5，其中，行星轮销轴4和行星轮轴5的数量均等于行星轮的个数N_p，梁单元的几何参数和材料参数定义为轴部件的真实设计参数，采用有限元方法计算轴部件的刚度矩阵K_s。

2)建立行星架的有限元模型：行星齿轮传动系统中的每个行星轮销轴4都需要通过行星架单元6与行星架轴2连接，从而使动力按照行星齿轮的传动比关系传递给行星架轴2，行星架单元6的空间结构如图3中的虚线所示，行星架单元6采用刚性梁单元模拟，刚度矩阵表示为K_r。

3)建立滚动轴承的有限元模型：行星轮轴5与行星轮销轴4之间通常采用滚动轴承7连接，本发明采用具有非线性耦合刚度特性的轴承单元模拟滚动轴承，滚动轴承刚的度矩阵K_b如式(2)所示：

$K_b=\left[\begin{array}{cccccc}\frac{\∂F_{bx}}{\∂{\mathrm{\δ}}_{bx}}& \frac{\∂F_{bx}}{\∂{\mathrm{\δ}}_{by}}& \frac{\∂F_{bx}}{\∂{\mathrm{\δ}}_{bz}}& \frac{\∂F_{bx}}{\∂{\mathrm{\θ}}_{bx}}& \frac{\∂F_{bx}}{\∂{\mathrm{\θ}}_{by}}& 0\\ \frac{\∂F_{by}}{\∂{\mathrm{\δ}}_{bx}}& \frac{\∂F_{by}}{\∂{\mathrm{\δ}}_{by}}& \frac{\∂F_{by}}{\∂{\mathrm{\δ}}_{bz}}& \frac{\∂F_{by}}{\∂{\mathrm{\θ}}_{bx}}& \frac{\∂F_{by}}{\∂{\mathrm{\θ}}_{by}}& 0\\ \frac{\∂F_{bz}}{\∂{\mathrm{\δ}}_{bx}}& \frac{\∂F_{bz}}{\∂{\mathrm{\δ}}_{by}}& \frac{\∂F_{bz}}{\∂{\mathrm{\δ}}_{bz}}& \frac{\∂F_{bz}}{\∂{\mathrm{\θ}}_{bx}}& \frac{\∂F_{bz}}{\∂{\mathrm{\θ}}_{by}}& 0\\ \frac{\∂M_{bx}}{\∂{\mathrm{\δ}}_{bx}}& \frac{\∂M_{bx}}{\∂{\mathrm{\δ}}_{by}}& \frac{\∂M_{bx}}{\∂{\mathrm{\δ}}_{bz}}& \frac{\∂M_{bx}}{\∂{\mathrm{\θ}}_{bx}}& \frac{\∂M_{bx}}{\∂{\mathrm{\θ}}_{by}}& 0\\ \frac{\∂M_{by}}{\∂{\mathrm{\δ}}_{bx}}& \frac{\∂M_{by}}{\∂{\mathrm{\δ}}_{by}}& \frac{\∂M_{by}}{\∂{\mathrm{\δ}}_{bz}}& \frac{\∂M_{by}}{\∂{\mathrm{\θ}}_{bx}}& \frac{\∂M_{by}}{\∂{\mathrm{\θ}}_{by}}& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]---\left(2\right)$

上式中，F_bx和F_by为轴承内外圈之间传递的径向力；F_bz为轴承所受轴向力；M_bx和M_by为轴承内外圈之间传递的径向弯矩；δ_bx和δ_by为轴承内外圈中心的相对径向变形；δ_bz为轴承内外圈中心的相对轴向变形；θ_bx和θ_by为轴承内外圈中心相对径向角变形。

4)建立行星齿轮有限元模型，包括太阳轮与行星轮的等效啮合模型以及齿圈与行星轮的等效啮合模型，建模过程如下：

①太阳轮与行星轮为外啮合齿轮副(如图4所示)，OXYZ为系统全局坐标系，Z轴正方向与X、Y轴满足右手定则，坐标原点O定义为太阳轮中心。o₂x₂y₂z₂为行星轮局部坐标系，行星轮方位角θ定义为x₂与X的夹角。o₁为太阳轮中心梁单元节点，o₂为行星轮中心梁单元节点，p₁为太阳轮等效啮合节点，p₂₁为行星轮在太阳轮一端的等效啮合节点，p₁和p₂₁均有6个自由度。o₁与p₁、o₂与p₂₁之间采用刚性梁单元连接，刚度矩阵表示为K_g21，p₁和p₂₁之间采用沿齿轮等效啮合力作用线方向的空间弹簧单元连接，坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₁表示为式(3)：

n₀₂₁＝[n_x21,n_y21,n_z21]^T (3)

上式中，n_x21、n_y21、n_z21分别为n₀₂₁在各坐标轴方向上的分量系数。

等效啮合节点p₁和p₂₁的坐标位置相同，可根据式(4)求得：

$\left\{\begin{array}{c}\stackrel{\‾}{o_1{p}_1}=\frac{N_1}{N_1+N_2}L\\ \stackrel{\‾}{o_2{p}_{21}}=\frac{N_2}{N_1+N_2}L\end{array}\right.---\left(4\right)$

上式中，N₁为太阳轮齿数；N₂为行星轮齿数；L为太阳轮与行星轮的中心距。

太阳轮与行星轮齿轮副的名义啮合分力如式(5)所示：

$\left\{\begin{array}{c}{F}_{t021}=2T_1/d_1\\ F_{r021}=F_{t021}{\mathrm{tan\α}}_n/cos\mathrm{\β}\\ F_{a021}=F_{t021}tan\mathrm{\β}\end{array}\right.---\left(5\right)$

上式中，F_t021为名义切向力；F_r021为名义径向力；F_a021为名义轴向力；T₁为行星轮分支输入转矩大小，设太阳轮总输入转矩T平均分配给各行星轮分支，则T₁＝T/N_p；d₁为太阳轮分度圆直径；α_n为法向压力角；β为螺旋角。

如果齿轮存在变位，齿轮副的真实啮合分力如式(6)所示：

$\left\{\begin{array}{c}{F}_{tw21}=T_1/\stackrel{\‾}{o_1{p}_1}\\ F_{rw21}=\sqrt{F_{t021}^2+F_{r021}^2-F_{tw21}^2}\\ F_{aw21}=F_{a021}\end{array}\right.---\left(6\right)$

上式中，F_tw21为真实切向力；F_rw21为真实径向力；F_aw21为真实轴向力

太阳轮与行星轮齿轮副的法向啮合力如式(7)所示：

$F_{n21}=\sqrt{F_{tw21}^2+F_{rw21}^2+F_{aw21}^2}---\left(7\right)$

图4太阳轮和行星轮在坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₁(如式(3)所示)的各分量系数如式(8)所示：

$\left\{\begin{array}{c}{n}_{x21}=-F_{rw21}/F_{n21}\\ n_{y21}=k_L{F}_{tw21}/F_{n21}\\ n_{z21}=-k_L{k}_R{F}_{aw21}/F_{n21}\end{array}\right.---\left(8\right)$

上式中，k_L为工况系数，当输入转矩T的方向为+Z时，k_L＝1，当输入转矩T的方向为-Z时，k_L＝-1；k_R为齿轮旋向系数，当太阳轮为右旋时，k_R＝1，当太阳轮为左旋时，k_R＝-1。

全局坐标系下太阳轮与行星轮齿轮副等效啮合力作用线方向矢量如式(9)所示：

$n_{21}^T=n_{021}^{T}H\left(\mathrm{\θ}\right)---\left(9\right)$

上式中，H(θ)为行星轮方位角θ对应的坐标变换矩阵，其可由式(10)得到：

$H\left(\mathrm{\θ}\right)=\left[\begin{array}{ccc}cos\mathrm{\θ}& -\sin \mathrm{\θ}& 0\\ sin\mathrm{\θ}& \cos \mathrm{\θ}& 0\\ 0& 0& 1\end{array}\right]---\left(10\right)$

太阳轮和行星轮等效啮合节点平动自由度之间的刚度耦合如式(11)所示：

$K_{m21}=k_{m21}{\left[\begin{array}{cc}{n}_{21}^T{n}_{21}& -n_{21}^T{n}_{21}\\ -n_{21}^T{n}_{21}& n_{21}^T{n}_{21}\end{array}\right]}_{6\×6}---\left(11\right)$

上式中，k_m21为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得。

太阳轮和行星轮等效啮合节点转动自由度之间的刚度耦合如式(12)所示：

$K_{m\mathrm{\θ}21}=k_{m\mathrm{\θ}21}{\left[\begin{array}{cc}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& -n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\\ -n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\end{array}\right]}_{6\×6}---\left(12\right)$

上式中，n_θ21为转动自由度耦合矢量，如式(13)所示：

n_θ21＝[-n_y21,n_x21,0]H(θ)； (13)

k_mθ21为太阳轮和行星轮啮合弯曲刚度系数，如式(14)所示：

$\begin{array}{c}{k}_{m\mathrm{\θ}21}=\frac{{\mathrm{\ΔM}}_{21}}{{\mathrm{\Δ\γ}}_{21}}=\frac{1}{{\mathrm{\Δ\γ}}_{21}}\[\underset{b/2}{\overset{b}{\∫}}{\mathrm{\Δ\γ}}_{21}z\frac{k_{m21}}{b}\left(z-\frac{b}{2}\right)dz-\underset{0}{\overset{b/2}{\∫}}{\mathrm{\Δ\γ}}_{21}z\frac{k_{m21}}{b}\left(\frac{b}{2}-z\right)dz\]\\ =\frac{b^2}{12}{k}_{m21}\end{array}---\left(14\right)$

上式中，ΔM₂₁表示太阳轮和行星轮接触线产生单位转角Δγ₂₁时，对齿宽中点的弯矩；b为有效齿宽，取太阳轮和行星轮齿宽的较小值；z为齿宽方向的局部坐标；dz为齿宽方向的微变量。

式(15)为太阳轮和行星轮等效啮合节点自由度对应的完整的齿轮等效啮合刚度矩阵：

$K_{m21}={\left[\begin{array}{cccc}{k}_{m21}{n}_{21}^T{n}_{21}& 0& -k_{m21}{n}_{21}^T{n}_{21}& 0\\ 0& k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& 0& -k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\\ -k_{m21}{n}_{21}^T{n}_{21}& 0& k_{m21}{n}_{21}^T{n}_{21}& 0\\ 0& -k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& 0& k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\end{array}\right]}_{12\×12}---\left(15\right)$

上式中，n₂₁为全局坐标系下的等效啮合力作用线方向矢量，如式(9)所示；n_θ21为转动自由度耦合矢量，如式(13)所示；k_m21为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得；k_mθ21为齿轮啮合弯曲刚度系数，如式(14)所示。

②齿圈与行星轮为内啮合齿轮副(如图4所示)，o₃为齿圈中心梁单元节点，其坐标位置与全局坐标系原点O、太阳轮中心梁单元节点o₁位置相同，o₂为行星轮中心梁单元节点，p₃为齿圈等效啮合节点，p₂₃为行星轮在齿圈一端的等效啮合节点，p₃和p₂₃均有6个自由度。o₃与p₃、o₂与p₂₃之间采用刚性梁单元连接，刚度矩阵表示为K_g23，p₃和p₂₃之间采用沿齿轮等效啮合力作用线方向的空间弹簧单元连接，坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₃表示为式(16)：

n₀₂₃＝[n_x23,n_y23,n_z23]^T (16)

上式中，n_x23、n_y23、n_z23分别为n₀₂₃在各坐标轴方向上的分量系数。

等效啮合节点p₃和p₂₃的坐标位置相同，可根据式(17)求得：

$\left\{\begin{array}{c}\stackrel{\‾}{o_3{p}_3}-\stackrel{\‾}{o_2{p}_{23}}=L\\ \frac{\stackrel{\‾}{o_2{p}_{23}}}{\stackrel{\‾}{o_3{p}_3}}=\frac{N_2}{N_3}\end{array}\right.---\left(17\right)$

上式中，N₂为行星轮齿数；N₃为齿圈齿数；L为齿圈与行星轮的中心距。

齿圈与行星轮齿轮副的名义啮合分力如式(18)所示：

$\left\{\begin{array}{c}{F}_{t023}=2T_3/d_3\\ F_{r023}=F_{t023}{\mathrm{tan\α}}_n/cos\mathrm{\β}\\ F_{a023}=F_{t023}tan\mathrm{\β}\end{array}\right.---\left(18\right)$

上式中，F_t023为名义切向力；F_r023为名义径向力；F_a023为名义轴向力；T₃为齿圈承受的总转矩在每个行星轮分支上的转矩分量；d₃为齿圈分度圆直径；α_n为法向压力角；β为螺旋角。

如果齿轮存在变位，齿轮副的真实啮合分力如式(19)所示：

$\left\{\begin{array}{c}{F}_{tw23}=T_3/\stackrel{\‾}{o_3{p}_3}\\ F_{w23}=\sqrt{F_{t023}^2+F_{r023}^2-F_{N23}^2}\\ F_{aw23}=F_{a023}\end{array}\right.---\left(19\right)$

上式中，F_tw23为真实切向力；F_rw23为真实径向力；F_aw23为真实轴向力

齿圈与行星轮齿轮副法向啮合力如式(20)所示：

$F_{n23}=\sqrt{F_{tw23}^2+F_{rw23}^2+F_{aw23}^2}---\left(20\right)$

图4齿圈和行星轮在坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₃的各分量系数如式(21)所示：

$\left\{\begin{array}{c}{n}_{x23}=-F_{rw23}/F_n\\ n_{y23}=k_L{F}_{tw23}/F_n\\ n_{z23}=-k_L{k}_R{F}_{aw23}/F_n\end{array}\right.---\left(21\right)$

上式中，工况系数k_L和齿轮旋向系数k_R的取值与式(8)相同。

全局坐标系下齿圈与行星轮齿轮副等效啮合力作用线方向矢量如式(22)所示：

$n_{23}^T=n_{023}^{T}H\left(\mathrm{\θ}\right)---\left(22\right)$

上式中，H(θ)为行星轮方位角θ对应的坐标变换矩阵，可由式(10)得到。

齿圈和行星轮等效啮合节点平动自由度之间的刚度耦合如式(23)所示：

$K_{m23}=k_{m23}{\left[\begin{array}{cc}{n}_{23}^T{n}_{23}& -n_{23}^T{n}_{23}\\ -n_{23}^T{n}_{23}& n_{23}^T{n}_{23}\end{array}\right]}_{6\×6}---\left(23\right)$

上式中，k_m23为齿圈和行星轮啮合刚度系数，采用ISO6336给出的计算方法求得。

齿圈和行星轮等效啮合节点转动自由度之间的刚度耦合如式(24)所示：

$K_{m\mathrm{\θ}23}=k_{m\mathrm{\θ}23}{\left[\begin{array}{cc}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& -n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\\ -n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\end{array}\right]}_{6\×6}---\left(24\right)$

上式中，n_θ23为转动自由度耦合矢量，如式(25)所示：

n_θ23＝[-n_y23,n_x23,0]H(θ) (25)

k_mθ23为齿轮啮合弯曲刚度系数，如式(26)所示：

$\begin{array}{c}{k}_{m\mathrm{\θ}23}=\frac{{\mathrm{\ΔM}}_{23}}{{\mathrm{\Δ\γ}}_{23}}=\frac{1}{{\mathrm{\Δ\γ}}_{23}}\[\underset{b/2}{\overset{b}{\∫}}{\mathrm{\Δ\γ}}_{23}z\frac{k_{m23}}{b}\left(z-\frac{b}{2}\right)dz-\underset{0}{\overset{b/2}{\∫}}{\mathrm{\Δ\γ}}_{23}z\frac{k_{m23}}{b}\left(\frac{b}{2}-z\right)dz\]\\ =\frac{b^2}{12}{k}_{m23}\end{array}---\left(26\right)$

上式中，ΔM₂₃表示齿圈和行星轮接触线产生单位转角Δγ₂₃时，对齿宽中点的弯矩；b为有效齿宽，取齿圈和行星轮齿宽的较小值；z为齿宽方向的局部坐标；dz为齿宽方向的微变量。

式(27)为齿圈和行星轮等效啮合节点自由度对应的完整的齿轮等效啮合刚度矩阵：

$K_{m23}={\left[\begin{array}{cccc}{k}_{m23}{n}_{23}^T{n}_{23}& 0& -k_{m23}{n}_{23}^T{n}_{23}& 0\\ 0& k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& 0& -k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\\ -k_{m23}{n}_{23}^T{n}_{23}& 0& k_{m23}{n}_{23}^T{n}_{23}& 0\\ 0& -k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& 0& k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\end{array}\right]}_{12\×12}---\left(27\right)$

上式中，n₂₃为全局坐标系下的等效啮合力作用线方向矢量，如式(22)所示；n_θ23为转动自由度耦合矢量，如式(25)所示；k_m23为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得；k_mθ23为齿轮啮合弯曲刚度系数，如式(26)所示。

5)建立行星齿轮传动系统静力学模型：根据轴部件、滚动轴承、行星架、齿轮之间的连接关系，采用有限元方法组集各部件的刚度矩阵，建立完整的行星齿轮传动系统静力学模型，系统静力学方程如式(28)示：

K(δ)δ＝F (28)

上式中，K(δ)为行星齿轮传动系统的非线性刚度矩阵，由轴模型刚度矩阵K_s、行星架单元刚度矩阵K_r、滚动轴承非线性刚度矩阵K_b、齿轮刚性梁单元刚度矩阵K_g21和K_g23、齿轮等效啮合单元刚度矩阵K_m21和K_m23组集而成；δ为系统模型节点自由度的静变形矢量；F为系统外载荷矢量，根据实际工况施加外载荷和边界条件。

6)传动系统静力学求解：采用牛顿-拉弗森方法求解系统非线性静力学方程，式(29)为迭代求解过程：

δ_k＝δ_k-1-K(δ_k-1)-¹[F_k-1-F] (29)

上式中，δ_k-1和δ_k分别表示第k-1次和第k次迭代后的系统节点自由度位移矢量；K(δ_k-1)为第k-1次迭代后系统节点自由度位移矢量δ_k-1对应的系统切线刚度矩阵；F_k-1为第k-1次迭代后的系统载荷矢量。

当相邻两次迭代的位移矢量结果之差的模满足式(30)时，计算收敛。

||δ_k-δ_k-1||＜ε (30)

上式中，ε为收敛容差，是一个很小的正数。

7)计算行星齿轮错位量，包括太阳轮和行星轮的错位量以及齿圈和行星轮的错位量：

如图5所示，当主动齿轮和从动齿轮接触线一端重合时，另一端沿端面啮合力作用线矢量方向n_t的位移即为齿轮错位量f_sh。

设静力学计算求得的太阳轮、行星轮和齿圈的等效啮合节点绕全局坐标系X轴和Y轴的转角变形分别为和根据图4所示的行星轮方位角θ，将全局坐标系下的转角变形按照式(31)变换到行星轮局部坐标系o₂x₂y₂z₂：

上式中，H(-θ)为行星轮方位角-θ对应的坐标变换矩阵，可由式(10)得到。

①太阳轮和行星轮各自在端面啮合力作用线方向上的错位量如式(32)所示：

$\left\{\begin{array}{c}{f}_{sh1}={\mathrm{\Φ}}_{x1}b{\mathrm{cos\α}}_{tw21}+{\mathrm{\Φ}}_{y1}b{\mathrm{sin\α}}_{tw21}\\ f_{sh21}={\mathrm{\Φ}}_{x2}b{\mathrm{cos\α}}_{tw21}+{\mathrm{\Φ}}_{y2}b{\mathrm{sin\α}}_{tw21}\end{array}\right.---\left(32\right)$

上式中，b为有效齿宽，取太阳轮和行星轮齿宽的较小值；α_tw21为考虑齿轮变位的真实端面压力角，如式(33)所示：

${\mathrm{\α}}_{tw21}=acos\left(\frac{d_1+d_2}{2L}{\mathrm{cos\α}}_{t21}\right)---\left(33\right)$

上式中，d₁为太阳轮分度圆直径；d₂为行星轮分度圆直径；L为齿轮副中心距；α_t21为分度圆端面压力角。

太阳轮和行星轮的总错位量如式(34)所示：

f_sh12＝|f_sh1-f_sh21| (34)

②齿圈和行星轮各自在端面啮合力作用线方向上的错位量如式(35)所示：

$\left\{\begin{array}{c}{f}_{sh3}={\mathrm{\Φ}}_{x3}b{\mathrm{cos\α}}_{tw23}+{\mathrm{\Φ}}_{y3}b{\mathrm{sin\α}}_{tw23}\\ f_{sh23}={\mathrm{\Φ}}_{x2}b{\mathrm{cos\α}}_{tw23}+{\mathrm{\Φ}}_{y2}b{\mathrm{sin\α}}_{tw23}\end{array}\right.---\left(35\right)$

上式中，b为有效齿宽，取齿圈和行星轮齿宽的较小值；α_tw23为考虑齿轮变位的真实端面压力角，如式(36)所示：

${\mathrm{\α}}_{tw23}=acos\left(\frac{d_3-d_2}{2L}{\mathrm{cos\α}}_{t23}\right)---\left(36\right)$

上式中，d₂为行星轮分度圆直径；d₃为齿圈分度圆直径；L为齿轮副中心距；α_t23为分度圆端面压力角。

齿圈和行星轮的总错位量如式(37)所示：

f_sh32＝|f_sh3-f_sh23| (37)

下面通过一个具体的实施例，用以说明本发明的效果。

以如图6所示的行星齿轮传动系统为例，行星轮初始方位如图7所示，行星轮个数N_p＝4，全局坐标系原点O定义为太阳轮中心，全局坐标系Z轴与输入轴平行，Z轴正方向与X、Y轴满足右手定则。

1)建立轴部件模型：依次建立图6轴部件的梁单元有限元模型，在与其他部件存在连接关系的位置、截面尺寸发生变化的位置建立梁单元节点，具体如下：

图8所示的输入轴有限元模型，包含3个梁单元和4个节点，其中2号节点为太阳轮中心节点，4号节点为输入转矩加载点。

图9所示的行星架轴有限元模型，包含5个梁单元和6个节点，其中4号节点为左侧行星架单元连接节点，5号节点为右侧行星架单元连接节点。

图10所示的齿圈轴有限元模型，包含5个梁单元和6个节点，其中2号节点为齿圈中心节点。

图11所示的行星轮销轴有限元模型，包含5个梁单元和6个节点，其中2号节点为左侧行星架单元连接节点，5号节点为右侧行星架单元连接节点，3号节点为左侧圆锥滚子轴承连接节点，4号节点为右侧圆锥滚子轴承连接节点。共建立N_p＝4个行星轮销轴模型。

图12所示的行星轮轴有限元模型，包含6个梁单元和7个节点，其中2号节点为左侧圆锥滚子轴承连接节点，6号节点为右侧圆锥滚子轴承连接节点，4号节点为行星轮中心。共建立N_p＝4个行星轮销轴模型。

梁单元的横截面形状均为圆或圆环，材料参数取钢材的参数，弹性模量为210GPa，泊松比为0.3，求得各轴部件的梁单元有限元模型的刚度矩阵K_s。

2)建立行星架模型：如图3和图6所示，每个行星轮销轴4均通过左、右两组行星架单元6与行星架轴2连接，则本实施例共需建立2N_p＝8个行星架单元，均建立为刚性梁单元，长度均为太阳轮轴与行星轮轴的中心距L＝64mm，求得行星架单元的刚度矩阵K_r。

3)建立滚动轴承模型：如图6所示，每个行星轮由一对朝向相反的圆锥滚子轴承7支撑，则本实施例中的行星齿轮传动系统共包含2N_p＝8个圆锥滚子轴承。圆锥滚子轴承的型号均为FAG30305A，轴承内径为25mm，外径为62mm，宽度为18.25mm，平均直径为43.5mm，滚子数为12，滚子直径为9.25mm，滚子有效长度为12.22mm，根据滚子轴承的非线性刚度计算公式，求得每个轴承各自的非线性刚度矩阵K_b。

4)建立行星齿轮模型：

本实施例中的行星齿轮参数如表1所示。

表1行星齿轮参数

行星轮的初始方位如图7所示，求得各分支行星轮中心在全局坐标系中的坐标分别为：(64,0,0)、(0,64,0)、(-64,0,0)、(0,-64,0)，太阳轮中心o₁、齿圈中心o₃与全局坐标系原点O重合，坐标均为(0,0,0)。

求得全局坐标系太阳轮与各分支行星轮等效啮合节点坐标分别为：(21.029,0,0)、(0,21.029,0)、(-21.029,0,0)、(0,-21.029,0)；齿圈与各分支行星轮的等效啮合节点坐标分别为：(106.971,0,0)、(0,106.971,0)、(-106.971,0,0)、(0,-106.971,0)。

求得全局坐标系太阳轮与各分支行星轮的等效啮合力作用线方向矢量分别为：[0.37357,0.89932,-0.22733]^T、[-0.89932,0.37357,-0.22733]^T、[-0.37357,-0.89932,-0.22733]^T、[0.89932,-0.37357,-0.22733]^T；齿圈与各分支行星轮的等效啮合力作用线方向矢量分别为：[-0.37357,0.89932,0.22733]^T、[-0.89932,-0.37357,0.22733]^T、[0.37357,-0.89932,0.22733]^T、[0.89932,0.37357,0.22733]^T。

采用ISO6336计算方法，求得太阳轮与行星轮之间的等效啮合刚度系数为748520N/mm，齿圈与行星轮之间的等效啮合刚度系数为890050N/mm。

进一步求得齿轮刚性梁单元刚度矩阵K_g21、K_g23和齿轮等效啮合单元刚度矩阵K_m21、K_m23。

5)建立行星齿轮传动系统静力学模型：采用有限元方法组集各部件的刚度矩阵，包括：轴模型刚度矩阵K_s、行星架单元刚度矩阵K_r、滚动轴承非线性刚度矩阵K_b、齿轮刚性梁单元刚度矩阵K_g21和K_g23、齿轮等效啮合单元刚度矩阵K_m21和K_m23，获得行星齿轮传动系统的非线性刚度矩阵K(δ)，所建立的行星齿轮传动系统静力学模型共包含84个节点，504个自由度。

如图6所示，输入转矩11施加在输入轴1上，对应图8中的4号节点，输入转矩T＝1000N·m；在行星架轴2上施加集中力F_X＝5000N，对应图9中的2号节点。

约束图9行星架轴1号节点的轴向转动自由度，并约束图10齿圈轴6号节点的所有自由度，以消除行星齿轮传动系统的轴向转动刚体位移，并实现齿圈固定；约束图9行星架轴1号节点、图8输入轴4号节点的平动自由度，以消除刚体位移。

6)传动系统静力学求解：采用牛顿-拉弗森方法迭代求解系统静力学方程，以相邻两次迭代的位移矢量结果之差的模小于收敛容差作为收敛准则，收敛容差取为10^-4mm，迭代8次计算收敛，耗时6秒，为了验证传动系统静力学求解的正确性，检查系统的载荷传递。

求得图9行星架轴的1号节点的轴向转动约束产生的支反转矩为6086.9Nm，该转矩即为输入转矩1000Nm在系统静平衡状态下对应的输出转矩，满足式(38)所示的行星齿轮理论传动比关系：

$i=1+\frac{N_3}{N_1}=1+\frac{117}{23}=6.087---\left(38\right)$

求得图10齿圈轴的5号节点的轴向转动约束产生的支反转矩为5086.9Nm，同样满足式(39)所示的齿圈与太阳轮之间的理论传动比关系：

$i^{\′}=\frac{N_3}{N_1}=\frac{117}{23}=5.087---\left(39\right)$

由上述结果可知，行星齿轮传动系统静力求解正确。

7)计算行星齿轮错位量：各行星轮分支齿轮错位量计算结果如表2和表3所示，由于在行星架轴6上施加了集中力F_X＝5000N，行星齿轮系统承受非轴对称载荷，导致各行星轮分支的齿轮错位量不同。

表2太阳轮-行星轮错位量

行星轮分支	太阳轮错位量/μm	行星轮错位量/μm	总错位量/μm
				1	0.15	-14.15	14.00
2	-6.57	-17.79	11.22
				3	-0.15	-15.51	15.66
4	6.57	-11.91	18.48

表3齿圈-行星轮错位量

行星轮分支	齿圈错位量/μm	行星轮错位量/μm	总错位量/μm
				1	0.06	-11.83	11.89
2	0.11	-10.89	11.00
				3	-0.06	-11.83	11.77
4	-0.11	-12.73	12.62

在施加T＝1000N·m和F_X＝5000N的工况下，连续改变图7所示的行星轮方位，使各行星轮分支绕+Z方向转动角度θ，θ由0度逐渐变化到360度，绘制各行星轮分支齿轮错位量f_sh随行星轮方位角θ变化曲线，如图13所示，各分支错位量均随方位角θ呈周期性变化。此外，由于本实施例中的齿圈固定，且齿圈与行星轮的啮合刚度比太阳轮与行星轮的啮合刚度大，图13所示的齿圈-行星轮错位量的幅值和变化量均要比太阳轮-行星轮小。

如果仅施加输入转矩T＝1000N·m，求得的齿轮错位量如表4和表5所示，因为系统承受轴对称载荷，所以各行星轮分支的受载和变形均相同。此外，该工况下输入轴和齿圈轴没有发生弯曲变形，所以太阳轮和齿圈错位量均为零。

表4太阳轮-行星轮错位量

行星轮分支	太阳轮错位量/μm	行星轮错位量/μm	总错位量/μm
				1	0	-14.83	14.83
2	0	-14.83	14.83
				3	0	-14.83	14.83
4	0	-14.83	14.83

表5齿圈-行星轮错位量

行星轮分支	齿圈错位量/μm	行星轮错位量/μm	总错位量/μm
				1	0	-11.84	11.84
2	0	-11.84	11.84
				3	0	-11.84	11.84
4	0	-11.84	11.84

在仅施加输入转矩的工况下，各行星轮分支的错位量相同，研究不同转矩大小对行星齿轮错位量的影响，如图14所示，由于滚动轴承具有非线性刚度特性，行星齿轮错位量随转矩大小的变化曲线呈现出非线性特性。

综上所述，本发明提出的建模和计算方法能够全面准确地模拟行星齿轮传动系统的结构特征，准确计算各行星轮分支齿轮副在任意工况和任意行星轮方位下的齿轮错位量，克服了ISO6336计算方法的不足，且具有较高的计算效率。本发明可广泛应用于行星齿轮传动系统和同类齿轮传动系统的设计建模和计算分析。

上述各实施例仅用于对本发明的目的、技术方案和有益效果进行了进一步详细说明，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种获取行星齿轮错位量的有限元方法，包括以下步骤：1)建立轴部件的有限元模型：采用欧拉-伯努利空间梁单元建立轴部件的有限元模型，获取轴部件的刚度矩阵；2)建立行星架的有限元模型：采用刚性梁单元对行星架单元进行模拟，获得行星架的刚度矩阵；3)建立滚动轴承的有限元模型：采用具有非线性耦合刚度特性的轴承单元对滚动轴承进行模拟，获得滚动轴承的刚度矩阵；4)建立行星齿轮有限元模型：分别建立太阳轮与行星轮的等效啮合模型和齿圈与行星轮的等效啮合模型，获得太阳轮与行星轮的等效啮合刚度矩阵和齿圈与行星轮的等效啮合刚度矩阵；5)建立行星齿轮传动系统静力学模型：根据轴部件、滚动轴承、行星架、齿轮之间的连接关系，采用有限元方法组集各部件的刚度矩阵，建立完整的行星齿轮传动系统静力学模型；6)传动系统静力学求解：采用牛顿-拉弗森方法求解系统非线性静力学方程；7)计算行星齿轮错位量：即太阳轮和行星轮的错位量以及齿圈和行星轮的错位量。

2.如权利要求1所述的一种获取行星齿轮错位量的有限元方法，其特征在于：所述步骤1)中的轴部件包括输入轴、行星架轴、齿圈轴、行星轮销轴和行星轮轴，其中，所述行星轮销轴和行星轮轴的数量均等于行星轮的个数。

3.如权利要求1所述的一种获取行星齿轮错位量的有限元方法，其特征在于：所述步骤4)中建立太阳轮与行星轮的等效啮合模型并获得太阳轮与行星轮的等效啮合刚度矩阵的过程如下：

太阳轮与行星轮为外啮合齿轮副，建立系统全局坐标系OXYZ和行星轮局部坐标系o₂x₂y₂z₂，其中，系统全局坐标系OXYZ的Z轴正方向与X、Y轴满足右手定则，坐标原点O定义为太阳轮中心；定义行星轮方位角θ为x₂与X的夹角，太阳轮中心梁单元节点为o₁，行星轮中心梁单元节点为o₂，太阳轮等效啮合节点为p₁，行星轮在太阳轮一端的等效啮合节点为p₂₁，p₁和p₂₁均有6个自由度；o₁与p₁、o₂与p₂₁之间采用刚性梁单元连接，刚度矩阵表示为K_g21，p₁和p₂₁之间采用沿齿轮等效啮合力作用线方向的空间弹簧单元连接，坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₁表示为式(3)：

n₀₂₁＝[n_x21,n_y21,n_z21]^T (3)

上式中，n_x21、n_y21、n_z21分别为n₀₂₁在各坐标轴方向上的分量系数；

等效啮合节点p₁和p₂₁的坐标位置相同，根据式(4)求得：

$\left\{\begin{array}{c}\stackrel{\‾}{o_1{p}_1}=\frac{N_1}{N_1+N_2}L\\ \stackrel{\‾}{o_2{p}_{21}}=\frac{N_2}{N_1+N_2}L\end{array}\right.---\left(4\right)$

上式中，N₁为太阳轮齿数；N₂为行星轮齿数；L为太阳轮与行星轮的中心距；

太阳轮与行星轮齿轮副的名义啮合分力如式(5)所示：

$\left\{\begin{array}{c}{F}_{t021}=2T_1/d_1\\ F_{r021}=F_{t021}{\mathrm{tan\α}}_n/cos\mathrm{\β}\\ F_{a021}=F_{t021}tan\mathrm{\β}\end{array}\right.---\left(5\right)$

上式中，F_t021为名义切向力；F_r021为名义径向力；F_a021为名义轴向力；T₁为行星轮分支输入转矩大小，设太阳轮总输入转矩T平均分配给各行星轮分支，则T₁＝T/N_p；d₁为太阳轮分度圆直径；α_n为法向压力角；β为螺旋角；

如果齿轮存在变位，齿轮副的真实啮合分力如式(6)所示：

$\left\{\begin{array}{c}{F}_{tw21}=T_1/\stackrel{\‾}{o_1{p}_1}\\ F_{rw21}=\sqrt{F_{t021}^2+F_{r021}^2-F_{tw21}^2}\\ F_{aw21}=F_{a021}\end{array}\right.---\left(6\right)$

上式中，F_tw21为真实切向力；F_rw21为真实径向力；F_aw21为真实轴向力；

太阳轮与行星轮齿轮副的法向啮合力如式(7)所示：

$F_{n21}=\sqrt{F_{tw21}^2+F_{rw21}^2+F_{aw21}^2}---\left(7\right)$

太阳轮和行星轮在坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₁的各分量系数如式(8)所示：

$\left\{\begin{array}{c}{n}_{x21}=F_{rw21}/F_{n21}\\ n_{y21}=k_L{F}_{tw21}/F_{n21}\\ n_{z21}=k_L{k}_R{F}_{aw21}/F_{n21}\end{array}\right.---\left(8\right)$

上式中，k_L为工况系数，当输入转矩T的方向为+Z时，k_L＝1，当输入转矩T的方向为-Z时，k_L＝-1；k_R为齿轮旋向系数，当太阳轮为右旋时，k_R＝1，当太阳轮为左旋时，k_R＝-1；

全局坐标系下太阳轮与行星轮齿轮副等效啮合力作用线方向矢量如式(9)所示：

$n_{21}^T=n_{021}^{T}H\left(\mathrm{\θ}\right)---\left(9\right)$

上式中，H(θ)为行星轮方位角θ对应的坐标变换矩阵，其由式(10)得到：

$H\left(\mathrm{\θ}\right)=\left[\begin{array}{ccc}cos\mathrm{\θ}& -sin\mathrm{\θ}& 0\\ sin\mathrm{\θ}& \cos \mathrm{\θ}& 0\\ 0& 0& 1\end{array}\right]---\left(10\right)$

太阳轮和行星轮等效啮合节点平动自由度之间的刚度耦合如式(11)所示：

$K_{m21}=k_{m21}{\left[\begin{array}{cc}{n}_{21}^T{n}_{21}& -n_{21}^T{n}_{21}\\ -n_{21}^T{n}_{21}& n_{21}^T{n}_{21}\end{array}\right]}_{6\×6}---\left(11\right)$

上式中，k_m21为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得；

太阳轮和行星轮等效啮合节点转动自由度之间的刚度耦合如式(12)所示：

$K_{m\mathrm{\θ}21}=k_{m\mathrm{\θ}21}{\left[\begin{array}{cc}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& -n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\\ -n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& n_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\end{array}\right]}_{6\×6}---\left(12\right)$

上式中，n_θ21为转动自由度耦合矢量，如式(13)所示：

n_θ21＝[-n_y21,n_x21,0]H(θ)； (13)k_mθ21为太阳轮和行星轮啮合弯曲刚度系数，如式(14)所示：

$\begin{array}{c}{k}_{m\mathrm{\θ}21}=\frac{{\mathrm{\ΔM}}_{21}}{{\mathrm{\Δ\γ}}_{21}}=\frac{1}{{\mathrm{\Δ\γ}}_{21}}\[\underset{b/2}{\overset{b}{\∫}}{\mathrm{\Δ\γ}}_{21}z\frac{k_{m21}}{b}(z-\frac{b}{2})dz-\underset{0}{\overset{b/2}{\∫}}{\mathrm{\Δ\γ}}_{21}z\frac{k_{m21}}{b}(\frac{b}{2}-z)dz\]\\ =\frac{b^2}{12}{k}_{m21}\end{array}---\left(14\right)$

上式中，ΔM₂₁表示太阳轮和行星轮接触线产生单位转角Δγ₂₁时，对齿宽中点的弯矩；b为有效齿宽，取太阳轮和行星轮齿宽的较小值；z为齿宽方向的局部坐标；dz为齿宽方向的微变量；

式(15)为太阳轮和行星轮等效啮合节点自由度对应的完整的齿轮等效啮合刚度矩阵：

$K_{m21}={\left[\begin{array}{cccc}{k}_{m21}{n}_{21}^T{n}_{21}& 0& -k_{m21}{n}_{21}^T{n}_{21}& 0\\ 0& k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& 0& -k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\\ -k_{m21}{n}_{21}^T{n}_{21}& 0& k_{m21}{n}_{21}^T{n}_{21}& 0\\ 0& -k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}& 0& k_{m\mathrm{\θ}21}{n}_{\mathrm{\θ}21}^T{n}_{\mathrm{\θ}21}\end{array}\right]}_{12\×12}---\left(15\right)$

上式中，n₂₁为全局坐标系下的等效啮合力作用线方向矢量，如式(9)所示；n_θ21为转动自由度耦合矢量，如式(13)所示；k_m21为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得；k_mθ21为齿轮啮合弯曲刚度系数，如式(14)所示。

4.如权利要求3所述的一种获取行星齿轮错位量的有限元方法，其特征在于：所述步骤4)中建立齿圈与行星轮的等效啮合模型并获得齿圈与行星轮的等效啮合刚度矩阵的过程如下：

齿圈与行星轮为内啮合齿轮副，定义齿圈中心梁单元节点为o₃，其坐标位置与全局坐标系原点O、太阳轮中心梁单元节点o₁位置相同，行星轮中心梁单元节点为o₂，齿圈等效啮合节点为p₃，行星轮在齿圈一端的等效啮合节点为p₂₃，p₃和p₂₃均有6个自由度；o₃与p₃、o₂与p₂₃之间采用刚性梁单元连接，刚度矩阵表示为K_g23，p₃和p₂₃之间采用沿齿轮等效啮合力作用线方向的空间弹簧单元连接，坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₃表示为式(16)：

n₀₂₃＝[n_x23,n_y23,n_z23]^T (16)

上式中，n_x23、n_y23、n_z23分别为n₀₂₃在各坐标轴方向上的分量系数；

等效啮合节点p₃和p₂₃的坐标位置相同，根据式(17)求得：

$\left\{\begin{array}{c}\stackrel{\‾}{o_3{p}_3}-\stackrel{\‾}{o_2{p}_{23}}=L\\ \frac{\stackrel{\‾}{o_2{p}_{23}}}{\stackrel{\‾}{o_3{p}_3}}=\frac{N_2}{N_3}\end{array}\right.---\left(17\right)$

上式中，N₂为行星轮齿数；N₃为齿圈齿数；L为齿圈与行星轮的中心距。

齿圈与行星轮齿轮副的名义啮合分力如式(18)所示：

$\left\{\begin{array}{c}{F}_{t023}=2T_3/d_3\\ F_{r023}=F_{t023}{\mathrm{tan\α}}_n/cos\mathrm{\β}\\ F_{a023}=F_{t023}tan\mathrm{\β}\end{array}\right.---\left(18\right)$

上式中，F_t023为名义切向力；F_r023为名义径向力；F_a023为名义轴向力；T₃为齿圈承受的总转矩在每个行星轮分支上的转矩分量；d₃为齿圈分度圆直径；α_n为法向压力角；β为螺旋角；

如果齿轮存在变位，齿轮副的真实啮合分力如式(19)所示：

$\left\{\begin{array}{c}{F}_{tw23}=T_3/\stackrel{\‾}{o_3{p}_3}\\ F_{rw23}=\sqrt{F_{t023}^2+F_{r023}^2-F_{tw23}^2}\\ F_{aw23}=F_{a023}\end{array}\right.---\left(19\right)$

上式中，F_tw23为真实切向力；F_rw23为真实径向力；F_aw23为真实轴向力；

齿圈与行星轮齿轮副法向啮合力如式(20)所示：

$F_{n23}=\sqrt{F_{tw23}^2+F_{rw23}^2+F_{aw23}^2}---\left(20\right)$

齿圈和行星轮在坐标系o₂x₂y₂z₂下的等效啮合力作用线方向矢量n₀₂₃的各分量系数如式(21)所示：

$\left\{\begin{array}{c}{n}_{x23}=-F_{rw23}/F_n\\ n_{y23}=k_L{F}_{tw23}/F_n\\ n_{z23}=-k_L{k}_R{F}_{aw23}/F_n\end{array}\right.---\left(21\right)$

上式中，工况系数k_L和齿轮旋向系数k_R的取值与式(8)相同；

全局坐标系下齿圈与行星轮齿轮副等效啮合力作用线方向矢量如式(22)所示：

$n_{23}^T=n_{023}^{T}H\left(\mathrm{\θ}\right)---\left(22\right)$

上式中，H(θ)为行星轮方位角θ对应的坐标变换矩阵，由式(10)得到。

齿圈和行星轮等效啮合节点平动自由度之间的刚度耦合如式(23)所示：

$K_{m23}=k_{m23}{\left[\begin{array}{cc}{n}_{23}^T{n}_{23}& -n_{23}^T{n}_{23}\\ -n_{23}^T{n}_{23}& n_{23}^T{n}_{23}\end{array}\right]}_{6\×6}---\left(23\right)$

上式中，k_m23为齿圈和行星轮啮合刚度系数，采用ISO6336给出的计算方法求得；

齿圈和行星轮等效啮合节点转动自由度之间的刚度耦合如式(24)所示：

$K_{m\mathrm{\θ}23}=k_{m\mathrm{\θ}23}{\left[\begin{array}{cc}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& -n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\\ -n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& n_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\end{array}\right]}_{6\×6}---\left(24\right)$

上式中，n_θ23为转动自由度耦合矢量，如式(25)所示：

n_θ23＝[-n_y23,n_x23,0]H(θ) (25)

k_mθ23为齿轮啮合弯曲刚度系数，如式(26)所示：

$\begin{array}{c}{k}_{m\mathrm{\θ}23}=\frac{{\mathrm{\ΔM}}_{23}}{{\mathrm{\Δ\γ}}_{23}}=\frac{1}{{\mathrm{\Δ\γ}}_{23}}\[\underset{b/2}{\overset{b}{\∫}}{\mathrm{\Δ\γ}}_{23}z\frac{k_{m23}}{b}(z-\frac{b}{2})dz-\underset{0}{\overset{b/2}{\∫}}{\mathrm{\Δ\γ}}_{23}z\frac{k_{m23}}{b}(\frac{b}{2}-z)dz\]\\ =\frac{b^2}{12}{k}_{m23}\end{array}---\left(26\right)$

上式中，ΔM₂₃表示齿圈和行星轮接触线产生单位转角Δγ₂₃时，对齿宽中点的弯矩；b为有效齿宽，取齿圈和行星轮齿宽的较小值；z为齿宽方向的局部坐标；dz为齿宽方向的微变量；

式(27)为齿圈和行星轮等效啮合节点自由度对应的完整的齿轮等效啮合刚度矩阵：

$K_{m23}={\left[\begin{array}{cccc}{k}_{m23}{n}_{23}^T{n}_{23}& 0& -k_{m23}{n}_{23}^T{n}_{23}& 0\\ 0& k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& 0& -k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\\ -k_{m23}{n}_{23}^T{n}_{23}& 0& k_{m23}{n}_{23}^T{n}_{23}& 0\\ 0& -k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}& 0& k_{m\mathrm{\θ}23}{n}_{\mathrm{\θ}23}^T{n}_{\mathrm{\θ}23}\end{array}\right]}_{12\×12}---\left(27\right)$

上式中，n₂₃为全局坐标系下的等效啮合力作用线方向矢量，如式(22)所示；n_θ23为转动自由度耦合矢量，如式(25)所示；k_m23为齿轮啮合刚度系数，采用ISO6336给出的计算方法求得；k_mθ23为齿轮啮合弯曲刚度系数，如式(26)所示。

5.如权利要求1所述的一种获取行星齿轮错位量的有限元方法，其特征在于：所述步骤5)中，行星齿轮传动系统静力学方程为：

K(δ)δ＝F (28)

上式中，K(δ)为行星齿轮传动系统的非线性刚度矩阵，由轴模型刚度矩阵K_s、行星架单元刚度矩阵K_r、滚动轴承非线性刚度矩阵K_b、齿轮刚性梁单元刚度矩阵K_g21和K_g23、齿轮等效啮合单元刚度矩阵K_m21和K_m23组集而成；δ为系统模型节点自由度的静变形矢量；F为系统外载荷矢量。

6.如权利要求4所述的一种获取行星齿轮错位量的有限元方法，其特征在于：所述步骤7)中太阳轮和行星轮的错位量的求解过程如下：

设静力学计算求得的太阳轮、行星轮和齿圈的等效啮合节点绕全局坐标系X轴和Y轴的转角变形分别为和根据行星轮方位角θ，将全局坐标系下的转角变形按照式(31)变换到行星轮局部坐标系o₂x₂y₂z₂：

上式中，H(-θ)为行星轮方位角-θ对应的坐标变换矩阵；

太阳轮和行星轮各自在端面啮合力作用线方向上的错位量如式(32)所示：

$\left\{\begin{array}{c}{f}_{sh1}={\mathrm{\Φ}}_{x1}b{\mathrm{cos\α}}_{tw21}+{\mathrm{\Φ}}_{y1}b{\mathrm{sin\α}}_{tw21}\\ f_{sh21}={\mathrm{\Φ}}_{x2}b{\mathrm{cos\α}}_{tw21}+{\mathrm{\Φ}}_{y2}b{\mathrm{sin\α}}_{tw21}\end{array}\right.---\left(32\right)$

上式中，b为有效齿宽，取太阳轮和行星轮齿宽的较小值；α_tw21为考虑齿轮变位的真实端面压力角，如式(33)所示：

${\mathrm{\α}}_{tw21}=acos\left(\frac{d_1+d_2}{2L}{\mathrm{cos\α}}_{t21}\right)---\left(33\right)$

上式中，d₁为太阳轮分度圆直径；d₂为行星轮分度圆直径；L为齿轮副中心距；α_t21为分度圆端面压力角。

太阳轮和行星轮的总错位量如式(34)所示：

f_sh12＝|f_sh1-f_sh21|。 (34)

7.如权利要求6所述的一种获取行星齿轮错位量的有限元方法，其特征在于：所述步骤7)中齿圈和行星轮的错位量的求解过程如下：

齿圈和行星轮各自在端面啮合力作用线方向上的错位量如式(35)所示：

$\left\{\begin{array}{c}{f}_{sh3}={\mathrm{\Φ}}_{x3}b{\mathrm{cos\α}}_{tw23}+{\mathrm{\Φ}}_{y3}b{\mathrm{sin\α}}_{tw23}\\ f_{sh23}={\mathrm{\Φ}}_{x2}b{\mathrm{cos\α}}_{tw23}+{\mathrm{\Φ}}_{y2}b{\mathrm{sin\α}}_{tw23}\end{array}\right.---\left(35\right)$

上式中，b为有效齿宽，取齿圈和行星轮齿宽的较小值；α_tw23为考虑齿轮变位的真实端面压力角，如式(36)所示：

${\mathrm{\α}}_{tw23}=acos\left(\frac{d_3-d_2}{2L}{\mathrm{cos\α}}_{t23}\right)---\left(36\right)$

上式中，d₂为行星轮分度圆直径；d₃为齿圈分度圆直径；L为齿轮副中心距；α_t23为分度圆端面压力角。

齿圈和行星轮的总错位量如式(37)所示：

f_sh32＝|f_sh3-f_sh23|。 (37)