基于差齒面拓撲的輪齒承載擬赫茲接觸分析

魏冰陽 李家琦 王文勝

1.河南科技大學機電工程學院，洛陽，4710032.河南科技大學工程力學系，洛陽，471023 3.大連理工大學工業裝備結構分析國家重點實驗室，大連，116023

0 引言

輪齒接觸仿真已成為現代齒輪傳動設計與性能控制不可或缺的一項技術。20世紀90代以來，輪齒接觸分析(tooth contact analysis，TCA)技術在復雜齒面的設計中逐漸得到應用。之后，發展出了ease-off差齒面拓撲分析方法[1-2]。TCA技術與ease-off拓撲方法的結合使輪齒接觸仿真分析更趨完善。WANG等[3]利用加工時齒面拓撲誤差的實時分布對誤差的變化趨勢進行了分析與預控。文獻[4-5]利用ease-off拓撲解決了弧齒錐齒輪高階傳動誤差拓撲，以及準雙曲面齒輪的修形設計與預控問題。文獻[6-7]發展了ease-off等距變換分析方法，據此提出了曲面綜合弧齒錐齒輪加工參數計算方法。利用TCA技術、ease-off曲面拓撲方法，蔣進科等[8]研究了高階傳動誤差斜齒輪的修形設計與實現方法。上述研究較好地解決了ease-off差齒面拓撲設計與齒面性能控制的問題，但對ease-off曲面的拓撲映射、嚙合信息的解析還不夠，更欠缺承載接觸分析(loaded TCA，LTCA)的求解計算，因此，有必要基于ease-off曲面拓撲，發展LTCA計算方法。有限元建模方法難以用于齒面拓撲修形的參數化、高效率設計計算。勢能解析法較好地解決了漸開線圓柱齒輪的嚙合剛度計算問題，與有限元方法相比，平均剛度誤差在5%以內[9-11]，且計算效率遠高于有限元方法。鑒于此，筆者利用ease-off曲面幾何分析和輪齒剛度勢能計算方法，建立了ease-off齒面幾何拓撲與LTCA計算模型，給出了從齒面幾何拓撲到力學分析的一體化計算流程。

1 拋物面拓撲修形的ease-off曲面

1.1 齒條產成漸開線圓柱齒輪模型

以u、v為曲面參數的拋物面方程如下：

w(u,v)=a0+a1u2+a2uv+a3v2

(1)

式中，u為齒長參數；v為齒廓切向參數；a0～a3為多項式系數；a1控制拋物面u方向的法曲率，對應齒向修形；a2控制拋物面的扭轉方向；a3控制v方向的法曲率，對應齒廓修形。

圖1所示坐標系S0(O0X0Y0Z0)中，曲面w的法線為n，π平面為產形齒條面。利用矢量旋轉變換，繞X0旋轉壓力角α、繞Y0旋轉螺旋角β后，斜齒條方程為

圖1 齒條面坐標系

[x0y0z0]T=T0[uvw]T

(2)

式中，T0為旋轉變換矩陣。

建立齒條與圓柱齒輪的產成坐標系，如圖2所示。通過坐標變換，可求出圓柱齒輪齒面方程與法矢量：

圖2 齒條-齒面產成坐標系

(3)

rd=[x0y0+r0z0+r0φ]T

n0=T0[2a1u+a1v2a3v+a2u-1]T

式中，下標1、s分別表示修形齒面r1和標準漸開線齒面rs；下標0、d分別表示坐標系S0和Sd;M1d為坐標系Sd到S1的坐標變換矩陣；φ為齒輪轉角，由嚙合方程得到；r0為齒輪節圓半徑。

w=0時，拋物面退化為平面，可作為標準齒條展成標準漸開線圓柱齒輪rs。漸開線齒面的修形量以r1與rs之間的法向誤差——離差zd表示。以u、v為水平坐標軸、zd為垂直坐標軸繪制誤差曲面，可得到ease-off差齒面。

1.2 基于齒條公切產形的ease-off曲面映射

齒條平面Σ0左右兩側分別與小輪齒廓Σ1、大輪齒廓Σ2相切(圖2)，Σ1、Σ0、Σ2在展成運動中一直保持公切嚙合運動狀態，Σ0分別與Σ1、Σ2構成曲面映射關系。以齒條平面Σ0為坐標映射平面(以u、v為參數)，得到大小輪修形齒面的差齒面。對大小輪ease-off曲面的對應點高度zdg和zdp(下標g、p分別表示大輪和小輪)做代數和，可得到大小輪差齒面ease-off的離差zd，如圖3所示。

圖3 差齒面離差圖

本文以一對漸開線齒輪為例建立ease-off曲面映射。齒輪參數如下：齒數比z1/z2=21/62，壓力角α=20°，螺旋角β=13°，法向模數mn=5 mm，齒寬B=80 mm，變位系數x1=0.1，x2=-0.1。拋物面參數為：a1=1.11×10-5mm-1，a2=6.83×10-4mm-1。

