兩齒輪正確嚙合和相互耦合接觸的SPH計算方法研究

熱合買提江·依明， 買買提明·艾尼

（1．新疆大學數學與系統科學學院，烏魯木齊 830046；2．新疆大學機械工程學院，烏魯木齊 830047）

兩齒輪正確嚙合和相互耦合接觸的SPH計算方法研究

熱合買提江·依明1， 買買提明·艾尼2

（1．新疆大學數學與系統科學學院，烏魯木齊 830046；2．新疆大學機械工程學院，烏魯木齊 830047）

基于在一對齒輪的接觸嚙合及其沖擊的動力學方程，通過調整核函數影響半徑提出了一對齒輪正確嚙合和相互耦合接觸的SPH計算方法，編程實現了正確嚙合一對齒輪模型離散成SPH粒子的SPH前處理程序建立了一對嚙合齒輪的SPH離散模型并改變工況下進行了嚙合沖擊過程的動態數值仿真，分析了一對齒輪沖擊嚙合過程中的最大動應力隨時間的變化規律和齒廓面上動應力的分布及變化過程，對比探討了嚙合沖擊過程中的動應力傳播及接觸區域附近的接觸動應力分布和變化過程。對齒輪傳動系統的動態接觸強度評價和優化設計提供新的SPH數值算法。

齒輪嚙合；SPH計算方法；嚙合沖擊；動應力

在承受載荷和傳遞動力的齒輪傳動的研究中，用解析方法和簡單的數值方法計算和仿真時，做了大量的簡化，不能正確的反映實際情況［1］。計算機技術發展和新的數值算法的出現，對齒輪傳動分析方面顯示了強大的優勢。比如，用彈性接觸有限元法可以分析齒輪接觸區域位移和應力的變化［2－3］。用解析法和有限元法能解決嚙合輪齒的齒面變形、齒面載荷分布，可以分析齒面變形對傳動的影響［4］。基于ANSYS／LS-DYNA對齒輪嚙入沖擊過程可以進行數值仿真［5］。但基于網格（有限差分，有限元、邊界元等）法分析時，用固定網格很難實時追蹤齒輪傳動過程。特別是，當面臨大變形和裂紋擴散時，出現網格纏結和扭曲，因此在求解過程中需要實時進行跟蹤和網格再劃分，這樣不僅計算費用昂貴、費時，而且會使計算精度受損。近年來發展較快的無網格方法之一：光滑粒子流體動力學（Smoothed Particle Hydrodynamics，SPH）方法引人注目，被廣泛應用在流體和固體力學領域。SPH法的基本思想是將齒輪模型離散成有限多個有密度、體積、位置和速度等屬性的粒子［6－7］，計算過程粒子之間不需要網格連接，而通過核函數將它們聯系起來，從而保證了物體的連續性，同時滿足了能量、動量和質量的守恒性。由于無網格方法可以完全拋開網格并用拉格朗日方法來描述，因此直接對物體的實時運動和變形過程的進行追蹤和坐標更新、容易編程實現，從而避免網格奇異性帶來的計算誤差、保證了動態變形和沖擊破碎的計算精度。雖然SPH法在一對物體高速相互沖擊接觸和破碎過程的數值分析方面有很強的優勢和能力并有諸多的研究成果［8－10］，但是用SPH法像齒輪傳動這樣有正確嚙合要求接觸體的模型建立和進行實時動態數值分析的算法等方面的研究很少。

本研究首先用SPH法對一對輪齒正確嚙合沖擊問題進行研究，在保證核函數一致性的前提下核函數影響半徑的選擇進行改進，提出了適用于兩個獨立連續體接觸和相互影響的SPH法，從而解決了一對齒輪正確嚙合用SPH法計算和模擬問題。然后，用在不同轉速條件下對一對輪齒正確嚙合沖擊過程進行了實時動態數值仿真，分析了齒廓面上動應力的實時動態分布及變化過程。研究結果表明，SPH法在研究齒輪正確嚙合和沖擊問題的實時動態分析方面有著較好的潛力。

1 一對齒輪正確嚙合和輪齒相互耦合接觸的SPH描述方法

1.1 SPH基本描述

在SPH法中，首先在插值理論的基礎上，采用核近似方法將函數及其梯度用核函數的積分表示，把偏微分形式的連續體的方程轉化為積分形式的方程，然后用粒子近似方法將連續形式的積分方程轉化成離散形式的方程。

對于任一個函數f（x），在一點x的值可近似為

用分部積分、高斯定理和核函數性質，梯度?f（x）的核近似為

函數f（x）及其梯度?f（x）積分形式表達式（1）和（2），用SPH法的粒子離散思路進行離散，有

