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

（1．新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，烏魯木齊 830046；2．新疆大學(xué)機械工程學(xué)院，烏魯木齊 830047）

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

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

（1．新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，烏魯木齊 830046；2．新疆大學(xué)機械工程學(xué)院，烏魯木齊 830047）

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

齒輪嚙合；SPH計算方法；嚙合沖擊；動應(yīng)力

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

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

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

1.1 SPH基本描述

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

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

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

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

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

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

在連續(xù)介質(zhì)力學(xué)和線彈性動力學(xué)理論的基礎(chǔ)上齒輪嚙合沖擊的動力學(xué)方程和本構(gòu)方程可寫為：

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

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

1.3 齒輪正確嚙合及相互動態(tài)耦合接觸的SPH算法

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

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

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

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

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

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

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

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

2.2 輪齒參數(shù)和材料屬性

輪齒參數(shù)：模數(shù)m＝40 mm，齒數(shù)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 初始和邊界條件

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

表示：

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

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

3 數(shù)值仿真和結(jié)果分析

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

圖5 最大動應(yīng)力隨時間的變化規(guī)律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軸方向的最大動正應(yīng)力σx，即σx＝max｛σix｝，i＝1…N（N是總粒子數(shù)）；最大動米塞斯等效［13］（Mises Equivalent Stress）應(yīng)力σe，即σe＝max｛σie｝，i＝1…N（N是總粒子數(shù)）；隨著時間的變化規(guī)律。從圖可知，剛嚙合時，最大動應(yīng)力隨著嚙合沖擊時間的延長而增大。當(dāng)沖擊嚙合時間2μs之前應(yīng)力波動比較大，之后應(yīng)力逐漸趨于平穩(wěn)增加狀態(tài)。動應(yīng)力隨著沖擊速度的增大而變大。

圖6是沖擊后不同時刻動等效應(yīng)力σie分布云圖，這些不同時刻的動應(yīng)力云圖顯示嚙合輪齒中彈性波傳播過程?？梢钥闯觯瑒拥刃?yīng)力從初始的嚙合區(qū)域向齒根方向和向齒廓法向方向擴展，最后傳播到輪齒整體。

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

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

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

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

4 結(jié) 論

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

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

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

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

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

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

［1］陳霞，夏巨諶，胡國安，等．直齒圓錐齒輪嚙合過程數(shù)值模擬［J］．機械設(shè)計，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］．計算力學(xué)學(xué)報，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］何建，唐平，王善，等．柱形桿高速碰撞薄靶板的數(shù)值仿真［J］．系統(tǒng)仿真學(xué)報，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，F(xiàn)rancis 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］薜守義．彈塑性力學(xué)［M］．北京：中國建材工業(yè)出版社，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

國家自然科學(xué)基金（51075346）；國家重點基礎(chǔ)研究發(fā)展計劃973項目（2011CB706600）

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

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

買買提明·艾尼 男，博士，教授，博士生導(dǎo)師，1958年生