2 Ease-off曲面拓撲的輪齒接觸分析

2.1 差曲線與接觸路徑

由式(3)可知小輪的轉角φi=φ(u(i,j),v(i,j))，其中，i為差曲線序號，j為第i條差曲線上點的序號。根據離差zd的定義可得zdi=zd(u(i,j),v(i,j))。函數zdi對應ease-off曲面上的曲線si——差曲線，如圖3所示。差曲線反映齒面修形后由線接觸變為點接觸而產生的差曲率，差曲率的大小和方向決定了齒面瞬時接觸橢圓長軸的大小和方向。按順序遍歷φi可得一族差曲線：

si=s(φi)

(4)

將式(4)中的序列差曲線置于與齒面相切的位置，差曲線與齒面的切點構成了齒面接觸路徑(圖4)，差曲線瀑布的垂直高度反映沿瞬時接觸線方向的兩齒面間隙zs。

圖4 齒面接觸路徑仿真

2.2 接觸路徑與傳動誤差

差曲線脊點(圖3、圖4中差曲線si的極值點)軌跡的映射為齒面的接觸路徑[12]。差曲線脊點的離差為該點的傳動誤差ETi(μm)。以小輪轉角φi(°)為橫坐標，傳動誤差為縱坐標，繪制無載傳動誤差EUT曲線，如圖5所示，3條拋物線表示前齒Ti-1、當前齒Ti與后齒Ti+1的順序嚙合過程。波浪線為不同載荷下的承載傳動誤差ELT。

1～9、14～22、27～35為三齒接觸區 9～14、22～27為二齒接觸區

綜合式(1)～式(4)，在一次完成齒面數字化后，通過ease-off曲面解析，可得齒面接觸路徑、傳動誤差、接觸線等輪齒嚙合特性信息。

3 嚙合剛度計算與承載接觸分析

3.1 輪齒嚙合剛度計算

齒輪副嚙合剛度kv的計算公式為

(5)

式中，kp1～kp4分別為小輪單齒的彎曲剛度、剪切剛度、徑向壓縮剛度和基體剛度[12]；kg1～kg4分別為大輪單齒的彎曲剛度、剪切剛度、徑向壓縮剛度和基體剛度[12]；kc為赫茲接觸剛度。

將斜齒輪沿縱向細分為有限個單位長度的直齒圓柱齒輪，將其看作僅受單位力作用的單元體(每一單元的端面齒廓相同)，如圖6所示。沿齒廓方向，以v為變量，將齒廓劃分為n個單元，每個齒廓單元的厚度為dy。圖6中，Sz、ly分別為點M到輪齒中線與齒根圓的距離。對dy進行數值積分，計算齒廓不同位置的剛度，遍歷vi(vmin≤vi≤vmax；i=1,2,…,n)得到齒廓時變的輪齒單位剛度矩陣Kv：

Kv=[kv1kv2…kvn]

(6)

注：Ei+1為后齒誤差；kvi為當前齒的單元剛度；si為當前齒的差曲線；di為當前齒的變形

(7)

3.2 承載接觸分析

齒面修形引起的點接觸效應導致接觸線上除接觸脊點外的位置均存在間隙zs(圖4中差曲線的高度坐標)。載荷F作用下，當前齒Ti接觸脊產生的法向變形量為δi，單元j的變形量為δi-zsj。假設載荷由m個單元承擔(見圖7)，則輪齒變形協調公式為

(8)

多齒接觸(圖7)要考慮同一時間前齒Ti-1和后齒Ti+1上各單元的接觸順序，后齒Ti+1的接觸發生在當前齒Ti的齒間間隙消除之后即δi≥ETi+1。接觸線上單元接觸的競爭條件始終是間隙小的先接觸，因此在載荷的作用下，各接觸線單元按照自身間隙量zsj由小到大的順序依次發生接觸，直至m個單元分擔的載荷F1～Fm滿足收斂條件。

按照齒面接觸路徑，遍歷接觸線序列，求出轉角φ對應的嚙合剛度Kφ、脊點變形量δφ、傳動誤差ETφ、接觸線載荷Φφ，以及載荷分擔比λφ=Φφ/Φ、承載傳動誤差LTEφ=δφ+ETφ。

式(1)到式(8)的計算過程就是由ease-off曲面到輪齒承載接觸分析的計算過程，計算無需求解非線性方程組，適宜數值化計算。

大輪扭矩MT為750 N·m、3000 N·m、6000 N·m時，對應的法向作用力F為5 kN、20 kN、40kN。沿齒廓劃分140個單元，沿齒向劃分200個單元，在嚙合區間內，小輪轉角φ依次取35個數值，代表35個嚙合位置。計算得到的承載傳動誤差如圖5所示，3種載荷下承載傳動誤差波動均在2 μm以內，3000 N·m的波動僅為0.5 μm。750 N·m、3000 N·m、6000 N·m這3種載荷下的齒面最大相對位移依次為9.13 μm、16.84 μm、24.80 μm；隨著載荷的增大，由2或3個齒分擔的區段增多，齒面實際重合度增大；3000 N·m時的重合度為1.86，6000 N·m時的重合度為2.32。

