圆弧齿线圆柱齿轮接触应力预测模型研究*

童 钦，张 祺，梁 君，侯 力

(1.绵阳师范学院 机电工程学院，四川 绵阳 621000;2.攀枝花学院 智能制造学院，四川 攀枝花 617000;3.四川大学 机械工程学院，四川 成都 610065)

0 引 言

齿轮被广泛应用于工业的各个行业[1]，其中，渐开线圆柱齿轮是应用最广泛的圆柱齿轮。但渐开线圆柱齿轮中，直齿圆柱齿轮的承载力有限；而斜齿轮则存在轴向力，人字齿轮的加工工艺又过于复杂，故渐开线齿轮都存在种种的不足[2]。

为解决渐开线齿轮的不足，日本学者长谷川吉三郎等人[3]提出了一种新型圆弧齿线齿轮传动装置。这种新型的齿轮传动具有啮合性能好、重合度高、无轴向力、传动平稳等优点。针对这种齿轮，学者们也进行了深入的研究。Tseng等人[4-7]先后对其数学模型、根切条件、接触等方面进行了研究；曹磊等人[8]对圆弧齿线齿轮的精确三维建模方式进行了研究；马登秋等人[9]研究了其接触分布；陈帅等人[10]研究了其动态接触；王虹等人[11-13]针对圆弧齿线齿轮的接触性能进行了研究。

为研究圆弧齿线齿轮设计参数对其接触应力的影响，笔者提出利用Kriging代理模型来建立圆弧齿线齿轮的设计参数(齿宽、模数、压力角、齿线半径)与圆弧齿轮接触应力的数学模型；为提高Kriging模型的建模精度，笔者采用鲸鱼优化算法(WOA)对传统Kriging模型变异函数的参数进行优化。

1 圆弧齿轮数学模型

根据其成形原理[14-16]，圆弧齿线圆柱齿轮的成形原理坐标系如图1所示。

图1 圆弧齿线圆柱齿轮成形原理坐标系R—加工齿坯的分度圆半径；RT—加工刀具的刀盘半径；Rn—加工刀具分度圆切线方向上内刃半径；Rw—加工刀具分度圆切线方向上外刃半径，Rw=RT+πm/4；m—齿轮模数；ω—加工刀具的刀盘旋转角速度；VT—加工刀具的刀盘移动速度

图1中，S(O—XYZ)为静止坐标系，S1(O1—X1Y1Z1)为齿坯固化坐标，ST(OT—XTYTZT)为刀具坐标；它以VT=R×ω的速度相对于S(O—XYZ)坐标移动。

在ST(OT—XTYTZT)坐标系统中，刀具切制过程中所形成的曲面参数方程的矢量表达式为：

(1)

式中：a—加工刀具切削刃的压力角；θ—齿坯当前包络点到中间截面的转角；q—刀盘坐标系中刀具切削面上的点沿切削锥面母线方向距离包络参考点的长度。

将ST(OT—XTYTZT)变换到S(O—XYZ)，其坐标转换关系为：

(2)

将S(O—XYZ)变换到S1(O1—X1Y1Z1)，其坐标转换关系为：

(3)

如果将ST(OT—XTYTZT)转换到S(O—XYZ)，则坐标变换矩阵为：T1T=T10T0T。

1.1 刀具曲面的单位法向矢量

刀具曲面方程为：

(4)

(5)

1.2 刀具和齿轮在啮合点的相对速度

由于：

(6)

则刀具与齿轮在啮合点处的相对速度为：

(7)

1.3 啮合函数Г

基于啮合原理，啮合函数表示如下：

(8)

因为ω1≠0，由上式可以得出:

(9)

1.4 共轭曲面

将式(9)代入式(1)中，则刀具与被加工齿轮接触线方程为：

(10)

将ST(OT—XTYTZT)中的坐标转换到S1(O1—X1Y1Z1)后，可以获得所切制齿轮的齿面方程为：

(11)

1.5 瞬时接触线

啮合函数Г=Г(q,θ,φ)是参数变量q,θ,φ的函数。因为φ1=ω1t，在此将φ1看作是某一瞬时代入Г的一个函数，从而可得到此刻齿条刀具和齿坯的瞬时接触线的函数表达式。

如果把任意时刻的φ1值代入式(7)，就能得到整个表面的接触线方程为：

(12)

1.6 齿廓方程

在齿轮轴向中间截面，由其展成坐标系可知：b=0，那么θ也为零。将这两个值代入式(11)中，得到中间截面的齿廓方程为：

(13)

根据式(13)可以看出，该齿轮轴向对称面的齿廓形状为渐开线。

