接觸應力對飛輪轉子動力特性的影響

任正義，楊立平，曹志好，曹 飛

(1.哈爾濱工程大學工程訓練國家級實驗教學示范中心，黑龍江 哈爾濱150001；2.哈爾濱工程大學機電工程學院，黑龍江 哈爾濱150001)

1 引言

環境惡化和資源枯竭是現如今人類面臨的巨大挑戰，人類在發展的道路上，正在尋找節能環保高效的能源來代替傳統能源。飛輪儲能系統作為新能源領域的代表，越來越受到人們的重視［1］。飛輪儲能系統是將電能轉化為機械能再到電能的一種能量轉化裝置，主要由轉子系統、高速永磁無刷電機、控制箱、支架和冷卻系統等組成［2］。如今，人們對飛輪儲能系統的研究日漸成熟。韓邦成等對工作轉速為5000rpm的飛輪轉子進行了模態分析，并通過試驗驗證有限元計算的結果［3］。陳強等對飛輪儲能系統有限元分析，求解其臨界轉速和模態振型［4］。劉佩等對1MW/60MJ軸系轉軸進行優化分析和模態分析，并未考慮過盈配合處應力的影響［5］。付雄新等提出了將安裝基礎運動轉化為對轉子激勵的動力學求解方法來分析飛輪電池轉子在路面輸入引起安裝基礎振動變化情況下的振幅變化［6］。以上論文都是在無接觸應力的基礎上分析的，實際情況是軸系結構主要部件之間都是過盈裝配，存在接觸應力。因此建立有接觸應力的高精度飛輪儲能系統的動力學模型顯得尤為重要。

2 飛輪儲能系統

2.1 飛輪儲能系統結構組成

該飛輪儲能系統采用立式左右對稱結構布局，主要包括機械部分和控制部分兩大模塊。機械部分的結構示意圖，如圖1所示。其中金屬飛輪轉動慣量0.59kg.m^2，最大額定轉速為16000rpm，儲能量為600Wh。電機采用電動/發電一體的永磁無刷高速直流電機，額定功率為7.4Kw，額定轉速為16000rpm。該系統主要由循環水冷系統進行冷卻。另外，電機的轉子安裝在轉軸上，使結構更加緊湊，但也會使轉子軸系結構過長。推力盤和軸向電磁軸承安裝在轉軸上部分，主要作用是使轉子系統懸浮在空中。上下兩端各有電磁軸承和保護軸承作為支撐，對轉子系統的動力學特性至關重要［7］。

圖1 600Wh飛輪儲能系統結構布局Fig.1 Structure Layout of 600Wh Flywheel Energy Storage System

2.2 飛輪轉自三維模型

飛輪轉子系統是整個儲能系統的關鍵部分，大約95%的儲能量儲存在這里。利用Soildworks將飛輪轉子零件裝配，得到的三維圖，如圖2所示。

圖2 飛輪轉子三維圖Fig.2 Graphic Model of Flywheel Rotor

轉子系統主要部分的材料及屬性，如表1所示。

表1 軸系部件材料屬性Tab.1 Material Properties of Shafting Components

3 接觸應力分析

3.1 靜態過盈量計算

該600Wh飛輪儲能系統設計的額定轉速為16000rpm，不平衡質量產生的離心力也會很大。為了保證軸上過盈配合處零件在高速旋轉的時候與轉軸之間不發生松脫，合適的過盈裝配量顯得尤為重要。軸上過盈配合處的結構簡圖，如圖3所示。a表示軸徑，b表示配合面的半徑，c表示與軸過盈配合處部件外徑，p表示過盈配合面之間的靜壓力，E1表示軸的彈性模量，E2表示配合件的彈性模量，x、y分別為轉軸和配合件的徑向位移，μ1表示軸的泊松比，μ2表示配合件的泊松比，r1表示軸的半徑，r2表示配合件半徑，Δ表示配合面的過盈量。

圖3 過盈配合處的結構簡圖Fig.3 Structural Sketch of Interference Fit

根據彈性力學理論［8］得：

根據上式三個方程可以求解出r1=r2=b處配合面的接觸壓力為：

由于該飛輪轉軸為實心，即a=0時，則配合處的接觸壓力為：

要滿足使用要求，需滿足接觸壓力值小于材料的屈服強度即：

根據關系式（4），以配合處最低屈服強度的材料屬性參數代入，求得各過盈配合處最大過盈量，如表2所示。

表2 各配合處最大靜態過盈量Tab.2 Maximum Static Interference at Each Mating Point

3.2 動態過盈量計算

當轉子以角速度ω旋轉時，由F=mrω2得到，飛輪轂、軸承轉子電機轉子和推力盤的離心力大于轉軸的離心力，會使上述各部分的內徑變形比轉軸的外徑變形大，兩者之間存在一個變形差值。隨著轉速的增大，變形差越大，當差值大于兩者之間的過盈量時，轉軸和各部分之間會發生脫動。轉動時與軸過盈配合處各部件的徑向位移為y1，軸的徑向位移為x1，配合面的過盈量減少值為Δ1，則［8］

