彭 龍，李光軍，崔亞?wèn)|，王 芳

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

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

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

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

1.1 Timoshenko 梁理論

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

本文采用基于Timoshenko 梁理論的有限元法計(jì)算轉(zhuǎn)子的臨界轉(zhuǎn)速及振動(dòng)模態(tài)。圖1 為Timoshenko模型，運(yùn)動(dòng)方程包含剪切變形、轉(zhuǎn)動(dòng)慣量、陀螺力矩，在復(fù)坐標(biāo)中可寫(xiě)成式(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)動(dòng)慣量及直徑轉(zhuǎn)動(dòng)慣量，κ為Timoshenko 剪切系數(shù)。剪切系數(shù)κ對(duì)Timoshenko 旋轉(zhuǎn)梁的振動(dòng)分析結(jié)果有很大影響，其大小與材料特性和截面幾何形狀有關(guān)[14]，對(duì)圓截面κ=6(1+μ)2/(7+12μ+4μ4)。

圖1 Timoshenko 梁Fig.1 Timoshenko beam

1.2 轉(zhuǎn)子單元運(yùn)動(dòng)方程

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

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

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

圖2 軸單元Fig.2 Shaft element

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

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

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

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

式中

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

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

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

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

圖4 Campbell 圖Fig.4 Campbell diagram

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

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

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

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

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

圖6 轉(zhuǎn)子擺動(dòng)模態(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）

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

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

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

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

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

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

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

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

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

3 飛輪測(cè)試

3.1 測(cè)試方法

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

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

3.2 試驗(yàn)結(jié)果分析

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

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

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

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

圖16 上徑向磁軸承支撐位置轉(zhuǎn)軸的振動(dòng)位移響應(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 的儲(chǔ)能飛輪坎貝爾圖，在不同支撐剛度下計(jì)算出轉(zhuǎn)子臨界轉(zhuǎn)速，結(jié)果得出隨著支撐的剛度增加，轉(zhuǎn)子的臨界轉(zhuǎn)速也在增大。可利用調(diào)節(jié)主動(dòng)磁軸承的剛度以調(diào)整轉(zhuǎn)子的固有振動(dòng)頻率，避開(kāi)有害共振。

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

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

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