荀 超, 吳海軒, 王云龍
(南京工程學院 機械工程學院,南京 211167)
齒輪廣泛應用于各類機械傳動中,在多種內外激勵的共同作用下,復雜齒輪機構的振動與噪聲問題突出,引起的齒輪失效風險嚴重制約了其應用。這就需要在設計階段,對復雜齒輪機構的動力學行為進行準確預測。齒輪傳動機構是具有參數激勵和“非光滑”特性的非線性系統。時變的嚙合剛度是齒輪傳動機構振動的主要內部激勵之一[1]。齒側間隙導致在振動劇烈時,輪齒嚙合會出現脫齒。脫齒使齒輪機構的動態特性具有明顯的間隙非線性特征。復雜齒輪傳動機構存在多處輪齒嚙合,機構內部有多處參數激勵與脫齒,其間復雜的耦合關系增加了對齒輪機構動態特性分析的難度。二級齒輪傳動機構是復雜齒輪傳動系統主要的基本組成形式,針對二級齒輪傳動機構開展非線性動態特性分析,可為復雜齒輪機構的傳動性能研究提供理論方法指導。
數值積分方法需對時間歷程進行分段,各時間段內的嚙合剛度按照定值處理。從而將一個嚙合周期內的時變嚙合剛度曲線,切分成分段定值曲線。通過依次求解各分段內的響應,從而獲取機構在整個時間歷程上的振動響應[2]。通過數值積分求解,能夠得到齒輪接觸特性與系統振動之間的耦合關系[3-4]、齒輪機構的非線性振動特性[5]和失穩特性[6]等。數值積分能夠提供較準確的計算結果,但其計算耗時較長,且無法明確解釋各種動力學現象與系統參數之間的定量關系。
諧波平衡法和多尺度法是常用的齒輪非線性動態特性解析方法。Kahraman等[7-9]使用諧波平衡法分析了軸-軸承-齒輪系統的非線性動態特性,得到了系統響應的幅頻曲線。Al-Shyyab等[10]借助諧波平衡法,對多級齒輪系統展開了周期1和次諧波振動分析。通過實例分析發現,脫齒引起的剛度“軟化”現象對多級齒輪系統動態特性具有顯著影響。在脫齒發生的頻率范圍內,系統會失穩甚至出現混沌。此后,有學者將諧波平衡法應用于分析行星齒輪系統[11-13]、非圓齒輪[14]的非線性動態特性。然而傳統的諧波平衡法無法處理超諧、次諧響應。Liu等[15-16]采用多尺度法,考慮時變嚙合剛度與脫齒現象,分析了多級齒輪系統和行星齒輪系統的非線性動力學特性。多尺度法依賴于攝動參數為小變量,假定齒輪系統的脫齒時間與嚙合周期相比(脫齒率)為小量。當齒輪出現輕微脫齒,多尺度法具有較高的計算精度,適用于分析系統的主共振、超諧共振、亞諧共振的響應和穩定性問題等。然而,在齒輪脫齒率較高時,多尺度法的計算精度明顯下降。當脫齒率達到30%時,多尺度法與數值積分之間的計算誤差將超過10%[17]。數值積分的計算結果表明,在主共振區域內,齒輪系統的脫齒率會達到40%,甚至更高。此時多尺度法的計算精度難以保證,需要尋求一種可適用于齒輪系統作強非線性振動時的理論分析方法。
廖世俊等[18]提出的同倫攝動方法,不同于多尺度方法,其分析過程不依賴于任何小變量,針對弱非線性和強非線性動力學問題均可適用。同倫方法提供了一個可以調整和控制系統近似級數解收斂域的途徑[19],已成功地應用于求解諸多類型的非線性問題中[20-23]。Wen等[24-26]已證實同倫方法可應用于分析單自由度的齒輪系統非線性動力學行為。多級齒輪系統是復雜的多自由度非線性系統,各自由度之間耦合關系密切。對多級齒輪系統非線性動態特性的分析難度遠大于單自由度齒輪系統。
采用同倫法推導二級齒輪傳動機構主共振、亞諧共振和超諧共振的幅頻響應的解析表達式,提高對其各類共振區域內非線性動態響應的精度,揭示時變嚙合剛度、嚙合脫齒等激勵對機構非線性動態特性的影響機理。
一種常見的二級齒輪傳動機構,左側的小齒輪為主動輪,右側為輸出輪,中間為惰輪,如圖1所示。齒輪機構的動態傳遞誤差主要來自齒輪扭轉方向的振動。為明確時變嚙合剛度的波動對齒輪機構扭轉振動和動態傳遞誤差的激勵機理,僅考慮齒輪扭轉方向的自由度,建立純扭轉動力學模型,如圖2所示。其中齒輪的主體部分視為質點,輪齒嚙合借助線性彈簧模擬。惰輪兩側的嚙合剛度分別km1、km2,u1、u2、u3分別代表各齒輪的扭轉位移量。

