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

（1北京泓慧國際能源技術(shù)發(fā)展有限公司，北京101300；2北京航空航天大學(xué)，北京100191）

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

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

1 轉(zhuǎn)子系統(tǒng)分析理論基礎(chǔ)

1.1 Timoshenko 梁理論

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

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

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

圖1 Timoshenko 梁Fig.1 Timoshenko beam

1.2 轉(zhuǎn)子單元運動方程

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

引入Timoshenko 梁的典型位移函數(shù)[15]：

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

圖2 軸單元Fig.2 Shaft element

軸段單元內(nèi)任一截面的位移用該單元的節(jié)點位移表示為

只考慮軸段的橫向振動，且轉(zhuǎn)子軸段截面為對稱圓截面，則在固定坐標(biāo)系中的動能和勢能表達(dá)式為

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

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

式中

2 轉(zhuǎn)子動力學(xué)有限元分析

2.1 臨界轉(zhuǎn)速計算

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

圖3 飛輪轉(zhuǎn)子有限元模型Fig.3 FE model of flywheel rotor

圖4 Campbell 圖Fig.4 Campbell diagram

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

表1 臨界轉(zhuǎn)速的計算Table 1 The calculation of the critical speed

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

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

圖5 轉(zhuǎn)子平動模態(tài)振型Fig.5 Rotor translation mode

圖6 轉(zhuǎn)子擺動模態(tài)振型Fig.6 Rotor vibration mode

圖7 轉(zhuǎn)子一階彎曲模態(tài)振型（oxz 平面）Fig.7 First order bending mode of rotor（oxz plane）

圖8 轉(zhuǎn)子一階彎曲模態(tài)（oyz 平面）Fig.8 First order bending mode of rotor（oyz plane）

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

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

2.2 不平衡響應(yīng)的計算

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

圖9 不平衡力激勵下的轉(zhuǎn)子軸的位移響應(yīng)圖Fig.9 Displacement response diagram of rotor bearing under excitation of unbalanced force

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

圖10 轉(zhuǎn)軸平動的軌跡圖Fig.10 Trajectory chart of translation of spindle

圖11 轉(zhuǎn)軸擺動的軌跡圖Fig.11 Trajectory chart of axis swing

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

3 飛輪測試

3.1 測試方法

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

圖13 測試的儲能飛輪產(chǎn)品Fig.13 The energy storage flywheel for testing

3.2 試驗結(jié)果分析

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

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

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

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

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

4 結(jié) 論

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

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

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

（4）試驗測試結(jié)果驗證了仿真分析結(jié)果的正確性。說明了對250 kW/3 kW·h 儲能飛輪的轉(zhuǎn)子動力學(xué)設(shè)計的參數(shù)是合理的，滿足產(chǎn)品的設(shè)計指標(biāo)要求。