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

Abstract

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

Description

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

技术领域

本发明属于齿轮测量技术领域，具体涉及一种故障齿轮啮合刚度仿真分析方法，特别是一种基于修正能量法和有限元法的含损伤性单齿故障的圆柱直齿轮时变啮合刚度的仿真分析方法。

背景技术

齿轮是机械传动设备中应用最广泛而且是最容易出现故障的零件之一。齿轮的啮合刚度对于齿轮传动的动力学性能有着显著影响。时变的啮合刚度是齿轮传动系统振动响应的主要动态激励源之一。轮齿变形与啮合刚度随啮合位置变化规律的研究是轮齿修形、动态特性、故障诊断以及寿命预测等研究的基础。所以，有必要深入地研究探讨故障齿轮的啮合刚度变化规律性及其快速有效地仿真分析方法。

能量法和有限元法是目前常用的齿轮啮合刚度计算方法。能量法基于弹性力学知识推导出刚度的积分公式，利用数值计算软件计算出时变啮合刚度的数值解，该方法计算效率较高但是精度和对于故障齿轮啮合刚度的计算不太理想。有限元法一般是利用有限元分析软件建立齿轮副接触有限元模型，通过非线性接触分析功能计算齿轮的时变啮合刚度，计算精度较高，但是计算量很大，效率不高。鉴于此，本发明利用国家标准方法的齿轮平均啮合刚度计算结果，先检验有限元法和能量法的计算误差，提出基于刚度修正系数的修正能量法，改进能量法的计算精度；再结合有限元法精确高效地计算损伤性故障齿轮的时变啮合刚度，本方法对研究损伤性故障圆柱直齿轮系统振动产生与扩展机理及有效的故障诊断技术有着非常重要的意义。

发明内容

本发明为了精确高效地求解含损伤性单齿故障齿轮副时变啮合刚度，提出了一种结合有限元法和修正能量法的齿轮啮合刚度仿真分析方法，采用此方法仿真计算损伤性单齿故障齿轮的时变啮合刚度达到了即保证精度，又提高计算效率的有益效果。

为实现上述目的，本发明的技术方案如下：

一种含损伤性故障圆柱直齿轮啮合刚度仿真分析方法，包括以下具体步骤：

步骤1正常齿轮啮合刚度仿真计算

步骤1.1用国家标准方法计算齿轮副平均啮合刚度

选定标准渐开线圆柱直齿轮的参数和材料特性：模数m，齿数z₁、z₂，齿轮宽度H，齿轮轴孔半径r，弹性模量E，泊松比ν，材料密度ρ；使用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)式中ε_α为齿轮端面重合度，有：

$\mathrm{\ϵ}=\frac{1}{2\mathrm{\π}}[z_1(\tan {\mathrm{\α}}_{a1}-\tan {\mathrm{\α}}_0)+z_2(\tan {\mathrm{\α}}_{a2}-\tan {\mathrm{\α}}_0)]---\left(4\right)$

其中，

α0＝20°为标准压力角。

步骤1.2用有限元法计算齿轮副时变啮合刚度

通过在齿轮副有限元模型的接触面上设定接触单元，将这些接触单元与齿轮副进行连接和组装，建立接触系统整体平衡方程为：

$\left[\begin{array}{cc}K+\mathrm{\α}{K}_p^T& B^T\\ B& 0\end{array}\right]\left\{\begin{array}{c}u\\ \mathrm{\λ}\end{array}\right\}=\left\{\begin{array}{c}R-\mathrm{\α}{\mathrm{\γ}}_p\\ -{\mathrm{\γ}}_L\end{array}\right\}---\left(5\right)$

其中

$B=\underset{e}{\mathrm{\Σ}}{\∫}_{\stackrel{\‾}{S_e}}{L}^T\mathrm{Nds},$ $K_p=\underset{e}{\mathrm{\Σ}}{\∫}_{\stackrel{\‾}{S_e}}{N}^T\mathrm{Nds},$ ${\mathrm{\γ}}_L=\underset{e}{\mathrm{\Σ}}{\∫}_{\stackrel{\‾}{S_e}}{L}^T{g}^0\mathrm{ds},$ ${\mathrm{\γ}}_p=\underset{e}{\mathrm{\Σ}}{\∫}_{\stackrel{\‾}{S_e}}{N}^T{g}^0\mathrm{ds};$

K为总体刚度矩阵，u为节点位移矩阵，R为对角矩阵，N为单元形函数，L为算子，λ为节点处接触内力矩阵。

采用迭代法求解平衡方程，提取出轮齿轴孔节点的切向位移量u_y，则齿轮副的啮合刚度：

$k_{\mathrm{\θ}}^1=\frac{{\mathrm{fN}}_1{r}^2}{r_b^2{u}_y}---\left(6\right)$

式中N₁为节点数，r_b为小齿轮基圆半径，r为小齿轮轴孔半径，f为施加在每个轴孔节点上的切向力。求解(6)式即可求得一个啮合位置齿轮副的啮合刚度值。以小齿轮为主动轮，且顺时针转动，初始啮合位置定义为齿轮副刚进入双齿啮合区的临界位置。在齿轮副一个啮合周期内，对于双齿啮合区间和单齿啮合区间分别计算n_ys和n_yd个啮合位置的刚度，第i个啮合位置对应的小齿轮转角位移为θ_i，啮合刚度为k_θi ¹，则一个啮合周期总的平均啮合刚度为：

