基于HyperMesh的磁懸浮儲能飛輪模態(tài)分析

查生凱， 唐景春

（合肥工業(yè)大學(xué) 機械與汽車工程學(xué)院，合肥 230009）

0 引 言

飛輪儲能作為一種先進的物理儲能技術(shù)已經(jīng)部分應(yīng)用于航空航天、電動汽車、電磁彈射、通信、醫(yī)療、電力等領(lǐng)域［1］。飛輪的儲能量與轉(zhuǎn)子轉(zhuǎn)動慣量和飛輪轉(zhuǎn)速的平方成正比，因此提高飛輪儲能密度有兩種途徑：一是增加飛輪轉(zhuǎn)子的輪緣質(zhì)量以增大轉(zhuǎn)動慣量，二是提高飛輪轉(zhuǎn)速。前者主要用于沖壓機械，由于其體積大不便推廣使用。后者應(yīng)用領(lǐng)域廣闊，但受到三個方面的制約：其一是飛輪在高速連續(xù)運轉(zhuǎn)時，軸承由于摩擦導(dǎo)致磨損、發(fā)熱、升溫、失效，潤滑油變質(zhì)，限制了飛輪轉(zhuǎn)速的提高；其二是要確保輪緣材料在離心力作用下不會因強度不足發(fā)生失效；其三是高速旋轉(zhuǎn)時，飛輪受到的空氣阻力和軸承的摩擦力矩會顯著增加，使得儲能效率大大減低。近年來磁懸浮軸承技術(shù)與復(fù)合材料技術(shù)的發(fā)展正在逐步解決上述難題。

美國、瑞士等工業(yè)發(fā)達國家長期從事電磁軸承儲能飛輪的研究與開發(fā)，同時進行了磁軸承相關(guān)技術(shù)研究［2］。文獻［3］分析了多層復(fù)合材料轉(zhuǎn)子的應(yīng)力分布與自由模態(tài)。文獻［4］根據(jù)平面應(yīng)力假設(shè)分析得到復(fù)合材料儲能飛輪的最大徑向應(yīng)力位置。

本文針對一種磁懸浮儲能飛輪，利用HyperMesh軟件對每個零件網(wǎng)格劃分，并由零件裝配生成磁懸浮儲能飛輪的有限元模型。利用彈簧單元springs模擬磁軸承。利用Rasioss求解器求解系統(tǒng)模態(tài)頻率與振型。根據(jù)分析結(jié)果可以有效預(yù)估系統(tǒng)的振動特性，為結(jié)構(gòu)優(yōu)化與性能改進提供理論依據(jù)。

1 儲能飛輪系統(tǒng)的結(jié)構(gòu)設(shè)計

根據(jù)儲能量的設(shè)計要求，利用三維設(shè)計軟件Pro/E建模和運動仿真，運用ANSYS軟件對主要部件進行力學(xué)分析，完成磁懸浮儲能飛輪樣機的初步結(jié)構(gòu)設(shè)計工作，其結(jié)構(gòu)如圖1所示，主要有飛輪轉(zhuǎn)子、電機、機械保護軸承、徑向磁軸承、軸向混合磁軸承、真空室（維持內(nèi)部真空，減小空氣阻力）等關(guān)鍵部件組成。圖2為儲能飛輪系統(tǒng)實物圖。

2 模態(tài)分析原理

工程中的彈性結(jié)構(gòu)可以離散化為有限個質(zhì)量元件、彈性元件和阻尼元件組成的N自由度線性系統(tǒng)，其動力學(xué)微分方程為［5］

其中：［M］為結(jié)構(gòu)質(zhì)量矩陣；［C］為結(jié)構(gòu)阻尼矩陣；［ K］為剛度矩陣；｛F｝為激勵向量；｛X｝、｛ X˙｝、｛ X¨｝分別為各離散質(zhì)量的N維位移、速度、加速度向量。

對上式作傅立葉變換，可得

