一种基于振动可靠性和遗传算法的齿轮齿廓修形方法

摘要

一种基于振动可靠性和遗传算法的齿轮齿廓修形方法，属于可靠性设计技术领域。本发明提供一种基于振动可靠性和遗传算法的齿轮齿廓修形方法，该方法不但可减小由弹性变形和制造误差引起的啮合冲击，而且还可减小齿轮啮合激励，使得齿轮系统传动平稳，减小了振动和噪声，提高了齿轮系统的可靠性和使用寿命。本发明包括以下几个步骤：步骤一：建立齿轮和修形齿轮啮合弹性变形虚拟样机；步骤二：进行渐开线齿廓和齿廓修形时齿轮传递误差的动态仿真；步骤三：分析齿轮副随机参数啮合传递误差的可靠性敏感度；步骤四：采用遗传算法确定齿轮齿廓修形的最优参数；步骤五：检验参数的正确性。

说明书

技术领域

本发明属于可靠性设计技术领域，特别是涉及一种基于振动可靠性和遗传算法的齿轮齿廓修形方法。 

背景技术

齿轮修形技术被列为齿轮16项关键技术中I类攻关问题。齿轮修形技术是高精度齿轮传动设计和制造的关键技术，它是降低振动和提高高速重载齿轮传动可靠性的一个最重要的手段。齿廓修形为齿轮修形技术的一种方式，齿廓修形为把原来的渐开线齿廓在齿顶或接近齿根圆角的部位修去一部分来缓和啮合刚度的变化，减少由于基节误差和受载变形所引起的啮入和啮出冲击，改善齿面润滑状态，防止胶合发生。国内外对齿轮齿廓修形技术做了大量的研究，唐增宝等将齿轮振动加速度均方根值最小作为优化目标，提出了齿轮动态性能优化设计方法，并利用该方法求得动态性能最佳的齿廓修形量和修形长度。Tavakoli等把齿廓修形作为减小轮齿啮合激励的有效手段，用于消除轮齿的啮入和啮出冲击并最大限度地减小齿轮传动误差的波动。这些研究方法虽然确定了修形量的大小及修形部位的选择等，但是由于轮齿啮合刚度和弹性变形的计算模型还不够精确，而修形参数一般都能精确到微米，所以确定修形参数的方法还无法准确的指导齿轮齿廓修形设计；而且由于考虑齿轮实际的动态特性来设计齿轮修形参数是非常困难的，所以目前大部分减少齿轮振动的方法是基于静态计算。这类方法多数研究都集中在轮齿的几何干涉、载荷突变等问题上，然而精确的确定齿廓修形参数的方法尚不成熟，难以满足大型复杂机械传动重载、高速、高精度的要求。 

随着计算机技术的不断发展，大型的有限元软件的计算速度已经可以接受。因此，将可靠性理论、优化算法与三维齿轮有限元模型结合起来确定齿轮齿廓修形的最优参数，这对研究齿轮齿廓修形设计具有重要的实际指导意义，为齿轮齿廓修形的理论研究工作开辟了新思路和新方法。 

发明内容

针对现有技术存在的上述不足，本发明提供一种基于振动可靠性和遗传算法的齿轮齿廓修形方法，该方法不但可减小由弹性变形和制造误差引起的啮合冲击，而且还可减小齿轮啮合激励，使得齿轮系统传动平稳，减小了振动和噪声，提高了齿轮系统的可靠性和使用寿命。 

为了实现上述目的，本发明采用如下技术方案，一种基于振动可靠性和遗传算法的齿轮齿廓修形方法，包括以下几个步骤： 

步骤一：建立齿轮和修形齿轮啮合弹性变形虚拟样机； 

步骤二：进行渐开线齿廓和齿廓修形时齿轮传递误差的动态仿真； 

步骤三：分析齿轮副随机参数啮合传递误差的可靠性敏感度； 

步骤四：采用遗传算法确定齿轮齿廓修形的最优参数； 

步骤五：检验参数的正确性。 

步骤一所述的建立齿轮和修形齿轮啮合弹性变形虚拟样机，其具体包括以下几个步骤： 

步骤A：建立齿廓曲线的参数方程及全齿模型； 

步骤B：划分齿轮有限元网格； 

步骤C：处理轮齿边界条件及载荷，建立齿轮有限元模型； 

步骤D：选择齿廓修形方式，若是直线修形，则建立直线修形有限元模型；若是抛物线修形，则建立抛物线修形有限元模型。 

步骤二所述的进行渐开线齿廓和齿廓修形时齿轮传递误差的动态仿真，其具体包括以下几个步骤： 

步骤A：进行渐开线齿廓齿轮传递误差的动态仿真； 

步骤B：进行齿廓修形齿轮传递误差的动态仿真； 

步骤C：将传递误差曲线进行对比分析。 

步骤三所述的分析齿轮副随机参数啮合传递误差的可靠性敏感度，其具体包括以下几个步骤： 

步骤A：确定随机输入变量的个数及其分布类型，并判断所需样本点的数目N； 

步骤B：随机抽取输入变量，指定输出变量，得到确定的有限元模型； 

步骤C：利用确定性有限元方法构造随机输入变量与随机输出变量之间的对应关系，并输出结果； 

步骤D：计算所有样本点的输出； 

步骤E：利用响应面法用二次多项式拟合样本点，确定随机输出变量与输入变量之间的响应面函数表达式和极限状态函数表达式； 

步骤F：对响应面函数进行可靠性灵敏度分析，输出结果。 

步骤四所述的采用遗传算法确定齿轮齿廓修形的最优参数，其具体包括以下几个步骤： 

步骤A：对修形参数进行基因编码，并随机初始化群体； 

步骤B：评价群体； 

步骤C：根据停止准则判断是否结束，若结束，则转去执行步骤E；否则，执行步骤D； 

步骤D：对群体应用染色体选择算子、交叉算子和变异算子；然后再评价群体，并转去执行步骤C； 

步骤E：优化修形参数，对直线修形的参数和抛物线修形的参数进行优化。 

步骤五所述的检验参数的正确性，其具体包括以下几个步骤： 

步骤A：进行齿轮非线性动力学建模； 

步骤B：对齿轮非线性动力学方程进行求解； 

步骤C：进行最优方案的减振论证。 

本发明的有益效果： 

1、本发明考虑了振动可靠性这一重要因素： 

