汪 建， 张 俊

(福州大学 机械工程及自动化学院， 福州 350116)

传递误差是引起齿轮传动系统振动与噪声的主要原因[1-2]。工程上可通过齿轮修形来改善齿面受力、降低传递误差波动，进而抑制传动系统的振动和噪声[3-5]。因此，明确轮齿修形参数和齿轮副传递误差波动量之间的映射关系，对于指导齿轮修形、改善传动性能具有重要的理论与工程意义。

为制定合理的齿轮修形方案，学术界和工业界开展了大量研究。王成等[6]建立了2自由度单级齿轮传动非线性动力学模型，对比研究了齿廓修形参数对齿轮动态特性的影响。张俊等[7]运用多体方法建立了风电齿轮箱的动力学模型，探讨了一种面向最小动态传动误差波动量的风电齿轮修形方法。Bonori等[8]尝试将遗传算法和齿轮有限元分析相结合，用来确定最优的修形参数。遵循类似思路，袁哲等[9]以减小齿轮传递误差波动为目标,对齿轮修形参数进行了优化设计。

然而，由于齿轮系统具有复杂的非线性，理论上难以给出修形参数和传递误差之间的解析表达。为此，部分学者采用数值拟合方法给出两者间的关系。Park[10]运用响应面法建立线性模型来预测修形对静态传递误差的影响。Zhang等[11]运用响应面法拟合出修形量与动态传递误差波动量间的函数关系，并运用蒙特卡洛法对齿轮修形的可靠性进行评估。方宗德等[12]利用回归分析方法，给出了最佳修形量和反映修形效果的传动误差幅值比的回归公式。Zhang等[13]综合运用Kriging法和遗传算法，开展了大型球磨机齿轮传动系统的可靠性设计。

文献检索表明，响应面法和Kriging法是目前较为常用的两种拟合方法，但二者在拟合特定场合关系函数上的优劣尚无定论。有鉴于此，本文以某兆瓦级风电齿轮箱中行星级传动的外啮合斜齿轮副为对象，分别采用响应面法和Kriging法构建给定工况下齿轮修形参数和齿轮副传递误差波动量间的映射函数。在此基础上，将两种拟合方法所得的响应函数结果与Romax数值仿真结果对比，以评估响应面法和Kriging法的拟合精度和修形参数对传递误差波动的灵敏度。希冀通过上述研究，探索出用于构建齿轮传递误差与修形参数间关系函数的高效拟合方法，据此优选修形参数，为后续的减振降噪和传动可靠性研究提供理论依据。

1 基准修形参数的确定

不失一般性，不妨以图1所示的某兆瓦级风电齿轮箱中行星级传动的外啮合斜齿轮副为例，研究齿轮修形量与齿轮副动态传递误差波动量之间的关系。

(a)行星级传动 (b) 太阳轮-行星轮啮合副图1 某兆瓦级风电齿轮箱行星级及其外啮合齿轮副Fig.1 S-P pair in the planetary stage of a MW wind turbine

图1所示的斜齿行星传动中，太阳轮制成一端浮动的齿轮轴，行星轮采用两端对称布置的圆柱滚子轴承支承，而内齿圈与箱体固结。其中，行星架作为输入构件(输入转速为n=29 r/min)，太阳轮为输出构件(输出转矩为To=18 990 N·m)。传动齿轮为ISO 5级精度，模数mn=10 mm、齿宽B=220 mm、法向压力角αn=20°、端面压力角αt=20.158°、螺旋角β=7.493°、变位系数分别为xs=0.186，xp=0.39。

不妨以第1对外啮合齿轮副(s-p1)为例，开展齿轮箱的负载轮齿接触分析(LTCA)。为叙述方便，下文略去表征行星轮的下标1。借助Romax软件的Designer模块，可仿真获知s-p外啮合副沿啮合线方向的综合错位量以及太阳轮齿面接触单位长度载荷分布，其结果如图2所示。