式中：｛ X（ ω）｝ 為位移｛ X｝ 的傅立葉變換；｛ F（ ω）｝ 為激勵｛F｝的傅立葉變換；［H（ω）］為系統(tǒng)頻率響應(yīng)函數(shù)矩陣。

［H（ω）］中元素H1p（ω）的含義是：當(dāng)其它點上激勵為零時，1點的響應(yīng)譜與p點激勵譜之比，即

通過系統(tǒng)測試獲得矩陣中元素，得到頻響函數(shù)矩陣，即可獲得系統(tǒng)的全部動力學(xué)特性。

在結(jié)構(gòu)阻尼、小阻尼或比例阻尼情況下，模態(tài)參數(shù)與頻響函數(shù)之間的關(guān)系為

式中：｛ φ｝i為第 i階振型；ki為第 i階模態(tài)剛度；mi為第 i階模態(tài)質(zhì)量；ci為第i階模態(tài)阻尼。

運用模態(tài)參數(shù)識別理論進行模態(tài)參數(shù)識別，可獲得各階模態(tài)參數(shù)。

3 建立有限元模型與施加邊界條件

3.1 零件有限元建模

將儲能飛輪簡化為底座、下支撐、中間圓筒、電機定子、上支撐、真空罩、轉(zhuǎn)子7大部分，利用三維設(shè)計軟件分別建立實體模型，在有限元軟件HyperMesh環(huán)境下劃分網(wǎng)格。轉(zhuǎn)子的有限元模型如圖3所示。

儲能飛輪有限元模型共包含97 620個節(jié)點，208 773個單元。根據(jù)零件的實際結(jié)構(gòu)，底座、下支撐與真空罩采用一階六面體單元Hex8單元劃分網(wǎng)格，其他零件采用一階四面體單元tetra4單元劃分網(wǎng)格。

3.2 儲能飛輪零件結(jié)合部建模

定子是由多個零件通過螺栓聯(lián)接組成的，轉(zhuǎn)子由磁軸承支撐。這些零件間關(guān)聯(lián)方式的模擬是否合理直接影響計算結(jié)果的準(zhǔn)確性。

3.2.1 螺栓聯(lián)接處理

螺栓聯(lián)接作為常用聯(lián)接在機械結(jié)構(gòu)中被廣泛使用，本系統(tǒng)共有36處螺栓聯(lián)接。采用connectors中的bolt單元模擬螺栓聯(lián)接，其主要思想是用2組剛性單元RBE2分別聯(lián)接2個被聯(lián)接件的內(nèi)孔節(jié)點，然后再用梁單元聯(lián)接2組RBE2單元的中心節(jié)點。

3.2.2 接觸處理

為了真實地反映系統(tǒng)的動態(tài)特性，不應(yīng)將組成定子的各零件作一體化處理，應(yīng)在零件間的接觸面上定義摩擦接觸單元。接觸力與位移之間的關(guān)系如圖4所示。

圖4 接觸力與位移之間的關(guān)系

圖4中，縱坐標(biāo)FX為接觸力，橫坐標(biāo)Us-Um為兩接觸面的相對位移量，U0為初始間隙量，STIFF為接觸面相對剛度。STIFF的值由接觸面周圍單元的剛度決定。由圖4可知：當(dāng)兩接觸面由初始位置壓縮時接觸剛度為定值，接觸力與壓縮位移量成正比；當(dāng)兩接觸面由初始位置分離時，為了避免出現(xiàn)奇異值取接觸面的剛度為10-14×STIFF。

接觸面間的摩擦力與滑動距離之間的關(guān)系如圖5所示。

圖5 摩擦力與滑動距離之間的關(guān)系

圖5中，橫坐標(biāo)d為兩接觸面相對滑動距離，F(xiàn)為接觸面摩擦力，μ1為靜摩擦因數(shù)，μ2為動摩擦因數(shù)，F(xiàn)x為法向接觸力。由圖5可知：隨著滑動距離d的增大接觸面摩擦力F逐漸增加到靜摩擦力 μ1×Fx，當(dāng)滑動距離 d 繼續(xù)增加時，F(xiàn) 為定值 μ2×Fx。