同理，在轴向非对称面上，令z1=b，根据式(11)中Z1和q的表达式，可以得到：

φ={[-b/tanθ+(RT±πm/4)cosθ]/sin2α-
(RT+πm/4)cosθ-RT}/R

(14)

从而可得到非中间截面的齿廓的表达式为：

(15)

2 基于有限元的圆弧齿线圆柱齿轮接触应力分析

基于有限元的圆弧齿线圆柱齿轮接触应力分析步骤如下：

(1)模型材料属性定义

在ABAQUS中建立材料信息，如弹性模量E=2.08 MPa，泊松比=0.298，材料17CrNiMo6。

(2)建立分析步和相互作用

建立分析步，主要包括：定义分析类型(static)、定义分析增量步、确定迭代方法，创建场变量与历史变量，并确定输出参数，开启非线性，定义相互作用为接触，设置圆弧齿线齿轮副有限元分析的接触类型为“无摩擦”；主动与从动以各自旋转中心建立MPC(multi-point constrain)约束。

(3)施加约束和力矩

在主动轮上施加扭矩，其大小为6.08×104N·mm；同时，对主动轮添加MPC约束，并将其旋转轴方向设置为自由，其他旋转和平移设置为固定。对于从动轮而言，MPC固定约束,全向固定。

(4)网格划分

笔者采用扫掠方式对齿轮副进行网格画分。网格类型采用C3D8I，齿轮整体的单元大小设置为2 mm,接触区域进行局部细分，大小设置为0.02 mm。如果不满足要求，则要继续调整参数。在后续计算中，可以进一步细化网格，提高划分质量，并通过多次试算对比分析结果；若计算结果变化不大，则该网格可作为最终分析网格。

分析过程中，应对接触区域进行局部细分，每对齿轮的网格数量大概在1.2×106个左右。

齿轮网格划分结果如图2所示。

图2 齿轮网格划分结果

(5)求解与可视化

考虑每对齿轮的网格数量大概在1.2×106个左右，故求解时采用并行计算。求解后的主动轮和从动轮的应力云图如图3所示。

由图3的仿真分析结果可以看出：接触区域为靠近分度圆附近；主动轮的接触区域该在分度圆以上，且最大应力值为503.3 MPa；从动轮的接触区域该在分度圆以下，且最大应力值为501 MPa，故接触应力的最大值为503.3 MPa；且主动轮与从动轮之间的差值为2.3 MPa，仅为接触应力最大值的0.457 0%，几乎一致，可以忽略。

同样，通过图3可以看出，在载荷的作用下，笔者所研究的齿轮其接触区域为椭圆，正确地印证了点接触齿轮在载荷的作用下接触区域为椭圆的事实。

图3 主动轮和从动轮的应力云图

3 基于代理模型的齿线圆柱接触应力预测模型

3.1 试验设计与响应结果

常用的试验设计方法主要有均匀试验设计、正交试验设计、拉丁方试验设计等方法，本文研究了齿宽、模数、压力角、齿线半径与接触力之间的关系，利用回归正交试验法原理设计仿真方案选择了4个因素3个水平,齿轮设计因素水平表如表1所示。

表1 齿轮设计因素水平

本文采用正交试验进行了样本抽样，然后利用有限元方法得到了不同样本数据下的齿轮的接应力，根据正交表可知，抽样的样本为9个样本,L9(34)正交试验圆弧齿轮接触应力仿真结果如表2所示。

表2 L9(34)正交试验圆弧齿轮接触应力仿真结果

3.2 改进Kriging模型

3.2.1 Kriging代理模型

在解决非线性程度较高的问题时，Kriging(克里金)模型可较容易地获得理想的拟合结果，其插值结果定义为已知样本函数响应值的线性加权，即：

(16)

式中：fj(x)—函数，一般为多项式；βj—相对应的系数；Z(x)—静态随机过程，其满足均值为0，方差为σ2。

且对于设计空间内，不同两点处所对应的随机变量之间的协方差为：

Cov[Z(xi),Z(xj)]=σ2R(xi,xj)

(17)

(18)

式中：R(xi,xj)—相关性函数，它表示不同位置处随机变量之间的相关性，常用的相关性函数为高斯型函数。

式(18)中，θ为Kriging模型的变差函数的参数，其大小通过极大似然估计法求解优化问题的方式来确定，即：

(19)

为保证Kriging预测值与真实函数值之间的均方根误差(RMSE)最小，Kriging模型的近似表达式为：

(20)

3.2.2 鲸鱼算法

鲸鱼算法的原理来自于座头鲸的“泡泡网”觅食行为，基于这一特殊捕食策略的数学表达式如下：

