基于彈性圓柱殼的柔性齒輪副建模及動力學分析

張 洋，關先磊，王青山

（中南大學 高性能復雜制造國家重點實驗室，長沙 410083）

齒輪傳動系統廣泛應用于汽車、飛機等各種場合，被認為是最重要的機械部件之一。在工作過程中，傳動齒輪存在橫向與面內的彈性振動，由此會引起結構失效等動力學現象。這種結構失效現象在高速輕質重載齒輪中經常出現，例如高速動車牽引齒輪傳動系統以及航空發動機齒輪系統。因此，建立考慮輪體柔性的齒輪傳動動力學模型具有十分重要的意義。

為了獲得高效齒輪，人們對齒輪系統的動態建模和振動特性分析進行了大量的研究。對于齒輪的輪體結構，通常有三種建模方法。其一是扭振模型即剛性模型，Han和Chu[1]研究了齒輪傳動轉子系統周期角角速度運動下的動力學特性。Chen等[2]研究了雙螺旋齒輪的旋轉動力學，考慮到了軸承、回轉效果，其中采用鐵木辛柯梁作為軸單元，齒輪則視為剛性體。He 等[3]提出了一種新穎的剛性齒輪嚙合模型，并采用改進的勢能法計算了齒輪偏心對時效嚙合剛度的影響。Shi等[4]完成了考慮時變參數和側隙的嚙合圓柱齒輪對的建模和分析。其二是以彈性環作為齒輪的輪體模型，Christopher和Robert[5]研究以耦合彈性環為模型，分析研究高速柔性齒輪副的振動。Chen等[6]研究了帶撓性正齒環行星齒輪組的動態仿真。Fan 等[7]根據殼理論和鐵木辛柯梁理論建立了行星齒輪內環和轉軸模型。其三是以柔性圓盤作為齒輪輪體模型。Vinayak 和Singh[8]建立了柔性齒輪的多體動力學模型，考慮了分布齒輪嚙合剛度和齒胚模型，求解了一系列的模態問題和響應問題。

綜合比較上述幾種模型可知，剛性模型具有模型簡單，計算量小的優點，但是剛性模型沒有考慮到齒輪齒面變形，和實際齒輪貼合度存在差距，尤其是齒輪輪體厚度較薄的齒輪，剛性模型和柔性模型的差距較大。相比較于剛性模型，柔性輪體模型則更加符合工程實際，因此柔性模型的建模和分析是非常有必要的?；趶椥原h的柔性齒輪副模型只適用于齒輪的齒寬較小的情況，而以柔性圓盤為齒輪輪體的齒輪模型雖然適用范圍較廣，但是存在模型收斂性不好的缺點?；诂F有柔性齒輪副模型研究的不足，本文基于彈性圓柱殼單元，建立一種新穎的柔性齒輪輪體模型，相比于上述兩種柔性模型，基于彈性圓柱殼的齒輪副模型，具有良好的收斂性，可適用于齒輪寬度較大的情況，還可以進行齒輪寬度方向的動力學特性研究。

綜上所述，本文首先基于微分求積有限元法和旋轉圓柱殼理論，推導旋轉圓柱殼的動能和勢能，完成柔性齒輪的輪體建模。兩圓柱齒輪由空間固定的具有時變嚙合剛度的彈簧耦合，表示齒輪之間的嚙合效應。并且每個圓柱殼用彈性基礎支承，來表示軸以及軸承的彈性。最后分別使用拉格朗日方程和Newmark 迭代來求解柔性齒輪副的模態和動態響應。本文以某航空航天齒輪副為例，通過收斂性分析，研究了本文模型的數值穩定性；然后與ANSYS的有限元模型進行模態特性的對比，驗證本文模型的正確性。接下來，數值計算了大轉速范圍下柔性齒輪副系統固有頻率隨轉速的變化規律。最后通過與傳統的基于齒輪體剛性假設的扭振模型進行對比，闡述齒輪柔性對系統動力學特性帶來的影響。

1 理論推導

1.1 柔性齒輪副的幾何模型