本发明对齿轮进行齿廓修形来消除理论渐开线齿轮的啮入和啮出冲击，减小振动和噪声，并考虑多种工程实际结构中的随机变量的影响，这是一种动态的计算方法，计算结果更精确，更符合现代高速、重载齿轮的实际要求。通过分析计算得出了各随机修形参数对传递误差的可靠性敏感度，得到敏感度图和散点图。敏感度图表明要改变结构的可靠性应该修改哪些设计变量，散点图则进一步表明应该怎样改动这些设计变量及改动的大致范围。减小各随机变量的离散范围并提高或降低相应的参数，可以有效提高齿轮的传递误差可靠性，可靠性敏感度对于指导工程实际的可靠性设计具有重要意义。 

2、本发明采用遗传算法进行齿轮齿廓修形参数的优化： 

优化设计的修形参数能大幅度的减小齿轮的振动，避免啮合瞬间冲击的产生。本发明通过分别对直线修形和抛物线修形时的齿轮修形参数进行优化设计，确定了这两种修形曲线所能达到的最优减振效果，得到最优的修形曲线，对齿轮齿廓修形设计有重要的意义。 

3、本发明对参数的正确性进行检验： 

本发明在齿轮副间隙扭转振动非线性动力学模型基础上，对方程进行求解并绘制了修形前和修形后的幅频曲线，通过对其前后振幅的对比，由振动的角度验证了优化后齿轮修形参数的正确性，具有理论依据。 

综上所述，本发明不但减小了由弹性变形和制造误差引起的啮合冲击，而且还减小了振动和噪声；有效的提高了齿轮的传递误差可靠性，使得齿轮系统传动平稳，提高了齿轮系统的可靠性和使用寿命，对于指导工程实际的可靠性优化设计具有重要意义。 

附图说明

图1为本发明的修形方法的程序流程图； 

图2为建立齿轮和修形齿轮啮合弹性变形虚拟样机的程序流程图； 

图3为进行渐开线齿廓和齿廓修形时齿轮传递误差的动态仿真的程序流程图； 

图4为分析齿轮副随机参数啮合传递误差的可靠性敏感度的程序流程图； 

图5为采用遗传算法确定齿轮齿廓修形的最优参数的程序流程图； 

图6为检验参数的正确性的程序流程图； 

图7为齿廓渐开线的生成原理图； 

图8为齿条型刀具结构图； 

图9为齿条型刀具加工齿轮的过程图； 

图10为齿轮啮合弹性变形虚拟样机的实体模型图； 

图11为齿轮啮合虚拟试验样机的有限元模型图； 

图12为齿轮修形参数图； 

图13为齿廓直线修形的模型图； 

图14为齿廓渐开线在b点的切线斜率图； 

图15为齿廓抛物线修形的模型图； 

图16为直线修形的齿轮实体模型图； 

图17为抛物线修形的齿轮实体模型图； 

图18为直线修形齿轮的啮合虚拟样机有限元模型图； 

图19为抛物线修形齿轮的啮合虚拟样机有限元模型图； 

图20为齿轮的传递误差动态曲线图； 

图21为直线修形的齿轮传递误差动态曲线图； 

图22为抛物线修形的齿轮传递误差动态曲线图； 

图23为抽样样本图； 

图24为频率分布直方图； 

图25为累积分布函数反映齿轮失效的概率图； 

图26为随机参数的柱状图和饼图； 

图27为极限状态函数值g(x)与各随机输入变量间的散点图； 

图28为基因编码图； 

图29为直线修形的齿轮传递误差动态曲线图； 

图30为抛物线修形的齿轮传递误差动态曲线图； 

图31为直齿圆柱齿轮副纯扭转振动的动力学模型图； 

图32为修形前后幅频曲线对比图。 

具体实施方式

下面以圆柱直齿轮为例，来进行基于振动可靠性和遗传算法的齿轮齿廓修形，齿轮的材料为16Cr3NiWMoVNbE，主动轮齿数z₁为26，从动轮齿数z₂为40，模数m为3mm，压力角α₀为20°，齿宽B为20mm，顶隙系数c^*为0.25，主、从动轮齿形偏差分别为Ef_p＝8μm， Ef_g＝8μm，综合齿距偏差ES＝10μm，主、从动轮的最大修形量为S_p和S_g，修形角度为α_p和α_g，S_p＝0.06mm，α_p＝63°，S_g＝0.06mm，α_g＝63°，本发明的程序流程图如图1所示。 

步骤一：建立齿轮和修形齿轮啮合弹性变形虚拟样机； 

齿轮的传递误差与综合变形量相关，精确分析需要通过建立三维实体模型并进行有限元计算。为了模拟齿轮原始加工误差的随机性，本发明选用ANSYS有限元软件来进行齿轮的参数化建模。 

步骤A：建立齿廓曲线的参数方程及全齿模型； 

齿廓渐开线的生成原理，如图7所示，在以齿轮回转中心为极点，以x轴为极轴的极坐标系中，齿廓渐开线方程可表示为 

$\left(\begin{array}{c}{r}_k=r_b/\cos {\alpha }_k\\ {\theta }_k=\tan {\alpha }_k-{\alpha }_k\end{array}\right)---\left(1\right)$

式中，r_k为渐开线上任意一点的矢径，r_b为基圆半径，α_k为该点的压力角，θ_k为该点的展角。 

通过坐标变换，得到直角坐标系下渐开线的参数方程，即齿廓曲线的参数方程为：

$\left(\begin{array}{c}{x}_k=\frac{r_b}{\cos {\alpha }_k}\sin (\frac{\phi }{2}+\mathrm{inv}{\alpha }_0-\mathrm{inv}{\alpha }_k)\\ y_k=\frac{r_b}{\cos {\alpha }_k}\cos (\frac{\phi }{2}+\mathrm{inv}{\alpha }_0-\mathrm{inv}{\alpha }_k)\end{array}\right)---\left(2\right)$

式中，x_k、y_k分别为渐开线上任意一点在直角坐标系下的横、纵坐标；α₀为齿轮分度圆上的压力角；φ为分度圆齿厚所对应的中心角，对于齿数为z的标准直齿圆柱齿轮，φ＝π/z。 

