基于差齒面拓?fù)涞妮嘄X承載擬赫茲接觸分析

魏冰陽(yáng) 李家琦 王文勝

1.河南科技大學(xué)機(jī)電工程學(xué)院，洛陽(yáng)，4710032.河南科技大學(xué)工程力學(xué)系，洛陽(yáng)，471023 3.大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連，116023

0 引言

輪齒接觸仿真已成為現(xiàn)代齒輪傳動(dòng)設(shè)計(jì)與性能控制不可或缺的一項(xiàng)技術(shù)。20世紀(jì)90代以來，輪齒接觸分析(tooth contact analysis，TCA)技術(shù)在復(fù)雜齒面的設(shè)計(jì)中逐漸得到應(yīng)用。之后，發(fā)展出了ease-off差齒面拓?fù)浞治龇椒╗1-2]。TCA技術(shù)與ease-off拓?fù)浞椒ǖ慕Y(jié)合使輪齒接觸仿真分析更趨完善。WANG等[3]利用加工時(shí)齒面拓?fù)湔`差的實(shí)時(shí)分布對(duì)誤差的變化趨勢(shì)進(jìn)行了分析與預(yù)控。文獻(xiàn)[4-5]利用ease-off拓?fù)浣鉀Q了弧齒錐齒輪高階傳動(dòng)誤差拓?fù)洌约皽?zhǔn)雙曲面齒輪的修形設(shè)計(jì)與預(yù)控問題。文獻(xiàn)[6-7]發(fā)展了ease-off等距變換分析方法，據(jù)此提出了曲面綜合弧齒錐齒輪加工參數(shù)計(jì)算方法。利用TCA技術(shù)、ease-off曲面拓?fù)浞椒?，蔣進(jìn)科等[8]研究了高階傳動(dòng)誤差斜齒輪的修形設(shè)計(jì)與實(shí)現(xiàn)方法。上述研究較好地解決了ease-off差齒面拓?fù)湓O(shè)計(jì)與齒面性能控制的問題，但對(duì)ease-off曲面的拓?fù)溆成?、嚙合信息的解析還不夠，更欠缺承載接觸分析(loaded TCA，LTCA)的求解計(jì)算，因此，有必要基于ease-off曲面拓?fù)?，發(fā)展LTCA計(jì)算方法。有限元建模方法難以用于齒面拓?fù)湫扌蔚膮?shù)化、高效率設(shè)計(jì)計(jì)算。勢(shì)能解析法較好地解決了漸開線圓柱齒輪的嚙合剛度計(jì)算問題，與有限元方法相比，平均剛度誤差在5%以內(nèi)[9-11]，且計(jì)算效率遠(yuǎn)高于有限元方法。鑒于此，筆者利用ease-off曲面幾何分析和輪齒剛度勢(shì)能計(jì)算方法，建立了ease-off齒面幾何拓?fù)渑cLTCA計(jì)算模型，給出了從齒面幾何拓?fù)涞搅W(xué)分析的一體化計(jì)算流程。

1 拋物面拓?fù)湫扌蔚膃ase-off曲面

1.1 齒條產(chǎn)成漸開線圓柱齒輪模型

以u(píng)、v為曲面參數(shù)的拋物面方程如下：

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

(1)

式中，u為齒長(zhǎng)參數(shù)；v為齒廓切向參數(shù)；a0～a3為多項(xiàng)式系數(shù)；a1控制拋物面u方向的法曲率，對(duì)應(yīng)齒向修形；a2控制拋物面的扭轉(zhuǎn)方向；a3控制v方向的法曲率，對(duì)應(yīng)齒廓修形。

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

圖1 齒條面坐標(biāo)系

[x0y0z0]T=T0[uvw]T

(2)

式中，T0為旋轉(zhuǎn)變換矩陣。

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

圖2 齒條-齒面產(chǎn)成坐標(biāo)系

(3)

rd=[x0y0+r0z0+r0φ]T

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

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

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

1.2 基于齒條公切產(chǎn)形的ease-off曲面映射

齒條平面Σ0左右兩側(cè)分別與小輪齒廓Σ1、大輪齒廓Σ2相切(圖2)，Σ1、Σ0、Σ2在展成運(yùn)動(dòng)中一直保持公切嚙合運(yùn)動(dòng)狀態(tài)，Σ0分別與Σ1、Σ2構(gòu)成曲面映射關(guān)系。以齒條平面Σ0為坐標(biāo)映射平面(以u(píng)、v為參數(shù))，得到大小輪修形齒面的差齒面。對(duì)大小輪ease-off曲面的對(duì)應(yīng)點(diǎn)高度zdg和zdp(下標(biāo)g、p分別表示大輪和小輪)做代數(shù)和，可得到大小輪差齒面ease-off的離差zd，如圖3所示。

圖3 差齒面離差圖

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

2 Ease-off曲面拓?fù)涞妮嘄X接觸分析

2.1 差曲線與接觸路徑

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

si=s(φi)

(4)

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

圖4 齒面接觸路徑仿真

2.2 接觸路徑與傳動(dòng)誤差

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

1～9、14～22、27～35為三齒接觸區(qū) 9～14、22～27為二齒接觸區(qū)

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

3 嚙合剛度計(jì)算與承載接觸分析

3.1 輪齒嚙合剛度計(jì)算

齒輪副嚙合剛度kv的計(jì)算公式為

(5)

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

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

Kv=[kv1kv2…kvn]

(6)

注：Ei+1為后齒誤差；kvi為當(dāng)前齒的單元?jiǎng)偠?；si為當(dāng)前齒的差曲線；di為當(dāng)前齒的變形

(7)

3.2 承載接觸分析

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

(8)

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

按照齒面接觸路徑，遍歷接觸線序列，求出轉(zhuǎn)角φ對(duì)應(yīng)的嚙合剛度Kφ、脊點(diǎn)變形量δφ、傳動(dòng)誤差ETφ、接觸線載荷Φφ，以及載荷分擔(dān)比λφ=Φφ/Φ、承載傳動(dòng)誤差LTEφ=δφ+ETφ。

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

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

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

圖8 齒間載荷分擔(dān)比

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

(a)載荷為3000 N·m

4 擬赫茲單元線接觸齒面應(yīng)力計(jì)算

如圖10所示，任意瞬時(shí)的齒面接觸應(yīng)變區(qū)域一般為橢圓形。邊緣發(fā)生接觸時(shí)將出現(xiàn)應(yīng)力集中，瞬時(shí)接觸區(qū)發(fā)生畸變，一端邊緣接觸導(dǎo)致橢圓一端畸變，兩端邊緣接觸導(dǎo)致橢圓兩端畸變。載荷較小時(shí)，接觸線上載荷分布于齒面中部，接觸區(qū)呈橢圓形，該接觸問題是無限邊界的赫茲接觸問題，接觸應(yīng)力可按赫茲點(diǎn)接觸彈性應(yīng)變模型計(jì)算；載荷較大時(shí)，齒面邊緣接觸，接觸橢圓出現(xiàn)畸變，該接觸問題變成有限邊界彈性接觸問題。

圖10 假設(shè)瞬時(shí)接觸形變區(qū)

接觸線上的載荷分布已知時(shí)，邊緣接觸問題可采用擬赫茲線接觸的辦法求解，即沿接觸線長(zhǎng)度方向?qū)⒔佑|區(qū)細(xì)分為m個(gè)長(zhǎng)度為Δl的線接觸單元，如圖11所示。根據(jù)ease-off曲面的拓?fù)浣Y(jié)構(gòu)，將每個(gè)單元都當(dāng)作圓柱同平面的接觸。由于Δl很小，單元j上的作用力可按均布載荷qj處理，即qj=Fj/Δl。

圖11 接觸區(qū)有限元模型

圓柱體與平面發(fā)生赫茲接觸，單元j接觸區(qū)中點(diǎn)處應(yīng)力最大，為

(9)

(10)

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

接觸區(qū)半寬為

(11)

沿接觸區(qū)y方向任意點(diǎn)的接觸應(yīng)力為

(12)

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

(a)擬赫茲接觸

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

(a)擬赫茲接觸LTCA

(a)擬赫茲接觸LTCA

5 結(jié)論

(1)利用曲面多項(xiàng)式、齒條-齒輪等切共軛產(chǎn)形原理能準(zhǔn)確構(gòu)建數(shù)值齒面和ease-off差齒面模型，實(shí)現(xiàn)對(duì)齒面拓?fù)浣Y(jié)構(gòu)形狀與性質(zhì)的控制。

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

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