柔性齒輪副的模型如圖1所示。齒輪簡化為圓柱殼模型，齒輪嚙合副簡化為彈簧模型。分別在兩個圓柱殼的中心層上定義正交圓柱坐標系來描述殼體的振動，坐標軸xi，θi，zi(i=1，2)分別沿著軸向，周向以及徑向。以該坐標系為參考，殼體的位移用ui，vi，wi(i=1，2)來表示。齒輪的轉速為Ωi，根據齒輪嚙合原理可知Ω2=-R1Ω1/R2，齒輪的基圓半徑為Ri，殼厚為hi，齒寬為li。假設齒輪的密度、泊松比、彈性模量以及剪切模量分別用ρi，μi，Ei，Gi表示。每個齒輪分別采用徑向彈簧（kri）和切向彈簧（kθi）來代表軸的彈性，嚙合彈簧的剛度為km。

圖1 柔性齒輪副模型

1.2 齒輪副系統的能量分析

根據1 階剪切變形理論，殼體的位移分量可由中間面的位移來表示[9]：

殼體上任意點的位置矢量和速度矢量分別為：

其中：i、j、k分別表示x、θ、z方向上的單位矢量。

旋轉圓柱殼的動能為：

將式（2）代入式（3）可得：

圓柱殼的彈性勢能公式為：

其中應力應變為：

其中：

E為彈性模量，μ為泊松比。

兩齒輪對應的嚙合點之間用彈簧連接，其彈簧的嚙合勢能為：

其中：v1和v2分別表示主、從動輪上嚙合點位置的切向位移，km(t)為時變嚙合彈簧的剛度。

在圓柱殼中間環內部，分別采用徑向彈簧和切向彈簧模擬軸的彈性。支撐彈簧勢能為：

kri和kθi分別表示徑向和切向單位長度上的彈簧剛度。但是在有限元方法中，彈簧只能添加在節點上，因此，采用一下等效方式，得到節點上的等效彈簧剛度：

其中：C表示圓環周長，N表示有限元網格的節點數，keri表示施加在節點上的等效彈簧剛度。

因此系統的總應變能為：

1.3 系統的離散及求解

將系統的的總應變能和總動能代入拉格朗日方程可得齒輪嚙合系統的動力學方程為：

即可得到圓柱殼齒輪嚙合結構的整體質量矩陣M、剛度矩陣K、阻尼矩陣C分別為：

其中：

其中：Mk是對稱矩陣。A和B為微分權系數矩陣；C為積分權系數矩陣[10]，E為單位矩陣。

最后使用Newmark 積分迭代方法進行動態響應求解[11]。Newmark積分迭代方法的算法流程圖如圖2所示。

2 算例分析

本文采用實際工程中的齒輪參數作為算例對象，主從動齒輪的尺寸參數如表1所示。材料參數為：彈性模量E=2.1×1011Pa；密度ρ=7 850 kg/m3；泊松比μ=0.3。需要特別說明支撐彈簧設置和載荷設置。在研究系統模態時，同時考慮徑向和周向的支撐彈簧，模擬轉軸和載荷作用下的工況。而在進行動力響應計算時，只在主動輪上施加周向載荷（即扭矩），同時去掉周向支撐彈簧，保留徑向支撐彈簧；而從動輪上則不施加載荷，繼續保留周向和徑向支撐彈簧。

表1 齒輪副尺寸參數

根據參考文獻得到的齒輪嚙合剛度的解析模型[12]，可得到如圖3所示的嚙合剛度曲線，取算術平均值可得平均嚙合剛度1.350 3×109N/m，這將作為后續進行模態特性分析的嚙合剛度。齒輪的嚙合剛度平均分配到對應的有限元嚙合節點上。

圖3 時變嚙合剛度曲線

2.1 收斂性和模型驗證

首先進行模態收斂性分析，通過計算主從動齒輪嚙合模型的模態，來驗證微分求積有限元法的收斂性。分別沿兩圓柱殼的軸向和周向劃分單元數為M和N，并且每個單元的微分求積節點數為n。由于計算量過大，先確定軸向單元數M=2；研究周向單元數和節點數對收斂性的影響。圖中偏差計算公式為(fi-fend)/fend%，其中fi為某單元數N和單元節點數n時的頻率，fend為N和n同時取最大值時得到的頻率，由圖4易得，模型的收斂性極好，根據收斂性結果，在后續的分析中選擇單個圓柱殼軸向單元數量為2，周向單元數量為2×6，單元節點為6進行計算。