HyperMesh11.0定義了三種標(biāo)準(zhǔn)接觸類型：SLIDE（滑動），STICK（黏結(jié)），F(xiàn)REEZE（ 固聯(lián)）［6］。 在分析過程中，軟件會自動計算這三種接觸的接觸剛度與摩擦因數(shù)。根據(jù)各零件的實際接觸情況，本文選擇的接觸類型是SLIDE。

3.2.3 磁軸承支撐處理

上下徑向磁軸承固定件通過螺釘固定于上下支撐，上下徑向磁軸承轉(zhuǎn)動件固定于轉(zhuǎn)軸上的軸套，軸向磁軸承轉(zhuǎn)動件嵌裝于飛輪體的下端面，軸向磁軸承固定件嵌裝于下支撐的上端面。為了簡化模型，采用了一體化處理。磁軸承定轉(zhuǎn)子之間通過彈簧阻尼單元springs聯(lián)接起來。磁軸承的剛度為［7］

式中：p為軸承的平均周長；J1、J2為磁化強度；μ0為真空磁導(dǎo)率；R為間隙量；S1、S2為磁軸承的橫截面積。

根據(jù)磁軸承的結(jié)構(gòu)特點，將徑向磁軸承簡化為8根彈簧［8］， 其 剛 度 為2.458×106N/m，軸向磁軸承簡化為8根彈簧，其剛度為0.162×106N/m。其中徑向磁軸承等效的彈簧單元聯(lián)接如圖6所示。

3.3 施加邊界條件

系統(tǒng)由4個水平調(diào)節(jié)螺栓固定支撐，給有限元模型施加如下邊界條件：用剛性單元RBE2聯(lián)接螺栓內(nèi)孔與中心節(jié)點，約束4個中心節(jié)點的6個自由度。其中一個約束如圖7所示。

圖7 邊界條件

4 模擬計算與結(jié)果分析

磁懸浮儲能飛輪的模態(tài)分析在HyperMesh環(huán)境下，采用Radioss求解器求解模態(tài)。Radioss提供了兩種蘭索斯法（Lanczos） 用于提取特征值與特征向量。如果EIGRL卡片中沒有定義上限且待求的特征值階數(shù)小于50，系統(tǒng)會自動選擇速度較快的蘭索斯法［9］。模態(tài)分析提取了前6階模態(tài)，其固有頻率與模態(tài)振型如表1所述，部分振型如圖8～圖10所示。模態(tài)振型是結(jié)構(gòu)在特定頻率的振動激勵作用下產(chǎn)生的變形形式，為了使結(jié)構(gòu)的變形形式更加直觀，選用5倍比例因子將其放大。

表1 磁懸浮儲能飛輪模態(tài)頻率與振型

計算分析結(jié)果表明，磁懸浮儲能飛輪在低頻段的振動模態(tài)主要為轉(zhuǎn)子的剛性偏轉(zhuǎn)，在高頻段的振動模態(tài)主要為定子彎曲與軸向上下振動、轉(zhuǎn)子彎曲振動。從表1可以看出，第1階與第2階的模態(tài)振型都是轉(zhuǎn)子相對定子偏轉(zhuǎn)，但是偏轉(zhuǎn)方向不同，第3階與第4階模態(tài)頻率接近且振型都是定子彎曲，但是彎曲的方向不同。低頻振動對系統(tǒng)的控制精度影響較大，應(yīng)采用模態(tài)控制方法加以控制。系統(tǒng)在高頻段的模態(tài)為定子的彎曲振動，可通過結(jié)構(gòu)優(yōu)化設(shè)計加以抑制。系統(tǒng)的更高階模態(tài)對系統(tǒng)的影響較小，可以忽略不計。操作實驗機時，應(yīng)盡量避開系統(tǒng)的固有頻率，以避免共振對系統(tǒng)性能的影響。