根据加工方式及加工刀具齿廓的不同，齿根过渡曲线有多种形式。采用齿条插刀或滚刀加工齿轮时，如果刀具齿廓的顶部只具有一个圆角，则过渡曲线为一整段延伸渐开线的等距曲线。齿条型刀具结构，如图8所示，图中，α₀为齿轮分度圆上的压力角，a为刀具圆角圆心C_ρ距中线的距离；b为刀具圆角圆心C_ρ距刀具齿槽中心线的距离；r_ρ为刀具圆角半径； 为齿高系数；c^*为径向间隙系数；m为模数。齿条型刀具加工齿轮的过程，如图9所示，图中，C为节点，nn为刀具圆角与过渡曲线接触点的公法线，a为刀具圆角圆心C_ρ距中线的距离，b为刀具圆角圆心C_ρ距刀具齿槽中心线的距离，r_ρ为刀具圆角半径，C_ρ为刀具圆角圆心，齿根过渡曲线的参数间具有如下关系： 

$\left(\begin{array}{c}a=h_a^{*}m+c^{*}m-r_{\rho }\\ b=\frac{\mathrm{\pi m}}{2}\\ r_{\rho }=\frac{\mathrm{\pi m}-4h_a^{*}m\tan {\alpha }_0}{4\cos {\alpha }_0}\\ c^{*}m=r_{\rho }(1-\sin {\alpha }_0)\end{array}\right)---\left(3\right)$

式中，α₀为齿轮分度圆上的压力角，m为模数， 为齿高系数，c^*为径向间隙系数。对于标准直齿圆柱齿轮，得出图9所示坐标系下延伸渐开线等距曲线的参数方程为： 

式中， xk、yk分别为渐开线上任意一点在直角坐标系下的横、纵坐标，α′_k为压力角，m为模数，z为齿数， 为圆心角，a为刀具圆角圆心C_ρ距中线的距离，b为刀具圆角圆心C_ρ距刀具齿槽中心线的距离，r_ρ为刀具圆角半径。 

在ANSYS中建立精确的齿轮齿廓比较困难，根据式(2)和式(4)，利用APDL循环建点命令建立关键点。对渐开线共建立了1000个点，齿根过渡曲线建立300个点，然后将各个关键点连线即可获得由标准渐开线和过渡曲线构成的齿廓。将单齿模型旋转镜像并拉伸成体，即可得到一个齿轮的三维模型。 

将主、从动齿轮分别沿其轴线旋转，由于分度圆上的齿厚和齿槽宽相等，所以主动齿轮旋转到节点啮合位置需旋转360/(4z₁)，从动齿轮旋转到节点啮合位置需旋转360/(4z₂)。建立齿轮啮合弹性变形虚拟样机的实体模型，即全齿模型，如图10所示。 

步骤B：划分齿轮有限元网格； 

有限元自动建模的关键是结构的网格自动生成。有限元网格生成技术主要可分成结构化网格生成和非结构化网格生成两大类。结构化网格生成法主要采用映射法，这种方法必须事先将结构分割成子区域，在子区域内进行网格自动生成，能够用于曲面网格的生成。非结构化网格生成法可以实现不同程度的自动化，按操作特点可分成三角化法、拓扑分解法、几何分解法和自由网格划分法等。网格的尺寸与有限元分析的效率和精度密切相关，网格越密，计算精度越高，但是计算效率越低。因此，可在齿轮接触区域采用较密的网格，在其他区域采用大尺寸的网格，本发明选择智能网格划分方式生成有限元模型。根据材料定义弹性模量 为2.06×10⁵MPa，泊松比为0.3，齿面摩擦系数为0.05。选择三维六面体单元SOLID45对样机的实体模型划分网格。 

步骤C：处理轮齿边界条件及载荷，建立齿轮有限元模型； 

在齿轮传动过程中，外载荷是从动齿轮的工作阻力和主动齿轮的驱动力矩。为了真实的模拟实际工作过程中齿轮的啮合行为，将齿轮做如下处理：假设在两齿轮轮齿接触的瞬间，从动轮固定不动，主动轮作定轴转动，其径向位移和轴向位移为零约束，即在以从动轮轴线为Z轴的局部圆柱坐标系下，对从动轮内圆柱面上的所有节点施加U_X、U_Y、U_Z约束，将其固定；再以主动轮轴线为Z轴的局部圆柱坐标系下，转换主动轮内圆柱面上所有节点的节点坐标系为圆柱坐标系，并对主动轮内圆柱面上的所有节点施加U_X、U_Z约束，并保持其旋转自由度U_Y；若主动轮内圆柱面的半径为r_z，主动轮输入转矩为T₁，则在主动轮内轴圆柱面的所有节点上平均施加沿旋转方向总大小为T₁/r_z的力。创建接触对、施加载荷及边界条件后，建立虚拟试验样机的有限元模型，如图11所示。 

步骤D：选择齿廓修形方式，若是直线修形，则建立直线修形有限元模型；若是抛物线修形，则建立抛物线修形有限元模型。 

齿轮修形参数，如图12所示，图中，S_p为主动轮最大修形量，S_g为从动轮最大修形量，α_p为主动轮修形角度，α_g为从动轮修形角度。齿廓修形方式包括直线修形和抛物线修形。 

直线修形：根据齿轮修形量和修形角度的值可以直接确定修形直线的具体位置，建立轮齿模型并且按直线修形，如图13所示，图中，a、b两点分别为直线与齿顶、渐开线齿廓的交点。 

抛物线修形：直线修形和抛物线修形的修形量和修形角度相同，均为S_p和S_g，α_p和α_g，即两种修形的轮廓曲线均过图13中的a、b两点。若设a、b两点的坐标分别为(x_a，y_a)和(x_b，y_b)，可得： 

$\left(\begin{array}{c}{y}_a=c_1{x}_a^2+c_2{x}_a+c_3\\ y_b=c_1{x}_b^2+c_2{x}_b+c_3\end{array}\right)---\left(5\right)$

修形后的齿廓曲线上，渐开线与抛物线连接处应光滑，以保证齿轮传动的稳定，即抛物线与渐开线在b点的切线斜率相等。齿廓渐开线在b点的切线斜率，如图14所示，图中，α′为啮合角，θ为展角，r_b为基圆半径，(x_b，y_b)为b点的坐标。 

由渐开线特性可知，切线cd与渐开线发生线be垂直，即可得cd平行于基圆圆心O与发 生线切点e的连线Oe。由此可得渐开线在b点处的切线斜率K_b。 