圖1 二級齒輪傳動機構Fig.1 Structure of multi-mesh gear set

圖2 二級齒輪機構的純扭轉集中質量動力學模型Fig.2 Lumped parameter model of amulti-mesh gear set
該齒輪機構的振動方程為
(1)
式中,Ft為負載扭矩在輸出齒輪分度圓上的等效阻力;M和x分別為機構的質量矩陣和位移向量。
(2)
Kmi分別為無量綱的嚙合剛度矩陣
(3)
時變嚙合剛度kmi(t)的傅里葉級數表達式分別為
(4)
式中,ωm為嚙合頻率。嚙合變形量δi及脫齒函數Θ(δi)分別為
(5)
脫齒函數的傅里葉級數表達式為
(6)
式中,βi為脫齒函數與時變嚙合剛度間的相位差。振動方程(1)的特征值方程為
(7)
式中,Km0為平均嚙合剛度矩陣
(8)
將x=Vz代入式(1),即可得到系統在模態空間的振動方程為
其中,
Cv=VTCV,Gi=VTKmiV,Fvt=VTFt
假定前兩階扭轉固有頻率分別為第k階和第h階,由式(9)可得第k階和第h階模態方程為
(11)
(12)
以表1中的二級齒輪傳動機構為研究對象,以數值積分(numerical integration, NI)計算結果為參考,分別對比同倫法(homotopy analysis method, HAM)和多尺度法(method of multi-scale, MMS)在主共振、亞諧共振和超諧共振附近的計算精度。

表1 示例齒輪機構的主要參數Tab.1 Parameters of an example gear set
通過有限元方法計算示例二級齒輪傳動機構的嚙合剛度,嚙合剛度曲線如圖3所示。當嚙合剛度取其平均值時,表1中齒輪機構的固有頻率分別為0(剛體模態)、4 049.1 Hz、5 744.5 Hz。

圖3 嚙合剛度曲線Fig.3 Mesh stiffness of the gear set
在ωm≈ωk主共振區域內,令Tk和Th分別為zk(t)和zh(t)的響應周期。令ak為zk(t)的最大值,即ak=max[zk(t)]。令

(13)
實際上,當ωμ≈ωκ時,第η階模態響應ζη(τ)為定值。經過式(14)的轉換,可得
τ=ωμτ,ζκ(τ)=?k+akvk(τ),zh(t)=?h
(14)
第k、h階模態方程進一步改寫為


(15)


(16)
第k階模態響應可寫為冪級數形式為
(17)
其中,cl和dl為各級冪級數系數。可令
(18)
由解表達式(17)和式(15)可設輔助線性算子為
(19)
由式(15)可令非線性算子為