(a)沿啮合线方向的综合错位量

(b)太阳轮齿面接触单位长度载荷图2 LTCA分析结果Fig.2 Results of LTCA

由图2可知，该齿轮副的最大综合啮合错位量为110.79 μm，其太阳轮最大齿面接触单位长度载荷为1 290.0 N/mm。考虑到太阳轮齿面接触单位长度载荷沿齿面距离呈线性分布，以及齿廓直线修形可能会产生的载荷突变情况，并结合前期研究，本文对太阳轮采取齿向线形修形和齿廓鼓形修形来改善齿面载荷分布并降低传递误差波动。

考虑到齿向载荷分布不均和齿轮副的啮合错位有很大关系，故齿向修形量初选为齿轮副的啮合错位量，即初定齿向线性修形量为110.79 μm。同时，齿廓修形量可初选为轮齿产生的变形量。该变形量与轮齿所受载荷的大小以及轮齿啮合刚度等因素有关，可由下式计算[14]

(1)

式中:δa为齿廓弹性变形量，μm；ωt为单位长度载荷，N/mm；cγ为轮齿平均啮合刚度，N/mm/μm，可根据式(2)计算得到

cγ=c′(0.75εa+0.25)

(2)

式中:c′为单对齿刚度；εa为端面重合度。根据文献[15]的计算公式，cγ最终可取为20.00 N/mm/μm。将由图2(b)所得的最大单位长度载荷1 290.0 N/mm代入式(1)可得δa=64.50 μm，故初选齿廓鼓形修形量为64.50 μm。

按照上述修形参数，在Romax中对太阳轮进行轮齿修形，其修形曲面是综合齿向修形与齿廓修形得到的，如图3(a)所示。修形后的齿面接触单位长度载荷分布如图3(b)所示。

由图3(b)可知，与修形前相比，修形后的齿廓载荷变化平缓，啮入啮出冲击较小。但就齿向载荷分布而言，仍存在一定的偏载。由此可见，按照经验修形公式拟定的修形方案，虽可有效改善齿面载荷状况，但并不能保证修形效果最佳。有鉴于此，下文将分别运用响应面法和Kriging法，以上述修形量为基准值，在其[-3σ1，3σ1]区间内对修形量进行寻优。其中，标准差σ1由齿轮的精度等级确定。

对于精度等级为5的标准齿轮，其齿距累积总公差Fp可由式(3)进行简单计算：

(3)

式中:mn为轮齿的法向模数；d为齿轮分度圆直径。

代入各项数据并根据公差计算圆整规则，得太阳轮的齿距累积总公差Fp=28 μm。将寻优范围的大小近似于齿距累积总公差值，故σ1=1/6·Fp，为便于后续处理，此处圆整取为σ1=4.50 μm。

(a)太阳轮修形曲面

(b)修形后单位长度载荷分布图3 修形曲面以及修形效果Fig.3 Modification surface and modification effect

2 基于响应面法和Kriging法的函数拟合

2.1 响应面法函数拟合

响应面函数的拟合形式一般分为一次函数拟合和二次函数描述。为保证拟合精度，本文选用含有交叉项的二次函数来拟合，其通式表达为：

(4)

为提高函数拟合的准确性和有效性，本文采用中心复合设计方法确定响应面法需要的样本点，其样本数为2k+2k+1[15]；自变量取值x可由式(5)求得：

x=μ+σPn

(5)

式中:μ表示随机变量分布的均值；σ为分布的标准差；Pn表示因素水平，根据中心复合设计方法分别取±α，±1，0。其中，α的取值由式(6)确定：

(6)

式中:k表示自变量个数，此处取k=2。

由式(5)可制定中心复合法的样本点数据，其结果如表1所示。将各样本数据输入Romax中构建的LTCA模型，可得s-p外啮合副的动态传递误差波动量Y。

表1 响应面法样本点及其响应值