$K_b=K_{\mathrm{oe}}=\tan ({\alpha }^{\prime }-\theta +\frac{\pi }{2})=-\cot ({\alpha }^{\prime }-\theta )---\left(6\right)$

由图中几何关系可得： 

$\theta =\arctan \left(\frac{x_b}{y_b}\right)---\left(7\right)$

${\alpha }^{\prime }=\arccos \left(\frac{L_{\mathrm{oe}}}{L_{\mathrm{ob}}}\right)---\left(8\right)$

式中，$L_{\mathrm{ob}}=\sqrt{x_b^2+y_b^2},$L_oe＝r_b。 

抛物线在b点处的切线斜率表达式为： 

$y_b^{\prime }=\frac{\mathrm{dy}}{\mathrm{dx}}=2c_1{x}_b+c_2---\left(9\right)$

根据抛物线和渐开线在b点切线斜率相等，再由式(5)可得 

$\left(\begin{array}{c}{y}_a=c_1{x}_a^2+c_2{x}_a+c_3\\ y_b=c_1{x}_b^2+c_2{x}_b+c_3\\ K_b=2c_1{x}_b+c_2\end{array}\right)---\left(10\right)$

式中，K_b为渐开线在b点处的切线斜率。 

可解得系数c₁、c₂、c₃，代入抛物线表达式即可得齿廓修形抛物线。编制参数化齿形曲线生成程序，得到齿廓抛物线修形的模型，如图15所示。 

分别对齿廓曲线进行直线修形和抛物线修形后，建立的实体模型如图16和图17所示。 

根据齿轮有限元网格划分法和边界条件及载荷处理方法，在齿轮接触区域采用较密的网格，在其他区域采用大尺寸的网格。创建接触对、施加载荷及边界条件后，建立的直线修形和抛物线修形齿轮的啮合虚拟样机有限元模型，如图18和图19所示。 

步骤二：渐开线齿廓和齿廓修形时齿轮传递误差的动态仿真； 

对步骤一得到的有限元模型进行传递误差的仿真。 

步骤A：进行渐开线齿廓齿轮传递误差的动态仿真； 

根据步骤一中的已知条件进行仿真计算，得到齿轮的传递误差曲线和传递误差的波动值，如图20所示，由图可见，齿轮在双齿啮合区的传递误差较小，而在单齿啮合区传递误差较大。所以，在单双齿交替区域产生了很大的波动，导致曲线呈阶梯形状，说明轮齿在啮入和啮出时会产生冲击，传递误差的最大波动为10μm。 

步骤B：进行齿廓修形齿轮传递误差的动态仿真； 

对齿轮采用有限元建模型进行仿真计算。修形参数按步骤一中选取，为了比较不同修形 曲线对传递误差的影响，分别对不同修形曲线进行仿真计算，图21、图22为分别用有限元模拟的按直线和抛物线修形时齿轮的传递误差曲线图。 

步骤C：将传递误差曲线进行对比分析； 

由图21可见，进行直线修形时，齿轮系统的传递误差的波动ΔTE＝7.98μm。对比修形前的传递误差曲线图20，这时的曲线已经没有明显的阶梯性，说明啮入啮出冲击已经不存在，而双齿啮合区的波峰说明单双齿交替是产生波动的主要原因，而在波峰处变化剧烈说明有冲击产生，振幅减小了大约20％。修形效果不好，需要对修形参数进行进一步的优化选取。 

由图22可见，进行抛物线修形时，齿轮系统的传递误差的波动ΔTE＝7.13μm。对比修形前的传递误差曲线图20，曲线没有明显的阶梯性，说明啮入啮出冲击已经不存在，在双齿啮合区存在较平缓的波峰，说明单双齿交替是产生波动的主要原因，但是曲线平缓，振幅减小了大约28％，效果比直线修形好。但是也需要进行修形参数的优化计算。 

虽然按抛物线修形时，齿轮系统的传递误差曲线波动小，而按直线修形时的波动大。但是由于对应不同类型的修形曲线，应取的修形量也是不同的。就是说在本实施方式中这两种修形曲线所能达到的最优效果并没有体现出来。所以不能说明抛物线修形一定比直线修形好。而本步骤的各种有限元模型及其传递误差的仿真正是为了后续的分析计算奠定基础。 

步骤三：分析齿轮副随机参数啮合传递误差的可靠性敏感度： 

由步骤二可知传递误差受各随机修形参数的影响，因此有必要采用振动可靠性的方法来对各个参数的随机性进行分析。 

步骤A：确定随机输入变量的个数及其分布类型，并判断所需样本点的数目N； 

为了考察修形参数对齿轮啮合传递误差可靠性的影响，设齿轮啮合传递误差波动ΔTE受主动轮齿顶最大修形量S_p、主动轮齿顶修形角α_p、从动轮齿顶最大修形量S_g、从动轮齿顶修形角α_g的影响，且S_p、α_p、S_g、α_g都是随机变量，最大许用传递误差波动ΔTE_max也是随机变量。假定以上随机变量均服从正态分布，其均值和标准差如表1所示，根据中心复合设计(k＝4)需要的取样点数为NS＝2^k+2k+1＝25。 

表1 

步骤B：随机抽取输入变量，指定输出变量，得到确定的有限元模型； 

随机抽取输入变量S_p、α_p、S_g、α_g，指定输出变量为ΔTE_max，即可建立确定的有限元 模型。 

步骤C：利用确定性有限元方法构造随机输入变量与随机输出变量之间的对应关系，并输出结果； 

根据确定的有限元模型进行啮合仿真计算，得到传递误差的波动ΔTE_max。 

步骤D：计算所有样本点的输出； 

概率水平分别为p₁＝0.01、p₂＝0.50、p₃＝0.99时，计算出样本点数值列入表2中。 

根据表2中的样本点经过25次有限元模拟，得到齿轮副传递误差波动的25个响应值，列入表2的最后一列。 

表2 

步骤E：利用响应面法用二次多项式拟合样本点，确定随机输出变量与输入变量之间的响应面函数表达式和极限状态函数表达式； 

响应面函数通常都选取含有交叉项的二次函数来描述，可表示为如下表达式： 

$\hat{Y}=C_0+\underset{i=1}{\overset{\mathrm{NR}}{\Sigma }}{C}_i{X}_i+\underset{i=1}{\overset{\mathrm{NR}}{\Sigma }}\underset{j=i}{\overset{\mathrm{NR}}{\Sigma }}{C}_{\mathrm{ij}}{X}_i{X}_j---\left(11\right)$

