人字齿轮最佳齿向修形设计的Kriging响应面法

余澍民

(中国电子科技集团第三十八研究所，安徽 合肥 230088)

0 引言

人字齿轮因其重合度高、结构紧凑等优点，广泛用于船舶和航空等各个传动领域[1]；但在人字齿轮传动中由于误差和变形不可避免，会使系统产生振动和噪声影响系统使用寿命和可靠性，而适当的修形可以有效地改善齿面载荷分布不均，降低受载产生的应力集中。因此修形在高速重载人字齿轮传动中显得尤为重要。

Wu等[2]通过有限元静力学分析确定最优齿廓修形类型及修形参数，通过动态接触试验证明齿廓修形在减振方面的作用，但未涉及到齿向修形理论。宋乐明等[3-4]展开了齿向修形的理论研究与实验验证，得出鼓形设计的齿轮即使在歪斜度较大的情况下也能保证接触端面相切不相交，避免棱边效应。但是并未考虑到由于制造、安装误差和变形的存在，两端斜齿轮传动时左右不完全对称造成两端传递扭矩的不同，导致两端受载变形的扭转角不相同，在对人字齿轮进行鼓形修形时，应对两侧斜齿修形量单独考虑。袁哲等[5]采用有限元计算时变啮合刚度，并结合遗传算法进行齿轮修形优化，取得了较好的减振优化效果。虽然基于有限元法的齿轮修形优化精度高，但因为修形所需要的大量样本通过有限元法计算庞大且操作烦琐，每一次样本需手动调节不同啮合位置，通过不同啮合位置结果拟合曲线然后由傅里叶展开得到时变啮合刚度函数[6-8]。为减少计算量，优化模型有必要对齿轮修形做进一步研究。

本文针对人字齿轮传动，对人字齿轮两侧鼓形修形量分开研究，采用BP神经网络与Kriging多响应预测模型结合，建立时变啮合刚度各参数的响应面，基于人字齿轮弯扭轴耦合动力学模型，以动载系数最小来确定两侧最佳鼓形修形参数，并与传统ISO修形方法对比，确定了本文最佳修形设计方法的有效性。

1 齿向修形理论

1.1 齿向修形方法

人字齿轮齿向修形主要包括鼓形修形、齿端修形和螺旋线修形。由于齿端修形只能进行端面附近的局部修形，而齿轮传动中由于误差和变形的存在，在端面容易造成应力集中，影响修形效果；螺旋线修形微量改变螺旋角的大小，其效果比齿端修形要好但是由于改变的角度很小，不易加工控制；鼓形修形是采用一种等半径圆弧来设计修形，可以有效地补偿齿轮制造误差和受载变形，如图1所示。本文主要研究的是鼓形修形[9]。

图1 鼓形修形示意

由于人字齿轮两侧传递扭矩不同导致变形的扭转角不一致，对人字齿轮两端修形量单独考虑，Δ1和Δ2分别为输出端斜齿和输入端斜齿鼓形修形量。

在齿轮手册[10]，鼓形修形量计算公式

Δ=0.5(fsh+1.5fHβ)+5～10

(1)

fsh,fHβ分别为综合变形产生的啮合齿向误差分量和螺旋线倾斜极限偏差;Δ为鼓形量。

对于高精度、高可靠度的高速齿轮，鼓形修形量Δ为上述计算值的60%～70%,即鼓形修形量为10≤Δ≤25 μm，加上制造误差5 μm。

2 预测啮合刚度的Kriging改进算法

为了探究人字齿轮修形参数与刚度的关系，建立Kriging预测模型。Kriging预测模型适用于强非线性空间估计问题，在工程优化领域得到广泛应用[11]。但是在求解过程中发现，当已知样本点达到75组以上时，Kriging预测模型才有实际意义。由于用有限元法获取样本复杂。为此对数据进行预处理，即采用BP神经网络通过20组训练样本数据预测出75组样本数据，为Kriging预测模型提供样本，大大减小了计算量。

2.1 有限元综合啮合刚度计算与数据处理

采用有限元法计算人字齿轮时变啮合刚度，人字齿轮参数如表1所示。

在1个齿轮啮合周期内均匀提取10个位置，可以转换为通过旋转相同的角度得到10个不同啮合位置的模型，分别导入ANSYS软件中给主动轮施加转矩T=5 094 N·m,从动轮约束6个自由度，求得从动轮啮合线上的变形量l，由此可求出人字齿轮综合啮合刚度。

