汪 建， 張 俊

(福州大學 機械工程及自動化學院， 福州 350116)

傳遞誤差是引起齒輪傳動系統(tǒng)振動與噪聲的主要原因[1-2]。工程上可通過齒輪修形來改善齒面受力、降低傳遞誤差波動，進而抑制傳動系統(tǒng)的振動和噪聲[3-5]。因此，明確輪齒修形參數(shù)和齒輪副傳遞誤差波動量之間的映射關系，對于指導齒輪修形、改善傳動性能具有重要的理論與工程意義。

為制定合理的齒輪修形方案，學術界和工業(yè)界開展了大量研究。王成等[6]建立了2自由度單級齒輪傳動非線性動力學模型，對比研究了齒廓修形參數(shù)對齒輪動態(tài)特性的影響。張俊等[7]運用多體方法建立了風電齒輪箱的動力學模型，探討了一種面向最小動態(tài)傳動誤差波動量的風電齒輪修形方法。Bonori等[8]嘗試將遺傳算法和齒輪有限元分析相結合，用來確定最優(yōu)的修形參數(shù)。遵循類似思路，袁哲等[9]以減小齒輪傳遞誤差波動為目標,對齒輪修形參數(shù)進行了優(yōu)化設計。

然而，由于齒輪系統(tǒng)具有復雜的非線性，理論上難以給出修形參數(shù)和傳遞誤差之間的解析表達。為此，部分學者采用數(shù)值擬合方法給出兩者間的關系。Park[10]運用響應面法建立線性模型來預測修形對靜態(tài)傳遞誤差的影響。Zhang等[11]運用響應面法擬合出修形量與動態(tài)傳遞誤差波動量間的函數(shù)關系，并運用蒙特卡洛法對齒輪修形的可靠性進行評估。方宗德等[12]利用回歸分析方法，給出了最佳修形量和反映修形效果的傳動誤差幅值比的回歸公式。Zhang等[13]綜合運用Kriging法和遺傳算法，開展了大型球磨機齒輪傳動系統(tǒng)的可靠性設計。

文獻檢索表明，響應面法和Kriging法是目前較為常用的兩種擬合方法，但二者在擬合特定場合關系函數(shù)上的優(yōu)劣尚無定論。有鑒于此，本文以某兆瓦級風電齒輪箱中行星級傳動的外嚙合斜齒輪副為對象，分別采用響應面法和Kriging法構建給定工況下齒輪修形參數(shù)和齒輪副傳遞誤差波動量間的映射函數(shù)。在此基礎上，將兩種擬合方法所得的響應函數(shù)結果與Romax數(shù)值仿真結果對比，以評估響應面法和Kriging法的擬合精度和修形參數(shù)對傳遞誤差波動的靈敏度。希冀通過上述研究，探索出用于構建齒輪傳遞誤差與修形參數(shù)間關系函數(shù)的高效擬合方法，據(jù)此優(yōu)選修形參數(shù)，為后續(xù)的減振降噪和傳動可靠性研究提供理論依據(jù)。

1 基準修形參數(shù)的確定

不失一般性，不妨以圖1所示的某兆瓦級風電齒輪箱中行星級傳動的外嚙合斜齒輪副為例，研究齒輪修形量與齒輪副動態(tài)傳遞誤差波動量之間的關系。

(a)行星級傳動 (b) 太陽輪-行星輪嚙合副圖1 某兆瓦級風電齒輪箱行星級及其外嚙合齒輪副Fig.1 S-P pair in the planetary stage of a MW wind turbine

圖1所示的斜齒行星傳動中，太陽輪制成一端浮動的齒輪軸，行星輪采用兩端對稱布置的圓柱滾子軸承支承，而內(nèi)齒圈與箱體固結。其中，行星架作為輸入構件(輸入轉(zhuǎn)速為n=29 r/min)，太陽輪為輸出構件(輸出轉(zhuǎn)矩為To=18 990 N·m)。傳動齒輪為ISO 5級精度，模數(shù)mn=10 mm、齒寬B=220 mm、法向壓力角αn=20°、端面壓力角αt=20.158°、螺旋角β=7.493°、變位系數(shù)分別為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?？紤]到太陽輪齒面接觸單位長度載荷沿齒面距離呈線性分布，以及齒廓直線修形可能會產(chǎn)生的載荷突變情況，并結合前期研究，本文對太陽輪采取齒向線形修形和齒廓鼓形修形來改善齒面載荷分布并降低傳遞誤差波動。