式中，C₀、C_i、C_ij(i＝1…NR；j＝i…NR)是待定系数。 

通过对随机参数向量的NS个样本点进行数值计算，得到NS个输出点(y₁，y₂，…，y_NS)，用最小二乘法对随机参数向量和结构响应进行回归分析 

$s=\underset{i=1}{\overset{\mathrm{NS}}{\Sigma }}{ϵ}^2=\underset{i=1}{\overset{\mathrm{NS}}{\Sigma }}{[y_i-(C_0+\underset{i=1}{\overset{\mathrm{NR}}{\Sigma }}{C}_i{x}_i+\underset{i=1}{\overset{\mathrm{NR}}{\Sigma }}\underset{j=i}{\overset{\mathrm{NR}}{\Sigma }}{C}_{\mathrm{ij}}{x}_i{x}_j)]}^2---\left(12\right)$

式中，ε为误差项，为使误差项最小，则有 

$\left(\begin{array}{c}\frac{\partial s}{{\partial C}_0}=0\\ \frac{\partial s}{{\partial C}_i}=0,i=\mathrm{1,2},L,\mathrm{NR}\\ \frac{\partial s}{{\partial C}_{\mathrm{ij}}}=0,i=\mathrm{1,2},L,\mathrm{NR};j=i,L,\mathrm{NR}\end{array}\right)---\left(13\right)$

根据表2中的数据和式(11)、(12)、(13)，得到齿轮副传递误差波动的响应面函数为： 

$\hat{Y}=173.086+5.11S_p-14.122{\alpha }_p-1.086S_g+5.325{\alpha }_g-0.017S_p^2$

$-0.195S_p{\alpha }_p+0.009S_p{S}_g+0.139S_p{\alpha }_g-1.587{\alpha }_p^2+0.164{\alpha }_p{S}_g---\left(14\right)$

$+3.399{\alpha }_p{\alpha }_g-0.006S_g^2-0.143S_g{\alpha }_g-1.728{\alpha }_g^2$

得到极限状态函数为： 

$g\left(X\right)={\mathrm{\Delta TE}}_{\max }-\hat{Y}={\mathrm{\Delta TE}}_{\max }-(173.086+5.11S_p-14.122{\alpha }_p-1.086S_g$

$+5.325{\alpha }_g-0.017S_p^2-0.195S_p{\alpha }_p+0.009S_p{S}_g+0.139S_p{\alpha }_g---\left(15\right)$

$-1.587{\alpha }_p^2+0.164{\alpha }_p{S}_g+3.399{\alpha }_p{\alpha }_g-0.006S_g^2-0.143S_g{\alpha }_g-1.728{\alpha }_g^2)$

步骤F：对响应面函数进行可靠性灵敏度分析，输出结果。 

指定齿轮修形参数S_p、α_p、S_g、α_g和最大许用传递误差波动ΔTE_max为随机输入变量，极限状态函数g(x)作为随机输出变量，并选择Monte-Carlo概率设计方法，进行5000次Monte-Carlo模拟，生成了g(x)的函数值的5000个样本，抽样样本如图23所示。对极限状态函数g(x)进行统计分析，得到频率分布直方图，如图24所示。 

图25所示的累积分布函数反映齿轮失效概率，在任一点的数值等于数据出现该点之下的概率。当g(x)＜0时齿轮失效，查看结果得不可靠度为0.32，说明齿轮可靠度为0.68。所以必须更改齿轮的修形参数来减小振动。 

若不考虑随机变量ΔTE_max，与齿轮传递误差相关的修形系数有4个，更改哪个参数以及该参数对齿轮传递误差可靠性的影响是增大还是减小，这都需要通过计算随机参数敏感度来 判断。概率敏感性可以采用统计显著性检验表示，给定显著性水平α，可将随机输入变量分为：对g(x)有显著性影响的和无显著影响的。显著性检验假定随机输入变量的敏感性为0，并计算其概率，当此概率值超过显著性水平α时，该随机输入变量的影响将被忽略；否则，就认为该随机输入变量对g(x)有显著性影响。取显著性水平α＝2.5％，用柱状图和饼图表示各随机输入变量对g(x)的影响，如图26所示。 

由饼图中可以看出各随机参数敏感性的比例对比。在柱状图中可看出各个参数的影响程度，敏感性最大的随机输入变量出现在最左边，并依照从大到小的顺序，依次向右排列；而且，敏感性有正负之分，随机输入变量出现在柱状图纵轴的正半轴说明敏感性为正，表示响应g(x)随该随机输入变量的增加而增加，敏感性为负表示响应g(x)随该随机输入变量的增加而减小。 

由图26可见，随机输入变量S_p、α_p、S_g、α_g均对g(x)有显著性影响，按照影响程度从大到小的次序排列依次为α_p、α_g、S_p、S_g，且g(x)随α_p的增大而增大，随S_p、α_g、S_g的增大而减小。表3列出了修形参数可靠性敏感度的值。 

表3 

敏感度图提示了要改变结构的可靠性应该修改哪些设计变量，而散点图则进一步提示了应该怎样改动这些设计变量及改动的大致范围。分别绘制极限状态函数值g(x)与各随机输入变量间的散点图，如图27所示。 

可见，减小各随机变量的离散范围和通过提高α_p及降低S_p、α_g、S_g的值，可以有效提高齿轮的传递误差可靠性。 

步骤四：采用遗传算法确定齿轮齿廓修形的最优参数： 

由步骤三得知了为提高齿轮系统的可靠性各随机参数的变化趋势和敏感度，这缩小了优化搜索参数的区间，并在相应的搜索区间中为步长的选取奠定了基础。 

步骤A：对修形参数进行基因编码，并随机初始化群体； 

1)对修形参数进行基因编码 

将修整高度换算成修形角度，得到修形角度和修形量的最大值分别为64.331°和0.06mm。为了更符合修形参数所在的区间，将搜索范围适当平移，结合遗传算法的编码方法，定义修形参数优化的取值范围，如表4所示。 

用遗传算法求解问题，必须先确定问题的编码方法，为了表示的方便，在进行编码时一般都采用二进制编码方式，二进制编码将问题空间的参数表示为0和1构成的染色体位串。由于16位的二进制数可以表示0至65535之间的值。根据表4的取值范围，若修形参数α_p、α_g、S_p、S_g各占16位，则它们的搜索跨度都是65535。可见精度足够，基因编码如图28所示。 

