周傳迪，柳亦兵，朱萬程，高峻澤，張昊隨

(華北電力大學(xué) 先進飛輪儲能技術(shù)研究中心，北京 102206)

0 引 言

飛輪儲能系統(tǒng)以其功率密度高、充放電速度快、無污染等優(yōu)點在解決新能源電力系統(tǒng)的不穩(wěn)定性和電網(wǎng)的變頻調(diào)速等方面?zhèn)涫荜P(guān)注[1-3]。由于儲能需求的增加使得飛輪儲能系統(tǒng)的設(shè)計儲能量大幅提高，進而導(dǎo)致飛輪的體積、重量和工作轉(zhuǎn)速有了非常大的提升，使得飛輪轉(zhuǎn)子的臨界轉(zhuǎn)速不可避免地落在工作轉(zhuǎn)速范圍內(nèi)[4]，導(dǎo)致飛輪產(chǎn)生較大的振動。因此，有必要準(zhǔn)確預(yù)測飛輪轉(zhuǎn)子-軸承系統(tǒng)的動態(tài)特性，以避免飛輪系統(tǒng)工作在臨界轉(zhuǎn)速附近，減少共振。

為了探索飛輪轉(zhuǎn)子-軸承系統(tǒng)的動力學(xué)特性，許多學(xué)者建立了飛輪轉(zhuǎn)子-軸承系統(tǒng)的動力學(xué)模型并進行了動力學(xué)特性分析。C.Tang[5，6]等建立了柔性支承的飛輪轉(zhuǎn)子-軸承-阻尼器的動力學(xué)模型，并計算了其在不同轉(zhuǎn)速下的動力學(xué)特性。C.peng[7]建立了磁懸浮飛輪的微振動動力學(xué)模型，包括轉(zhuǎn)子動力學(xué)、磁懸浮軸承控制系統(tǒng)，并對系統(tǒng)的微振動進行了研究。H.Wang[8]利用拉格朗日方程建立了具有混合軸承支承的儲能飛輪系統(tǒng)的動力學(xué)模型，計算了飛輪系統(tǒng)的動態(tài)特性。研究表明，飛輪轉(zhuǎn)子由于轉(zhuǎn)速較高，在通過臨界轉(zhuǎn)速時會產(chǎn)生較大的振動，為了解決這一問題，利用阻尼抑制飛輪轉(zhuǎn)子振動的方法越來越受到重視。阻尼器的類型包括擠壓油膜阻尼器[9-13]、擺錘調(diào)諧質(zhì)量阻尼器[14，15]和阻尼動力吸振器[16]等。通過相關(guān)研究表明，阻尼器能有效地降低飛輪轉(zhuǎn)子的振動。使轉(zhuǎn)子更好地通過臨界轉(zhuǎn)速。對于電磁軸承支承的飛輪轉(zhuǎn)子，通過調(diào)整剛度、阻尼等控制參數(shù)[15,17，18]和控制算法[19，20]同樣也可以抑制飛輪振動，效果顯著。

上述研究中，飛輪轉(zhuǎn)子普遍質(zhì)量較小，多數(shù)都在百公斤級以下，以小型剛性飛輪轉(zhuǎn)子為主。隨著儲能需求的增加，飛輪轉(zhuǎn)子的重量和轉(zhuǎn)速都逐漸增加，由于其重量大、轉(zhuǎn)速高，結(jié)構(gòu)軸向尺寸大，導(dǎo)致飛輪轉(zhuǎn)子成為柔性轉(zhuǎn)子，而且其電機轉(zhuǎn)子部分的質(zhì)量和剛度效應(yīng)往往不可忽略，對動力學(xué)模型的精確性要求較高，其振動特性還需進行探究。本文針對某立式重載儲能飛輪轉(zhuǎn)子-軸承系統(tǒng)的動力學(xué)建模問題進行研究，加入電機轉(zhuǎn)子部分的質(zhì)量，建立飛輪轉(zhuǎn)子-電機轉(zhuǎn)子-軸承系統(tǒng)的動力學(xué)模型，并分析軸承支撐剛度、阻尼以及飛輪升速速率對飛輪轉(zhuǎn)子振動特性的影響，為飛輪轉(zhuǎn)子的動特性分析和振動抑制提供一定參考。

