基于LiToSim的圓柱齒輪參數(shù)化建模與接觸分析

李少奇，李九一，馮志強(qiáng),2

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

0 引言

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

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

1 齒輪模型的數(shù)學(xué)表達(dá)式

圖1 圓柱齒輪輪廓曲線

1.1 齒根過渡曲線方程

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

(1)

(2)

1.2 齒頂圓圓弧方程

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

(3)

(4)

1.3 漸開線方程

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

(5)

(6)

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

(7)

上式為超越方程，可用計(jì)算機(jī)程序迭代求解，其中xB、yB、xC、yC為B、C兩點(diǎn)的坐標(biāo)。

1.4 齒根圓圓弧方程

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

(8)

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

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

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

圖2 圓柱齒輪截面輪廓圖

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

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

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

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

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

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

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

步驟1Q9為頂部曲線的中點(diǎn)。

步驟2Q10是P4點(diǎn)的垂線和齒輪Y軸的交點(diǎn)。

步驟3Q11和Q12是左右圓角曲線的中點(diǎn)。

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

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

步驟6Q16為Y軸與底部曲線的交點(diǎn)。

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

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

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

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

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

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

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

設(shè)xi(u)、yi(u)、xj(u)、yj(u)、xm(v)、ym(v)及xn(v)、yn(v)分別代表計(jì)算平面邊界對(duì)應(yīng)的物理平面的邊界的節(jié)點(diǎn)坐標(biāo)，則(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í)針數(shù)4個(gè)點(diǎn)的坐標(biāo)。

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

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

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

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

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

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

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

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

2.3 單元節(jié)點(diǎn)編號(hào)

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

三維單元?jiǎng)t采用了C3D8單元類型即八節(jié)點(diǎn)線性六面體單元，與CPE4單元編號(hào)方式相反，C3D8單元由上下2個(gè)平行面上的節(jié)點(diǎn)順時(shí)針方式編號(hào)決定，如圖7(b)所示。

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

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

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

3 程序化實(shí)現(xiàn)

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

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

4 齒輪接觸驗(yàn)證

4.1 赫茲理論解

圓柱齒輪赫茲接觸應(yīng)力的理論解[6]的計(jì)算過程如下所示：

(10)

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

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

(11)

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

(12)

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

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

(13)

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

表1 齒輪參數(shù)

續(xù)表(表1)

表2 材料參數(shù)

表3 赫茲理論解

4.2 模型前處理及求解

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

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

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

為了驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確性，將同樣的模型 在Abaqus中計(jì)算并對(duì)比其計(jì)算結(jié)果。

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

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

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

圖11 Mises應(yīng)力云圖

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

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

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

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

圖14 齒輪接觸應(yīng)力云圖

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

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

5 結(jié)論

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