表1 人字齿轮参数

(2)

(3)

(4)

θ为从动轮变形角度；ra为从动轮齿顶圆半径；rb为从动轮基圆半径；T为主动轮施加的转矩;αt为端面压力角;βb为基圆螺旋角。

由于齿轮刚度激励通常采用傅里叶变换，取其一阶分量[12]，即

(5)

A为综合啮合刚度的半均值;ε为比例系数；φ为啮合相位，只要确定了这3个参数，齿轮振动微分方程的刚度激励也就确定。

由1.1节知，两侧鼓形量Δ1和Δ2的取值范围都在[10,25] μm.通过拉丁超立抽样方法抽样，在输出端鼓形量Δ1区间[10,25] μm、输入端鼓形量Δ2区间[10,25] μm进行抽样，抽样75个点，其中前20组为BP神经网络训练样本，后55组为BP预测样本，其中5组为检验样本。

利用训练好的网络进行预测，并抽取5个预测的样本，用有限元方法求解，求解后的数据如表2所示。

对BP神经网络预测结果与ANSYS有限元计算结果进行Pearson相关性分析，分析结果如表3所示。

由表3可见，Sig<0.01，A,ε,φ相关系数分别为0.918，0.929，0.924，查表3可得结论为极强相关，故BP神经网络预测数据是合理的。

表2 BP预测样本与有限元计算比较

表3 Pearson相关性分析

2.2 Kriging预测模型

Kriging预测模型采用Gaussion随机过程函数，在某一点的预测是通过空间上已知所有点的数据加权求和得来的[13]。真实响应值是由回归多项式和随机过程2个部分组成，即

f(x)Tβi,l+zl(x),l=1,…,q

(6)

f(x)=[f1(x),f2(x),…,fp(x)]T为回归多项式基函数向量;{βi,l}为多项式参数;zl(x)为第l个响应分量对应的Gaussian随机过程函数协方差矩阵，表示为

(7)

w=[w1,w2,…,wn]T和x=[x1,x2,…,xn]T分别为2个不同的随机变量；R(θ,w,x)为带有参数θ的相关函数，Gaussian相关函数R(θ,w,x)可以表示为

(8)

考虑Kriging线性预测，未知点x的预测值可以表示为

(9)

c=c(x)=[c1(x),c2(x),…,cm(x)]T为待求权系数向量;Yi,j为试验设计点S=[s1,s2,…,sm]T处对应的响应值。

为保证预测值无偏和预测均方误差最小，即采用Lagrangian求最小值，可得

(10)

将BP得到的样本数据导入Kriging预测模型中，分别得到A,ε,φ的预测面，利用拉丁抽样抽取5组数据，并采用ANSYS有限元法计算再与其对比，如图2所示。

图2 各响应值Kriging预测面与有限元结果对比

表4为ANSYS有限元法和Kriging多响应预测法的结果，Pearson相关性分析结果如表5所示，发现在Sig<0.01的情况下，A,ε,φ相关性达到0.988，0.948，0.998查表为极强相关，证明采用Kriging多响应预测模型的精确性。

表4 Kriging多响应预测与ANSYS有限元计算

表5 Pearson相关性分析

3 人字齿轮动力学模型求解

3.1 人字齿轮动力学方程建立

本文建立一对人字齿轮传动系统的三维离散模型，如图3所示。对于单个人字齿来说，将其分为左右斜齿和中间退刀槽部分，每一部分具有径向、轴向、轴线扭转共4个自由度，加入输入盘和输出盘2个自由度，利用梁单元法将左右两边斜齿耦合起来获得人字齿轮传动系统共26个自由度的动力学方程。

图3 人字齿轮动力学模型

(11)

系统的阻尼由比例阻尼矩阵表示[14]，为了推导的方便，在公式中不列举。推导出的人字齿轮副中单侧斜齿轮啮合副j(j=L,R)的弯扭轴耦合线性动力学方程为

(12)

Jin,Jout分别为输入和输出盘的转动惯量;θin,θout分别为输入输出盘的扭转位移;θpM,θgM分别为主、从动齿轮中间连接部分的扭转位移量;kin,kout分别为输入盘与主动轮、输出盘与从动轮之间的扭转刚度，cg=cos(g),sg=sin(g)，齿轮副δj相对啮合位移为

