含损伤性单齿故障圆柱直齿轮啮合刚度仿真分析方法

摘要

本发明涉及一种含损伤性单齿故障圆柱直齿轮啮合刚度仿真分析方法。该方法首先基于有限元法和国家标准方法的平均刚度计算结果，提出单、双齿啮合区间啮合刚度修正系数，用于改进能量法计算正常齿轮啮合刚度的精度。其次，针对齿轮故障部位，结合三维建模软件和有限元分析软件建立损伤性故障齿轮的有限元模型，采用计算机语言编制仿真计算程序，计算其时变啮合刚度；最后整合两部分的计算结果求解出故障齿轮的完整啮合刚度。该方法充分综合了修正能量法和有限元法的优势，既保证了计算精度，又提高了计算效率。应用本方法仿真计算的损伤性单齿故障齿轮啮合刚度可有效地用于齿轮系统的振动响应机理研究。

主权项

1.一种含损伤性单齿故障圆柱直齿轮啮合刚度仿真分析方法，其特征在于，包括以下步骤：步骤1正常齿轮啮合刚度计算步骤1.1用国家标准方法计算齿轮副平均啮合刚度选定标准渐开线圆柱直齿轮的参数和材料特性：模数m，齿数z₁、z₂，齿轮宽度H，齿轮轴孔半径r，弹性模量E，泊松比v，材料密度ρ；使用GB-T3480提供的计算公式：c′＝c_thC_MC_RC_B                        (1)计算单对齿刚度，其中C_M＝0.8为理论修正系数，C_R＝1为轮柸结构系数，C_B＝1为基本齿廓系数，单对齿刚度的理论值：$c_{\mathrm{th}}=\frac{1}{q},q=0.04723+\frac{0.15551}{z_1}+\frac{0.25791}{z_2}---\left(2\right)$再由公式：c_γ＝(0.75ε_α+0.25)H·c′×10⁶        (3)即可计算出齿轮平均啮合刚度c_γ，(3)式中ε_α为齿轮端面重合度，有：$ϵ=\frac{1}{2\pi }[z_1(\tan {\alpha }_{a1}-{\tan \alpha }_0)+z_2({\tan \alpha }_{a2}-{\tan \alpha }_0)]---\left(4\right)$其中，α₀＝20°为标准压力角；步骤1.2用有限元法计算齿轮副时变啮合刚度通过在齿轮副有限元模型的接触面上设定接触单元，将这些接触单元与齿轮副进行连接和组装，建立接触系统整体平衡方程为：$\left[\begin{array}{cc}K+\alpha K_p^T& B^T\\ B& 0\end{array}\right]\left\{\begin{array}{c}u\\ \lambda \end{array}\right\}=\left\{\begin{array}{c}R-{\mathrm{\alpha \gamma }}_p\\ -{\gamma }_L\end{array}\right\}---\left(5\right)$其中$B=\underset{e}{\Sigma }{\int }_{\overline{S_e}}{L}^T\mathrm{Nds},K_p=\underset{e}{\Sigma }{\int }_{\overline{S_e}}{N}^T\mathrm{Nds},{\gamma }_L=\underset{e}{\Sigma }{\int }_{\overline{S_e}}{L}^T{g}^0\mathrm{ds},{\gamma }_p=\underset{e}{\Sigma }{\int }_{\overline{S_e}}{N}^T{g}^0\mathrm{ds};$K为总体刚度矩阵，u为节点位移矩阵，R为对角矩阵，N为单元形函数，L为算子，λ为节点处接触内力矩阵；采用迭代法求解平衡方程，提取出轮齿轴孔节点的切向位移量u_y，则齿轮副的啮合刚度：$k_{\theta }^1=\frac{{\mathrm{fN}}_1{r}^2}{r_b^2{u}_y}---\left(6\right)$式中N₁为节点数，r_b为小齿轮基圆半径，f为施加在每个轴孔节点上的切向力；以小齿轮为主动轮，且顺时针转动，初始啮合位置定义为齿轮副刚进入双齿啮合区的临界位置；在齿轮副一个啮合周期内，对于双齿啮合区间和单齿啮合区间分别计算n_ys和n_yd个啮合位置的刚度，第i个啮合位置对应的小齿轮转角位移为θ_i，啮合刚度为k_θi¹；一个啮合周期总的平均啮合刚度为：$c_y=\frac{\underset{i=1}{\overset{n_{\mathrm{ys}}+n_{\mathrm{yd}}}{\Sigma }}{k}_{\mathrm{\theta i}}^1}{n_{\mathrm{ys}}+n_{\mathrm{yd}}}---\left(7\right)$为了校准有限元法的计算精度，做如下比较：$\lambda =\frac{|c_y-c_{\gamma }|}{c_{\gamma }}---\left(8\right)$根据λ的取值不同，以小齿轮基于初始啮合位置的角位移θ为变量的有限元法最终啮合刚度仿真结果可表示为：则双、单齿啮合区间的啮合刚度平均值c_ys和c_yd可分别表示为：$c_{\mathrm{ys}}=\frac{\underset{i=1}{\overset{n_{\mathrm{ys}}}{\Sigma }}{k}_{\mathrm{\theta i}}^y}{n_{\mathrm{ys}}}---\left(10\right),$$c_{\mathrm{yd}}=\frac{\underset{i=1}{\overset{n_{\mathrm{yd}}}{\Sigma }}{k}_{{\theta }_i}^y}{n_{\mathrm{yd}}}---\left(11\right)$步骤1.3用能量法计算齿轮副时变啮合刚度齿轮的啮合刚度由赫兹刚度k_h、弯曲刚度k_b、径向压缩刚度k_a和剪切刚度k_s组成；分别可由下列公式求得：$k_{h,i}=\frac{\mathrm{\pi EH}}{4(1-v^2)}---\left(12\right)$$\frac{1}{k_{b1,i}}={\int }_{{\alpha }_{1,i}}^{{\alpha }_2}\frac{3{\{1+{\cos \alpha }_{1,i}[({\alpha }_2-\alpha )\sin \alpha -\cos \alpha \left]\right\}}^2({\alpha }_2-\alpha )\cos \alpha }{2\mathrm{EH}{[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}^3}\mathrm{d\alpha }---\left(13\right)$$\frac{1}{k_{s1,i}}={\int }_{-{\alpha }_{1,i}}^{{\alpha }_2}\frac{1.2(1+v)({\alpha }_2-\alpha )\cos \alpha {\cos }^2{\alpha }_{1,i}}{\mathrm{EL}[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}\mathrm{d\alpha }---\left(14\right)$$\frac{1}{k_{a1,i}}={\int }_{{-\alpha }_{1,i}}^{{\alpha }_2}\frac{({\alpha }_2-\alpha )\cos \alpha {\sin }^2{\alpha }_{1,i}}{2\mathrm{EH}[\sin \alpha +({\alpha }_2-\alpha )\cos \alpha ]}\mathrm{d\alpha }---\left(15\right)$$\frac{1}{k_{b2,i}}={\int }_{{\alpha }_{1,i}^{\prime }}^{{\alpha }_2^{\prime }}\frac{3{\{1+{\cos \alpha }_{1,i}^{\prime }[({\alpha }_2^{\prime }-\alpha )\sin \alpha -\cos \alpha \left]\right\}}^2({\alpha }_2^{\prime }-\alpha )\cos \alpha }{2\mathrm{EH}{[\sin \alpha +({\alpha }_2^{\prime }-\alpha )\cos \alpha ]}^3}\mathrm{d\alpha }---\left(16\right)$$\frac{1}{k_{s2,i}}={\int }_{-{\alpha }_{1,i}^{\prime }}^{{\alpha }_2^{\prime }}\frac{1.2(1+v)({\alpha }_2^{\prime }-\alpha )\cos {\cos }^2{\alpha }_{1,i}^{\prime }}{\mathrm{EH}[\sin \alpha +({\alpha }_2^{\prime }-\alpha )\cos \alpha ]}\mathrm{d\alpha }---\left(17\right)$$\frac{1}{k_{a2i}}={\int }_{-{\alpha }_{1,i}^{\prime }}^{{\alpha }_2^{\prime }}\frac{({\alpha }_2^{\prime }-\alpha )\cos \alpha {\sin }^2{\alpha }_{1,i}^{\prime }}{2\mathrm{EH}[\sin \alpha +({\alpha }_2^{\prime }-\alpha )\cos \alpha ]}\mathrm{d\alpha }---\left(18\right)$其中，i＝1，2分别对应于齿轮啮合时的第一对轮齿和第二对轮齿，α₂、α′₂分别是小、大齿轮的齿基半角；${\alpha }_{\mathrm{1,1}}={\theta }_1+{\alpha }_{\mathrm{1,1}}^0=\theta -\frac{\pi }{2z_1}-\mathrm{inv}{\alpha }_0+$$\tan \left[\arccos \frac{z_1\cos {\alpha }_0}{\sqrt{{(z_2+2)}^2+(z_1+z_2{)}^2-2(z_2+2)(z_1+z_2)\cos (\arccos \frac{z_2\cos {\alpha }_0}{z_2+1}-{\alpha }_0)}}\right]---\left(19\right)$${\alpha }_{\mathrm{1,1}}^{\prime }{\alpha }_{\mathrm{1,1}}^{O\prime }-{\theta }_2=\tan \left(\arccos \frac{z_2\cos {\alpha }_0}{z_2+2}\right)-\frac{\pi }{{2z}_2}-{\mathrm{inv\alpha }}_0-\frac{z_1}{z_2}\theta ---\left(20\right)$${\alpha }_{\mathrm{1,2}}={\alpha }_{\mathrm{1,1}}+\frac{2\pi }{z_1}=\theta -\frac{3\pi }{{2z}_1}-{\mathrm{inv\alpha }}_0+$$\tan \left[\arccos \frac{z_1\cos {\alpha }_0}{\sqrt{{(z_2+2)}^2+(z_1+z_2{)}^2-2(z_2+2)(z_1+z_2)\cos (\arccos \frac{z_2\cos {\alpha }_0}{z_2+2}-{\alpha }_0)}}\right]---\left(21\right)$${\alpha }_{\mathrm{1,2}}^{\prime }{=\alpha }_{\mathrm{1,1}}^{\prime }-\frac{2\pi }{z_2}=\tan \left(\arccos \frac{z_2\cos {\alpha }_0}{z_2+2}\right)-\frac{5\pi }{{2z}_2}-{\mathrm{inv\alpha }}_0-\frac{z_1}{z_2}\theta ---\left(22\right)$invα₀＝tanα₀-α                        (23)以上各式中，θ为小齿轮的基于初始啮合位置的转角位移；则在小齿轮角位移为θ时齿轮副的总啮合刚度：$k_{\theta }^2=k_1+k_2=\underset{i=1}{\overset{2}{\Sigma }}\frac{1}{\frac{1}{k_{h,i}}+\frac{1}{k_{b1,i}}+\frac{1}{k_{s1,i}}+\frac{1}{k_{a1,i}}+\frac{1}{k_{b2,i}}+\frac{1}{k_{s2,i}}+\frac{1}{k_{a2,i}}}---\left(24\right)$在齿轮副一个啮合周期内，双齿啮合区间和单齿啮合区间分别计算n_ns和n_nd个啮合位置的刚度，第i个啮合位置对应的小齿轮转角为θ_i，啮合刚度为k_θi²，则双、单齿啮合区间的啮合刚度平均值c_ns和c_nd可分别表示为：$c_{\mathrm{ns}}=\frac{\underset{i=1}{\overset{n_{\mathrm{ns}}}{\Sigma }}{k}_{{\theta }_i}^2}{n_{\mathrm{ns}}}---\left(25\right),$$c_{\mathrm{nd}}=\frac{\underset{i=1}{\overset{n_{\mathrm{nd}}}{\Sigma }}{k}_{{\theta }_i}^2}{n_{\mathrm{nd}}}---\left(26\right)$一个啮合周期总的平均啮合刚度为：$c_n=\frac{\underset{i=1}{\overset{n_{\mathrm{ns}}+n_{\mathrm{nd}}}{\Sigma }}{k}_{\mathrm{\theta i}}^2}{n_{\mathrm{ns}}+n_{\mathrm{nd}}}---\left(27\right)$步骤2修正能量法的计算结果基于步骤1中所得国家标准方法和有限元法的计算结果，定义并计算出用于修正能量法的单齿啮合区间啮合刚度修正系数μ_d和双齿啮合区间啮合刚度修正系数μ_s如下：${\mu }_s=\frac{c_{\mathrm{ns}}}{c_{\mathrm{ys}}}---\left(28\right),$${\mu }_d=\frac{c_{\mathrm{nd}}}{c_{\mathrm{yd}}}---\left(29\right)$修正后的正常齿轮任意啮合位置啮合刚度为：步骤3齿轮副故障部位的啮合刚度仿真计算选定齿轮的损伤性单齿故障的类型和故障特性，建立故障齿轮副的有限元模型，之后采用步骤1中的有限元法仿真计算故障齿轮副在啮合位置θ的啮合刚度k_θ³；根据步骤1中λ的取值不同，故障部位的啮合刚度最终仿真结果可表示为：因为损伤性单齿故障会影响齿轮副的两个啮合周期的啮合刚度，所以这个步骤中需要计算故障轮齿处于啮合时的两个连续地啮合周期的啮合刚度；步骤4损伤性单齿故障齿轮完整啮合刚度仿真计算以初始啮合的两对轮齿中左边的轮齿为基准，若故障轮齿是逆时针第a个轮齿；小齿轮旋转一周时，齿轮副有z₁个啮合周期；其中在[1，a-1]和[a+2，z₁]啮合周期的啮合刚度与无故障齿轮的啮合刚度相同，即采用步骤2中修正能量法所仿真计算的结果；在第a个和第a+1个啮合周期的啮合刚度则因故障轮齿处于啮合状态，所以齿轮副的啮合刚度采用步骤3中有限元法所仿真计算的结果；至此，即可得到含损伤性单齿故障齿轮的一个完整周期的啮合刚度。