齒輪齒條動態(tài)嚙合及傳動分析的SPH數(shù)值算法研究

唐宇峰，胡光忠，曹修全，陽明君

(1.四川輕化工大學(xué) 機(jī)械工程學(xué)院，四川 宜賓 644000;2.重大危險源測控四川省重點實驗室，成都 610031)

用傳統(tǒng)的解析法進(jìn)行齒輪齒條嚙合及傳動力學(xué)研究時，由于大量的簡化難以正確反映實時的動態(tài)力學(xué)行為。而近年來發(fā)展起來的以網(wǎng)格為主的數(shù)值方法(如有限元法、有限差分法等)則可以對嚙合接觸區(qū)域的力學(xué)狀態(tài)及沖擊過程進(jìn)行計算，如：Hwang等[1]采用有限元發(fā)法通過齒圓柱齒輪和斜齒圓柱齒輪兩個例子對一對齒輪在旋轉(zhuǎn)過程中的接觸應(yīng)力進(jìn)行了分析，并指出考慮一對配合齒輪的接觸應(yīng)力的齒輪設(shè)計比美國AGMA標(biāo)準(zhǔn)更加準(zhǔn)確；賈海濤等[2]采用有限元法分析了圓柱斜齒輪齒面接觸動應(yīng)力；閆曉青等[3]采用有限元法討論了齒條內(nèi)夾雜物、齒面蝕坑等缺陷對齒輪齒條運行的影響。然而，固定的網(wǎng)格在實時追蹤輪齒傳動，特別是當(dāng)面臨大變形及裂紋擴(kuò)展時容易產(chǎn)生網(wǎng)格畸變而造成計算誤差甚至計算中止，因此限制了其應(yīng)用。

光滑粒子流體動力學(xué)(smoothed particle hydrodynamics，SPH)方法是近年來逐漸發(fā)展起來的一種純拉格朗日無網(wǎng)格計算方法，目前已應(yīng)在許多計算力學(xué)領(lǐng)域，如爆炸力學(xué)、流體力學(xué)和固體力學(xué)等領(lǐng)域得到了應(yīng)用，如：強洪夫等[4]將SPH方法應(yīng)用于寬速域巖石侵徹問題，得到了花崗巖靶板受碰撞侵徹的大應(yīng)變、高應(yīng)變率變形；陳佩佩等[5]采用SPH方法對非飽和巖土介質(zhì)滲流問題進(jìn)行了討論；王維國等[6]采用SPH-FEM(smoothed particle hydrodynamics-finite element method)耦合的方法，研究了小當(dāng)量集團(tuán)裝藥在不同埋置深度時的土中爆炸效應(yīng)；駱釗等[7]針對復(fù)雜固壁邊界提出了改進(jìn)，并將其應(yīng)用于沙土滑坡大變形的模擬；Frissane等[8]采用SPH方法研究了射彈撞擊對人體的損害；Becker等[9]通過模擬和試驗的方法對比了SPH方法、光滑粒子伽遼金法和有限元法在進(jìn)行高速沖擊計算中的優(yōu)劣。近年來，學(xué)者們逐步將SPH方法引入到機(jī)械工程領(lǐng)域，并取得了一些可喜的成果。如：熱合買提江等[10]討論了SPH分析中的齒輪建模方法；Imin等[11-12]研究了兩齒輪嚙合和相互耦合接觸的SPH計算方法，將嚙合粒子之間的作用力視為一種系統(tǒng)內(nèi)力，分析了一對齒輪沖擊嚙合過程中最大動應(yīng)力隨時間的變化規(guī)律；Keller等[13]研究了旋轉(zhuǎn)直齒圓柱齒輪在噴油沖擊過程中的復(fù)雜兩相流動，并預(yù)測了3個不同射流傾角的直齒圓柱齒輪上噴油沖擊過程中的流動現(xiàn)象。然而總體來說，SPH方法在機(jī)械領(lǐng)域內(nèi)的研究尚屬于起步階段，特別是對類似齒輪齒條等具有復(fù)雜幾何形狀的動態(tài)嚙合及傳統(tǒng)分析研究還鮮見報道。因此，構(gòu)建基于SPH的新型齒輪齒條動態(tài)嚙合及傳動的數(shù)值計算方法：一方面能夠為齒輪、齒條、軸承等復(fù)雜接觸界面問題的動態(tài)力學(xué)分析提供科學(xué)依據(jù)及新的計算方法；另一方面也能夠推動SPH方法在該領(lǐng)域內(nèi)的進(jìn)一步發(fā)展。