δj=[(xpj-xgj)sψ+(ypj-ygj)cψ+rpjθpj+

rgjθgj]cβ+(zpj-zgj)sβ-ej

(13)

ej为齿轮副的静态传递误差，这里取5 μm。

人字齿轮左右两侧斜齿轮通过中间轴段固连在一起，利用有限元理论，采用Euler梁单元等效人字齿轮中间轴段，每个节点具有4个自由度，Euler梁质量矩阵Me和刚度矩阵Ke如下，具体的矩阵见文献[15]。

(14)

(15)

图4为人字齿轮三片式模型，按照L、M、R节点组装，中间节点通过轴承或花键刚度矩阵，将齿轮和支撑结构相连。

图4 人字齿轮三片式模型

连接结构的刚度和质量矩阵为

(16)

(17)

ζ=p,g，下标e1,e2分别表示梁单元1和梁单元2。

系统的整体方程为

(18)

M为质量矩阵;q为位移矢量;K为整个系统的刚度矩阵同，由啮合刚度阵和支承刚度矩阵组成;F为力矩阵;C为阻尼矩阵，定义为

C=μK+vM

(19)

μ和v分别为比例系数和阻尼系数。

3.2 不同修形参数下的人字齿轮动载系数

由3.1节可求出不同修形下人字齿轮副的动态响应，即可得到不同修形下的动载系数，其中动载系数定义为动态啮合力与静态啮合力之比，即Kpg=Fpg/Fn，取每次稳态下最大动载系数建立响应面如图5所示。

图5 不同修形参数下的人字齿轮动载系数

通过响应面可以得到在某些修形参数下，如输出端修形量Δ1=20.253 2 μm，输入端修形量Δ2=24.6 203 μm时，动载系数达到最大值为1.654 3，这是因为人字齿轮传动时由于制造和误差所引起的齿面载荷从动力输入端至输出端逐渐降低，且随着齿宽的增加，两侧载荷相差增大，如果两侧修形量都偏大，会造成输出端斜齿修形过大，导致齿面接触相交，产生应力集中，动载系数变大；而当输出端修形量Δ1=12.468 4 μm，输入端修形量Δ2=20.063 3 μm时，动载系数最小为1.140 4，此时修形量能够有效补偿输入端到输出端由于变形产生的不同啮合歪斜度，使得齿轮传动时接触面相切而不相割，大大减少了动态啮合力。为了验证本文修形参数的正确性，与传统的ISO修形方法对比，结果如表6所示。

表6 两种修形方法动载系数对比

由表6可知，未修形下动载系数幅值为1.942 3，采用传统的ISO修形方法得到的动载系数幅值为1.342 3，而本文修形方法的动载系数幅值为1.140 4。

本文修形方法相比于未修形下动载系数下降了(1.942 3-1.140 4)/1.942 3×100%=41.29%。相比于ISO修形下动载系数下降了(1.342 3-1.140 4)/1.342 3×100%=15.04%。从而证明本文修形方法大大减少了未修形下的动载系数幅值，且比ISO方法效果更好。

4 结束语

利用BP神经网络快速预测人字齿轮啮合刚度；Kriging多响应面建立全局各参数响应面保证能够在全范围内得到最佳齿向鼓形量，实现人字齿轮齿向修形优化。

采用人字齿轮弯扭轴耦合动力学模型求解不同鼓形下的动载系数并建立响应面，结果表明两侧最佳鼓形量并不一致，输入端斜齿鼓形量为20.063 3 μm，输出端斜齿鼓形量为12.468 4 μm，两侧不同的修形量能够有效地补偿两端由于变形产生的不同啮合歪斜度，得到最优修形效果。

本文所得修形参数带入动力学求解动载系数为1.140 4，相比较ISO修形下动载系数为1.342 3和未修形参数下动载系数1.942 3，本文所得最佳修形参数更为优越，在人字齿轮减振，降噪等方面有重要意义。

文中采取的响应面优化算法，可以在全局参数下得到最优修形参数，但是响应面会随着工况、齿轮副参数改变而改变，不能一概而论，且提高动力学求解效率是响应面得到的关键一步，有待进一步优化。