李少奇，李九一，馮志強,2

(1.西南交通大學 力學與工程學院， 成都 610031；2.巴黎-薩克雷大學 力學與熱能實驗室， 法國 埃夫里 91020)

0 引言

漸開線圓柱齒輪包括直齒輪、斜齒輪和人字齒輪，由于其結(jié)構(gòu)形式簡單，設(shè)計制造相對容易，因此被廣泛應用于各種機械裝置中。自Wilcox[1]于1973年開始使用有限元來確定齒輪在載荷下的應力以來，有限單元方法(FEM)已經(jīng)成為齒輪傳動研究和開發(fā)中的一項重要技術(shù)。在過去的幾年中，它被廣泛用于確定齒輪在不同載荷[2]下的接觸和彎曲應力。由于圓柱齒輪齒廓的特殊性，其網(wǎng)格模型的建立往往是齒輪傳動有限元分析中的一大難題。國內(nèi)外學者進行齒輪傳動分析用到的齒輪有限元模型一般通過以下3種方法建立：第1種是通過第三方CAD建模軟件建立齒輪模型，并導入商業(yè)有限元軟件劃分齒輪網(wǎng)格進行分析，此方法較為繁瑣，且模型兼容性較差，例如宋立偉等[3]應用Pro/E軟件建立齒輪的模型，并于Abaqus軟件中進行有限元分析。第2種是在有限元軟件中參數(shù)化建模，利用軟件自身的程序化語言來參數(shù)化控制模型的建立，此方法較為便捷，但使用范圍有限，例如劉洋等[4]借用Ansys提供的APDL語言開發(fā)了直齒輪的參數(shù)化建模模塊。第3種是通過編寫程序直接生成齒輪單元和節(jié)點網(wǎng)格模型，杜金成[5]利用VB編程語言開發(fā)了漸開線直齒圓柱齒輪參數(shù)化修形應用軟件，實現(xiàn)了對Ansys的直接調(diào)用。該方法較為便捷，但是由于其借用了APDL語言，其網(wǎng)格的劃分仍然依靠Ansys軟件自動劃分，在網(wǎng)格精度方面仍略有欠缺，且其每次生成模型都需要重新調(diào)用Ansys軟件，步驟略顯繁瑣。因此一種嵌入在CAE軟件內(nèi)部，并且可以直接生成齒輪網(wǎng)格模型的方法具有十分重要的意義。LiToSim作為一款國內(nèi)自主研發(fā)的CAE軟件擁有強大的前后處理模塊及固體力學求解器，可以實現(xiàn)幾何建模，復雜幾何體導入，網(wǎng)格劃分及其后處理顯示等功能，且其求解精度可媲美國際主流商用CAE軟件。目前LiToSim可以實現(xiàn)結(jié)構(gòu)靜態(tài)線性/非線性分析、動態(tài)分析、模態(tài)分析、熱力學分析等。

為了解決齒輪網(wǎng)格模型建立復雜這一難題，本文基于齒輪輪廓方程，利用超限映射法劃分齒輪六面體網(wǎng)格，生成參數(shù)化齒輪網(wǎng)格劃分程序并集成于國產(chǎn)CAE軟件LiToSim中。通過在LiToSim軟件中嵌入網(wǎng)格模型自動生成技術(shù)，可以快速建立齒輪六面體網(wǎng)格模型，極大加快有限元求解的進程。

1 齒輪模型的數(shù)學表達式

圖1 圓柱齒輪輪廓曲線

1.1 齒根過渡曲線方程

CD段曲線是齒根過渡曲線，這段曲線由加工齒輪的刀具來決定。齒條形加工刀具齒廓頂部具有2個圓角，所以其齒根過渡曲線是一段延伸漸開線的等距曲線。由于BC和DE線段要用到CD段曲線的首尾點，因此首先建立其參數(shù)方程如下：

(1)

(2)

1.2 齒頂圓圓弧方程

AB曲線是一段齒頂圓圓弧，其參數(shù)方程為：

(3)

(4)

1.3 漸開線方程

BC段曲線為一段漸開線，其參數(shù)方程為：

(5)

(6)