根据表1的数据，采用最小二乘法可得s-p外啮合齿轮副的动态传递误差波动量拟合式：

(7)

在拟合出响应函数后，可用其对自变量的偏导来表示响应对某变量的敏感度。由式(7)对x1求偏导，可得动态传递误差波动量对齿向修形的敏感度：

(8)

同样，由式(7)对x2求偏导，可得动态传递误差波动量对齿廓修形的敏感度：

(9)

2.2 Kriging法函数拟合

与响应面法相比，Kriging法是一种更加具有统计学特性的函数拟合方法。其半参数化模型由一个线性回归参数模型和一个非参数随机过程联合构成。模型的线性部分由多项式形式表达，非参数部分一般为随机分布，其通式表达为[16]：

(10)

式中:f(x)=[f1(x), …,fm(x)]T为变量x的多项式，提供模拟的全局近似，即响应的数学期望；线性回归系数β=[β1, …,βm]T；z(x)提供的是拟合函数的局部偏差的近似，即响应的局部变化，通常要求服从正态分布N(0,σ2)，其协方差为非零值，即z(x)同分布且非独立。

cov[z(xi),z(xj)]=σ2R[R(xi,xj)],i,j=1,…,n

(11)

式中:R(xi,xj)是任意两个样本点之间的相关方程，对拟合的精确程度起决定性作用，一般用高斯相关方程表示：

(12)

式中:θ的值是当式(13)取最大时的值。

(13)

经拟合值的无偏估计和最小化误差均方差，拟合方程(10)可以转换为：

(14)

式中:rT(x)为含有自变量x的矩阵。

(15)

β=(FTR-1F)-1FTR-1Y

(16)

式中:Y为样本点处的响应值列阵。

在Kriging法中，常用拉丁超立方法抽样。为方便与响应面法比较，此处样本区间选为[-3σ1,3σ1]，样本点个数也取为9。表2为Kriging法的样本点及其响应值。

表2 Kriging法样本点及其响应值

根据表2的样本点数据，运用Kriging法可得动态传递误差波动量的拟合函数为：

(17)

(18)

与响应面法类似，对Kriging法获得的拟合函数求偏导，可得响应对各自变量的参数敏感度。对式(17)中各元素(即式(12))求导，可得：

(19)

(20)

将式(19)、(20)代入式(17)，可得s-p外啮合齿轮副动态传递误差波动量对齿向修形和齿廓修形的敏感度。

3 拟合结果比较及修形参数优化

3.1 响应面法与Kriging法的比较

由式(7)、(17)可分别绘出响应面法及Kriging法拟合的动态传递误差波动量随齿向、齿廓修形量的变化，其结果如图4所示。

(a)响应面法 (b)Kriging法图4 动态传递误差波动量随修形量变化Fig.4 ΔDTE with respect to modification amounts

(a)响应面法 (b)Kriging法图5 动态传递误差波动量对修形量敏感度Fig.5 The sensitivity of ΔDTE with respect to modification amounts

由图4可知，两种方法拟合出的动态传递误差波动量随修形量的变化趋势基本一致。采用齿向线性修形和齿廓鼓形修形的综合修形策略后，可以有效降低齿轮副的动态传递误差波动量。LTCA分析显示，未修形时其动态传递误差波动量为26.19 μm，而修形后，动态传递误差波动量下降至8 μm。故，修形量在区间[97.29,124.29]、[51,78]中，当齿向和齿廓修形分别取97.29 μm和51.00 μm时，动态传递误差的波动量可降低至区间内最小。图5为进一步分析齿轮传递误差波动量对修形量的敏感度。

由图5可知，响应面法和Kriging法拟合的自变量敏感度变化趋势略有不同。采用响应面法拟合获得的动态传动误差波动量，其对齿向修形参数、齿廓修形参数的敏感度随修形量的增大单调递增。而采用Kriging法拟合获得的动态传动误差波动量，其对齿向修形量的敏感度随修形量的增大先递增后递减，但其对齿廓修形量的敏感度随修形量的增大先递增而后呈现平缓的保持。