$c_y=\frac{\underset{i=1}{\overset{n_{\mathrm{ys}}+\mathrm{nyd}}{\mathrm{\Σ}}}{k}_{\mathrm{\θi}}^1}{n_{\mathrm{ys}}+n_{\mathrm{yd}}}---\left(7\right)$

为了验证有限元法的计算精度，做如下比较：

$\mathrm{\λ}=\frac{|c_y-c_{\mathrm{\γ}}|}{c_{\mathrm{\γ}}}---\left(8\right)$

根据λ的取值不同，以小齿轮基于初始啮合位置的角位移θ为变量的有限元法最终啮合刚度仿真结果可表示为：

则双、单齿啮合区间的啮合刚度平均值c_ys和c_yd可分别表示为：

$c_{\mathrm{ys}}=\frac{\underset{i=1}{\overset{n_{\mathrm{ys}}}{\mathrm{\Σ}}}{k}_{\mathrm{\θi}}^y}{n_{\mathrm{ys}}}---\left(10\right),$ $c_{\mathrm{yd}}=\frac{\underset{i=1}{\overset{n_{\mathrm{yd}}}{\mathrm{\Σ}}}{k}_{{\mathrm{\θ}}_i}^y}{n_{\mathrm{yd}}}---\left(11\right)$

步骤1.3用能量法计算齿轮副时变啮合刚度

齿轮的啮合刚度由赫兹刚度k_h、弯曲刚度k_b、径向压缩刚度k_a和剪切刚度k_s组成。分别可由下列公式求得：

$k_{h,i}=\frac{\mathrm{\πEH}}{4(1-v^2)}---\left(12\right)$

$\frac{1}{k_{b1,i}}={\∫}_{{\mathrm{\α}}_{1,i}}^{{\mathrm{\α}}_2}\frac{3{\{1+\cos {\mathrm{\α}}_{1,i}[({\mathrm{\α}}_2-\mathrm{\α})\sin \mathrm{\α}-\cos \mathrm{\α}\left]\right\}}^2({\mathrm{\α}}_2-\mathrm{\α})\cos \mathrm{\α}}{2\mathrm{EH}{[\sin \mathrm{\α}+({\mathrm{\α}}_2-\mathrm{\α})\cos \mathrm{\α}]}^3}\mathrm{d\α}---\left(13\right)$

$\frac{1}{k_{s1,i}}={\∫}_{-{\mathrm{\α}}_{1,i}}^{{\mathrm{\α}}_2}\frac{1.2(1+v)({\mathrm{\α}}_2-\mathrm{\α})\cos \mathrm{\α}{\cos }^2{\mathrm{\α}}_{1,i}}{\mathrm{EH}[\sin \mathrm{\α}+({\mathrm{\α}}_2-\mathrm{\α})\cos \mathrm{\α}]}\mathrm{d\α}---\left(14\right)$

$\frac{1}{k_{a1,i}}={\∫}_{-{\mathrm{\α}}_{1,i}}^{{\mathrm{\α}}_2}\frac{({\mathrm{\α}}_2-\mathrm{\α})\cos \mathrm{\α}{\sin }^2{\mathrm{\α}}_{1,i}}{2\mathrm{EH}[\sin \mathrm{\α}+({\mathrm{\α}}_2-\mathrm{\α})\cos \mathrm{\α}]}\mathrm{d\α}---\left(15\right)$

$\frac{1}{k_{b2,i}}={\∫}_{{\mathrm{\α}}_{1,i}^{\′}}^{{\mathrm{\α}}_2^{\′}}\frac{3\{1+\cos {\mathrm{\α}}_{1,i}^{\′}{[({\mathrm{\α}}_2^{\′}-\mathrm{\α})\sin \mathrm{\α}-\cos \mathrm{\α}]\}}^2({\mathrm{\α}}_2^{\′}-\mathrm{\α})\cos \mathrm{\α}}{2\mathrm{EH}{[\sin \mathrm{\α}+({\mathrm{\α}}_2^{\′}-\mathrm{\α})\cos \mathrm{\α}]}^3}\mathrm{d\α}---\left(16\right)$

$\frac{1}{k_{s2,i}}={\∫}_{-{\mathrm{\α}}_{1,i}^{\′}}^{{\mathrm{\α}}_2^{\′}}\frac{1.2(1+v)({\mathrm{\α}}_2^{\′}-\mathrm{\α})\cos \mathrm{\α}{\cos }^2{\mathrm{\α}}_{1,i}^{\′}}{\mathrm{EH}[\sin \mathrm{\α}+({\mathrm{\α}}_2^{\′}-\mathrm{\α})\cos \mathrm{\α}]}\mathrm{d\α}---\left(17\right)$

$\frac{1}{k_{a2,i}}={\∫}_{-{\mathrm{\α}}_{1,i}^{\′}}^{{\mathrm{\α}}_2^{\′}}\frac{({\mathrm{\α}}_2^{\′}-\mathrm{\α})\cos \mathrm{\α}{\sin }^2{\mathrm{\α}}_{1,i}^{\′}}{2\mathrm{EH}[\sin \mathrm{\α}+({\mathrm{\α}}_2^{\′}-\mathrm{\α})\cos \mathrm{\α}]}\mathrm{d\α}---\left(18\right)$