采用SPH進(jìn)行齒輪齒條嚙合及傳動分析時，其基本思想是將齒輪及齒條模型離散為有限個無網(wǎng)格連接的離散點，而齒輪及齒條的物理力學(xué)參數(shù)及動態(tài)屬性，如密度、質(zhì)量、體積、位置、速度、加速度等都有離散點攜帶，并通過光滑核函數(shù)及質(zhì)量守恒、動量守恒和狀態(tài)方程建立離散點間的聯(lián)系。相對于傳統(tǒng)的基于網(wǎng)格的方法，SPH方法由于離散點間無網(wǎng)格連接，所以不會產(chǎn)生單元扭曲、網(wǎng)格畸變以及網(wǎng)格劃分奇異性帶來的計算誤差，特別適合求解高速傳遞的動態(tài)力學(xué)行為。本文在傳統(tǒng)SPH方法的基礎(chǔ)上，針對齒輪齒條嚙合及傳動過程的特點，通過對嚙合方式、受力特點等改進(jìn)，建立了適用于齒輪齒條嚙合及傳動分析的SPH方法。最后用某簡化的齒輪齒條動態(tài)傳動過程為例，對齒輪齒條動態(tài)嚙合及傳動過程進(jìn)行了分析。結(jié)果表明，該SPH方法可以較為準(zhǔn)確地計算齒輪齒條嚙合及傳動的全過程，揭示齒輪齒條機(jī)構(gòu)發(fā)生齒面點蝕及齒根破壞的力學(xué)原理，能夠為該領(lǐng)域的計算提供一種新的計算方法，另外也推動了SPH在復(fù)雜嚙合及傳動機(jī)械領(lǐng)域內(nèi)的進(jìn)一步發(fā)展。

1 齒輪齒條嚙合及傳動分析的SPH計算基礎(chǔ)

1.1 SPH理論基礎(chǔ)

在SPH中，其方程構(gòu)造主要有兩個關(guān)鍵步驟組成。其一就是將任意宏觀變量函f(x)數(shù)通過積分進(jìn)行表示，而函數(shù)f(x)的積分表示可通過對支持域內(nèi)相鄰粒子的值進(jìn)行累加求和近似，可用方程[14]表示為

(1)

(2)

而粒子近似就是將相關(guān)的連續(xù)積分表達(dá)式轉(zhuǎn)化為支持域內(nèi)所有與粒子疊加求和的離散化方程，即

(3)

式中：j為粒子編號；mj，ρj分別為粒子的質(zhì)量及密度；N為其支持域內(nèi)臨近粒子的總數(shù)。即任一函數(shù)值可用其支持域內(nèi)所有粒子相對應(yīng)的函數(shù)值進(jìn)行加權(quán)平均而近似表達(dá)。同樣的，函數(shù)f(x)的空間導(dǎo)數(shù)?·f(x)的粒子近似可表達(dá)為

(4)

式(1)～式(4)中，光滑函數(shù)Wij的選擇對離散近似起到了重要的作用，本文中選取的是一種新型的四次光滑核函數(shù)[15]，可表達(dá)為

W(R,h)=αd×

(5)

其光滑函數(shù)及其一階、二階導(dǎo)數(shù)，如圖1所示。

圖1 新的分段四次核函數(shù)及其一階、二階導(dǎo)數(shù)Fig.1 The NQS kernel and its one and two-order derivative

1.2 基于理想彈性模型的齒輪齒條嚙合及傳動分析方程的SPH離散

利用1.1節(jié)所述SPH的基本理論，SPH形式下的連續(xù)性方程和動量方程可以表示為

(6)

(7)

式中：v為速度；α和β為坐標(biāo)方向；σαβ為總應(yīng)力張量；Fα為外力引起的加速度，如重力、齒輪齒條相互作用力及固壁邊界作用力等；Πij為Monaghan[16]提出的人工黏度，可以防止粒子相互接近時產(chǎn)生的非物理穿透，其表達(dá)式為

(8)

(9)

(10)

同時，引入了XSPH的速度平均方法公式

式中，ε為修正因子，本文中選0.3。

根據(jù)線彈性動力學(xué)理論，SPH形式的本構(gòu)方程可以表示為

(11)