圖1 SPH粒子核近似示意圖Fig．1 SPH particle kernel approximation

1.2 基于彈性動力學理論的齒輪嚙合沖擊方程的的SPH離散

在連續介質力學和線彈性動力學理論的基礎上齒輪嚙合沖擊的動力學方程和本構方程可寫為：

式中：v為速度矢量；ρ為密度；σ為應力張量；f為外力；ε為應變張量；E為彈性模量；υ為泊松比，當α＝β時，δαβ＝1；當α≠β時，δαβ＝0；希臘字母上標α和β表示坐標方向。

用式（3）和式（4）并考慮人工粘度項Пij項，可得如式（9）～式（11）所示的連續方程、運動方程和幾何方程的SPH離散方程。

1.3 齒輪正確嚙合及相互動態耦合接觸的SPH算法

在用SPH法描述齒輪正確嚙合和沖擊過程的動態接觸行為時，兩個齒輪的相互作用可通過核估算來進行，因此原則上不需要定義特殊的接觸特性。但是，由于常規的SPH法可用一種影響半徑來保證了物體的連續性，而多體接觸問題中同時存在連續體和離散體問題，即各個接觸體自身內為連續體、接觸體之間為離散體。如果按常規的方法進行計算，將一對齒輪之間的關系視為連續體，這就帶來過早的相互影響、很難保證正確接觸問題。因此在SPH法中影響半徑的選擇直接影響正確接觸和計算精度。當然對于連續體來說，如果影響半徑選的太小，影響區域中沒有足夠的粒子，導致連續性降低、很難保證計算精度；如果影響半徑選的太大，保證了物體的連續性，但在計算區域無法反映變形的局域性，同樣也會影響計算精度。因此，通過調整核函數影響半徑可得到連續體和離散體的特征。對于兩個齒輪接觸面附近不同齒輪中，粒子間的相互作用應在兩個齒輪進入接觸時才可以參與計算。如果所有粒子選擇同樣影響半徑，兩個接觸面還未正式接觸就在兩齒輪接觸區域表面粒子間就產生相互作用。

本研究為了解決這個問題，對每一個粒子定義不同的影響半徑來處理兩個連續體之間的SPH接觸計算問題并用式（12）來定義影響半徑。式中：比較大的核半徑khii（i＝1，2）來計算和評價第i齒輪自身內粒子間的相互作用，一般核函數光滑長度h等于粒子間距的1．2倍左右，k是計算中有所使用的核函數來確定，比較小的核半徑khij（i，j＝1，2 i≠j）來計算第i和第j輪齒接觸面上粒子間的相互作用，一般選擇光滑長度h使得影響半徑等于粒子間距一倍左右即可（見圖2）。用修正SPH法［12］彌補影響半徑的選擇帶來的誤差。這樣可以有效防止沒有接觸的部位產生相互干涉問題并提高計算精度。

圖2 影響半徑的選擇Fig．2 Selection of influence radius

2 兩個齒輪的SPH離散粒子模型的建立和其他條件的設定

2.1 兩齒輪的SPH離散粒子模型的建立

齒輪離散SPH粒子模型的建立是齒輪離散成球形“粒子”確定每個粒子中心的坐標。本研究通過分析漸開線直齒輪每個部分的曲線，確定齒廓曲線和齒根過度曲線的方程，編程實現了齒輪的三維SPH離散粒子建模程序。將模型數據用免費公開可視化軟件RASMOL來顯示處理（見圖3）。

圖3 齒輪離散粒子模型Fig．3 Discretemodel of gear contact

實際齒輪傳動中，由于制造、安裝等誤差的原因兩個相互嚙合齒輪輪齒之間必然存在一定的間隙。這間隙在高速旋轉條件下的齒輪機構帶來一定的沖擊，產生噪音，嚴重時會產生輪齒的折斷、表面損傷等。因此，分析一對齒輪沖擊接觸所帶來的動態力學行為，具有一定的意義。由于一對齒輪整體SPH離散后，整體的扭矩載荷（包括阻力矩等）和定軸轉動等相關的邊界條件的計入難度很大，此外計算工作量也非常多。因此，本研究為了節省建模和計算工作量，只考慮了一對齒輪進入嚙合時的一對輪齒的沖擊問題（見圖4），并假設主動輪與從動輪處于完全接觸狀態。此外實際齒輪傳動中，從動輪受較大的工作載荷并扭矩也很大，因此從動輪輪齒可設為瞬間固定并兩輪齒正好在節圓處接觸嚙合。

2.2 輪齒參數和材料屬性