考慮到齒向載荷分布不均和齒輪副的嚙合錯位有很大關系，故齒向修形量初選為齒輪副的嚙合錯位量，即初定齒向線性修形量為110.79 μm。同時，齒廓修形量可初選為輪齒產(chǎn)生的變形量。該變形量與輪齒所受載荷的大小以及輪齒嚙合剛度等因素有關，可由下式計算[14]

(1)

式中:δa為齒廓彈性變形量，μm；ωt為單位長度載荷，N/mm；cγ為輪齒平均嚙合剛度，N/mm/μm，可根據(jù)式(2)計算得到

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

(2)

式中:c′為單對齒剛度；εa為端面重合度。根據(jù)文獻[15]的計算公式，cγ最終可取為20.00 N/mm/μm。將由圖2(b)所得的最大單位長度載荷1 290.0 N/mm代入式(1)可得δa=64.50 μm，故初選齒廓鼓形修形量為64.50 μm。

按照上述修形參數(shù)，在Romax中對太陽輪進行輪齒修形，其修形曲面是綜合齒向修形與齒廓修形得到的，如圖3(a)所示。修形后的齒面接觸單位長度載荷分布如圖3(b)所示。

由圖3(b)可知，與修形前相比，修形后的齒廓載荷變化平緩，嚙入嚙出沖擊較小。但就齒向載荷分布而言，仍存在一定的偏載。由此可見，按照經(jīng)驗修形公式擬定的修形方案，雖可有效改善齒面載荷狀況，但并不能保證修形效果最佳。有鑒于此，下文將分別運用響應面法和Kriging法，以上述修形量為基準值，在其[-3σ1，3σ1]區(qū)間內(nèi)對修形量進行尋優(yōu)。其中，標準差σ1由齒輪的精度等級確定。

對于精度等級為5的標準齒輪，其齒距累積總公差Fp可由式(3)進行簡單計算：

(3)

式中:mn為輪齒的法向模數(shù)；d為齒輪分度圓直徑。

代入各項數(shù)據(jù)并根據(jù)公差計算圓整規(guī)則，得太陽輪的齒距累積總公差Fp=28 μm。將尋優(yōu)范圍的大小近似于齒距累積總公差值，故σ1=1/6·Fp，為便于后續(xù)處理，此處圓整取為σ1=4.50 μm。

(a)太陽輪修形曲面

(b)修形后單位長度載荷分布圖3 修形曲面以及修形效果Fig.3 Modification surface and modification effect

2 基于響應面法和Kriging法的函數(shù)擬合

2.1 響應面法函數(shù)擬合

響應面函數(shù)的擬合形式一般分為一次函數(shù)擬合和二次函數(shù)描述。為保證擬合精度，本文選用含有交叉項的二次函數(shù)來擬合，其通式表達為：

(4)

為提高函數(shù)擬合的準確性和有效性，本文采用中心復合設計方法確定響應面法需要的樣本點，其樣本數(shù)為2k+2k+1[15]；自變量取值x可由式(5)求得：

x=μ+σPn

(5)

式中:μ表示隨機變量分布的均值；σ為分布的標準差；Pn表示因素水平，根據(jù)中心復合設計方法分別取±α，±1，0。其中，α的取值由式(6)確定：

(6)

式中:k表示自變量個數(shù)，此處取k=2。

由式(5)可制定中心復合法的樣本點數(shù)據(jù)，其結果如表1所示。將各樣本數(shù)據(jù)輸入Romax中構建的LTCA模型，可得s-p外嚙合副的動態(tài)傳遞誤差波動量Y。

表1 響應面法樣本點及其響應值

根據(jù)表1的數(shù)據(jù)，采用最小二乘法可得s-p外嚙合齒輪副的動態(tài)傳遞誤差波動量擬合式：

(7)

在擬合出響應函數(shù)后，可用其對自變量的偏導來表示響應對某變量的敏感度。由式(7)對x1求偏導，可得動態(tài)傳遞誤差波動量對齒向修形的敏感度：

(8)

同樣，由式(7)對x2求偏導，可得動態(tài)傳遞誤差波動量對齒廓修形的敏感度：

