人字齒輪最佳齒向修形設計的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，本文所得最佳修形參數更為優越，在人字齒輪減振，降噪等方面有重要意義。

文中采取的響應面優化算法，可以在全局參數下得到最優修形參數，但是響應面會隨著工況、齒輪副參數改變而改變，不能一概而論，且提高動力學求解效率是響應面得到的關鍵一步，有待進一步優化。