式中:φ為分度圓齒厚所對應的中心角，φ=(0.5πm+2xbmtanα)/0.5mz；αi是變參數(shù)，取值范圍為αb≥αi≥αc，αb，αc分別由下式求得：

(7)

上式為超越方程，可用計算機程序迭代求解，其中xB、yB、xC、yC為B、C兩點的坐標。

1.4 齒根圓圓弧方程

DE段曲線為一段齒根圓圓弧，其參數(shù)方程為：

(8)

1.5 漸開線圓柱齒輪截面輪廓

由式(1)(3)(5)(8)即可獲得漸開線圓柱齒輪輪廓曲線。由于齒輪輪廓線方程是以弧度為自變量，因此以均勻的弧度為自變量獲得的BD曲線段節(jié)點坐標是不理想的。為了解決這一問題，應加大在輪廓線方程第一次取點的密度，獲得一系列致密的節(jié)點坐標，并對這些節(jié)點坐標進行判斷，即當節(jié)點之間累加的距離達到程序設(shè)定的某一個值時(該值由節(jié)點密度和偏離率確定)，保存這個節(jié)點坐標到矩陣當中，最終獲得所需的曲線節(jié)點。

將獲得的輪廓線對稱并加以補充得到如圖2所示的齒輪截面輪廓。輪廓由8條邊界曲線組成，其中P1P2曲線為齒頂圓弧，P1P3和P2P4為漸開線，P3P5和P4P6為齒根過渡曲線和齒根圓圓弧(由于齒根圓圓弧曲線長度過短，不再將其作為獨立曲線段)，P5P7和P6P8為齒根圓半徑去除齒輪內(nèi)徑后的線段，P7P8為齒輪內(nèi)徑圓弧。

圖2 圓柱齒輪截面輪廓圖

2 有限元模型構(gòu)建

下文將以圖2所示的齒輪模型為例，介紹漸開線圓柱齒輪有限元網(wǎng)格模型的構(gòu)建方法。

2.1 齒輪截面區(qū)域劃分

由圖2可以看出，齒輪截面輪廓是一個具有復雜邊界的非凸域，因此在其上生成網(wǎng)格是一個非常棘手的問題。解決這種問題通常的辦法就是多塊網(wǎng)格劃分，即通過將復雜域分解為簡單的子域，然后在子域進行局部網(wǎng)格劃分，從而獲得整個齒輪截面網(wǎng)格。

為了提高生成齒輪網(wǎng)格的質(zhì)量，在進行齒輪截面區(qū)域劃分時需要滿足如下幾個要求：一是子域在其拓撲上必須等價為一個正方形。二是為了避免自折疊和變形的網(wǎng)格，子域必須盡可能地接近凸區(qū)域，并避免邊界曲線斜率的不連續(xù)。三是為了簡化邊界曲線的離散化和保證網(wǎng)格跨越邊界時的連續(xù)性，任何邊界曲線都不能被2個以上的子域共享。

考慮到上述要求，最終獲得了如圖3所示的區(qū)域劃分。該方法通過邊界點Pi{i=1,2,…,8}與輔助點Qi{i=9,10,…,16}的連線，將齒輪截面區(qū)域劃分為8個子域，輔助點由以下步驟確定：

圖3 單齒區(qū)域劃分圖

步驟1Q9為頂部曲線的中點。

步驟2Q10是P4點的垂線和齒輪Y軸的交點。

步驟3Q11和Q12是左右圓角曲線的中點。

步驟4Q13是左右圓角曲線中點的垂線與Y軸的交點。

步驟5令P3Q10和P4Q10曲線恒等于P5Q14和P6Q15曲線，從而獲得Q14和Q15點。

步驟6Q16為Y軸與底部曲線的交點。

通過以上步驟得到的子域滿足了本節(jié)開頭所提到的所有要求。此外，它最小化了子域?qū)吳€長度上的差異，有助于提高網(wǎng)格單元的長寬比，同時也最大化了邊界曲線在子域頂點上的正交性，有助于減少網(wǎng)格單元的偏斜度。

2.2 網(wǎng)格節(jié)點信息的獲取

2.2.1利用超限映射法獲得網(wǎng)格節(jié)點信息