式中：G為剪切模量；K為體積模量；εi,γγ為主應(yīng)變之和；δαβ為狄克拉函數(shù)；ωαβ為Jaumann應(yīng)力率，其表達(dá)式為

(12)

2 齒輪齒條嚙合及傳動作用力的SPH實現(xiàn)

在齒輪齒條沖擊及傳動過程中同時存在連續(xù)性及非連續(xù)性，即齒輪、齒條構(gòu)件內(nèi)部為連續(xù)體，而齒輪齒條構(gòu)件之間接觸具有離散性。雖然當(dāng)齒輪齒條粒子處于嚙合狀態(tài)時，其嚙合部分可采用連續(xù)體方法、核估算進(jìn)行計算，然而由于齒輪齒條間不斷處于嚙合及脫離嚙合狀態(tài)，當(dāng)其作用關(guān)系從嚙合向脫離嚙合狀態(tài)轉(zhuǎn)變時，采用核估算方法會在齒輪齒條間會產(chǎn)生與實際情況不符的拉伸力，與實際嚙合狀態(tài)不符。在本文中，對齒輪齒條嚙合及傳動過程進(jìn)行了以下考慮：

(1)如上所述，齒輪和齒條在各自構(gòu)件內(nèi)部為連續(xù)體，因此其構(gòu)件內(nèi)部的作用力可以通過核估算進(jìn)行計算。而在齒輪及齒條發(fā)生相互嚙合作用時，也可通過核估算方法來提高接觸力的作用的計算精度。然而，為正確表達(dá)齒輪齒條粒子嚙合狀態(tài)，避免齒輪齒條的提前作用，交界處粒子在計算嚙合沖擊力時，其影響域之和應(yīng)等于齒輪齒條嚙合粒子的平均間距，即

(13)

(2)當(dāng)齒輪齒條粒子從嚙合狀態(tài)向非嚙合狀態(tài)轉(zhuǎn)化時，雖然此時齒輪齒條粒子仍處于嚙合狀態(tài)，但如果仍按本章方法進(jìn)行處理，會在齒輪齒條粒子間產(chǎn)生與實際情況不符的拉應(yīng)力。因此，本文在此狀態(tài)下將連續(xù)方程中的粒子間相互嚙合作用力定義為0，并且在齒輪齒條粒子處于嚙合狀態(tài)時，給予齒輪及齒條粒子一個具有彈簧特性的相互作用的外力，該力與齒輪齒條粒子的變形及彈性模量有關(guān)，即距離越近說明變形越大，則外力隨著變形及彈性模量的增大而越大，而距離越遠(yuǎn)說明變形越小，則外力隨著變形及彈性模量的減小而越小。在二維條件下，此外力可表達(dá)為

(14)

式中：Feij,α為齒條粒子j對相嚙合的齒輪粒子i的外力；rr0為齒輪齒條嚙合粒子初始的平均間距；rr為齒輪齒條嚙合粒子某時刻的實際間距；sign為符號函數(shù)，即括號內(nèi)若為正則返回1，為負(fù)則返回-1；D為常數(shù)，取相對速度平方相等的量級；E為彈性模量。

基于本章理論，編制了適用于齒輪、齒條等高速嚙合及傳動分析的SPH計算程序，其計算流程如圖2所示。

圖2 采用SPH方法進(jìn)行齒輪齒條嚙合及傳動分析的計算流程Fig.2 Calculation flow of gear rack meshing and transmission analysis using SPH method

3 數(shù)值算例

本文以一對簡化的齒輪齒條算例為例進(jìn)行分析。算例中齒輪模數(shù)為4 mm，壓力角為20°，齒寬為20 mm，齒輪沿順時針旋轉(zhuǎn)，轉(zhuǎn)速為15 rad/s，為避免齒輪齒條過度沖擊，給予齒條沿X負(fù)方向，大小為0.353 m/s的初速度。齒輪齒條的材料均為球磨鑄鐵，彈性模量為173 GPa，泊松比為0.3，密度為7 800 kg/m3。為節(jié)約計算成本，對齒輪齒條結(jié)構(gòu)進(jìn)行了一定的簡化，粒子初始間距為0.001 m，齒輪齒條共離散1 053個粒子，其中齒輪粒子為395個，齒條粒子為658個，另外固壁邊界虛粒子322個，時間步長為4×10-8s，計算1×106個時間步，共計0.04 s，粒子初始裝配狀態(tài)如圖3所示。粒子在0.008 s，0.024 s，0.04 s時刻的運動狀態(tài)如圖4～圖6所示。