3.2 两种拟合方法与Romax仿真结果对比

为考察上述两拟合方法获得的齿轮动态传递误差波动量响应函数与真实函数之间的贴近程度，下文在自变量区间中按拉丁抽样方式选取16个样本点，分别运用Romax软件和RSM法、Kriging法计算各样本点的响应值，计算结果如表3所示。

由表3可得，由响应面法和Kriging法拟合而得的响应值与Romax仿真结果非常接近，表明两种方法在拟合动态传递误差波动量与修形参数间关系时具有良好的精度。

表3 两种拟合方法响应值与Romax仿真结果对比

进一步分析响应面法和Kriging法在拟合动态传递误差波动量与修形参数间关系的逼近程度。为直观计，图6给出了各样本点处拟合函数的响应值与仿真结果的相对误差。

图6 拟合相对误差Fig.6 The relative error of fitting

由图6可知，由Kriging法拟合的函数响应值与真实响应值间的计算误差大多小于5%，而由响应面法拟合的函数响应值的计算误差约为10%。两种拟合方法仅在第1、3号样本点处存在较大的相对误差，然进一步分析可知该样本点处动态传递误差波动量的绝对值非常小，分别为1.47 μm和1.13 μm，从而导致样本1、3的相对误差较大。因此，可以认为Kriging法较响应面法在拟合动态传递误差波动量时具有更高的精度，下文将基于Kriging法拟合的响应函数对修形参数进行优选。

3.3 修形量优选

根据式(17)，以动态传递误差波动量最小为优化目标，寻求最佳修形参数。计算可知，在寻优区间内，当齿向线性修形量97.29 μm、齿廓鼓形修形量51.00 μm时，动态传递误差波动量在区间内最小，为0.67 μm。

采用上述修形量对s-p齿轮副进行综合修形后，再对其进行LTCA分析，可知齿面载荷分布，结果如图7所示。

图7 优化修形后齿面单位长度载荷分布Fig.7 Load of per unit length after optimized modification

对比图3(b)可知，采用优化修形参数，可以有效消除齿面偏载现象。齿面最大单位长度载荷从1 290.0 N/mm降至923.1 N/mm，降幅达28.45%。此时，动态传递误差波动量降幅达97.44%。

4 结 论

(1)分别采用响应面法和Kriging法，拟合出了斜齿轮副动态传递误差波动量与齿轮修形参数的映射函数。

(2)当修形参数在基准值的[-3σ,3σ]区间内变化时，齿轮动态传递误差波动量对齿廓和齿向修形参数的敏感度随拟合方法的不同而略有差异。

(3)两种拟合方法获得的函数响应值与真实值均较为接近，能满足一般工程精度要求；但相比而言，Kriging法拟合函数的精度更高。

(4)基于Kriging法获得的拟合函数，以动态传递误差波动量在修形区间内最小为目标对修形参数进行了优选。结果表明，优选后的修形量可有效降低传递误差波动，改善齿面均载。

[ 1 ] 王旭东，林腾蛟，李润方，等．风力发电机组齿轮系统内部动态激励和响应分析[J]. 机械设计与研究, 2006, 22(3): 47-49.

WANG Xudong, LIN Tengjiao, LI Runfang, et al. The internal dynamic excitation and response analysis of the gear system for the wind-driven generator[J]. Mechanical Design and Research, 2006, 22(3): 47-49.

[ 2 ] BEACHAM M R, BELL D J, POWELL N N, et al. Development of transmission whine prediction tools[J]. SAE Technical Paper, 1999, 108(6): 1767-1780.

[ 3 ] 于蓬，章桐，冷毅,等. 电动车减/差速器振动特性分析及改进[J]. 振动与冲击，2015, 34(7): 85-92.