2)随机初始化群体 

遗传算法中，种群的初始化工作对算法的运行效果有很大的影响。第i代群体可表示为： 

P(i)＝(s₁，s₂...s_n)(16) 

式中，P(i)表示第i代群体，S表示染色体，n为群体规模。 

遗传算法中，种群的初始化工作对算法的运行效果有很大的影响。由于修形参数可根据经验公式确定大致的取值范围，所以将修形参数的取值区间作为边界条件，选择n作为群体的规模，随机产生n个符合边界条件的，位数为64的2进制编码的染色体作为初始种群，这里取n＝30。 

步骤B：评价群体； 

根据编码方式将染色体位串解码得修形量S_p和S_g，修形角度α_p和α_g。 

遗传算法对一个解的好坏用适应度函数来评价，适应度越大，解的质量越好。在应用遗传算法的实现之前应先定义适应度函数。总的说来，有两种方式可供选择：一种从轮齿所受冲击最小方面考虑；另一种是从齿轮的传递误差波动最小来考虑。由于齿轮传递误差恒定时齿轮无振动，所以将齿轮的传递误差波动倒数作为适应值来评价群体比较合理。定义适应值函数如下： 

$f=\frac{1}{|\max \left(\mathrm{TE}\right)-\min \left(\mathrm{TE}\right)|}---\left(17\right)$

步骤C：适应值计算，并根据停止准则判断是否结束，若结束，则转去执行步骤E；否则，执行步骤D； 

将上面第一步的解码结果代入有限元模型，计算出传递误差的最大和最小值，带入适应值函数式(17)就能算出该染色体的适应值。对每一个染色体计算其适应值并按照适应值大小将染色体排序。当计算到第K代时，若群体收敛则停止，转去执行步骤E，否则，执行步骤D。 

步骤D：对群体应用染色体选择算子、交叉算子和变异算子；然后再评价群体，并转去执行步骤C； 

1)对群体应用染色体选择算子： 

选择操作根据个体的适应度决定遗传到下一代的个体，首先在群体中删除不符合条件的个体，再补上符合条件的个体组成新的群体。这时就可以采用适应值比例选择来进行个体的选择，这是一种基本的选择方法，其中新群体中每个个体被选择的期望数值与其适应值和群体平均适应值的比例有关，可以采用轮赌盘的方式实现。对于给定规模为n的群体，P＝(s₁，s₂...s_n)，个体适应值为f(s_j)，该个体被选择的概率为： 

$p\left(s_j\right)=\frac{f\left(s_j\right)}{\underset{i=1}{\overset{n}{\Sigma }}f\left(s_i\right)},$j＝1，2，L，n    (18) 

即计算每个个体适应值在群体适应值总和中所占的比例，来表示该个体在选择过程中被选中的概率。 

2)对群体应用染色体交叉算子： 

交叉操作随机搭配群体内的个体以某个概率交换他们的部分染色体，得到新的个体。将父辈群体中的个体复制并随机两两配对，进行交叉操作。采用多点交叉方法，选出要交配的一对个体，随机选取三个交叉位置。在交叉点之间的变量间续的相互交换，但在第一个交叉点之前的不做交换。即第一个交叉点和第二个交叉点间交换，第四个交叉点后交换，其余的不交换。交叉后产生两个新的后代，最后产生n个个体构成子辈群体。 

3)对群体应用染色体变异算子： 

变异操作模拟自然界生物体进行中染色体上面某位基因发生的突变现象，从而改变染色体的结构和物理性状。变异操作就是对于群体内的每个个体，以某个概率改变基因座上的基因值为其他等位基因，从而产生新的个体。变异概率p_m影响着群体的多样性，其取值范围设定为： 

p_m＝(e^x-1)/(e^0.7-1)，x∈{0.0，L，0.7}(19) 

这里取p_m＝0.01。在遗传算法中，变异算子通过按变异概率p_m随机反转某位等位基因的 二进制字符来实现。一般包括两个步骤： 

以原始的变异概率p_m为基础，可以计算出群体中个体发生变异的概率P_N

P_N＝1-(1-p_m)^w    (20) 

式中，w为染色体的位数，这里w＝64。计算得P_N＝0.4744，取随机变量为random，若random≤p_N，则对该个体进行变异，否则表示不发生变异。 

若已经选中了个体，则将被选中的个体每一位都按变异概率p_m变异，插入到下一代中去即可。 

通过以上步骤产生新一代群体P(i+1)，这些算子的目的在于扩展有限个体的覆盖面，体现全局搜索的思想。 

步骤E：优化修形参数，对直线修形的参数和抛物线修形的参数进行优化。 

在不同的修形参数情况下，齿轮的传递误差波动也会相应的变化。为搜索到使传递误差波动最小的修形参数值，可以在有限元模型的基础上，利用遗传算法搜索修形参数的解空间，找到最优的修形参数。设一对齿轮的参数如表5所示。此例中忽略轴承的间隙，忽略输入转矩的波动。 

表5 

1)直线修形时修形参数的优化： 

采用齿廓直线修形有限元模型进行仿真计算和用遗传算法进行搜寻，将最佳染色体解码，得到齿轮修形参数如表6所示。 

表6 

按照表6修形数据进行仿真计算，得到直线修形的齿轮传递误差动态曲线，如图29所示。图中，虚线为修形前的传递误差曲线，可见在单双齿交替时的齿轮传递误差波动较大，说明啮入啮出冲击是传递误差波动的主要原因，其最大波动为10μm，实线为按照表6定参数修形后的传递误差曲线，这时曲线已无明显的阶梯性，说明啮入啮出冲击已经不再是传递误差波动的主要原因，其最大波动减小为3.897μm，减小了61％。 

2)抛物线修形时修形参数的优化： 

采用齿廓抛物线修形有限元模型进行仿真计算和用遗传算法进行搜寻，将最佳染色体解码，得到齿轮修形参数如表7所示。 

表7 

按照表7修形数据进行仿真计算，得到抛物线修形的齿轮传递误差动态曲线，如图30所示。图中虚线为修形前的传递误差曲线，可见在单双齿交替时的齿轮传递误差波动较大，说明啮入啮出冲击是传递误差波动的主要原因，其最大波动为10μm，实线为按照表7定参数抛物线修形后的传递误差曲线，这时曲线已无明显的阶梯性，说明啮入啮出冲击已经不存在，而双齿啮合区的波峰说明单双齿交替是产生波动的主要原因。其最大波动减小为8.017μm，减小了20％，传递误差减小的幅度不大，说明直线修形较好。 