動態過盈量為：

式中：ρ—材料密度。

以額定轉速16000rpm為例，由表1查得各個值代入式（5）和式（6），計算出各過盈配合處不發生脫動的最小動態過盈量，如表3所示。

表3 各配合處最小動態過盈量Tab.3 Minimum Static Interference at Each Mating Points

4 轉子系統有限元仿真與分析

4.1 轉子系統有限元模型的建立

過盈配合是一種高度的非線性接觸，利用Ansys Workbench進行模態分析時，Modal模塊會自動的將非線性接觸默認為線性接觸，導致計算結果不準確。考慮軸系裝配過盈量，需要先在Static Structural模塊進行靜應力分析，再將靜應力分析結果導入到Modal模塊進行模態求解。在進行有接觸應力的仿真過程中將默認的線性接觸改為非線性接觸的方法是在過盈配合處創建接觸對，在選擇目標面時需要遵循的原則是：當凸面與平面或凹面接觸時，應選擇平面或凹面為目標面［9］。因此在靜力學分析設置時，設置轉軸的表面為接觸面，設置飛輪轂、推力盤、電機轉子和磁軸承轉子的內表面為目標面。

為了精確計算模型，必須在設置接觸對的過程中設置合理的FKN（接觸剛度因子）值。對于綁定和不分離的線性接觸類型的問題，默認FKN=10；對于以體積為主的非線性接觸問題，默認FKN=1或選擇“Program Controlled”；對于以彎曲問題為主的非線性接觸，可以設置為FKN=0.01～0.1［9-10］。飛輪儲能系統轉子主要以彎曲為主，因此在設置FKN時大小設置為0.07。得到的600Wh飛輪轉子網格劃分結果，如圖4所示。

圖4 轉子有限元模型Fig.4 Finite Element Model of Rotor

4.2 轉子系統的有限元仿真結果分析

由于飛輪轂承載主要的能量，產生的力學變形主要發生在飛輪轂處，因此仿真時保持其他配合處過盈量不變，只改變飛輪轂和軸配合處的過盈量，得到的最大等效應力和前兩階臨界轉速結果，如表4所示。

表4 有限元仿真結果Tab.4 Finite Element Simulation Results

從表4可以看出，過盈量每增大0.01mm，等效應力大約增大16%左右，一階和二階臨界轉速也隨之增大。其中對第一階臨界轉速影響程度大于對第二階臨界轉速的影響程度，過盈量每增加0.1mm，一階臨界轉速增大12%左右。整體趨勢是，過盈量增大，等效應力和前兩階臨界轉速呈增大趨勢。當過盈量超過0.2mm時，最大等效應力超出了飛輪轂材料的屈服強度，飛輪轂開始變形，因此飛輪轂和轉軸過盈配合處過盈量最佳取值范圍是0.07mm～0.19mm。

4.3 應力理論值與仿真值對比分析

根據應力計算公式求得軸轂配合處各過盈量下應力的理論值，得到軸轂處應力的理論值和仿真值對比表，如表5所示。從表5中可以看出理論值與仿真值基本一致。

表5 應力理論值與仿真值對比分析Tab.5 Contrastive Analysis of Stress Theory Value and Simulated Value

4.4 轉子系統的模態分析

在轉軸與飛輪轂配合處設置過盈量為0.09mm，求解得到的等效應力圖和模態振型圖，如圖5、6、7所示。

由等效應力圖5可知，軸系的最大等效應力出現在飛輪轂與軸配合出最大值為187.41MPa，小于材料的屈服強度。由振型云圖6、7可知，轉子系統的一階固有頻率為64.62Hz，其振幅表現為從中心飛輪轂處向兩端減小，最大相對位移為6.2436mm。這是一階平動狀態。轉子系統的二階固有頻率為189.31Hz，其振幅表現為兩端開始變大，中間飛輪轂處出現分層現象，最大相對位移為4.0498mm。這是二階錐動狀態。

圖5 等效應力圖Fig.5 Equivalent Stress Diagram

圖6 一階模態振型Fig.6 First Mode Shapes

圖7 二階模態振型Fig.7 Second Mode Shapes

5 結論

通過理論計算公式計算出軸上各主要部件與轉軸過盈配合處過盈量的范圍，并且建立了考慮接觸應力的飛輪儲能系統有限元模型，然后進行有限元分析。分析發現，接觸應力對轉子系統的臨界轉速有很大的影響。通過改變飛輪轂和軸的過盈配合量，得出該配合處過盈量的大致范圍為0.07mm～0.19mm，為相同類型的飛輪儲能系統設計提供了參考依據。