圖4 柔性嚙合副模型的前4階頻率收斂性

為進一步驗證模型結果的正確性，接下來進行模態對比驗證。在ANSYS中建立有限元模型，圓柱殼采用Shell單元，每個圓柱殼40×80個單元；嚙合彈簧和支撐彈簧使用COMBIN14 單元有限元模型如圖5所示。

圖5 ANSYS模型

表2展示了本模型和ANSYS 模型的前12 階自然頻率對比，可以發現，兩者結果吻合很好，圖6展示了齒輪嚙合模型的前4 階模態振型，其結果也與ANSYS 結果吻合。上述結論充分驗證了模型的正確性。

圖6 柔性齒輪副前4階振型結果對比

表2 兩種模型自然頻率對比

2.2 坎貝爾圖分析

在進行動態響應分析之前，首先計算系統的自然頻率并得坎貝爾圖。在嚙合剛度為1.350 3×107N/m的情況下，計算不同轉速時系統的自然頻率，結果如圖7所示。在零速時，由于嚙合剛度的存在破壞了圓柱殼的軸對稱性，導致固有頻率不同。

圖7 柔性齒輪副坎貝爾圖

當速度從零開始增加時，一些固有頻率增加而另一些固有頻率減少。當自然頻率減小到0 時，出現臨界轉速。由于此模型的尺寸較小，導致在轉速增加到3 000 rad/s時，還未出現第一臨界轉速。

2.3 時間步長收斂性研究

從1.3的模型求解部分，可以知道進行動態響應計算時，需要選擇合適的時間步長，時間步長Δt的選擇會影響到動態響應結果的收斂性。時間步長的選取依據為：所關注頻段內的結構主要貢獻的若干振型的最小周期（最大頻率）的1/10～1/20[11]。設取第20 階頻率為最大頻率，取其倒數的1/c2的值作為一種時間步。對應的時間步長區間為：Δt=4.518×10-6～9.036×10-6s；

同時，還需考慮齒輪嚙合的頻率，不同轉速對應的時間步長為：

需要注意的是，取兩者時間步長的較小值?，F在轉速為100 rad/s 的工況下，選取主動輪上的中間嚙合點作為研究對象，c1取200，研究c2取值對收斂性的影響。表3為c2的取值對應的時間步長。

表3 時間步長取值

圖8是時間步長收斂性結果示意圖，其中局部放大圖的橫縱坐標單位與整體圖的單位一致。從圖中可以發現，在合理區間內，不同時間步長對模型的動態響應曲線影響很小，所得響應曲線幾乎完全重合。說明在時間步長足夠小的情況下，時間步長的選擇已經不影響動態響應分析結果，從而驗證了本模型在計算動態響應時時間步的收斂性很好。

圖8 動態響應時間步長收斂性

2.4 剛柔模型動態特性對比

接下來，通過和傳統的剛性扭振模型[13]進行動態特性結果對比，通過兩種模型的差異性，闡述輪體柔性對齒輪副系統的影響。兩種模型處于相同的工況之下，即載荷以及支承彈簧，時變嚙合剛度，模型尺寸都相同。圖9為不同轉速下剛柔模型的穩定之后的動態響應曲線對比，其中柔性模型的響應曲線為主動輪上中間嚙合點處的響應曲線，剛性模型的響應曲線為主動輪基圓半徑上的位移響應曲線。

由圖可知，兩種模型的動態響應周期完全一致，這是由于兩者的時變嚙合剛度一致。剛性模型的動態響應曲線波動范圍明顯大于柔性模型，振動平衡位置也大于柔性模型。這是由于柔性模型中圓柱殼模型的柔性引起。在剛性模型中，僅考慮了嚙合彈簧和支撐彈簧的柔性，因此在相同載荷下，外力做功全部轉化成彈簧的勢能；而在柔性模型中，外力做功一部分轉化成為彈簧的勢能，還有一部分轉化為柔性圓柱殼的應變能。換而言之，圓柱殼的柔性對于齒輪副系統有減振吸能效果，因此柔性模型中嚙合彈簧的勢能小于剛性模型中嚙合彈簧的勢能，從而表現為剛性結果的平衡位置和振動波動范圍都大于柔性模型。