圖3 齒輪齒條粒子初始狀態(tài)Fig.3 The initial state of particles of gear and rack

圖4 齒輪齒條在不同計算時刻時的運動狀態(tài)Fig.4 The motion state of gear and rack at different calculation times

為觀察齒輪齒條運動狀態(tài)，在齒輪齒條上隨機(jī)各選取了一個粒子觀察其速度隨時間的變化情況，圖5(a)為齒輪齒條粒子在X方向的速度變化曲線(vx)，圖5(b)為齒輪齒條粒子總速度變化曲線。由圖5可知，由于本算例中齒輪齒條初始速度不同，在運動初期產(chǎn)生了速度波動，而隨著計算時刻的增加，齒輪齒條速度趨近恒定，并與實際傳動速度相符。

由圖5可知，在約0.072 s時由于齒輪齒條間相互嚙合作用發(fā)生速度值的較大變化，選取該點查看其Mises應(yīng)力分布，如圖6所示?？芍藭r齒輪和齒條在此時產(chǎn)生嚙合額，且在嚙合點處應(yīng)力最大。隨后，應(yīng)力向四周進(jìn)行傳播，并在齒根處產(chǎn)生應(yīng)力集中，其應(yīng)力大小和發(fā)生時刻略小于嚙合點處應(yīng)力。分析結(jié)果表明，齒面點蝕及齒根破壞是齒輪齒條機(jī)構(gòu)破壞的主要形式。

圖5 齒輪齒條粒子在傳動過程中的速度變化規(guī)律Fig.5 The law of the velocity change of the gear and rack particles during transmission

圖6 齒輪齒條嚙合過程Mises應(yīng)力分布及傳遞Fig .6 Stress distribution and transmission of mises during gear bar meshing process

為了驗證結(jié)果的合理性，采用有限元方法進(jìn)行了對比，圖7為有限元計算的內(nèi)力結(jié)果及傳遞過程，可見其應(yīng)力大小、應(yīng)力傳遞過程與SPH方法的結(jié)果基本一致，證明了本文算法可以較好地模擬齒輪齒條嚙合及傳遞過程。

圖7 有限元方法的Mises應(yīng)力分布及傳遞Fig.7 The Mises stress distribution of the gear and rack using finite element method

如本章分析結(jié)果可知，通過本文提出的SPH方法，可以動態(tài)獲取齒輪齒條嚙合在任意時刻的運動狀態(tài)、應(yīng)力分布、速度分布、運動過程中最大應(yīng)力等參數(shù)，揭示結(jié)構(gòu)破壞的原因，可為齒輪的優(yōu)化設(shè)計、強度校核、判定損傷影響等提供可靠依據(jù)。然而應(yīng)該指出的是，為了減少計算成本，在本算例中對齒輪進(jìn)行了簡化，若要對完整的齒輪齒條進(jìn)行分析，則應(yīng)建立更加完善的齒輪齒條模型，同時研究更為高效的并行算法提高計算效率。

4 結(jié) 論

(1)本文提出了一種基于SPH的齒輪齒條動態(tài)嚙合及傳動分析的SPH算法，該算法可以解決以往SPH算法僅能分析準(zhǔn)靜態(tài)嚙合問題的不足，能夠有效地模擬齒輪齒條動態(tài)嚙合及傳動數(shù)值分析問題，推動了SPH在該領(lǐng)域內(nèi)的進(jìn)一步發(fā)展。

(2)對一個齒輪齒條動態(tài)嚙合及傳動的數(shù)值算例進(jìn)行了分析，分析結(jié)果有效揭示了齒輪齒條機(jī)構(gòu)發(fā)生齒面點蝕及齒根破壞的力學(xué)原理，并與有限元方法分析結(jié)果對比驗證了本文提出的SPH算法的正確性。該法可廣泛應(yīng)用于齒輪、齒條、軸承等復(fù)雜動態(tài)接觸機(jī)構(gòu)的問題分析，為該領(lǐng)域的數(shù)值分析提供了一種新的思路。

(3)由于本文基于串行SPH程序進(jìn)行計算，完整的齒輪齒條模型會導(dǎo)致計算成本過高而難以進(jìn)行。因此，發(fā)展高效的并行算法是未來需要進(jìn)一步做的工作，從而進(jìn)一步推動SPH在該領(lǐng)域內(nèi)的進(jìn)一步發(fā)展。