1 儲能飛輪系統(tǒng)結(jié)構(gòu)簡介

圖1所示為某立式重載飛輪轉(zhuǎn)子系統(tǒng)結(jié)構(gòu)示意圖，其主要結(jié)構(gòu)有上機械軸承、永磁軸承、飛輪轉(zhuǎn)子、電機轉(zhuǎn)子、電機、下機械軸承及飛輪外殼等。飛輪轉(zhuǎn)子由永磁軸承在軸向進行承載，保持轉(zhuǎn)子處于穩(wěn)定懸浮狀態(tài)。轉(zhuǎn)子上下兩端通過機械軸承支撐，承擔(dān)徑向載荷。雙向電機位于飛輪下部，電機轉(zhuǎn)子嵌套于飛輪轉(zhuǎn)子軸段上。飛輪系統(tǒng)的部件全部置于飛輪外殼中，由真空泵將內(nèi)部空氣抽出形成真空環(huán)境以降低損耗。

圖1中所示的飛輪轉(zhuǎn)子重量集中在飛輪轉(zhuǎn)子本體以及電機轉(zhuǎn)子兩部分，且飛輪轉(zhuǎn)子的重量與電機轉(zhuǎn)子的重量相差較大。由于電機轉(zhuǎn)子嵌套于飛輪轉(zhuǎn)子細長軸段上，與飛輪轉(zhuǎn)子下部軸段緊密連接，因此在建模時需要考慮電機轉(zhuǎn)子質(zhì)量及其剛度帶來的影響[21]。

圖1 飛輪結(jié)構(gòu)示意圖Fig.1 Schematic diagram of flywheel structure

2 飛輪轉(zhuǎn)子-軸承系統(tǒng)動力學(xué)建模

2.1 飛輪轉(zhuǎn)子-軸承系統(tǒng)動力學(xué)模型

采用集中質(zhì)量法，基于以下假設(shè)建立了如圖2所示的飛輪轉(zhuǎn)子-電機轉(zhuǎn)子-軸承系統(tǒng)的動力學(xué)模型：

(1)只考慮了轉(zhuǎn)子的橫向振動。

(2)飛輪轉(zhuǎn)子本體部分視為具有集中質(zhì)量的剛體。

(3)飛輪轉(zhuǎn)子的細長軸段是柔性的，將其視為無質(zhì)量的彈簧。

(4)考慮了電機轉(zhuǎn)子部分的質(zhì)量和剛度。

(5)軸承簡化為具有集中質(zhì)量、剛度和阻尼系數(shù)的線性結(jié)構(gòu)，且是各向同性的。

(6)轉(zhuǎn)子不平衡量假設(shè)均勻分布于轉(zhuǎn)子上下表面。

飛輪轉(zhuǎn)子、電機轉(zhuǎn)子及上下機械軸承均簡化為相應(yīng)的集中質(zhì)量，圖中m1、m5分別為上機械軸承和下機械軸承的質(zhì)量，m、m4分別為飛輪轉(zhuǎn)子和電機轉(zhuǎn)子的質(zhì)量。k1、c1分別上機械軸承處的剛度和阻尼，k5、c5分別下機械軸承處的剛度和阻尼。k2為飛輪轉(zhuǎn)子上部軸段剛度，k3為電機轉(zhuǎn)子上部軸段剛度，k4為電機轉(zhuǎn)子下部軸段剛度。

轉(zhuǎn)子轉(zhuǎn)動過程中由于受到不平衡力和力矩作用，轉(zhuǎn)子會發(fā)生偏轉(zhuǎn)，其各點的運動位置關(guān)系如圖2所示。圖中x1、x5分別為上軸承和下軸承質(zhì)量點位置坐標(biāo)，x2、x3分別為飛輪轉(zhuǎn)子本體部分上下端中心點坐標(biāo)，x4為電機質(zhì)心位置坐標(biāo)，xc為轉(zhuǎn)子形心的位置坐標(biāo)，xe為轉(zhuǎn)子質(zhì)心的位置坐標(biāo)，偏心距為e，θy為轉(zhuǎn)子在xoz平面的轉(zhuǎn)角。轉(zhuǎn)子為軸對稱結(jié)構(gòu)，轉(zhuǎn)速為ω，轉(zhuǎn)子的赤道轉(zhuǎn)動慣量為Jd，極轉(zhuǎn)動慣量為Jp。轉(zhuǎn)子在yoz上的各坐標(biāo)位置可參考xoz平面的坐標(biāo)位置，在此處不再列出。