其中，i＝1，2分别对应于齿轮啮合时的第一对轮齿和第二对轮齿，α₂、α′₂分别是小、大齿轮的齿基半角；

${\mathrm{\α}}_{\mathrm{1,1}}={\mathrm{\θ}}_1+{\mathrm{\α}}_{\mathrm{1,1}}^0=\mathrm{\θ}-\frac{\mathrm{\π}}{2z_1}-{\mathrm{inv\α}}_0+$

$\tan \left[\arccos \frac{z_1\cos {\mathrm{\α}}_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 {\mathrm{\α}}_0}{z_2+2}-{\mathrm{\α}}_0)}}\right]---\left(19\right)$

${\mathrm{\α}}_{\mathrm{1,1}}^{\′}={\mathrm{\α}}_{\mathrm{1,1}}^{O^{\′}}-{\mathrm{\θ}}_2=\tan \left(\arccos \frac{z_2\cos {\mathrm{\α}}_0}{z_2+2}\right)-\frac{\mathrm{\π}}{2z_2}-{\mathrm{inv\α}}_0-\frac{z_1}{z_2}\mathrm{\θ}---\left(20\right)$

${\mathrm{\α}}_{\mathrm{1,2}}={\mathrm{\α}}_{\mathrm{1,1}}+\frac{2\mathrm{\π}}{z_1}=\mathrm{\θ}-\frac{3\mathrm{\π}}{2z_1}-{\mathrm{inv\α}}_0+$

$\tan \left[\arccos \frac{z_1\cos {\mathrm{\α}}_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 {\mathrm{\α}}_0}{z_2+2}-{\mathrm{\α}}_0)}}\right]---\left(21\right)$

${\mathrm{\α}}_{1,2}^{\′}={\mathrm{\α}}_{\mathrm{1,1}}^{\′}-\frac{2\mathrm{\π}}{z_2}=\tan \left(\arccos \frac{z_2\cos {\mathrm{\α}}_0}{z_2+2}\right)-\frac{5\mathrm{\π}}{2z_2}-{\mathrm{inv\α}}_0-\frac{z_1}{z_2}\mathrm{\θ}---\left(22\right)$

invα₀＝tanα₀-α                   (23)

以上各式中，θ为小齿轮的基于初始啮合位置的转角位移；则在小齿轮角位移为θ时齿轮副的总啮合刚度：

$k_{\mathrm{\θ}}^2=k_1+k_2=\underset{i=1}{\overset{2}{\mathrm{\Σ}}}\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)$

求解(24)式即可求得齿轮一个啮合位置的啮合刚度。在齿轮副一个啮合周期内，双齿啮合区间和单齿啮合区间分别计算n_ns和n_nd个啮合位置的刚度，第i个啮合位置对应的小齿轮转角为θ_i，啮合刚度为k_θi ²，则双、单齿啮合区间的啮合刚度平均值c_ns和c_nd可分别表示为：

$c_{\mathrm{ns}}=\frac{\underset{i=1}{\overset{n_{\mathrm{ns}}}{\mathrm{\Σ}}}{k}_{\mathrm{\θi}}^2}{n_{\mathrm{ns}}}---\left(25\right),$ $c_{\mathrm{nd}}=\frac{\underset{i=1}{\overset{n_{\mathrm{nd}}}{\mathrm{\Σ}}}{k}_{\mathrm{\θi}}^2}{n_{\mathrm{nd}}}---\left(26\right)$

一个啮合周期总的平均啮合刚度为：

$c_n=\frac{\underset{i=1}{\overset{n_{\mathrm{ns}}+n_{\mathrm{nd}}}{\mathrm{\Σ}}}{k}_{\mathrm{\θi}}^2}{n_{\mathrm{ns}}+n_{\mathrm{nd}}}---\left(27\right)$

步骤2修正能量法的计算结果

基于步骤1中所得国家标准方法和有限元法的计算结果，定义并计算出用于修正能量法的单齿啮合区间啮合刚度修正系数μ_d和双齿啮合区间啮合刚度修正系数μ_s如下：

${\mathrm{\μ}}_s=\frac{c_{\mathrm{ns}}}{c_{\mathrm{ys}}}---\left(28\right),$ ${\mathrm{\μ}}_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中有限元法所仿真计算的结果；至此，即可得到含损伤性单齿故障齿轮的一个完整周期的啮合刚度。

本发明的有益效果是：利用国家标准方法的正常齿轮平均啮合刚度计算结果，检验有限元法和能量法的计算精度；并提出了针对能量法的刚度修正系数，改进了能量法的计算结果精度；再结合有限元法计算的故障部分齿轮啮合刚度，进而求解出损伤性故障齿轮完整的时变啮合刚度。该方法综合了修正能量法的高效性和有限元法的精确性，仿真分析的损伤性故障齿轮时变啮合刚度为齿轮系统故障的机理研究提供了准确、可靠的理论基础。

下面具体结合附图与实例对本发明作进一步的说明。

附图说明

图1是本发明的工作流程图；