圖10展示了不同轉速下剛柔模型的頻譜成分，圖10(a)為剛性模型，可以發現圖中頻率成分以一倍嚙合頻率和多倍嚙合頻率成分為主，當一倍頻率和多倍頻率的幅值線與一條平行于轉速軸的直線相交時，頻譜圖的幅值將會達到峰值，該條直線對應的頻率為該模型下的第1階固有頻率f1（扭轉模型的第1、2 階固有頻率分別為1 156.5 Hz 和15 368 Hz），這種現象叫做參數共振，這是由于嚙合頻率是該模型的激振頻率，當激勵頻率接近固有頻率或接近固有頻率的兩倍時，就會發生這種共振。由于頻譜圖的關注范圍為0～10 000 Hz，所以只出現了一個共振區域。除去參數共振區域附近，一倍頻和高倍頻對應的峰值隨轉速增加，總體上也呈現線性增加的趨勢，且一倍頻率的峰值明顯大于多倍頻率的峰值。

圖10(b)為柔性模型，可以發現圖中頻率成分仍然以一倍嚙合頻率和高倍嚙合頻率成分為主，但是峰值隨轉速的變化波動很大。和剛性模型類似，柔性模型的頻譜圖的幅值也會出現參數共振現象，但是由于柔性模型使用的圓柱殼單元，該模型的第一階固有頻率為1 230.5 Hz。除去峰值附近區域，一倍頻和高倍頻對應的峰值隨轉速增加呈現波動上升的趨勢，但是在高轉速下，一倍頻率的峰值并不占有明顯優勢。

圖10 剛柔模型的頻譜圖

從兩種模型頻譜圖的整體幅值量級上看，剛性模型的頻譜幅值大于柔性模型的頻譜幅值。圖9的動態響應曲線也可以很好說明這個現象，在轉速為10 rad/s和500 rad/s的轉速下，剛性的振動范圍都遠大于柔性的波動范圍，因此將時域曲線轉換為頻譜圖時，就會出現剛性頻譜幅值大于柔性頻譜幅值。

3 結語

(1)本文基于圓柱殼理論，考慮時變嚙合剛度以及彈性支撐的影響，建立了柔性齒輪副的等效簡化模型。該模型具有適用范圍廣以及收斂性好的優點，豐富了現有的理論知識。

(2)基于該模型，開展了柔性齒輪副動力學特性驗證和研究。推導了基于微分求積有限單元法的圓柱殼齒輪嚙合模型的線性微分方程組。在確定軸向單元數的前提下，分別研究不同周向單元數以及單元節點數下模型的自然頻率，討論了模型的數值收斂性。接下來，與ANSYS的有限元模型對比自然頻率，來驗證本文模型的正確性。

(3)建立了和圓柱殼齒輪嚙合副模型等效的集中質量剛性模型。通過兩種模型的動態響應曲線、位移響應頻譜圖的對比，來研究齒輪柔性的影響。

(4)在相同的工況下，柔性模型和剛性模型的動態響應曲線周期相同，并且兩種模型在同一轉速下的響應曲線形狀相似，但是柔性模型結果的平衡位置和振動波動范圍都小于剛性模型。這是由于柔性模型中圓柱殼模型的柔性引起。在剛性模型中，僅考慮了嚙合彈簧和支撐彈簧的柔性，因此在相同載荷下，外力做功全部轉化成彈簧的勢能；而在柔性模型中，外力做功一部分轉化成為彈簧的勢能，還有一部分轉化為柔性圓柱殼的應變能。

(5)從位移響應頻譜圖可以看出，兩種模型的頻率成分都以一倍嚙合頻率和高倍嚙合頻率成分為主。但是由于圓柱殼模型的柔性，導致了其固有頻率與剛性模型不一致，從而導致了參數共振區域的不一致，進而影響了一倍嚙合頻率和高倍嚙合頻率的幅值。并且由于齒輪柔性的影響，柔性模型頻譜圖的幅值整體上小于剛性模型的幅值。