相較于非結(jié)構(gòu)化網(wǎng)格，結(jié)構(gòu)化網(wǎng)格[7]具有網(wǎng)格速度生成快、網(wǎng)格質(zhì)量高、數(shù)據(jù)結(jié)構(gòu)簡單，可以更好地擬合幾何區(qū)域等特點，其缺點是僅適應于規(guī)則形狀的區(qū)域。然而在上節(jié)中，齒輪已經(jīng)被均勻地劃分為規(guī)則區(qū)域，這個缺點也就被克服了。

在結(jié)構(gòu)化網(wǎng)格劃分方法中，超限映射法[8]因其簡單的算法脫穎而出，即在由任意4條曲線圍成的區(qū)域之內(nèi)，通過建立物理平面[如圖4(a)]與計算平面[如圖4(b)]的映射函數(shù)來實現(xiàn)計算平面對物理平面的逼近，這就是超限映射法。理論上講，該方法比其他網(wǎng)格劃分方法對曲線邊界的逼近更為精確。

對于圖4(a)所示的需要劃分網(wǎng)格的區(qū)域ijmn，其中i，j分別代表上下2條邊，m，n分別代表左右2條邊，其4條邊界的坐標是已知的。例如曲邊m，由5個節(jié)點坐標來確定，則計算平面m邊也應由5個節(jié)點坐標確定，且其節(jié)點之間距離占m邊總長的比例與物理平面相等。同理其他3條邊也做類似處理。由于計算平面區(qū)域是一個邊長為1的正方形，因此將物理平面邊界的節(jié)點坐標(x,y)映射至計算平面就會得到計算平面邊界的節(jié)點坐標(u,v)的信息。

圖4 物理平面(a)與計算平面(b)示意圖

設(shè)xi(u)、yi(u)、xj(u)、yj(u)、xm(v)、ym(v)及xn(v)、yn(v)分別代表計算平面邊界對應的物理平面的邊界的節(jié)點坐標，則(x,y)與(u,v)的映射變換關(guān)系如下：

(9)

式中： 0≤u≤1，0≤v≤1，(xi(1)、yi(1))，(xj(1)、yj(1))，(xj(0)、yj(0))及(xi(0)、yi(0))分別代表圖4(a)中右上角順時針數(shù)4個點的坐標。

通過式(19)的映射關(guān)系，在保持節(jié)點間的拓撲關(guān)系不變的情況下，將計算平面的節(jié)點映射至物理平面，從而得到如圖4(a)所示的形狀規(guī)則的結(jié)構(gòu)化網(wǎng)格。

應用超限映射法將圖3中8個子域劃分網(wǎng)格，最終得到如圖5所示的網(wǎng)格節(jié)點信息，從圖中可以看出由于偏離率系數(shù)的存在，使得漸開線部分的節(jié)點是由密到疏，這使齒輪網(wǎng)格的長寬比更加合理。

圖5 單齒網(wǎng)格節(jié)點信息示意圖

2.2.2齒輪網(wǎng)格節(jié)點的旋轉(zhuǎn)與拉伸

由于齒輪三維模型是在二維的基礎(chǔ)上拉伸得到的，因此首先建立二維齒輪網(wǎng)格節(jié)點。對于上節(jié)的單齒網(wǎng)格節(jié)點，為了防止其旋轉(zhuǎn)后節(jié)點產(chǎn)生重復，刪除圖5單齒網(wǎng)格中加黑的節(jié)點后再進行旋轉(zhuǎn)。齒輪嚙合往往只涉及到幾個齒，所以既可以旋轉(zhuǎn)n-1次(n為齒數(shù))得到整個二維齒輪網(wǎng)格節(jié)點(如圖6所示，圖6模型由LiToSim軟件參數(shù)化建模生成)，也可旋轉(zhuǎn)指定次數(shù)得到想要顯示齒輪的齒的個數(shù)。

圖6 漸開線圓柱齒輪網(wǎng)格節(jié)點模型界面