YU Peng, ZHANG Tong, LENG Yi, et al. Vibration characteristics analysis and improvement of differential/ reducer of electric vehicle [J]. Journal of Vibration and Shork, 2015, 34(7): 85-92.

[ 4 ] 孙福兆. 基于Romax的拖拉机变速器齿轮修形分析[J]. 山东农业大学学报(自然科学版)，2014, 45(5): 702-705.

SUN Fuzhao. Analysis of gear modification of a tractor transmission based on Romax software [J]. Journal of Shandong Agriculture University (Natural Science Edition), 2014, 45(5): 702-705.

[ 5 ] 张策. 机床噪声原理及控制[M]. 天津：天津科学技术出版社，1982.

[ 6 ] 王成，刘辉，项昌乐. 波动转矩下渐开线直齿轮传动齿廓修形研究[J]. 振动与冲击，2014, 33(24): 32-38.

WANG Cheng, LIU Hui, XIANG Changle. Influences of profile modification on dynamic characteristics of involute spur gears under a fluctuating torque[J]. Journal of Vibration and Shork, 2014, 33(24): 32-38.

[ 7 ] 张俊，郭凡. 面向最小动态传递误差波动量的风电齿轮修形研究[J]. 机械传动，2014, 38(2): 16-20.

ZHANG Jun, GUO Fan. Research on gear modification of wind turbine gearbox based on minimum dynamic transmission error fluctuating quantity[J]. Journal of Mechanical Transmission, 2014, 38(2): 16-20.

[ 8 ] BONORI G, BARBIERI M, PELLIEANO M. Optimum profile modifications of spur gears by means of genetic algorithms [J]. Journal of Sound and Vibration, 2008, 313(3): 603-616.

[ 9 ] 袁哲，孙志礼，郭瑜. 直齿圆柱齿轮齿廓修形曲线优化设计[J]. 机械传动，2010, 34(5): 1-4.

YUAN Zhe, SUN Zhili, GUO Yu. Optimal design of profile modification curves for spur gears [J]. Journal of Mechanical Transmission, 2010, 34(5): 1-4.

[10] PARK C I L. Multi-objective optimization of the tooth surface in helical gears using design of experiment and the response surface method[J]. Journal of Mechanical Science and Technology, 2010, 24(3): 823-829.

[11] ZHANG Jun, GUO Fan. Statistical modification analysis of helical planetary gears based on Response Surface Method and Monte Carlo simulation[J]. Chinese Journal of Mechanical Engineering, 2015, 28(6): 1194-1230.

[12] 高扬，方宗德，蒋进科, 等. 基于回归分析的直齿齿廓修形计算研究[J]. 机械传动，2014, 38(7): 56-59.

GAO Yang, FANG Zongde, JIANG Jinke, et al. Research of straigh tooth profile modification calculation basedon the regression analysis[J]. Journal of Mechanical Transmission, 2014, 38(7): 56-59.

[13] ZHANG Guangyu, WANG Guoqiang, LI Xuefei, et al. Global optimization of reliability design for large ball mill gear transmission based on the Kriging model and genetic algorithm[J]. Mechanism and Machine Theory, 2013, 69(6): 321-336.

[14] 朱孝录. 齿轮传动设计手册[M]. 2版. 北京: 化学工业出版社, 2010.

[15] 张尚荣，谭平，杜永峰, 等. 基于响应面法的层间隔震结构地震易损性分析[J]. 振动与冲击，2014, 33(15): 42-48.

ZHANG Shangrong, TAN Ping, DU Yongfeng, et al. Seismic fragility analysis of Inter-story isolation structures based on response surface method[J]. Journal of Vibration and Shork, 2014, 33(15): 42-48.

[16] 张崎，李兴斯. 结构可靠性分析的模拟重要抽样方法[J]. 工程力学, 2007, 24(1): 33-36.

ZHANG Qi, LI Xingsi. Importance sampling approach in structural reliability analysis based on Kriging simulation [J]. Engineering Mechanics, 2007, 24(1): 33-36.