(9)

2.2 Kriging法函數(shù)擬合

與響應面法相比，Kriging法是一種更加具有統(tǒng)計學特性的函數(shù)擬合方法。其半?yún)?shù)化模型由一個線性回歸參數(shù)模型和一個非參數(shù)隨機過程聯(lián)合構成。模型的線性部分由多項式形式表達，非參數(shù)部分一般為隨機分布，其通式表達為[16]：

(10)

式中:f(x)=[f1(x), …,fm(x)]T為變量x的多項式，提供模擬的全局近似，即響應的數(shù)學期望；線性回歸系數(shù)β=[β1, …,βm]T；z(x)提供的是擬合函數(shù)的局部偏差的近似，即響應的局部變化，通常要求服從正態(tài)分布N(0,σ2)，其協(xié)方差為非零值，即z(x)同分布且非獨立。

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

(11)

式中:R(xi,xj)是任意兩個樣本點之間的相關方程，對擬合的精確程度起決定性作用，一般用高斯相關方程表示：

(12)

式中:θ的值是當式(13)取最大時的值。

(13)

經(jīng)擬合值的無偏估計和最小化誤差均方差，擬合方程(10)可以轉(zhuǎn)換為：

(14)

式中:rT(x)為含有自變量x的矩陣。

(15)

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

(16)

式中:Y為樣本點處的響應值列陣。

在Kriging法中，常用拉丁超立方法抽樣。為方便與響應面法比較，此處樣本區(qū)間選為[-3σ1,3σ1]，樣本點個數(shù)也取為9。表2為Kriging法的樣本點及其響應值。

表2 Kriging法樣本點及其響應值

根據(jù)表2的樣本點數(shù)據(jù)，運用Kriging法可得動態(tài)傳遞誤差波動量的擬合函數(shù)為：

(17)

(18)

與響應面法類似，對Kriging法獲得的擬合函數(shù)求偏導，可得響應對各自變量的參數(shù)敏感度。對式(17)中各元素(即式(12))求導，可得：

(19)

(20)

將式(19)、(20)代入式(17)，可得s-p外嚙合齒輪副動態(tài)傳遞誤差波動量對齒向修形和齒廓修形的敏感度。

3 擬合結果比較及修形參數(shù)優(yōu)化

3.1 響應面法與Kriging法的比較

由式(7)、(17)可分別繪出響應面法及Kriging法擬合的動態(tài)傳遞誤差波動量隨齒向、齒廓修形量的變化，其結果如圖4所示。

(a)響應面法 (b)Kriging法圖4 動態(tài)傳遞誤差波動量隨修形量變化Fig.4 ΔDTE with respect to modification amounts

(a)響應面法 (b)Kriging法圖5 動態(tài)傳遞誤差波動量對修形量敏感度Fig.5 The sensitivity of ΔDTE with respect to modification amounts

由圖4可知，兩種方法擬合出的動態(tài)傳遞誤差波動量隨修形量的變化趨勢基本一致。采用齒向線性修形和齒廓鼓形修形的綜合修形策略后，可以有效降低齒輪副的動態(tài)傳遞誤差波動量。LTCA分析顯示，未修形時其動態(tài)傳遞誤差波動量為26.19 μm，而修形后，動態(tài)傳遞誤差波動量下降至8 μm。故，修形量在區(qū)間[97.29,124.29]、[51,78]中，當齒向和齒廓修形分別取97.29 μm和51.00 μm時，動態(tài)傳遞誤差的波動量可降低至區(qū)間內(nèi)最小。圖5為進一步分析齒輪傳遞誤差波動量對修形量的敏感度。

由圖5可知，響應面法和Kriging法擬合的自變量敏感度變化趨勢略有不同。采用響應面法擬合獲得的動態(tài)傳動誤差波動量，其對齒向修形參數(shù)、齒廓修形參數(shù)的敏感度隨修形量的增大單調(diào)遞增。而采用Kriging法擬合獲得的動態(tài)傳動誤差波動量，其對齒向修形量的敏感度隨修形量的增大先遞增后遞減，但其對齒廓修形量的敏感度隨修形量的增大先遞增而后呈現(xiàn)平緩的保持。

3.2 兩種擬合方法與Romax仿真結果對比