輪齒參數：模數m＝40 mm，齒數z1＝21，z2＝42，壓力角a＝20°。取粒子之間的距離為1 mm，模型離散成69 205個粒子，其中主動輪齒含33 765個粒子，從動輪齒含35 440個粒子

材料屬性：彈性模量E＝206 GPa，泊松比ν＝0．3，密度ρ＝7 870 kg／m3。

2.3 初始和邊界條件

剛開始計算時主動輪齒回轉中心C為定軸以初始角速度ω轉動，被動輪齒保持靜止（見圖4），用下式

表示：

整個計算過程從動輪齒的齒根上粒子保持不動，用下式表示：

圖4 輪齒嚙合沖擊的SPH離散模型Fig．4 Discretemodel of the gear teeth meshing impact

3 數值仿真和結果分析

本文主動輪齒分別以ω＝60 r／min、ω＝120 r／min、ω＝240 r／min三種初始角速度轉動嚙合沖擊從動輪進行了計算和分析。

圖5 最大動應力隨時間的變化規律t＝0至20μsFig．5 Maximum Stress change with time From t＝0 to 20μs

圖5（a）和（b）分別表示兩輪齒沖擊嚙合開始t＝0μs至t＝20μs之間x軸方向的最大動正應力σx，即σx＝max｛σix｝，i＝1…N（N是總粒子數）；最大動米塞斯等效［13］（Mises Equivalent Stress）應力σe，即σe＝max｛σie｝，i＝1…N（N是總粒子數）；隨著時間的變化規律。從圖可知，剛嚙合時，最大動應力隨著嚙合沖擊時間的延長而增大。當沖擊嚙合時間2μs之前應力波動比較大，之后應力逐漸趨于平穩增加狀態。動應力隨著沖擊速度的增大而變大。

圖6是沖擊后不同時刻動等效應力σie分布云圖，這些不同時刻的動應力云圖顯示嚙合輪齒中彈性波傳播過程。可以看出，動等效應力從初始的嚙合區域向齒根方向和向齒廓法向方向擴展，最后傳播到輪齒整體。

圖6 等效動應力sie的傳播過程Fig．6 Propagation of elastic waves represented by equivalent stress distribution

圖7是沖擊后t＝10μs時刻動應力σx和動等效應力σe云圖的比較。圖顯示，不論是動應力σix還是等效動應力σie隨著沖擊速度的增加而增大，主動輪初時沖擊速度的增加，將引起輪齒動應力傳播的提前。

圖7 動應力云圖的比較t＝10μsFig．7 Comparison of stress at t＝10μs

圖8是t＝20μs時刻主動輪齒廓面正中間線上沿齒廓法向方向應力σix和等效應力σie的計算結果比較。結果顯示，最大動應力σix和最大等效動應力σie發生在節圓附近而且沖擊速度越大動應力值也越大。

4 結 論

圖8 主動輪齒廓面正中間曲線上動應力的分布t＝10μsFig．8 Stress on themiddle line along the contact faces of gear's tooth at t＝10μs

（1）本研究提出了選擇不同影響半徑的SPH算法，并通過考慮齒輪各種幾何參數和物理參數實現了齒輪SPH建模方法和數值模擬，解決了基于SPH的一對輪齒正確接觸嚙合的數值分析問題。

（2）通過一對齒輪傳動的一對輪齒進行簡單的SPH離散建模并給定不同工況和邊界條件進行了SPH數值分析，所得到的結果比較接近實際嚙合規律，驗證了本方法的等效性。這將進一步分析齒輪和軸承等復雜接觸界面問題的實時動態數值分析提供科學依據，同時對復雜接觸界面系統的數值建模和分析提供新的算法。

（3）通過一對輪齒的數值模擬，詳細分析了輪齒嚙合沖擊過程中接觸動應力的變化和傳播過程。結果表明，當從動輪瞬間固定時，主動輪轉速越高嚙合動應力也越高并隨時間的延長都趨于平穩增加的趨勢。嚙合動應力以應力波的形式從接觸處擴展到輪齒各部位。

由于目前考慮齒輪整體，進行SPH計算的難度大，因此本研究考慮到從動輪的扭矩設為很大的極端過載工況并只對一對輪齒進行瞬態狀態的SPH數值建模分析。然而，這不能完全接近一對齒輪傳動的真實嚙合過程，這部分將考慮到后續的研究之中并將與有限元等傳統方法進行對比分析，更進一步驗證本方法的等效性。

此外，在建立嚙合沖擊分析模型時，為了考慮數值計算的收斂性，只考慮了沖擊傳動比較平穩的兩輪齒節圓嚙合處開始嚙合的模型。但是，實際齒輪嚙合過程非常復雜，嚙合位置不斷變化，沖擊力也隨著嚙合位置的不同而發生變化，這部分也將后續的研究中進一步考慮計入。

