儲能飛輪的轉子動力學分析與測試試驗

彭 龍，李光軍，崔亞東，王 芳

（1北京泓慧國際能源技術發展有限公司，北京101300；2北京航空航天大學，北京100191）

飛輪儲能技術是通過電動/發電機將飛輪機械能與電能相互轉化的技術，具有充電時間短、響應速度快、能量密度較高、使用壽命長、環保無污染等優點[1-2]。隨著高強度復合材料、大功率電動/發電機、磁懸浮、真空和電力電子等技術的發展，飛輪儲能技術目前已經在電網調頻、電能質量控制、衛星儲能/姿態控制、車輛制動能量回收、不間斷電源和高功率脈沖電源等領域獲得了廣泛的應用[3-7]。近年來我國在飛輪儲能關鍵技術方面有所突破，部分大學和科研機構已經完成飛輪儲能系統試驗裝置或工程樣機的研制[8-10]。

飛輪轉子技術與磁軸承技術一直是飛輪儲能的關鍵技術[9]。本文以某公司的250 kW/3 kW·h 儲能飛輪產品為研究對象，利用轉子動力學理論與有限元仿真的方法分析，對儲能飛輪的轉子運動特性進行研究，通過對儲能飛輪轉子運動過程中的振動測試試驗來驗證分析結果的正確性。

1 轉子系統分析理論基礎

1.1 Timoshenko 梁理論

轉子動力學有限元法能對大型復雜轉子系統直接列出運動方程，考慮更復雜的邊界條件和更多的內部影響因素。有限單元法常將軸段用連續梁來描述，不同的梁模型對高速旋轉機械的動力學響應分析結果會有較大的影響。Timoshenko 梁理論考慮了剪切變形和轉動慣量的影響，能準確分析轉軸的高階橫向振動模態，被廣泛應用在旋轉軸的分析計算中[11-12]。

本文采用基于Timoshenko 梁理論的有限元法計算轉子的臨界轉速及振動模態。圖1 為Timoshenko模型，運動方程包含剪切變形、轉動慣量、陀螺力矩，在復坐標中可寫成式(1)形式[13]

式中，q、φ為復橫向位移和截面偏轉角位移，q=x+iy，φ=θx+iθy；F和M為復剪切力和彎矩，F=Fx+iFy，M=Mx+iMy；ρ、E及G分別材料的密度、彈性模量及剪切模量；A、Iρ及Id分別為截面積、極轉動慣量及直徑轉動慣量，κ為Timoshenko 剪切系數。剪切系數κ對Timoshenko 旋轉梁的振動分析結果有很大影響，其大小與材料特性和截面幾何形狀有關[14]，對圓截面κ=6(1+μ)2/(7+12μ+4μ4)。

圖1 Timoshenko 梁Fig.1 Timoshenko beam

1.2 轉子單元運動方程

采用如圖2 所示的彈性軸段單元來離散化轉子，該單元的廣義坐標是兩端節點的位移，僅考慮橫向振動時，單元位移向量包括4 個位移和4 個轉角

引入Timoshenko 梁的典型位移函數[15]：

式中，EI為梁的彎曲剛度，κGA剪切剛度。

圖2 軸單元Fig.2 Shaft element

軸段單元內任一截面的位移用該單元的節點位移表示為

只考慮軸段的橫向振動，且轉子軸段截面為對稱圓截面，則在固定坐標系中的動能和勢能表達式為

將式(5)代入上邊兩式的

忽略軸段內阻尼將(8)式代入Lagrange 方程，可得軸段單元運動方程[16]

式中

2 轉子動力學有限元分析

2.1 臨界轉速計算

250 kW/3 kW·h 儲能飛輪轉子組件，轉子長1400 mm，最大直徑600 mm，轉子為單圓盤偏置結構，兩端由金屬疊片組成的磁懸浮軸承支撐部位。飛輪的額定工作轉速為5000～9000 r/min。利用ANSYS 有限元分析軟件，建立轉子的有限元模型，采用軸對稱實體單元SOLID272 劃分網格，設置圓周方向上的傅里葉節點數為12，兩端用COMBIN14單元支撐，如圖3 所示。飛輪轉子兩端磁軸承支撐位置，通過利用彈簧單元等效磁軸承的彈性支撐形式，支撐剛度為5×105N/m。通過APDL 編程進行轉子臨界轉速的求解。通過對轉子施加0～20000 r/min 的轉速。打開柯氏效應（coriolis effect），采用reduced damped 方法提取前4 階模態頻率與振型，并繪制坎貝爾圖（Campbell diagram）如圖4所示，進而求得飛輪轉子的臨界轉速與振型。渦動頻率曲線與轉頻的直線（斜率為1）的交點為臨界轉速。計算不同剛度下轉子的一階、二階、三階臨界轉速。見表1。