齒間載荷分擔比如圖8所示，750 N·m扭矩下，0<λφ<1的區間為雙齒接觸區，對應轉角范圍[-11.42°，-5.15°)、(5.73°，11.99°]，λφ=1的區間為單齒接觸區，對應轉角范圍是[-5.15°，5.73°]，轉角φ<-11.42°與φ>11.99°的區間里，齒面沒有參與嚙合，齒面實際重合度為1.38，與圖5反映的結果一致。

圖8 齒間載荷分擔比

圖9給出了扭矩3000 N·m、6000 N·m下16個嚙合位置(以不同顏色區分)接觸線上的載荷分布。圖9a中，齒面中部承載，最大單元載荷為108.4 N；圖9b中，齒面幾乎全部受載，而在進入嚙合與退出嚙合的齒端位置有少部分沒有接觸，最大載荷為150.8 N，齒面邊緣發生了接觸。

(a)載荷為3000 N·m

4 擬赫茲單元線接觸齒面應力計算

如圖10所示，任意瞬時的齒面接觸應變區域一般為橢圓形。邊緣發生接觸時將出現應力集中，瞬時接觸區發生畸變，一端邊緣接觸導致橢圓一端畸變，兩端邊緣接觸導致橢圓兩端畸變。載荷較小時，接觸線上載荷分布于齒面中部，接觸區呈橢圓形，該接觸問題是無限邊界的赫茲接觸問題，接觸應力可按赫茲點接觸彈性應變模型計算；載荷較大時，齒面邊緣接觸，接觸橢圓出現畸變，該接觸問題變成有限邊界彈性接觸問題。

圖10 假設瞬時接觸形變區

接觸線上的載荷分布已知時，邊緣接觸問題可采用擬赫茲線接觸的辦法求解，即沿接觸線長度方向將接觸區細分為m個長度為Δl的線接觸單元，如圖11所示。根據ease-off曲面的拓撲結構，將每個單元都當作圓柱同平面的接觸。由于Δl很小，單元j上的作用力可按均布載荷qj處理，即qj=Fj/Δl。

圖11 接觸區有限元模型

圓柱體與平面發生赫茲接觸，單元j接觸區中點處應力最大，為

(9)

(10)

式中，Rg,j、Rp,j分別為接觸單元j在大小輪垂直于接觸線方向上的法曲率半徑；Rj為單元j的綜合法曲率半徑。

接觸區半寬為

(11)

沿接觸區y方向任意點的接觸應力為

(12)

將計算出的各單元的應力列陣合并，得到任意接觸線的應力分布矩陣σφ=[σ(j,y)]。對矩陣進行插值加密，繪制接觸應力云圖(圖12a)，其中，P5對應一端邊緣接觸，瞬時接觸橢圓在一端畸變，P16對應一完整的瞬時接觸橢圓，P27對應兩端邊緣接觸，瞬時接觸橢圓兩端畸變。采用MASTA軟件對齒面進行承載接觸分析。如圖12所示，應力云圖的梯度形狀大致相同，兩者最大接觸應力分別為1211.0 MPa、1172.04 MPa，說明擬赫茲接觸計算方法能真實反映輪齒邊緣接觸狀況。

(a)擬赫茲接觸

通過對每根接觸線的求解，可得齒面各點的應力分布矩陣，在保證插值精度的情況下繪制全齒面接觸應力變化云圖(圖13、圖14)，3000 N·m下，擬赫茲接觸與MASTA的計算結果接近，齒面最大接觸應力分別為1014.0 MPa、967.5 MPa；齒面應變均集中在中部，齒面邊緣無接觸。6000 N·m下，齒面邊緣發生了接觸應變，存在應力集中，這與圖9反映的載荷分布情況一致。由于齒面修形梯度控制方式不同，兩者應力云圖形狀存在一些差異。

(a)擬赫茲接觸LTCA

(a)擬赫茲接觸LTCA

5 結論

(1)利用曲面多項式、齒條-齒輪等切共軛產形原理能準確構建數值齒面和ease-off差齒面模型，實現對齒面拓撲結構形狀與性質的控制。

(2)提出的擬赫茲線接觸計算方法較好地解決了齒面邊緣接觸的應力應變求解問題，將該方法的計算結果與第三方軟件MASTA的計算結果進行對比，兩者有較好的一致性。

(3)利用ease-off差齒面對輪齒剛度進行解析，能便捷地獲得輪齒的時變嚙合剛度、承載傳動誤差與齒面載荷分布圖，實現參數化輪齒接觸分析，避免了非線性方程組求解、有限元建模和巨量運算問題。