图2是本发明正常齿轮副有限元模型；

图3是本发明仿真计算啮合刚度的APDL程序编制流程图；

图4是本发明有限元法仿真计算正常齿轮啮合刚度结果；

图5是本发明能量法仿真计算正常齿轮啮合刚度结果；

图6是本发明修正能量法和有限元法仿真计算正常齿轮啮合刚度结果；

图7是本发明有限元法仿真计算齿轮齿根裂纹故障部位啮合刚度结果；

图8是本发明齿根裂纹故障齿轮时变啮合刚度仿真结果；

具体实施方式

如图1所示，是本发明的一种含损伤性单齿故障圆柱直齿轮啮合刚度仿真分析方法的工作流程图。具体实施过程如下：

步骤1正常齿轮啮合刚度仿真计算

步骤1.1用国家标准方法计算齿轮副平均啮合刚度

选定标准渐开线圆柱直齿轮的参数和材料特性：模数m＝3.175，齿数z₁/z₂＝19/48，齿轮厚度H＝16mm，齿轮轴孔半径r＝20mm，弹性模量E＝6.028Gpa，泊松比v＝0.3，材料密度ρ＝7800kg/m³。利用GB-T 3480提供的计算公式：

c′＝c_thC_MC_RC_Bcosβ

计算单对齿刚度，其中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)式中ε_α为齿轮端面重合度，有：

$\mathrm{\ϵ}=\frac{1}{2\mathrm{\π}}[z_1(\tan {\mathrm{\α}}_{a1}-\tan {\mathrm{\α}}_0)+z_2(\tan {\mathrm{\α}}_{a2}-\tan {\mathrm{\α}}_0)]---\left(4\right)$

其中，

α0＝20°为标准压力角；

本例中计算结果c_γ＝3.12×10⁸N/m。

步骤1.2用有限元法计算齿轮副时变啮合刚度

利用三维建模软件SolidWorks建立上述步骤中选定齿轮的三维实体模型。将其保存为parasolid格式文件。在有限元分析软件Ansys中通过导入命令导入parasolid格式的齿轮副模型，经过材料属性定义，网格划分，设置接触对，定义约束，添加载荷后建立起齿轮副有限元仿真分析模型，最后设定分析类型，完成后的齿轮副有限元模型如图2。接触系统整体平衡方程为：

$\left[\begin{array}{cc}K+\mathrm{\α}{K}_p^T& B^T\\ B& 0\end{array}\right]\left\{\begin{array}{c}u\\ \mathrm{\λ}\end{array}\right\}=\left\{\begin{array}{c}R-\mathrm{\α}{\mathrm{\γ}}_p\\ -{\mathrm{\γ}}_L\end{array}\right\}---\left(5\right)$

其中

$B=\underset{e}{\mathrm{\Σ}}{\∫}_{\stackrel{\‾}{S_e}}{L}^T\mathrm{Nds},$ $K_p=\underset{e}{\mathrm{\Σ}}{\∫}_{\stackrel{\‾}{S_e}}{N}^T\mathrm{Nds},$ ${\mathrm{\γ}}_L=\underset{e}{\mathrm{\Σ}}{\∫}_{\stackrel{\‾}{S_e}}{L}^T{g}^0\mathrm{ds},$ ${\mathrm{\γ}}_p=\underset{e}{\mathrm{\Σ}}{\∫}_{\stackrel{\‾}{S_e}}{N}^T{g}^0\mathrm{ds};$

K为总体刚度矩阵，u为节点位移矩阵，R为对角矩阵，N为单元形函数，L为算子，λ为节点处接触内力矩阵；

采用迭代法求解平衡方程，提取出轮齿轴孔节点的切向位移量u_y，则齿轮副的啮合刚度：

$k_{\mathrm{\θ}}^1=\frac{{\mathrm{fN}}_1{r}^2}{r_b^2{u}_y}---\left(6\right)$

式中N_l为节点数，r_b为小齿轮基圆半径，r为小齿轮轴孔半径，f为施加在每个轴孔节点上的切向力；求解(6)式即可求得一个啮合位置齿轮副的啮合刚度值；以小齿轮为主动轮，且顺时针转动，初始啮合位置定义为齿轮副刚进入双齿啮合区的临界位置；在齿轮副一个啮合周期内，对于双齿啮合区间和单齿啮合区间分别计算n_ys和n_yd个啮合位置的刚度，第i个啮合位置对应的小齿轮转角位移为θ_i，啮合刚度为k_θi ¹，则一个啮合周期总的平均啮合刚度为：

$c_y=\frac{\underset{i=1}{\overset{n_{\mathrm{ys}}+\mathrm{nyd}}{\mathrm{\Σ}}}{k}_{\mathrm{\θi}}^1}{n_{\mathrm{ys}}+n_{\mathrm{yd}}}---\left(7\right)$

根据编制上述分析过程编制ansys的APDL程序，并加入循环结构和后处理计算分析命令，形成APDL计算分析程序，程序流程如图3所示。运行该程序可计算出c_y＝3.15×10⁸N/m。

为了验证有限元法的计算精度，做如下比较：