三維圓柱齒輪包括直齒輪、斜齒輪、人字形齒輪，其中直齒輪網(wǎng)格節(jié)點生成只需在圖6節(jié)點信息的基礎(chǔ)上在z軸方向作簡單復制即可。而斜齒輪與人字形齒輪網(wǎng)格節(jié)點的生成需要在z軸復制的同時旋轉(zhuǎn)一定的角度即可。

至此，圓柱齒輪的所有網(wǎng)格節(jié)點信息已全部生成。

2.3 單元節(jié)點編號

最后一步是單元節(jié)點編號，二維單元采用CPE4單元即四節(jié)點雙線性平面應變四邊形單元。其單元內(nèi)每個節(jié)點都按照逆時針的順序排列，如圖7(a)所示。

三維單元則采用了C3D8單元類型即八節(jié)點線性六面體單元，與CPE4單元編號方式相反，C3D8單元由上下2個平行面上的節(jié)點順時針方式編號決定，如圖7(b)所示。

圖7 CPE4單元和C3D8單元節(jié)點示意圖

至此，完整的齒輪CPE4四邊形網(wǎng)格和C3D8六面體有限元網(wǎng)格便建立完成，其中三維直齒輪、斜齒輪、人字形齒輪模型如圖8所示。

圖8 三維漸開線圓柱齒輪網(wǎng)格模型示意圖

3 程序化實現(xiàn)

以上所有操作均以C++語言編寫，在Qt平臺上運行并生成軟件[9]，嵌入LiToSim平臺。程序生成齒輪有限元網(wǎng)格模型的過程可以由圖9程序框圖概括表達。圖9展示了整個圓柱齒輪有限元網(wǎng)格模型構(gòu)建的基本流程，程序開始時，第一步需要輸入漸開線圓柱齒輪基本參數(shù)和網(wǎng)格密度等信息。第二步計算齒輪的節(jié)點坐標，并將節(jié)點xyz坐標信息儲存到事先定義好的VCORG坐標矩陣中。第三步是單元節(jié)點編號，構(gòu)建C3D8六面體單元或者CPE4四邊形平面應變單元，并將單元節(jié)點編號和單元編號存儲至事先定義好的ELEMENT矩陣空間中。第四步將齒輪單元及節(jié)點矩陣導至LiToSim中顯示。至此，漸開線圓柱齒輪有限元網(wǎng)格模型參數(shù)化建模模塊已經(jīng)建立完成。

圖9 模型構(gòu)建基本流程框圖

4 齒輪接觸驗證

4.1 赫茲理論解

圓柱齒輪赫茲接觸應力的理論解[6]的計算過程如下所示：

(10)

式中:σH為齒輪嚙合處的接觸應力；d1為小齒輪分度圓直徑；b為齒寬；u為齒數(shù)比，u=z2/z1。

ZH為節(jié)點區(qū)域系數(shù)，考慮到接觸節(jié)點處齒廓曲率的變化對接觸應力的影響，該系數(shù)將分度圓上切向力折算為節(jié)圓的法向力。其值可以查表獲得，也可通過下式計算：

(11)

ZE為彈性系數(shù)，該系數(shù)考慮了泊松比ν和彈性模量E對赫茲應力的的影響。其值可由下式計算：

(12)

Zε為重合度系數(shù)，該系數(shù)考慮了齒輪多齒嚙合時對單位齒寬載荷的影響。由于文中實例只考慮單齒嚙合時的狀況，因此該系數(shù)不在考慮范圍內(nèi)，此處其值為1。

Zβ為螺旋角系數(shù)，該系數(shù)考慮了斜齒輪接觸線的傾斜對接觸力的影響，若為直齒輪，其值為1。Ft為名義切向力，可由下式計算：

(13)

實驗所需的齒輪參數(shù)如表1所示，材料參數(shù)如表2所示，通過上述參數(shù)及式子可得到漸開線直齒圓柱齒輪接觸應力的理論解，如表3所示。

表1 齒輪參數(shù)

續(xù)表(表1)

表2 材料參數(shù)

表3 赫茲理論解

4.2 模型前處理及求解