圖2 飛輪轉(zhuǎn)子動力學(xué)簡化模型Fig.2 Simplified model of flywheel rotor dynamics

2.2 飛輪轉(zhuǎn)子-軸承系統(tǒng)運動方程

對圖2所示的飛輪轉(zhuǎn)子采用拉格朗日法建立轉(zhuǎn)子系統(tǒng)的運動方程，第二類拉格朗日方程如下：

(1)

飛輪轉(zhuǎn)子在轉(zhuǎn)速ω為定值狀態(tài)下，動能，勢能以及損耗能分析如下。

(1)動能T：

由圖2可得轉(zhuǎn)子形心坐標(biāo)(xc,yc)與質(zhì)心位移坐標(biāo)(xe,ye)的關(guān)系如下：

(2)

在小振動情況下，由圖2(b)可得

(3)

假設(shè)轉(zhuǎn)子本體的動能為Tz，軸承系統(tǒng)動能為Tb，電機轉(zhuǎn)子動能為Tr。則飛輪轉(zhuǎn)子-軸承系統(tǒng)的總動能為

T=Tz+Tb+Tr=

(2)勢能U：

飛輪轉(zhuǎn)子-軸承系統(tǒng)的總勢能為

(5)

(3)耗散能Z：

飛輪轉(zhuǎn)子-軸承系統(tǒng)總耗散能為

(6)

(4)外力Q：

轉(zhuǎn)子穩(wěn)態(tài)運動狀態(tài)，轉(zhuǎn)子上的外力為零，即

Q=0

(7)

將式(4)、(5)、(6)、(7)代入式(1)，得到飛輪轉(zhuǎn)子-軸承系統(tǒng)的穩(wěn)態(tài)渦動方程為

寫成矩陣形式：

(9)

(10)

m2=m/4+Jd/l22,m3=m/4-Jd/l22

(11)

式(9)為飛輪轉(zhuǎn)子-軸承系統(tǒng)的穩(wěn)態(tài)渦動方程，可以求解得到飛輪轉(zhuǎn)子-軸承系統(tǒng)的穩(wěn)態(tài)響應(yīng)。

在飛輪充放電過程中，由于轉(zhuǎn)子轉(zhuǎn)速會隨時間發(fā)生變化，此時ω是隨時間變化的變量，不再是定值，因此需要推導(dǎo)轉(zhuǎn)子的瞬態(tài)渦動響應(yīng)方程。

(12)

將式(12) 代入式(4)、(5)、(6)，并將(4)、(5)、(6)代入式(1)中，整理可得飛輪轉(zhuǎn)子-軸承系統(tǒng)的瞬態(tài)渦動方程為

(13)

其中：

(14)

由式(13)可以求解得到飛輪轉(zhuǎn)子-軸承系統(tǒng)的瞬態(tài)響應(yīng)。

3 飛輪轉(zhuǎn)子-軸承系統(tǒng)動力學(xué)特性分析

以圖1所示的大容量飛輪儲能轉(zhuǎn)子系統(tǒng)為對象進行動力學(xué)特性建模分析，根據(jù)轉(zhuǎn)子結(jié)構(gòu)尺寸，利用三維建模軟件Solidworks建立轉(zhuǎn)子模型，計算轉(zhuǎn)子的質(zhì)量、轉(zhuǎn)動慣量、軸段剛度等精確參數(shù)，參考文獻[22]中的方法，先設(shè)定了軸承的支撐剛度和阻尼的取值范圍，并根據(jù)實測數(shù)據(jù)進行了調(diào)整，最終得到各參數(shù)取值列于表1。

表1 飛輪轉(zhuǎn)子-軸承系統(tǒng)各參數(shù)值Tab.1 Parameters of flywheel rotor-bearing system

3.1 臨界轉(zhuǎn)速分析