$\mathrm{\&lambda;}=\frac{|c_y-c_{\mathrm{\&gamma;}}|}{c_{\mathrm{\&gamma;}}}\frac{|3.15\&times;{10}^8-3.12\&times;{10}^8|}{3.12\&times;{10}^8}\&ap;0.9\%<1\%---\left(8\right)$

所以以小齿轮基于初始啮合位置的角位移θ为变量的有限元法最终啮合刚度仿真结果可表示为：

$k_{\mathrm{\&theta;}}^y=k_{\mathrm{\&theta;}}^1---\left(9\right)$

有限元法仿真求解的正常齿轮时变啮合刚度，如图4。双、单齿啮合区间的啮合刚度平均值c_ys和c_yd可分别表示为：

$c_{\mathrm{ys}}=\frac{\underset{i=1}{\overset{n_{\mathrm{ys}}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{\&theta;i}}^y}{n_{\mathrm{ys}}}---\left(10\right),$ $c_{\mathrm{yd}}=\frac{\underset{i=1}{\overset{n_{\mathrm{yd}}}{\mathrm{\&Sigma;}}}{k}_{{\mathrm{\&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}}={\&Integral;}_{{\mathrm{\&alpha;}}_{1,i}}^{{\mathrm{\&alpha;}}_2}\frac{3{\{1+\cos {\mathrm{\&alpha;}}_{1,i}[({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\sin \mathrm{\&alpha;}-\cos \mathrm{\&alpha;}\left]\right\}}^2({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}}{2\mathrm{EH}{[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}^3}\mathrm{d\&alpha;}---\left(13\right)$

$\frac{1}{k_{s1,i}}={\&Integral;}_{-{\mathrm{\&alpha;}}_{1,i}}^{{\mathrm{\&alpha;}}_2}\frac{1.2(1+v)({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}{\cos }^2{\mathrm{\&alpha;}}_{1,i}}{\mathrm{EH}[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}\mathrm{d\&alpha;}---\left(14\right)$

$\frac{1}{k_{a1,i}}={\&Integral;}_{-{\mathrm{\&alpha;}}_{1,i}}^{{\mathrm{\&alpha;}}_2}\frac{({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}{\sin }^2{\mathrm{\&alpha;}}_{1,i}}{2\mathrm{EH}[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}\mathrm{d\&alpha;}---\left(15\right)$

$\frac{1}{k_{b2,i}}={\&Integral;}_{{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}^{{\mathrm{\&alpha;}}_2^{\&prime;}}\frac{3\{1+\cos {\mathrm{\&alpha;}}_{1,i}^{\&prime;}{[({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\sin \mathrm{\&alpha;}-\cos \mathrm{\&alpha;}]\}}^2({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}}{2\mathrm{EH}{[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}^3}\mathrm{d\&alpha;}---\left(16\right)$

$\frac{1}{k_{s2,i}}={\&Integral;}_{-{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}^{{\mathrm{\&alpha;}}_2^{\&prime;}}\frac{1.2(1+v)({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}{\cos }^2{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}{\mathrm{EH}[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}\mathrm{d\&alpha;}---\left(17\right)$

$\frac{1}{k_{a2,i}}={\&Integral;}_{-{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}^{{\mathrm{\&alpha;}}_2^{\&prime;}}\frac{({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}{\sin }^2{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}{2\mathrm{EH}[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}\mathrm{d\&alpha;}---\left(18\right)$

其中，i＝1，2分别对应于齿轮啮合时的第一对轮齿和第二对轮齿，α₂、α’₂分别是小、大齿轮的齿基半角；

${\mathrm{\&alpha;}}_{\mathrm{1,1}}={\mathrm{\&theta;}}_1+{\mathrm{\&alpha;}}_{\mathrm{1,1}}^0=\mathrm{\&theta;}-\frac{\mathrm{\&pi;}}{2z_1}-{\mathrm{inv\&alpha;}}_0+$

$\tan \left[\arccos \frac{z_1\cos {\mathrm{\&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 {\mathrm{\&alpha;}}_0}{z_2+2}-{\mathrm{\&alpha;}}_0)}}\right]---\left(19\right)$

${\mathrm{\&alpha;}}_{\mathrm{1,1}}^{\&prime;}={\mathrm{\&alpha;}}_{\mathrm{1,1}}^{O^{\&prime;}}-{\mathrm{\&theta;}}_2=\tan \left(\arccos \frac{z_2\cos {\mathrm{\&alpha;}}_0}{z_2+2}\right)-\frac{\mathrm{\&pi;}}{2z_2}-{\mathrm{inv\&alpha;}}_0-\frac{z_1}{z_2}\mathrm{\&theta;}---\left(20\right)$

${\mathrm{\&alpha;}}_{\mathrm{1,2}}={\mathrm{\&alpha;}}_{\mathrm{1,1}}+\frac{2\mathrm{\&pi;}}{z_1}=\mathrm{\&theta;}-\frac{3\mathrm{\&pi;}}{2z_1}-{\mathrm{inv\&alpha;}}_0+$

$\tan \left[\arccos \frac{z_1\cos {\mathrm{\&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 {\mathrm{\&alpha;}}_0}{z_2+2}-{\mathrm{\&alpha;}}_0)}}\right]---\left(21\right)$

${\mathrm{\&alpha;}}_{1,2}^{\&prime;}={\mathrm{\&alpha;}}_{\mathrm{1,1}}^{\&prime;}-\frac{2\mathrm{\&pi;}}{z_2}=\tan \left(\arccos \frac{z_2\cos {\mathrm{\&alpha;}}_0}{z_2+2}\right)-\frac{5\mathrm{\&pi;}}{2z_2}-{\mathrm{inv\&alpha;}}_0-\frac{z_1}{z_2}\mathrm{\&theta;}---\left(22\right)$

invα₀＝tanα₀-α                   (23)

以上各式中，θ为小齿轮的基于初始啮合位置的转角位移；则在小齿轮角位移为θ时齿轮副的总啮合刚度：

$k_{\mathrm{\&theta;}}^2=k_1+k_2=\underset{i=1}{\overset{2}{\mathrm{\&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)$

求解(24)式即可求得齿轮一个啮合位置的啮合刚度；在齿轮副一个啮合周期内，双齿啮合区间和单齿啮合区间分别计算n_ns和n_nd个啮合位置的刚度，第i个啮合位置对应的小齿轮转角为θ_i，啮合刚度为k_θi ²，则双、单齿啮合区间的啮合刚度平均值c_ns和c_nd可分别表示为：

$c_{\mathrm{ns}}=\frac{\underset{i=1}{\overset{n_{\mathrm{ns}}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{\&theta;i}}^2}{n_{\mathrm{ns}}}---\left(25\right),$ $c_{\mathrm{nd}}=\frac{\underset{i=1}{\overset{n_{\mathrm{nd}}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{\&theta;i}}^2}{n_{\mathrm{nd}}}---\left(26\right)$

一个啮合周期总的平均啮合刚度为：

$c_n=\frac{\underset{i=1}{\overset{n_{\mathrm{ns}}+n_{\mathrm{nd}}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{\&theta;i}}^2}{n_{\mathrm{ns}}+n_{\mathrm{nd}}}---\left(27\right)$

本实例中计算结果c_n＝8.84×10⁸N/m。

本实例中三种方法计算的正常齿轮平均啮合刚度结果比较如表1所示。从表中可以明显知道，有限元法的计算结果较能量法精确，其与标准方法的误差为0.9％。而能量法的计算结果偏大，但因能量法的结果也能反映时变啮合刚度的主要特征，而且在实际实验中发现其计算速度是有限元法的十倍左右。

(将下表从说明书附图移到了说明书中)好的。

表1正常齿轮啮合刚度c(10⁸N/m)

方法	最大值	最小值	平均值
				有限元法	3.66	2.38	3.15
能量法	10.55	5.78	8.84
				国家标准方法	------	------	3.12

步骤2修正能量法的计算结果

基于步骤1中所得国家标准方法和有限元法的计算结果，定义并计算出用于修正能量法的单齿啮合区间啮合刚度修正系数μ_d和双齿啮合区间啮合刚度修正系数μ_s如下：

${\mathrm{\&mu;}}_s=\frac{c_{\mathrm{ns}}}{c_{\mathrm{ys}}}---\left(28\right),$ ${\mathrm{\&mu;}}_d=\frac{c_{\mathrm{nd}}}{c_{\mathrm{yd}}}---\left(29\right)$

修正后的正常齿轮任意啮合位置啮合刚度为：

利用此修正能量法计算的正常齿轮啮合刚度和有限元法的计算结果如图6所示。从图中可以看出，两种方法的计算结果已经非常的接近，这充分地说明了修正能量法的有益效果：改进计算精度的同时，又保持了其计算效率。

步骤3齿轮副故障部位的啮合刚度仿真计算

本实例采用了实际情况中常见的齿根裂纹故障(其它损伤性故障亦适用)，裂纹深度为3mm。利用SolidWorks建立故障齿轮高精度实体模型，导入ansys，步骤1中的有限元法，划分网格时采用混合分网的方式划分有限元网格：用切割命令将齿轮的正常部分和故障部分分离，正常部位用体扫掠划分网格，而有故障部位则用自由网格划分方法。这样即能快速简易地建立故障齿轮副的有限元模型，计算出故障齿轮副在啮合位置θ的啮合刚度k_θ ³；因为步骤1中λ＜1％，所以故障部位的啮合刚度最终仿真结果为：

$k_{\mathrm{\&theta;}}^g=k_{\mathrm{\&theta;}}^3---\left(31\right)$

因为损伤性单齿故障会影响齿轮副的两个啮合周期的啮合刚度，所以这个步骤中需要计算故障轮齿处于啮合时的两个连续地啮合周期的啮合刚度。本实例有限元法计算齿轮故障部位啮合刚度结果如图7所示。

步骤4损伤性单齿故障齿轮完整啮合刚度仿真计算

以初始啮合的两对轮齿中左边的轮齿为基准，故障轮齿是逆时针第2个轮齿；小齿轮旋转一周时，齿轮副有19个啮合周期；其中在1和[4，19]啮合周期的啮合刚度与无故障齿轮的啮合刚度相同，即采用步骤2中修正能量法所仿真计算的结果；在第2个和第3个啮合周期的啮合刚度则因故障轮齿处于啮合状态，所以齿轮副的啮合刚度采用步骤3中有限元法所仿真计算的结果；至此，即可得到含损伤性单齿故障齿轮的一个完整周期的啮合刚度，如图8。

本发明修正了能量法计算齿轮正常部位啮合刚度的结果精度，用有限元法计算齿轮故障部位啮合刚度的结果，将两者整合成故障齿轮的完整啮合刚度，达到了高效精确的目的。仿真分析的损伤性单齿故障齿轮时变啮合刚度为齿轮系统故障的机理研究提供了准确、可靠的理论基础。

Claims (1)

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)式中ε_α为齿轮端面重合度，有：

$\mathrm{\&epsiv;}=\frac{1}{2\mathrm{\&pi;}}[z_1(\tan {\mathrm{\&alpha;}}_{a1}-{\tan \mathrm{\&alpha;}}_0)+z_2({\tan \mathrm{\&alpha;}}_{a2}-{\tan \mathrm{\&alpha;}}_0)]---\left(4\right)$

其中，

α₀＝20°为标准压力角；

步骤1.2用有限元法计算齿轮副时变啮合刚度

通过在齿轮副有限元模型的接触面上设定接触单元，将这些接触单元与齿轮副进行连接和组装，建立接触系统整体平衡方程为：

$\left[\begin{array}{cc}K+\mathrm{\&alpha;}{K}_p^T& B^T\\ B& 0\end{array}\right]\left\{\begin{array}{c}u\\ \mathrm{\&lambda;}\end{array}\right\}=\left\{\begin{array}{c}R-{\mathrm{\&alpha;\&gamma;}}_p\\ -{\mathrm{\&gamma;}}_L\end{array}\right\}---\left(5\right)$

其中

$B=\underset{e}{\mathrm{\&Sigma;}}{\&Integral;}_{\stackrel{\&OverBar;}{S_e}}{L}^T\mathrm{Nds},K_p=\underset{e}{\mathrm{\&Sigma;}}{\&Integral;}_{\stackrel{\&OverBar;}{S_e}}{N}^T\mathrm{Nds},{\mathrm{\&gamma;}}_L=\underset{e}{\mathrm{\&Sigma;}}{\&Integral;}_{\stackrel{\&OverBar;}{S_e}}{L}^T{g}^0\mathrm{ds},{\mathrm{\&gamma;}}_p=\underset{e}{\mathrm{\&Sigma;}}{\&Integral;}_{\stackrel{\&OverBar;}{S_e}}{N}^T{g}^0\mathrm{ds};$

K为总体刚度矩阵，u为节点位移矩阵，R为对角矩阵，N为单元形函数，L为算子，λ为节点处接触内力矩阵；

采用迭代法求解平衡方程，提取出轮齿轴孔节点的切向位移量u_y，则齿轮副的啮合刚度：

$k_{\mathrm{\&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}}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{\&theta;i}}^1}{n_{\mathrm{ys}}+n_{\mathrm{yd}}}---\left(7\right)$

为了校准有限元法的计算精度，做如下比较：

$\mathrm{\&lambda;}=\frac{|c_y-c_{\mathrm{\&gamma;}}|}{c_{\mathrm{\&gamma;}}}---\left(8\right)$

根据λ的取值不同，以小齿轮基于初始啮合位置的角位移θ为变量的有限元法最终啮合刚度仿真结果可表示为：

则双、单齿啮合区间的啮合刚度平均值c_ys和c_yd可分别表示为：

$c_{\mathrm{ys}}=\frac{\underset{i=1}{\overset{n_{\mathrm{ys}}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{\&theta;i}}^y}{n_{\mathrm{ys}}}---\left(10\right),$ $c_{\mathrm{yd}}=\frac{\underset{i=1}{\overset{n_{\mathrm{yd}}}{\mathrm{\&Sigma;}}}{k}_{{\mathrm{\&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}}={\&Integral;}_{{\mathrm{\&alpha;}}_{1,i}}^{{\mathrm{\&alpha;}}_2}\frac{3{\{1+{\cos \mathrm{\&alpha;}}_{1,i}[({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\sin \mathrm{\&alpha;}-\cos \mathrm{\&alpha;}\left]\right\}}^2({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}}{2\mathrm{EH}{[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}^3}\mathrm{d\&alpha;}---\left(13\right)$

$\frac{1}{k_{s1,i}}={\&Integral;}_{-{\mathrm{\&alpha;}}_{1,i}}^{{\mathrm{\&alpha;}}_2}\frac{1.2(1+v)({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}{\cos }^2{\mathrm{\&alpha;}}_{1,i}}{\mathrm{EL}[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}\mathrm{d\&alpha;}---\left(14\right)$

$\frac{1}{k_{a1,i}}={\&Integral;}_{{-\mathrm{\&alpha;}}_{1,i}}^{{\mathrm{\&alpha;}}_2}\frac{({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}{\sin }^2{\mathrm{\&alpha;}}_{1,i}}{2\mathrm{EH}[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}\mathrm{d\&alpha;}---\left(15\right)$

$\frac{1}{k_{b2,i}}={\&Integral;}_{{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}^{{\mathrm{\&alpha;}}_2^{\&prime;}}\frac{3{\{1+{\cos \mathrm{\&alpha;}}_{1,i}^{\&prime;}[({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\sin \mathrm{\&alpha;}-\cos \mathrm{\&alpha;}\left]\right\}}^2({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}}{2\mathrm{EH}{[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}^3}\mathrm{d\&alpha;}---\left(16\right)$

$\frac{1}{k_{s2,i}}={\&Integral;}_{-{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}^{{\mathrm{\&alpha;}}_2^{\&prime;}}\frac{1.2(1+v)({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos {\cos }^2{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}{\mathrm{EH}[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}\mathrm{d\&alpha;}---\left(17\right)$

$\frac{1}{k_{a2i}}={\&Integral;}_{-{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}^{{\mathrm{\&alpha;}}_2^{\&prime;}}\frac{({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}{\sin }^2{\mathrm{\&alpha;}}_{1,i}^{\&prime;}}{2\mathrm{EH}[\sin \mathrm{\&alpha;}+({\mathrm{\&alpha;}}_2^{\&prime;}-\mathrm{\&alpha;})\cos \mathrm{\&alpha;}]}\mathrm{d\&alpha;}---\left(18\right)$

其中，i＝1，2分别对应于齿轮啮合时的第一对轮齿和第二对轮齿，α₂、α′₂分别是小、大齿轮的齿基半角；

${\mathrm{\&alpha;}}_{\mathrm{1,1}}={\mathrm{\&theta;}}_1+{\mathrm{\&alpha;}}_{\mathrm{1,1}}^0=\mathrm{\&theta;}-\frac{\mathrm{\&pi;}}{2z_1}-\mathrm{inv}{\mathrm{\&alpha;}}_0+$

$\tan \left[\arccos \frac{z_1\cos {\mathrm{\&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 {\mathrm{\&alpha;}}_0}{z_2+1}-{\mathrm{\&alpha;}}_0)}}\right]---\left(19\right)$

${\mathrm{\&alpha;}}_{\mathrm{1,1}}^{\&prime;}{\mathrm{\&alpha;}}_{\mathrm{1,1}}^{O\&prime;}-{\mathrm{\&theta;}}_2=\tan \left(\arccos \frac{z_2\cos {\mathrm{\&alpha;}}_0}{z_2+2}\right)-\frac{\mathrm{\&pi;}}{{2z}_2}-{\mathrm{inv\&alpha;}}_0-\frac{z_1}{z_2}\mathrm{\&theta;}---\left(20\right)$

${\mathrm{\&alpha;}}_{\mathrm{1,2}}={\mathrm{\&alpha;}}_{\mathrm{1,1}}+\frac{2\mathrm{\&pi;}}{z_1}=\mathrm{\&theta;}-\frac{3\mathrm{\&pi;}}{{2z}_1}-{\mathrm{inv\&alpha;}}_0+$

$\tan \left[\arccos \frac{z_1\cos {\mathrm{\&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 {\mathrm{\&alpha;}}_0}{z_2+2}-{\mathrm{\&alpha;}}_0)}}\right]---\left(21\right)$

${\mathrm{\&alpha;}}_{\mathrm{1,2}}^{\&prime;}{=\mathrm{\&alpha;}}_{\mathrm{1,1}}^{\&prime;}-\frac{2\mathrm{\&pi;}}{z_2}=\tan \left(\arccos \frac{z_2\cos {\mathrm{\&alpha;}}_0}{z_2+2}\right)-\frac{5\mathrm{\&pi;}}{{2z}_2}-{\mathrm{inv\&alpha;}}_0-\frac{z_1}{z_2}\mathrm{\&theta;}---\left(22\right)$

invα₀＝tanα₀-α                        (23)

以上各式中，θ为小齿轮的基于初始啮合位置的转角位移；则在小齿轮角位移为θ时齿轮副的总啮合刚度：

$k_{\mathrm{\&theta;}}^2=k_1+k_2=\underset{i=1}{\overset{2}{\mathrm{\&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}}}{\mathrm{\&Sigma;}}}{k}_{{\mathrm{\&theta;}}_i}^2}{n_{\mathrm{ns}}}---\left(25\right),$ $c_{\mathrm{nd}}=\frac{\underset{i=1}{\overset{n_{\mathrm{nd}}}{\mathrm{\&Sigma;}}}{k}_{{\mathrm{\&theta;}}_i}^2}{n_{\mathrm{nd}}}---\left(26\right)$

一个啮合周期总的平均啮合刚度为：

$c_n=\frac{\underset{i=1}{\overset{n_{\mathrm{ns}}+n_{\mathrm{nd}}}{\mathrm{\&Sigma;}}}{k}_{\mathrm{\&theta;i}}^2}{n_{\mathrm{ns}}+n_{\mathrm{nd}}}---\left(27\right)$

步骤2修正能量法的计算结果

基于步骤1中所得国家标准方法和有限元法的计算结果，定义并计算出用于修正能量法的单齿啮合区间啮合刚度修正系数μ_d和双齿啮合区间啮合刚度修正系数μ_s如下：

${\mathrm{\&mu;}}_s=\frac{c_{\mathrm{ns}}}{c_{\mathrm{ys}}}---\left(28\right),$ ${\mathrm{\&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中有限元法所仿真计算的结果；至此，即可得到含损伤性单齿故障齿轮的一个完整周期的啮合刚度。

Images