(20)
式中:q∈[0,1]為嵌入參數;Vk(τ;q)為τ和q的實函數;Ak(q)、Λk(q)和Λh(q)則為q的實函數。令?k為非零輔助參數,而Hk(τ)表示非零輔助函數。則0階變形方程可構造為
(1-q)Lk[Vk(τ;q)-vk,0(τ)]=q?kHk(τ)Nk
(21)
當q=0時,顯然有
Vk(τ;0)=vk,0(τ),Ak(0)=ak,0,
Λk(0)=?k,0,Λh(0)=?h,0
(22)
當q=1時,則
Vk(τ;1)=vk(τ),Ak(1)=ak,
Λk(1)=?k,Λh(1)=?h
(23)
由此可見,當嵌入參數q由0連續增至1的過程中,Vk(τ;q)由初始猜測解vk,0(τ)逼近精確解vk(τ)。同樣,Ak(q)、Λk(q)和Λh(q)則分別由初始猜測解ak,0,?k,0和?h,0向相應的精確解ak,?k和?h靠近。
0階變形方程式中包含輔助參數?k和輔助函數Hk(τ),假定?k和Hk(τ)都構造合理,使n>1時存在

(24)
由泰勒理論和式(20),能夠以θ的級數形式將各式展開為

(25)
假定?k和Hk(τ)都構造合理,使得式(25)在q=1處收斂。借助式(25),解表達式為

(26)
定義如下幾個向量

(27)
令0階變形方程對q求導n次,并除以n!,最終令q=0,即可得到高階變形方程為

(28)
其中,

(29)
且
(30)
簡便起見,令
Hk(τ)=1
(31)
當n=1時


(32)
Rk,1可簡化為
Rk,1=c1,0+c1,1ej(τ+βk)+d1,1e-j(τ+βk)+
c2,1ej2(τ+βk)+d2,1e-j2(τ+βk)+c3,1ej3(τ+βk)+
d3,1e-j3(τ+βk)
(33)
若c11≠0,d11≠0,則式(33)中包含永年項,因此必須強制令

(34)
將實部與虛部分離,即可得到模態振幅與相位方程組為

(35)
由式(35)化簡可得

(36)
式中,ψ為Ξ2的相位。消除式(36)中的三角函數,即可得到系統在主共振區域內的幅頻關系式為
(37)
由此可求得第k階模態振動幅值為
(38)
由幅頻響應表達式可見,系統主共振響應與阻尼λk和Ξ1,2,3關系密切。Ξ1中體現著平均嚙合剛度與脫齒函數Θ第0階諧波成分的作用;Ξ2則體現出時變嚙合剛度的作用;而Ξ3則體現出平均嚙合剛度、脫齒函數Θ第0階和第2階諧波成分的作用。式(15)和式(16)中,第k階和第h階模態響應的平均值滿足方程

(39)
即可得到?k,0和?h,0為


(40)
當ν=1,由式(6)可得脫齒函數的傅里葉級數表達式為
(41)
第κ階模態振幅脫齒邊界取決于?k,0、?h,0和系統模態振型。由式(40)中?k,0和?h,0的表達式可見,傳動轉矩、平均嚙合剛度、脫齒函數共同決定著模態振幅的脫齒邊界。

(42)
利用式(38)的幅頻響應表達式和式(36)、式(40)、式(41),給定一個模態響應幅值ak,0,即可得到響應的嚙合頻率ωm,進而得到機構在第k階模態主共振區域內的幅頻響應曲線。圖4(a)為同倫法所得從動齒輪幅頻曲線。與多尺度法和數值積分結果對比可見,在嚙合頻率遠離固有頻率、未發生脫齒的頻率范圍中,同倫法和多尺度法所得曲線與數值積分的曲線之間的誤差很小,且同倫法得到了與多尺度法、數值計算方法一致的脫齒幅值邊界。

圖4 第1階模態主共振,HAM計算結果對比Fig.4 Comparisons between the HAM and NI for primary resonance in the first mode
對照圖4(b)、圖4(c)中的脫齒率曲線可見,多尺度法所得幅頻曲線與數值積分結果之間的偏差,隨著脫齒率的增大而增大。將某一嚙合頻率處的模態振幅間的差值定義為不同算法間的誤差。如圖4(a)所示,多尺度法、同倫法所得幅頻響應曲線的最大幅值分別為67.8 μm、65.3 μm,對應嚙合頻率處的數值積分結果分別為56.7 μm和62.3 μm。以數值積分的結果為參照,多尺度法和同倫法對應的相對誤差分別為19.6%和4.8%。由于多尺度法依賴于脫齒率較小的假設,此處主動輪嚙合處的脫齒率超過了40%,導致多尺度法的誤差較大,而同倫法仍能具有較高的計算精度。
圖5(a)中為各方法所得中間齒輪第二階扭轉模態主共振幅頻曲線對比,圖5(b)、圖5(c)展示了相應的脫齒率曲線。與圖4中結果所得結論相同,以數值計算結果為參考,隨著脫齒率的增大,多尺度方法的計算精度降低,而同倫法仍能保持較高的計算精度。