［1］陳霞，夏巨諶，胡國安，等．直齒圓錐齒輪嚙合過程數值模擬［J］．機械設計，2006，23（4）：21－23．

CHEN Xia，XIA Ju-chen，HU Guo-an，etal．Numerical value simulation on matching course of spur bevel gears［J］．Journal of Machine Design，2006；23（4）：21－23．

［2］Hwang SC，Lee JH，etal．Contact stress analysis for a pair of mating gears［J］．Mathematical and Computer Modeling，2013，57：40－49．

［3］Vijayakar SM，Houser D R．Contact analysis of gear using a combined finite element and structure integral method［R］．AGMA Technical Paper，No．16，1991．

［4］Chen JS，Litvin F L，Shabana A A．Computerized simulation of meshing and contact of loaded gear drives［J］．Proc．Inter．Gearing'94 1994：161－166．

［5］唐進元，劉欣，戴進．基于ANSYS／LS-DYNA的齒輪傳動線外嚙合沖擊研究［J］．振動與沖擊，2007，26（9）：40－44．

TANG Jin-yuan，LIU Xin，DAI Jin．A study based on ANSYS／LS-DYNA gear transmission meshing impact［J］．Journal of Vibration and Impact，2007，26（9）：40－44．

［6］Monaghan JJ．Smoothed particle hydrodynamics［J］．Reports on Progress in Physics，2005，68：1703－1759．

［7］丁樺，龍麗平，伍彥峰．SPH方法在模擬線彈性波傳播中的運用［J］．計算力學學報，2005，22（3）：320－325．

DING Hua，LONG Li-ping，WU Yan-feng．Application of smoothed particle hydrodynamicsmethod to the simulation of elastic wave propagation in solid［J］．Chinese Journal of Computational Mechanics，2005，22（3）：320－325．

［8］Shintate K，Sekine H．Numerical simulation of hypervelocity impacts of a projectile on laminated composite plate targets by means of improved SPH method［J］．Composites：Part A，2004，35：683－692．

［9］何建，唐平，王善，等．柱形桿高速碰撞薄靶板的數值仿真［J］．系統仿真學報，2005，17（9）：2107－2110．

HE Jian，TANG Ping，WANG Shan，et al．Numerical simulation fo columnar bar high velocity impact on thin plate［J］．Journal of System Simulation，2005，17（9）：2107－2110．

［10］Liu G R，Liu MB．Smoothed particle hydrodynamics-a meshless particlemethod［M］．Singapore，World Scientific Publishing，2003．

［11］Watkins S J，Bhattal A S，Francis N，et al．A new prescription for viscosity in smoothed particle hydrodynamics［J］．Astronomy and Astrophysics Supplement Series，1996，119（1）：177－187．

［12］Chen JK，Berauni J E，Carney T C．A corrective smoothed particle method for boundary value problems in heat conduction［J］．International Journal for NumericalMethod in Engineering，1999，46：231－252．

［13］薜守義．彈塑性力學［M］．北京：中國建材工業出版社，2006．

SPH algorithm for proper meshing and coup ling contact of gears

Rahmatjan Imin1，2，Mamtimin Geni1
（1．Schools of Mathematics and System Science，Xinjiang University，Urumqi830046，China；2．School of Mechanical Engineering，Xinjiang University，Urumqi830047，China）

Based on the kinetic equations of a pair of gears under contactmeshing and corresponding impact，a smoothed particle hydrodynamics（SPH）method was developed to simulate the propermeshing and coupling contact of gears by adjusting the influence radius of kernel function．SPH pre-processor program was compiled to build a model of correctmeshing of gears，and the dynamic numerical simulation was carried out with considering the meshing impact process of gears in different conditions．The maximum dynamic stress changing with time under the impact between two meshing gears and the distribution and changing process of dynamic stresses on tooth-profile surface were analyzed．The dynamics stress propagation in the meshing impact process，the dynamic contact stresses distribution and its changing process in the vicinty of contact areawere discussed．The paper provides the an effective new SPH numerical algorithm to dynamic contact strength evaluation and optimization design of gear transmission system．

gearmeshing；SPH method；meshing impact；dynamic stress

TH212；TH213．3

A

10．13465／j．cnki．jvs．2015．12．012

國家自然科學基金（51075346）；國家重點基礎研究發展計劃973項目（2011CB706600）

2014－04－04 修改稿收到日期：2014－11－24

熱合買提江·依明 男，博士，講師，1974年生

買買提明·艾尼 男，博士，教授，博士生導師，1958年生