圖8 轉(zhuǎn)子第2階振型

圖9 轉(zhuǎn)子第6階振型

圖10 定子第4階振型

5 實驗?zāi)B(tài)分析

按照實驗?zāi)B(tài)分析原理，首先根據(jù)儲能飛輪的實際結(jié)構(gòu)，確定激勵點與測試點，建立測量模型，在激勵點施加激勵并測量振動響應(yīng)，將時間響應(yīng)轉(zhuǎn)換到頻域或直接利用時間響應(yīng)進行模態(tài)參數(shù)估計。一般來說，一個系統(tǒng)的動態(tài)響應(yīng)是它的若干模態(tài)振型的疊加。在一個給定的頻帶范圍內(nèi)只有一個模態(tài)是重要的，其模態(tài)參數(shù)可以單獨確定。采用單自由度法對系統(tǒng)進行峰值檢測與模態(tài)阻尼計算［10］。識別出來的固有頻率與模態(tài)阻尼如表2所述。

表2 實驗?zāi)B(tài)頻率與模態(tài)阻尼

6結(jié) 論

1）本文運用HyperMesh軟件對各個零件分別建模，然后利用聯(lián)接單元bolt、springs裝配。相比定子一體化建模方式，模型真實合理，結(jié)果更加準(zhǔn)確可靠。

2）由模態(tài)分析得到了系統(tǒng)前6階模態(tài)的固有頻率與振型，找出系統(tǒng)設(shè)計的薄弱環(huán)節(jié)（磁軸承安裝位置），這些部位最容易因振動而產(chǎn)生碰撞破環(huán)。

3）對比表1與表2的前4階模態(tài)頻率，可以看出實驗?zāi)B(tài)頻率值比有限元分析的結(jié)果略低一點，這是因為采用有限元方法進行模態(tài)分析時往往忽略系統(tǒng)阻尼的影響，而實際系統(tǒng)中存在阻尼，導(dǎo)致實驗?zāi)B(tài)頻率較有限元分析結(jié)果低。

［1］ 湯雙清.飛輪儲能技術(shù)及應(yīng)用［M］.武漢：華中科技大學(xué)出版社，2007：1-18.

［2］ Schweitzer G，Maslen E H.Magnetic bearings∶Theory,Design，and Application to Rotating Machinery［M］.Berlin：Springerverlag Berlin Heidelberg，2009∶1-195.

［3］ 孔繁鑫，黃勤.高能飛輪轉(zhuǎn)子動態(tài)特性及應(yīng)力分析［J］.機械設(shè)計與制造，2014（5）：28-30.

［4］ Wen Shaobo.Analysis of maximum radial stress location of composite energy storage flywheel rotor［J］.Archive of Applied Mechanics，2014，84（7）∶1007-1010.

［5］ 郭昌進，楊喜，王金麗.基于ANSYSWorkbench的甘蔗葉粉碎機機架模態(tài)分析［J］.農(nóng)機化研究，2014（8）：23-26.

［6］ 王鈺棟，金磊，洪清泉.HyperMesh&HyperView應(yīng)用技巧與高級實例［M］.北京：機械工業(yè)出版社，2012：257-408.

［7］ Yonnet J P.Permanent Magnet Bearings and Couplings［J］.IEEE Transactions on Magnetics，1981，17（1）∶1168-1173.

［8］ 小颯工作室.最新經(jīng)典ANSYS及Workbench教程［M］.北京：電子工業(yè)出版社，2004：378-395.

［9］ 歐賀國，方獻軍，洪清泉.Radioss理論基礎(chǔ)與工程應(yīng)用［M］.北京：機械工業(yè)出版社，2013：121-123.

［10］ Heylen W，Lammens S，et al.Modal Analysis Theory and Testing［D］.Katholieke Universteit Leuven，Departement Werktuigkunde，2006：50-72.