有限元模擬實例來自于文獻[10]，采用厚度為0.1 mm的三維直齒輪，齒輪各項參數(shù)上節(jié)中已給出，使用LiToSim軟件建立模型并設(shè)置模型參數(shù)。模型整體結(jié)構(gòu)如圖10所示，由左側(cè)的主動齒輪和右側(cè)的從動齒輪組成。齒輪傳動模型一般由3個或者多個齒構(gòu)成，此處考慮齒輪嚙合的一般情況，即只有單對齒嚙合的情況，因此為了保證模擬的質(zhì)量和計算結(jié)果的準確性，采用三齒齒輪模型。

圖10 通用齒輪傳動的有限元模型整體結(jié)構(gòu)示意圖

為了更精確地模擬齒輪嚙合過程，將每個齒輪幾何形狀的邊緣耦合到了與齒輪旋轉(zhuǎn)中心一致的參考點O1和O2上。如圖10所示。2個齒輪的參考點用于指定有限元模型的載荷和邊界條件。在2個參考點處沿x和y方向的位移受到限制。在主動齒輪的參考點處，允許繞z方向旋轉(zhuǎn)，其中施加逆時針扭矩T，而在從動齒輪的參考點處圍繞z方向的旋轉(zhuǎn)處則要被限制?？紤]到齒輪厚度較薄，因此限制了2個齒輪z向的位移。齒輪的接觸區(qū)域如圖10圓內(nèi)所示，接觸方法采用增廣拉格朗日法，偵測算法采用主從面，剛度因子設(shè)為1，摩擦因數(shù)為0.1。

為了驗證計算結(jié)果的準確性，將同樣的模型 在Abaqus中計算并對比其計算結(jié)果。

4.3 有限元分析結(jié)果

圖11為在Abaqus和LiToSim上進行計算后的Mises應力云圖，在接觸區(qū)域表面以下一定距離處，Mises應力達到了最大值，此外在右齒根圓角區(qū)域也產(chǎn)生較大的Mises應力[11]。這些應力達到最大值的區(qū)域如圖12(a)和圖13(a)所示。

圖12(b)為經(jīng)過LiToSim軟件計算獲得的圖11(a)區(qū)域A的Mises應力云圖。在網(wǎng)格模型相同的情況下，Misese應力最大值與Abaqus計算的結(jié)果相差1.84%，且其應力云圖分布也近乎一致。

圖11 Mises應力云圖

圖12 區(qū)域A的Mises應力云圖

圖13(b)為LiToSim獲得的圖11(a)區(qū)域B處的Mises應力云圖。從圖中可以觀察到該區(qū)域的最大Mises應力位于齒根表面，其值達到了238.2 MPa，與Abaqus計算結(jié)果相差2.1%，且其應力分布也幾乎完全一致，由以上結(jié)果充分證明了LiToSim計算結(jié)果的準確性。

圖13 區(qū)域B的Mises應力云圖

圖14為由Abaqus計算的齒輪相互嚙合時的接觸壓力，可以看出，齒輪的最大接觸壓力接近 1 248 MPa，由表3可知，齒輪的赫茲理論解為 1 236 MPa,其誤差僅有0.9%，由此證明了文章所提出的齒輪網(wǎng)格模型的準確性。

圖14 齒輪接觸應力云圖

表4展示了分別由Abaqus與LiToSim計算的A、B區(qū)域的Mises應力最大值和Abaqus計算的最大接觸壓力與赫茲理論解。由上述結(jié)果可知，通過在LiToSim中嵌入齒輪參數(shù)化建模的程序，不但免去了建立模型，劃分網(wǎng)格的繁瑣步驟，且模型的計算結(jié)果也能夠得到保證，證明了本文所提方法的可行性。

表4 最大應力與接觸壓力結(jié)果

5 結(jié)論

介紹了一種漸開線圓柱齒輪有限元網(wǎng)格模型精確化建模的方法，并以三維直齒輪為例通過Abaqus和LiToSim驗證了圓柱齒輪參數(shù)化建模方法的可行性與準確性。在實際的工程分析和優(yōu)化設(shè)計中，通過參數(shù)化建立漸開線圓柱齒輪網(wǎng)格模型，既避免了設(shè)計人員繁瑣的勞動，又提高了工作效率和設(shè)計精度，滿足了實際工程需要。