圖3所示為飛輪轉(zhuǎn)子系統(tǒng)的坎貝爾圖。圖中所示曲線分別為飛輪轉(zhuǎn)子前兩階正反進動固有頻率隨轉(zhuǎn)速變化曲線以及轉(zhuǎn)子轉(zhuǎn)速曲線。由于轉(zhuǎn)子陀螺效應(yīng)，導(dǎo)致轉(zhuǎn)子固有頻率隨轉(zhuǎn)速而變化，其中正進動曲線隨轉(zhuǎn)速增加而增大，反進動曲線隨轉(zhuǎn)速增加而減小。飛輪轉(zhuǎn)速曲線在飛輪轉(zhuǎn)子工作范圍內(nèi)與4條進動曲線共有4個交點，然而由于轉(zhuǎn)子在實際轉(zhuǎn)動過程中，不平衡量始終與轉(zhuǎn)子同向旋轉(zhuǎn)，因此轉(zhuǎn)子的實際臨界轉(zhuǎn)速應(yīng)為轉(zhuǎn)速曲線與飛輪轉(zhuǎn)子正進動曲線的交點。從圖3所示的坎貝爾圖中可以得到飛輪轉(zhuǎn)子的一階臨界轉(zhuǎn)速為1 900 r/min，二階臨界轉(zhuǎn)速為7 300 r/min。飛輪轉(zhuǎn)子在充放電過程中，轉(zhuǎn)速在最低(3 000 r/min)和最高工作轉(zhuǎn)速(10 000 r/min)之間不斷變化，需要經(jīng)常經(jīng)過第二階臨界轉(zhuǎn)速，因此第二階臨界轉(zhuǎn)速對轉(zhuǎn)子動力學(xué)特性影響較大，需要重點關(guān)注。

圖3 飛輪轉(zhuǎn)子坎貝爾圖Fig.3 Flywheel rotor campbell diagram

3.2 軸承支撐剛度對轉(zhuǎn)子振動特性的影響

為探究軸承剛度對飛輪轉(zhuǎn)子動力學(xué)特性的影響，將上下軸承剛度在0.5e9～3e9 N/m范圍內(nèi)取不同值，計算轉(zhuǎn)子穩(wěn)態(tài)振動響應(yīng)。圖4和圖5以瀑布圖形式分別示出飛輪轉(zhuǎn)子上下軸承處的頻響函數(shù)曲線隨軸承剛度變化情況?？梢钥闯?，飛輪轉(zhuǎn)子頻響函數(shù)曲線有兩個峰值，分別對應(yīng)飛輪轉(zhuǎn)子的第一和第二階臨界轉(zhuǎn)速。

圖4 上軸承處振動響應(yīng)隨剛度變化情況Fig.4 Vibration response at the upper bearing changes with stiffness

圖5 下軸承處振動響應(yīng)隨剛度變化情況Fig.5 Vibration response at the lower bearing changes with stiffness

圖6示出飛輪轉(zhuǎn)子兩階臨界轉(zhuǎn)速隨軸承剛度變化曲線，可以看出，隨著軸承剛度增加，飛輪轉(zhuǎn)子臨界轉(zhuǎn)速逐漸增大，當(dāng)剛度較小時，兩個臨界轉(zhuǎn)速增加較為明顯，當(dāng)剛度增大到一定程度以后，臨界轉(zhuǎn)速的增加速度變緩。

圖6 飛輪轉(zhuǎn)子臨界轉(zhuǎn)速隨剛度變化曲線Fig.6 Critical speed of flywheel rotor changes with stiffness

圖7示出飛輪轉(zhuǎn)子過臨界時，上下軸承處的共振幅值隨軸承支撐剛度的變化情況，可以看到，隨著剛度增加，上下軸承處共振幅值出現(xiàn)高低交錯的波動現(xiàn)象，在某些剛度下，幅值達到最大值，在某些剛度降到最小，而且波動周期不是定值。這種現(xiàn)象表明，軸承剛度在一定范圍的調(diào)整，可以改變轉(zhuǎn)子的動剛度，降低轉(zhuǎn)子系統(tǒng)振動。

圖7 飛輪轉(zhuǎn)子臨界轉(zhuǎn)速對應(yīng)峰值隨剛度變化曲線Fig.7 Critical speed of flywheel rotor corresponds to peak value versus stiffness change curve

3.3 軸承阻尼對轉(zhuǎn)子振動特性的影響