D=|C·X*(t)-X(t)|

(21)

X(t+1)=X*(t)-A·D

(22)

式中：t—当前迭代次数；X(t)—当前一座头鲸的坐标向量；X(t+1)—下一次迭代后的目标坐标向量；X*(t)—到目前得到的最佳位置向量，它将随时间不断更新；D—当前这条座头鲸和最佳位置之间的距离。

上式(21～22)中，A和C是系数，其表达式分别表示为:

A=2a·r-a

(23)

C=2r

(24)

式中：a—在值域[0,2]上并随迭代时间线性递减的参数；r—区间[0,1]内的随机向量。

当|A|>1时，对应鲸鱼群的游走觅食行为。利用种群的随机个体坐标Xrand来定位导航寻找食物,其数学表达式如下：

X(t+1)=Xrand(t)-A·D

(25)

当|A|<1时，对应鲸鱼群的包围捕食和攻击猎物这两种行为，其数学模型描述如下：

X(t+1)=X*(t)-D·ebl·cos(2πl)

(26)

式中：b—与螺旋形状的常数；l—区间[-1,1]上的随机数。

由于鲸鱼的收缩包围机制和螺旋更新位置是一种同步行为，笔者在数学上选取概率相同方式来对其进行位置更新，于是可以得到以下表达式：

(27)

根据基于Kriging和WOA算法，基于WOA算法的Kriging代理模型改进流程如图4所示。

图4 基于WOA算法的Kriging代理模型改进流程

4 数值仿真

笔者利用Matlab数字仿真平台，分别利用Kriging和基于WOA算法改进的Kriging算法，建立了圆弧齿轮接触应力的预测模型。

优化前后，相关系数(R2)均方根误差(RMSE)以及相对最大绝对误差(RMAE)评价指标情况如表3所示。

表3 优化前后评价指标情况

从表3中可以看出：

(1)相关系数(R2)从0.992 2提高到了0.997 4，提高了0.52%，且更加接近于1；优化后的Kriging模型全局近似能力更好(RMSE用于表示估计值与真实值之间减值大小)；

(2)优化后，RMSE从2.856 9降低到1.654 0，降低了42.11%，则说明优化后的Kriging模型能够更好地对样本进行估计；

(3)优化后，RMAE从0.132 2降低到0.075 4，降低了42.97%，且更加接近于0，则说明其优化后减少了局部误差。

综合来看，采用WOA算法对Kriging算法进行改进，提升了Kriging算法的拟合能力和精度。

利用Kriging和基于WOA算法改进的Kriging算法，笔者分别建立了圆弧齿轮接触应力预测模型的测试集残差图，分别如图5、图6所示。

图5 基于Kriging的圆弧齿轮接触应力预测模型残差图

图6 基于WOA-Kriging的圆弧齿轮接触应力预测模型残差图

从图5、图6中可以看出：

(1)基于传统Kriging算法和基于WOA算法改进的Kriging算法均能很好地建立圆弧齿轮接触应力预测模型，精度都在可使用的范围之内；

(2)Kriging的误差范围在[-2,4]，而基于WOA算法改进的Kriging算法误差范围在[0,3]，由此可见，采用改进的Kriging算法，其精度得到了明显的提高。

5 结束语

为研究圆弧齿线圆柱齿轮设计参数对其接触应力的影响，笔者利用Kriging代理模型建立了齿轮的齿宽、模数、压力角、齿线半径与圆弧齿轮接触应力的数学模型，同时，提出了一种基于鲸鱼优化算法(WOA)的Kriging模型建模方法，通过鲸鱼优化算优化传统Kriging模型变异函数的参数，提高Kriging模型的建模精度。

研究所得到的结论如下：

(1)基于WOA算法改进的Kriging算法，相关系数(R2)从0.992 2提高到了0.997 4，MSE从2.856 9降低到1.654 0，RMAE从0.132 2降低到0.075 4，改进后的算法其相关系数(R2),均方根误差(RMSE)以及相对最大绝对误差(RMAE)均得到了不同程度的改良，提高了传统Kriging算法的全局近似能力，减少了局部误差，提升了拟合精度；

(2)基于传统的Kriging算法和基于WOA算法改进的Kriging算法均能很好地建立圆弧齿轮接触应力预测模型，精度都在可使用的范围之内，但是Kriging的误差范围在[-2,4]，而基于WOA算法改进的Kriging算法误差范围在[0,3]；精度得到了明显的提高，能够为其优化设计参数和各参数对其接触应力影响的显著性分析，提供更加精确的数学模型。