圖3 飛輪轉子有限元模型Fig.3 FE model of flywheel rotor

圖4 Campbell 圖Fig.4 Campbell diagram

表1為磁軸承剛度對臨界轉速的影響，隨著磁軸承的剛度增加，轉子的臨界轉速也在增大?？衫谜{節磁軸承的剛度以調整轉子的固有振動頻率，避開有害共振。

表1 臨界轉速的計算Table 1 The calculation of the critical speed

實際表明，磁懸浮軸承控制柔性轉子的困難很大[16-17]。這是因為一般磁懸浮軸承的剛度要比滾動軸承小1～2 個數量級，而且控制參數的穩定域限制了磁懸浮軸承的剛度阻尼的變化范圍。因此磁懸浮控制是有限的，為了避免這個問題，要求磁軸承支撐的為剛性轉子結構，即額定轉速下無彎曲模態臨界轉速。

通過轉子的一階臨界轉速對應的振型（圖5）可以看出轉子在做徑向的平動，左端位移小，右端位移大，這是由于轉子的平動模態產生。轉子的二階臨界轉速對應的振型（圖6）是由于轉子的擺動模態產生的徑向的擺動，基本以圓盤為中心，左端偏離距離大于右端。圖7、圖8 分別表示轉子在兩個不同方向的一階彎曲模態的振型，在坎貝爾圖中，分別對應的是由一階彎曲模態產生的轉子的反向渦動頻率與正向渦動頻率。如圖4 所示，在20000 r/min的轉速范圍內，轉頻與一階彎曲模態產生的反向渦動頻率曲線相交，為三階臨界轉速，振型如圖7所示。

圖5 轉子平動模態振型Fig.5 Rotor translation mode

圖6 轉子擺動模態振型Fig.6 Rotor vibration mode

圖7 轉子一階彎曲模態振型（oxz 平面）Fig.7 First order bending mode of rotor（oxz plane）

圖8 轉子一階彎曲模態（oyz 平面）Fig.8 First order bending mode of rotor（oyz plane）

通過臨界轉速所對應的振型可以看出轉子的一階、二階臨界轉速為轉子剛性振動，三階臨界轉速為轉子彎曲的柔性振動。儲能飛輪設計時，工作轉速范圍應避開一階、二階剛性臨界轉速，且不超過第三階彎曲臨界轉速。由于在臨界轉速附近一定范圍內也會引起振動，所以在實際設計中轉子的最高工作轉速應小于0.7 倍的彎曲臨界轉速，才能保證轉子為剛性轉子，避免受到臨界轉速的干擾。

儲能飛輪的最高工作轉速為9000 r/min，小于計算得出的第三階彎曲臨界轉速13238 r/min 的0.7倍，因此飛輪的轉速設定是合理的。

2.2 不平衡響應的計算

轉子模型采用了BEAM188 單元建模，使用COMBIN214 單元支承轉子，支承剛度為1×106N/m。在轉子中心施加1.88×10-4N 的不平衡力。利用完全法對轉子不平衡響應進行諧響應分析。設置不平衡力的激勵與轉子轉頻同步，指定執行加載步驟的子步驟數為500，定義了諧波分析的頻率范圍為0～500 Hz，通過命令打開科里奧利效應，然后進行求解。繪制出在0～500 Hz 轉頻內不平衡量激勵下的轉軸上三點（上徑向磁軸承支撐位置、轉子中心位置、下徑向磁軸承支撐位置）的位移響應圖，如圖9 所示。

圖9 不平衡力激勵下的轉子軸的位移響應圖Fig.9 Displacement response diagram of rotor bearing under excitation of unbalanced force

在不平衡力激勵下的轉軸出現三次明顯波動，對應的頻率分別為10.4 Hz、35.2 Hz、280.0 Hz，對應轉軸的軌跡圖分別如圖10、圖11、圖12 所示。通過轉軸的軌跡圖可以看出，轉子在不平衡力激勵下產生的三次振動，與2.1 節中分析的三階臨界轉速相對應。由表1 可知，在支承剛度為1×106N/m時，轉子的臨界轉速分別為534.60 r/min、1788.95 r/min、13263.92 r/min（即8.91 Hz、29.8 Hz、221.06 Hz），與在不平衡力激勵下的振動頻率存在約20%的誤差，主要原因是由于采用了BEAM188 單元求轉動軌跡簡化模型造成的。綜上分析得出，在頻率280.0 Hz 時是由一階彎曲模態反向渦動引起的臨界轉速，在不平衡激勵力作用下，在其臨界轉速附近具有較大的振動幅值，如圖9 所示。通過對轉子進行不平衡響應分析，驗證了臨界轉速求解的正確性，并說明了一階彎曲模態引起的臨界轉速在一定范圍內都會引起轉子的振動，在臨界轉速時達到峰值。

圖10 轉軸平動的軌跡圖Fig.10 Trajectory chart of translation of spindle

圖11 轉軸擺動的軌跡圖Fig.11 Trajectory chart of axis swing

圖12 轉軸一階彎曲的軌跡圖Fig.12 Trajectory chart of first order bending rotation axis

3 飛輪測試

3.1 測試方法

250 kW/3 kW·h 的儲能飛輪產品如圖13 所示，它主要由儲能飛輪本體、變頻器、磁軸承控制系統、真空系統及主控系統組成。測試方法是在儲能飛輪從0～9000 r/min 的升速過程中，利用徑向磁軸承的電渦流位移傳感器采集轉子徑向位移信號，通過磁軸承的控制器輸出到示波器，每間隔200 r/min 存儲一次位移信號的波形和頻譜數據。然后利用MATLAB 繪制頻譜的三維瀑布圖，與仿真分析的Campbell 圖進行對比分析。利用波形數據繪制出最大振動位移響應圖與不平衡力下的轉子位移響應圖進行對比分析，通過對比分析結果來驗證仿真分析的準確性。

圖13 測試的儲能飛輪產品Fig.13 The energy storage flywheel for testing

3.2 試驗結果分析

由于飛輪轉子是立式旋轉，主要儲能的圓盤位于轉子偏下部位，在轉子擺動時，上徑向磁軸承位置的擺動幅值比較大，通過圖6 能明顯看出，所以采集上徑向磁軸承Y向的轉子振動位移數據進行分析，利用MATLAB 軟件編程繪制了上徑向磁軸承位置振動位移的頻譜瀑布圖，如圖14 所示。

從圖14中，明顯能看出轉子的轉頻及3倍轉頻，它主要是由不平衡量激勵引起，還有由一階彎曲模態而產生的正向渦動與反向渦動。與圖4 的Campbell 圖對比發現，由于儲能飛輪結構部件模態響應及其他部件的影響，對采集的振動位移信號產生了許多低頻干擾信號，致使轉子的平動和擺動頻率無法明顯從圖14 中明顯看出，此點不足在下文中能通過對比振動位移曲線來進行分析驗證。通過對比圖4與圖14能發現由一階彎曲模態引起正向渦動頻率和反向渦動頻率曲線是完全一致的，證明了臨界轉速的分析結果與實際測試結果一致。

圖14 頻譜瀑布圖Fig.14 Spectral waterfall chart

通過實測的上徑向磁軸承Y向的轉子振動位移數據，繪制出最大振動位移隨轉速變化的曲線，如圖15 所示。從圖9 中提取出在不平衡力的激勵下上徑向磁軸承支撐位置轉軸的振動響應，如圖16 所示。通過圖15 和圖16 的對比得出，在0～9000 r/min（0～150 Hz）范圍內，實測與仿真分析出的曲線變化規律完全相同。在低轉速下，出現了兩次峰值，第一次峰值明顯高于第二次峰值，然后曲線降低，趨近于平緩。通過不平衡響應仿真分析結果，可知實測曲線的兩個峰值分別對應轉子平動模態的臨界轉速與擺動模態的臨界轉速，驗證了分析結果前兩階剛性臨界轉速的正確性。

圖16 上徑向磁軸承支撐位置轉軸的振動位移響應圖（仿真）Fig.16 Vibration displacement response diagram of the support position of the upper radial magnetic bearing(simulation)

4 結 論

（1）利用ANSYS 軟件仿真分析，繪制了250 kW/3 kW·h 的儲能飛輪坎貝爾圖，在不同支撐剛度下計算出轉子臨界轉速，結果得出隨著支撐的剛度增加，轉子的臨界轉速也在增大?？衫谜{節主動磁軸承的剛度以調整轉子的固有振動頻率，避開有害共振。

（2）通過對轉子進行不平衡響應分析，計算出了轉子上不同位置的振動位移曲線，和每個臨界轉速對應的運動軌跡圖，形象的展示了轉子在臨界轉速時的運動形態，并且充分驗證了臨界轉速求解的正確性。

（3）通過分析得出，儲能飛輪的工作轉速范圍應避開前兩階剛性臨界轉速，且不能超過彎曲臨界轉速的0.7 倍。

（4）試驗測試結果驗證了仿真分析結果的正確性。說明了對250 kW/3 kW·h 儲能飛輪的轉子動力學設計的參數是合理的，滿足產品的設計指標要求。