圖5 第2階扭轉模態主共振,HAM結果對比Fig.5 Comparisons between the HAM and NI for primary resonance in the second mode
對于ωm≈2ωk,參照主共振中的同倫推導過程,即可得到亞諧共振的幅頻關系式為

(43)
由幅頻響應表達式可見,亞諧共振響應同樣與阻尼λk和Ξ1,2,3關系密切。Ξ1中體現著平均嚙合剛度與脫齒函數Θ第1階諧波成分的作用;Ξ2則體現出時變嚙合剛度的作用;Ξ3則體現出平均嚙合剛度、脫齒函數Θ第0階和第2階諧波成分的作用。
圖6、圖7分別對比了同倫法、多尺度法和數值積分前兩階扭轉模態亞諧共振的幅頻響應曲線。其中,虛線為不穩定解分支,兩條垂直線標定著系統的不穩定邊界。幅頻響應曲線在出現脫齒前保持垂直。出現脫齒現象后,原本垂直的兩條幅頻響應曲線向左傾斜。與主共振情況相似,較低的分支為不穩定解,而另一分支為穩定解。

圖6 第1階扭轉模態亞諧共振,HAM計算結果對比Fig.6 Comparisons between the HAM and NI for sub-harmonic resonance in the first mode

圖7 第二階扭轉模態亞諧共振,HAM計算結果對比Fig.7 Comparisons between the HAM and NI for sub-harmonic resonance in the second mode
由圖6和圖7可知,以數值計算結果為參照,同倫法比多尺度法可以在更高的脫齒率范圍內保持計算精度,這是因為同倫法不依賴于任何小變量。
考慮時變嚙合剛度的前兩階諧波成分,依照與主共振相似的同倫法推導過程,即可得到超諧共振ωm≈1/2ωk時的幅頻關系式
(44)
其中,

系統亞諧共振響應同樣與阻尼λk和Ξ1,2,3關系密切。Ξ1中體現著平均嚙合剛度與脫齒函數Θ第1階諧波成分的作用;Ξ2則體現出時變嚙合剛度第二階諧波成分的作用;Ξ3則體現出平均嚙合剛度、脫齒函數Θ第0和第2階諧波成分的作用。
圖8中對比了超諧共振ωm≈1/2ωk同倫法與數值積分所得幅頻響應曲線。圖中實線為數值積分所得曲線,圈線為同倫法所得曲線。

圖8 超諧共振,HAM計算結果對比Fig.8 Comparisons between the HAM and NI for super-harmonic resonance
由圖8可知,當嚙合頻率接近前兩階扭轉固有頻率1/2時,幅頻響應曲線的頂部向左輕微傾斜,說明系統發生了輕微的脫齒現象。同倫法與數值積分所得幅頻響應曲線達到了很好的一致,同倫法同樣適用于計算齒輪系統超諧共振響應。
(1) 考慮時變嚙合剛度、脫齒、各級齒輪嚙合間的耦合關系,建立了二級直齒傳動機構的剛柔混合模型,通過同倫法一階變形公式,推導得到了二級齒輪系統主共振、亞諧共振和超諧共振的幅頻關系表達式。
(2) 同倫法可適用于求解含有間隙非線性因素的,多自由度齒輪機構的共振、亞諧共振和超諧共振等非線性響應,能夠準確預測脫齒發生的嚙合頻率、剛度“軟化”等。
(3) 與多尺度法相比,同倫法推導過程不依賴于任何小變量,因而在脫齒率較高時,以數值積分結果為參照,同倫法具有比多尺度法更高的計算精度。本研究為復雜齒輪傳動機構的非線性動態特性分析提供了更為精確的理論分析方法。