為探究軸承阻尼對飛輪轉(zhuǎn)子動力學(xué)特性的影響，將上下軸承阻尼在1 500～1e6 Ns/m范圍內(nèi)取不同值，計算飛輪轉(zhuǎn)子穩(wěn)態(tài)振動響應(yīng)。圖8和圖9示出飛輪轉(zhuǎn)子上下軸承處的頻響函數(shù)曲線隨軸承阻尼變化的瀑布圖，圖10和圖11為兩個臨界轉(zhuǎn)速及其共振幅值隨軸承阻尼變化關(guān)系曲線，可以看到，在文中所計算的阻尼范圍內(nèi)，兩階臨界轉(zhuǎn)速隨阻尼的變化而有小幅度改變，在阻尼范圍為1 500～2e5 Ns/m區(qū)間內(nèi)，阻尼對系統(tǒng)臨界轉(zhuǎn)速影響較小，兩階臨界轉(zhuǎn)速改變不明顯，在阻尼范圍為2e5～1e6 Ns/m區(qū)間內(nèi)，第二階臨界轉(zhuǎn)速相較于第一階臨界轉(zhuǎn)速變化明顯。此外，共振幅值隨著阻尼的增大而呈明顯的下降趨勢，說明軸承阻尼對飛輪轉(zhuǎn)子軸承處的振動具有抑制作用，且當(dāng)阻尼較小時，振動幅值下降較為明顯，當(dāng)阻尼增大到一定程度時振動幅值下降變緩。

圖8 上軸承處振動響應(yīng)隨阻尼變化情況Fig.8 Vibration response at the upper bearing changes with damping

圖9 下軸承處振動響應(yīng)隨阻尼變化情況Fig.9 Vibration response at the lower bearing changes with damping

圖10 飛輪轉(zhuǎn)子臨界轉(zhuǎn)速隨阻尼變化曲線Fig.10 Curve of the critical speed of flywheel rotor with damping

圖11 飛輪轉(zhuǎn)子臨界轉(zhuǎn)速峰值隨阻尼變化曲線Fig.11 Curve of the critical speed peak of flywheel rotor with damping

3.4 升速速率對升速瞬態(tài)響應(yīng)的影響

飛輪儲能系統(tǒng)通過飛輪轉(zhuǎn)子升降轉(zhuǎn)速實現(xiàn)動能的存儲和釋放，因此在運行過程中，轉(zhuǎn)速可能反復(fù)通過臨界轉(zhuǎn)速，為了避免轉(zhuǎn)子在通過臨界轉(zhuǎn)速時產(chǎn)生較大振動，可以通過使飛輪轉(zhuǎn)子快速通過其臨界轉(zhuǎn)速來達到降低振動的目的。本節(jié)分析飛輪轉(zhuǎn)子在不同升速速率下的振動響應(yīng)的變化情況，為飛輪轉(zhuǎn)子如何通過臨界轉(zhuǎn)速提供理論依據(jù)。

對飛輪轉(zhuǎn)子的瞬態(tài)渦動方程進行求解，參數(shù)取表1所列參數(shù)，計算得到飛輪轉(zhuǎn)子在不同升速速率下的振動響應(yīng)情況。共計算了4種升速速率下的飛輪轉(zhuǎn)子振動響應(yīng)，升速速率分別取2π rad/s、5π rad/s、10π rad/s以及20π rad/s。圖12示出不同升速速率下飛輪轉(zhuǎn)子上軸承處的振動響應(yīng)幅值曲線，圖13示出了飛輪轉(zhuǎn)子通過臨界轉(zhuǎn)速時的共振峰值大小變化趨勢，圖14示出了飛輪轉(zhuǎn)子臨界轉(zhuǎn)速變化趨勢，可以看出，當(dāng)轉(zhuǎn)子以低速率升速時，過臨界的振動幅度較大，當(dāng)轉(zhuǎn)子升速速率提高時，轉(zhuǎn)子過臨界振動幅值明顯減小，而臨界轉(zhuǎn)速略有增大。分析表明，通過提升轉(zhuǎn)子升速速率，可以在一定程度上降低轉(zhuǎn)子過臨界的振動幅度。

圖12 飛輪轉(zhuǎn)子不同升速速率下的振動響應(yīng)Fig.12 Vibration response of flywheel rotor at different speed-up rates