為考察上述兩擬合方法獲得的齒輪動態(tài)傳遞誤差波動量響應函數(shù)與真實函數(shù)之間的貼近程度，下文在自變量區(qū)間中按拉丁抽樣方式選取16個樣本點，分別運用Romax軟件和RSM法、Kriging法計算各樣本點的響應值，計算結果如表3所示。

由表3可得，由響應面法和Kriging法擬合而得的響應值與Romax仿真結果非常接近，表明兩種方法在擬合動態(tài)傳遞誤差波動量與修形參數(shù)間關系時具有良好的精度。

表3 兩種擬合方法響應值與Romax仿真結果對比

進一步分析響應面法和Kriging法在擬合動態(tài)傳遞誤差波動量與修形參數(shù)間關系的逼近程度。為直觀計，圖6給出了各樣本點處擬合函數(shù)的響應值與仿真結果的相對誤差。

圖6 擬合相對誤差Fig.6 The relative error of fitting

由圖6可知，由Kriging法擬合的函數(shù)響應值與真實響應值間的計算誤差大多小于5%，而由響應面法擬合的函數(shù)響應值的計算誤差約為10%。兩種擬合方法僅在第1、3號樣本點處存在較大的相對誤差，然進一步分析可知該樣本點處動態(tài)傳遞誤差波動量的絕對值非常小，分別為1.47 μm和1.13 μm，從而導致樣本1、3的相對誤差較大。因此，可以認為Kriging法較響應面法在擬合動態(tài)傳遞誤差波動量時具有更高的精度，下文將基于Kriging法擬合的響應函數(shù)對修形參數(shù)進行優(yōu)選。

3.3 修形量優(yōu)選

根據(jù)式(17)，以動態(tài)傳遞誤差波動量最小為優(yōu)化目標，尋求最佳修形參數(shù)。計算可知，在尋優(yōu)區(qū)間內(nèi)，當齒向線性修形量97.29 μm、齒廓鼓形修形量51.00 μm時，動態(tài)傳遞誤差波動量在區(qū)間內(nèi)最小，為0.67 μm。

采用上述修形量對s-p齒輪副進行綜合修形后，再對其進行LTCA分析，可知齒面載荷分布，結果如圖7所示。

圖7 優(yōu)化修形后齒面單位長度載荷分布Fig.7 Load of per unit length after optimized modification

對比圖3(b)可知，采用優(yōu)化修形參數(shù)，可以有效消除齒面偏載現(xiàn)象。齒面最大單位長度載荷從1 290.0 N/mm降至923.1 N/mm，降幅達28.45%。此時，動態(tài)傳遞誤差波動量降幅達97.44%。

4 結 論

(1)分別采用響應面法和Kriging法，擬合出了斜齒輪副動態(tài)傳遞誤差波動量與齒輪修形參數(shù)的映射函數(shù)。

(2)當修形參數(shù)在基準值的[-3σ,3σ]區(qū)間內(nèi)變化時，齒輪動態(tài)傳遞誤差波動量對齒廓和齒向修形參數(shù)的敏感度隨擬合方法的不同而略有差異。

(3)兩種擬合方法獲得的函數(shù)響應值與真實值均較為接近，能滿足一般工程精度要求；但相比而言，Kriging法擬合函數(shù)的精度更高。

(4)基于Kriging法獲得的擬合函數(shù)，以動態(tài)傳遞誤差波動量在修形區(qū)間內(nèi)最小為目標對修形參數(shù)進行了優(yōu)選。結果表明，優(yōu)選后的修形量可有效降低傳遞誤差波動，改善齒面均載。

[ 1 ] 王旭東，林騰蛟，李潤方，等．風力發(fā)電機組齒輪系統(tǒng)內(nèi)部動態(tài)激勵和響應分析[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]. 山東農(nóng)業(yè)大學學報(自然科學版)，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 ] 王成，劉輝，項昌樂. 波動轉(zhuǎn)矩下漸開線直齒輪傳動齒廓修形研究[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 ] 張俊，郭凡. 面向最小動態(tài)傳遞誤差波動量的風電齒輪修形研究[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 ] 袁哲，孫志禮，郭瑜. 直齒圓柱齒輪齒廓修形曲線優(yōu)化設計[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版. 北京: 化學工業(yè)出版社, 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.