需要指出的是，虽然在本例中直线修形优于抛物线修形，但是这并不适用于所有范围，有些情况下抛物线修形会比直线好。进行修形曲线设计时，必须根据实际工况来确定修形曲线类型和曲线参数。 

步骤五：检验参数的正确性： 

由步骤四可知在本例中齿廓直线修形对齿轮传递误差波动的减弱作用要优于齿廓抛物线修形。但是这个直线修形方案只是代表了对传递误差波动的减弱效果较好。该修形方案对齿轮系统振动的减弱作用，还需要通过齿轮非线性振动模型来对其进行减振的验证。 

步骤A：进行齿轮非线性动力学建模； 

采用集中质量法建立直齿圆柱齿轮副纯扭转振动的动力学模型，如图31所示。其中，时变啮合刚度k_m和啮合阻尼c_m是沿着啮合线方向的；θ_p和θ_g是齿轮的角位移， 和 是作用在齿轮上的力矩；I_p和I_g为转动惯量，R_p和R_g为基圆半径； 为齿轮的传递误差，齿侧间隙为2b；则齿轮副的扭转振动分析模型为： 

$\left(\begin{array}{c}{I}_p\frac{d^2{\theta }_p}{d{\bar{t}}^2}+R_p{c}_m(R_p\frac{{\mathrm{d\theta }}_p}{d\bar{t}}-R_g\frac{{\mathrm{d\theta }}_g}{d\bar{t}}-\frac{d\bar{e}}{d\bar{t}})+k_m{R}_p\bar{f}(R_p{\theta }_p-R_g{\theta }_g-\bar{e})={\bar{T}}_p\left(\bar{t}\right)\\ I_g\frac{d^2{\theta }_g}{d{\bar{t}}^2}+R_g{c}_m(R_p\frac{{\mathrm{d\theta }}_p}{d\bar{t}}-R_g\frac{{\mathrm{d\theta }}_g}{d\bar{t}}-\frac{d\bar{e}}{d\bar{t}})+k_m{R}_g\bar{f}(R_p{\theta }_p-R_g{\theta }_g-\bar{e})=-{\bar{T}}_g\left(\bar{t}\right)\end{array}\right)---\left(21\right)$

式中，I_p、I_g-主、被动齿轮的转动惯量； 

-啮合齿轮的传递误差； 

R_p、R_g-基圆半径； 

c_m-啮合阻尼； 

θ_p、θ_g-主、被动齿轮的扭转振动位移； 

k_m-时变啮合刚度； 

-作用在主、被动齿轮上的外力矩； 

-具有齿侧间隙时轮齿啮合力的非线性函数。 

为无量纲化的间隙非线性描述函数，应表示为： 

$f\left(\bar{x}\left(\bar{t}\right)\right)=\frac{\bar{f}\left(\bar{x}\left(\bar{t}\right)\right)}{k_m}=\left(\begin{array}{cc}\bar{x}\left(\bar{t}\right)-b& \bar{x}\left(\bar{t}\right)>b\\ 0& -b<\bar{x}\left(\bar{t}\right)<b\\ \bar{x}\left(\bar{t}\right)+b& \bar{x}\left(\bar{t}\right)<-b\end{array}\right)---\left(22\right)$

式中， -传递误差。 

为了求解的方便，可以将式(21)简化为单自由度动力学模型。在没有加工、安装误差的情况下，齿轮处于理想啮合状态，齿轮在传动过程中的弹性变形产生的相对变形称为动态传递误差。由于主、被动齿轮沿啮合线方向的振动位移分别为R_pθ_p和R_gθ_g，则系统的动态传递误差 为： 

${\bar{x}}_d\left(\bar{t}\right)=R_p{\theta }_p\left(\bar{t}\right)-R_g{\theta }_g\left(\bar{t}\right)---\left(23\right)$

实际上，齿轮副实际传递误差是动态传递误差和传递误差之差，定义为 则： 

$\bar{x}\left(\bar{t}\right)={\bar{x}}_d\left(\bar{t}\right)-\bar{e}\left(\bar{t}\right)---\left(24\right)$

设 m_p和m_g称为齿轮的当量质量，将式(23)、(24)代入(21)整理得： 

$\frac{m_p{m}_g}{m_p+m_g}\frac{d^2\bar{x}}{d{\bar{t}}^2}+c_m\frac{d\bar{x}}{d\bar{t}}+k_m\bar{f}\left(\bar{x}\right)=\frac{m_p{m}_g}{m_p+m_g}(\frac{{\bar{T}}_p\left(\bar{t}\right)}{R_p{m}_p}+\frac{{\bar{T}}_g\left(\bar{t}\right)}{R_g{m}_g}-\frac{d^2\bar{e}}{d{\bar{t}}^2})---\left(25\right)$

设齿轮的等效质量$m_e=\frac{m_p{m}_g}{m_p+m_g},$主动齿轮上作用的啮合力${\bar{F}}_p\left(\bar{t}\right)=\frac{{\bar{T}}_p\left(\bar{t}\right)}{R_p},$从动齿轮上作用的啮合力${\bar{F}}_g\left(\bar{t}\right)=\frac{{\bar{T}}_g\left(\bar{t}\right)}{R_g},$则上式变为： 

$m_e\frac{d^2\bar{x}}{d{\bar{t}}^2}+c_m\frac{d\bar{x}}{d\bar{t}}+k_m\bar{f}\left(\bar{x}\right)=\frac{m_e}{m_p}{\bar{F}}_p\left(\bar{t}\right)+\frac{m_e}{m_g}{\bar{F}}_g\left(\bar{t}\right)-m_e\frac{d^2\bar{e}}{d{\bar{t}}^2}---\left(26\right)$

齿轮传动系统的激励分为两类，对于由旋转质量的不平衡、几何偏心、原动机和负载力矩波动引起的频率较低的激励，一般处理成外部激励，而将由加工误差、轮齿弹性变形等引起的频率较高的激励处理成内部激励，令 

