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

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

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

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

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

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

1.1 SPH理論基礎

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

(1)

(2)

而粒子近似就是將相關的連續積分表達式轉化為支持域內所有與粒子疊加求和的離散化方程，即

(3)

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

(4)

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

W(R,h)=αd×

(5)

其光滑函數及其一階、二階導數，如圖1所示。

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

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

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

(6)

(7)

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

(8)

(9)

(10)

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

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

根據線彈性動力學理論，SPH形式的本構方程可以表示為

(11)

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

(12)

2 齒輪齒條嚙合及傳動作用力的SPH實現

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

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

(13)

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

(14)

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

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

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

3 數值算例

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

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

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

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

由圖5可知，在約0.072 s時由于齒輪齒條間相互嚙合作用發生速度值的較大變化，選取該點查看其Mises應力分布，如圖6所示。可知此時齒輪和齒條在此時產生嚙合額，且在嚙合點處應力最大。隨后，應力向四周進行傳播，并在齒根處產生應力集中，其應力大小和發生時刻略小于嚙合點處應力。分析結果表明，齒面點蝕及齒根破壞是齒輪齒條機構破壞的主要形式。

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

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

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

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

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

4 結 論

(1)本文提出了一種基于SPH的齒輪齒條動態嚙合及傳動分析的SPH算法，該算法可以解決以往SPH算法僅能分析準靜態嚙合問題的不足，能夠有效地模擬齒輪齒條動態嚙合及傳動數值分析問題，推動了SPH在該領域內的進一步發展。

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

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