圖13 飛輪轉(zhuǎn)子共振峰值隨升速速率變化曲線Fig.13 Curve of resonant peak value of flywheel rotor with speed-up rate

圖14 飛輪轉(zhuǎn)子臨界轉(zhuǎn)速隨升速速率變化曲線Fig.14 Critical Speed of flywheel rotor variation curve with speed up rate

3.5 實驗驗證

為掌握飛輪系統(tǒng)的運行特性，并驗證文中建立的飛輪轉(zhuǎn)子動力學(xué)模型的合理性，對大容量儲能飛輪轉(zhuǎn)子系統(tǒng)進行升速試驗。在飛輪轉(zhuǎn)子外殼的上軸承和下軸承位置各放置一個振動傳感器，用于測試升速過程中飛輪轉(zhuǎn)子上下軸承處的振動情況。升速過程從轉(zhuǎn)速為0開始，直至升速到10 000 r/min結(jié)束。

圖15示出升速過程中的軸承振動響應(yīng)幅值變化趨勢，從圖中可以看出，在升速過程中，上下軸承處均有三個較為明顯的共振峰，對應(yīng)轉(zhuǎn)速分別為1 712 r/min，2 661 r/min以及7 336 r/min。與圖12數(shù)值計算結(jié)果相對比，數(shù)值計算結(jié)果與飛輪轉(zhuǎn)子實測結(jié)果基本符合，略有差異。飛輪數(shù)值計算結(jié)果在轉(zhuǎn)速4 000 r/min以內(nèi)范圍內(nèi)有一個共振峰，而實際測試結(jié)果有多個共振峰，分析原因可能是飛輪升速實測過程中傳感器安裝于飛輪外殼上，測得的數(shù)據(jù)受到飛輪外殼和其他部件的影響，從而在該轉(zhuǎn)速范圍內(nèi)產(chǎn)生了共振所致。

圖15 飛輪轉(zhuǎn)子實測升速振動響應(yīng)曲線Fig.15 Flywheel rotor actual measured speed-up vibration response curve

4 結(jié) 論

本文建立了立式重載飛輪轉(zhuǎn)子的動力學(xué)模型，對飛輪轉(zhuǎn)子動力學(xué)方程進行了動特性參數(shù)求解，獲得了飛輪轉(zhuǎn)子的固有頻率、臨界轉(zhuǎn)速以及升速過程的響應(yīng)曲線。通過分析軸承支撐剛度、阻尼以及飛輪轉(zhuǎn)子升速速率對飛輪轉(zhuǎn)子動力學(xué)特性的影響，可以得出以下結(jié)論：

(1)本文研究的大容量儲能飛輪轉(zhuǎn)子在0～10 000 r/min轉(zhuǎn)速范圍內(nèi)存在兩階臨界轉(zhuǎn)速，與飛輪轉(zhuǎn)子實際振動測試所獲取的結(jié)果相一致，表明建立的飛輪轉(zhuǎn)子-軸承系統(tǒng)動力學(xué)模型較為合理。

(2)飛輪轉(zhuǎn)子的臨界轉(zhuǎn)速隨著軸承支撐剛度的增大而增大，當(dāng)支撐剛度增大到一定程度時，臨界轉(zhuǎn)速增加緩慢；飛輪轉(zhuǎn)子過臨界時的共振幅值隨剛度增加出現(xiàn)波動現(xiàn)象，在某些剛度下出現(xiàn)極值；飛輪轉(zhuǎn)子臨界轉(zhuǎn)速隨軸承阻尼的改變而略有變化，但是飛輪轉(zhuǎn)子振動幅值隨阻尼的增大呈減小的趨勢。

(3)通過改變飛輪轉(zhuǎn)子的升速速率可以改變飛輪轉(zhuǎn)子振動響應(yīng)的幅值大小，升速速率越快則飛輪轉(zhuǎn)子的振動響應(yīng)幅值越小，因此在飛輪轉(zhuǎn)子通過臨界轉(zhuǎn)速時可考慮通過控制策略提升轉(zhuǎn)子的升速速率，進而降低飛輪轉(zhuǎn)子過臨界時產(chǎn)生的振動，使飛輪轉(zhuǎn)子穩(wěn)定運行。