$\bar{F}\left(\bar{t}\right)=\frac{m_e}{m_p}{\bar{F}}_p\left(\bar{t}\right)+\frac{m_e}{m_g}{\bar{F}}_g\left(\bar{t}\right)-m_e\frac{d^2\bar{e}\left(\bar{t}\right)}{d{\bar{t}}^2}{\bar{F}}_a\left(\bar{t}\right)+{\bar{F}}_h\left(\bar{t}\right)---\left(27\right)$

式中， 为总激励； 

为外部激励力，${\bar{F}}_a\left(\bar{t}\right)=\frac{m_e}{m_p}{\bar{F}}_p\left(\bar{t}\right)+\frac{m_e}{m_g}{\bar{F}}_g\left(\bar{t}\right);$

为内部激励力，${\bar{F}}_h\left(\bar{t}\right)=-m_e\frac{d^2\bar{e}\left(\bar{t}\right)}{d{\bar{t}}^2}.$

则式(26)表示为： 

$m_e\frac{d^2\bar{x}}{d{\bar{t}}^2}+c_m\frac{d\bar{x}}{d\bar{t}}+k_m\bar{f}\left(\bar{x}\right)=\bar{F}\left(\bar{t}\right)---\left(28\right)$

式(28)即为齿轮副含齿侧间隙和时变啮合刚度的非线性振动微分方程，它是单自由度、变参数的强非线性二阶微分方程。在国际单位制下，微分方程中各物理量之间的数量级相差很大，如刚度系数的数量级在10⁷～10⁹之间，阻力系数的数量级在10²～10⁴之间，而振动响应位移的数量级则在μm级，这使得对方程求解计算带来很大困难；因此，一般都要对像式(28)这类微分方程进行归一化无量纲处理。再者，由于强非线性微分方程，一般都要用数值方法求解，在进行数值求解时，如果同一方程中各量的数量级相差很大，则使误差控制值和步长值难以选择。由式(28)可知，齿轮副动力学模型的固有频率为 式中，k₀为时变啮合刚度的平均值，m_e为等效质量，定义无量纲时间为 则： 

$\frac{d\bar{x}\left(\bar{t}\right)}{d\bar{t}}=\frac{d\bar{x}\left(\bar{t}\right)}{\mathrm{dt}}·\frac{\mathrm{dt}}{d\bar{t}}=\frac{d\bar{x}\left(\bar{t}\right)}{\mathrm{dt}}{\omega }_m---\left(29\right)$

$\frac{d^2\bar{x}\left(\bar{t}\right)}{d{\bar{t}}^2}=\frac{d}{d\bar{t}}\left[\frac{d\bar{x}\left(\bar{t}\right)}{d\bar{t}}\right]=\frac{d}{\mathrm{dt}}{\omega }_n\left[\frac{d\bar{x}\left(\bar{t}\right)}{\mathrm{dt}}{\omega }_n\right]=\frac{d^2\bar{x}\left(\bar{t}\right)}{d{t}^2}{\omega }_n^2---\left(30\right)$

将式(29)、(30)带入式(28)，则有： 

$m_e{\omega }_n^2\frac{d^2\bar{x}}{d{t}^2}+c_m{\omega }_n\frac{d\bar{x}}{\mathrm{dt}}+k_m\left(\bar{t}\right)\bar{f}\left(\bar{x}\right)=\bar{F}\left(\bar{t}\right)---\left(31\right)$

设齿侧间隙的一半b为位移标称尺度，将 代入式(31)，得： 

$\frac{d^{2}x}{d{t}^2}+\frac{c_m}{m_e{\omega }_n}\frac{\mathrm{dx}}{\mathrm{dt}}+\frac{k_m\left(\bar{t}\right)\bar{f}\left(\bar{x}\right)}{m_e{\omega }_n^{2}b}=\frac{\bar{F}\left(\bar{t}\right)}{m_e{\omega }_n^{2}b}---\left(32\right)$

设各变量相应的无量纲量分别为： 

$\zeta =\frac{c_m}{2m_e{\omega }_n},k\left(t\right)=\frac{k_m\left(\bar{t}\right)}{m_e{\omega }_n^2},f\left(x\right)=\frac{\bar{f}\left(\bar{x}\right)}{b},F\left(t\right)=\frac{\bar{F}\left(\bar{t}\right)}{m_{e}b{\omega }_n^2}---\left(33\right)$

则可得到无量纲形式的齿轮副分析模型为： 

$\stackrel{··}{x}\left(t\right)+2\zeta \stackrel{·}{x}\left(t\right)+k\left(t\right)f\left(x\left(t\right)\right)=F\left(t\right)---\left(34\right)$

式中，f(x(t))为无量纲化的齿侧间隙非线性函数。 

$f\left(x\left(t\right)\right)=\frac{\bar{f}\left(\bar{x}\left(\bar{t}\right)\right)}{b}=\left(\begin{array}{cc}x\left(t\right)-1& x\left(t\right)>1\\ 0& -1\le x\left(t\right)\le 1\\ x\left(t\right)+1& x\left(t\right)<-1\end{array}\right)---\left(35\right)$

步骤B：对齿轮非线性动力学方程进行求解； 

齿轮副含间隙和时变啮合刚度的非线性动力学方程的求解方法分为解析求解法和数值仿真法。解析法主要是谐波平衡法、多尺度法、渐进法等，由于受到所假设的激励和响应形式的限制，解析法的计算精度不太理想，而且用解析法虽然能在弱非线性情形下求解非线性微分方程，但当系统处于强非线性时，则很难求解，所以采用数值法求解。 

步骤C：进行最优方案的减振论证。 

无量纲频率 说明无量纲频率就是频率比，计算得修形前的参数为ω_e＝0.6929，ξ＝0.02，ε＝0.283，φ_k＝0，φ_e＝π，F_m＝0.1，e₁/b＝0.05；修形后的参数为ω_e＝0.7，ξ＝0.02，ε＝0.0962，φ_k＝0，φ_e＝π，F_m＝0.1，e₁/b＝0.015。为绘制幅频曲线奠定了基础。 

若初始条件为x(0)＝0， 以频率比ω_e作为横坐标，分别计算修形前和修形后系 统的幅频曲线，如图32所示，其中，虚线为修形前的幅频曲线，实线为按照表6所示数据修形后的幅频曲线。由图可见，修形前系统在ω_e＝1处的振幅达到了2.8×10^-3。修形后，振幅减小到0.3×10^-3左右，可见得到的修形参数有效的减小了齿轮扭转振动振幅，可见，修形参数的正确性得到了有力的证明。