基于等效磁荷法用蒙特卡洛法計(jì)算永磁軸承磁力

張鋼，張堅(jiān)，張海龍，孟慶濤，樊曼

(上海大學(xué) 機(jī)電工程與自動(dòng)化學(xué)院，上海 200072)

永磁懸浮軸承(以下簡稱永磁軸承)是利用永磁體在氣隙中產(chǎn)生的磁場實(shí)現(xiàn)轉(zhuǎn)軸穩(wěn)定懸浮于空間的一種新型高技術(shù)軸承。由于其具有壽命長、無磨損、無需潤滑等顯著優(yōu)點(diǎn)，因此得到了廣泛關(guān)注和研究[1-7]。承載能力和剛度計(jì)算是永磁軸承研究與設(shè)計(jì)過程中的重要環(huán)節(jié)。而永磁軸承的磁力解析計(jì)算，目前常用的方法是利用等效磁荷法或分子電流法建立磁力模型，通過求解一個(gè)復(fù)雜的四重積分得到[4-7]。由于四重積分求解比較復(fù)雜，有學(xué)者直接用Matlab中的符號積分函數(shù)求解，不僅計(jì)算周期長，而且往往得不到正確的數(shù)值解；也有學(xué)者用C++進(jìn)行編程求解，但是編程過程較為復(fù)雜[7]。

基于以上原因，提出應(yīng)用蒙特卡洛法對基于等效磁荷法建立的永磁軸承磁力模型的四重積分方程組進(jìn)行求解；并與ANSYS分析得到的磁力值進(jìn)行比對。

1 永磁軸承磁力模型

永磁軸承的基本結(jié)構(gòu)如圖1所示，其中A為動(dòng)磁環(huán)，B為靜磁環(huán)，c0為兩磁環(huán)間的初始?xì)庀丁．?dāng)與轉(zhuǎn)子固定的動(dòng)磁環(huán)受到z軸正方向外力Fz作用時(shí)，將在沿z軸正方向上產(chǎn)生一軸向位移z(0

圖1 永磁軸承結(jié)構(gòu)圖

圖2 動(dòng)磁環(huán)A受力示意圖

依據(jù)等效磁荷理論[1]，可建立面1和面3之間的磁力F13,z的解析式

dr1dr3dαdβ/{[(c0-z+2L)2+(r3cosβ-r1cosα)2+

(r3sinβ-r1sinα)2]3/2}；

(1)

同理可建立面2和面4、面1和面4以及面2和面3的磁力F24,z，F(xiàn)14,z，F(xiàn)23,z的解析式

dαdβ/{[(c0-z)2+ (r4cosβ-r2cosα)2+

(r4sinβ-r2sinα)2]3/2}，

(2)

dαdβ/{[(c0-z+L)2+ (r4cosβ-r1cosα)2+(r4sinβ-r1sinα)2]3/2}，

(3)

dr2dr3dαdβ/{[(c0-z+L)2+(r3cosβ-r2cosα)2+ (r3sinβ-r2sinα)2]3/2}，

(4)

式中：Br為磁環(huán)剩磁(假設(shè)兩磁環(huán)剩磁相同)；μ0為真空磁導(dǎo)率；(r1,α)，(r2,α)，(r3,β),(r4,β)分別為面1，2，3，4上任一點(diǎn)的極坐標(biāo)。

規(guī)定磁力正方向?yàn)樽鴺?biāo)軸正方向，則軸向外載荷Fz為

Fz=F13,z+F24,z-F14,z-F23,z。

(5)

永磁軸承的剛度是指動(dòng)磁環(huán)在受到外擾力時(shí)，抵抗偏離平衡位置的能力。因此可以定義永磁軸承軸向剛度Kz為

(6)

2 蒙特卡洛法對磁力的求解

蒙特卡洛法求解多重積分的思想為任意一個(gè)積分都可看作是某個(gè)隨機(jī)變量的數(shù)學(xué)期望[8]。因此，在利用蒙特卡洛法計(jì)算多重積分時(shí)，采用隨機(jī)變量的算術(shù)平均值作為其近似值。

2.1 基本蒙特卡洛積分法

設(shè)D為n維空間Rn的一個(gè)區(qū)域，f(x)∈D?Rn→R，區(qū)域D上的n重積分用下式表示

(7)

I可以被認(rèn)為等于區(qū)域D的測度乘以函數(shù)f(p)的期望。計(jì)算時(shí)首先要選擇概率密度函數(shù)g(p1,p2,…，pn)，并且在區(qū)域D內(nèi)滿足g(p1,p2,…，pn)>0。

則(7)式可改寫為

dp1dp2…dpn。

(8)

無論區(qū)域D的形狀如何，根據(jù)概率密度函數(shù)g(p1,p2,…,pn)進(jìn)行抽樣時(shí)(抽樣次數(shù)為N)，均可以給出積分的估計(jì)為

(9)

標(biāo)準(zhǔn)差為

(10)

(11)

積分誤差為

(12)

式中：Xα為與置信度1-α有關(guān)的參數(shù)，當(dāng)α已知，可以查到相應(yīng)的Xα；σ為標(biāo)準(zhǔn)差。

基本的蒙特卡洛法就是找一個(gè)包含區(qū)域D的超立方體(測度已知，為Mc)，在D內(nèi)隨機(jī)生成N(N一般足夠大)個(gè)均勻分布的點(diǎn)，統(tǒng)計(jì)落入?yún)^(qū)域D的點(diǎn)，假設(shè)有m個(gè)，則區(qū)域D的測度為

(13)

函數(shù)f的期望為

(14)

從而有

(15)

將(15)式應(yīng)用到四重積分，可得

(16)

式中：bj為四重積分的4個(gè)積分上限；aj為對應(yīng)的4個(gè)積分下限；j=1,2,3,4。

2.2 磁力的求解流程

如圖3所示，按照圖示流程即可將(16)式用程序?qū)崿F(xiàn)，求解磁力非常方便。

圖3 蒙特卡洛法計(jì)算磁力流程圖

采用Matlab進(jìn)行編程，生成超立方體內(nèi)的隨機(jī)數(shù)可以調(diào)用函數(shù)unifrnd，調(diào)用格式為[9]：unifrnd(積分下限，積分上限，1，N)。

值得注意的是，蒙特卡洛法是一種近似方法，是概率統(tǒng)計(jì)的方法，定義不同的N值，積分值不會(huì)相同。下面將通過定義不同的N值和多次求解磁力來探究N值對結(jié)果精度以及可重復(fù)性的影響。

3 永磁軸承磁力計(jì)算分析

對照圖1的軸承結(jié)構(gòu)，文中計(jì)算所采用的尺寸為R1=10 mm，R2=20 mm，R3=10 mm，R4=20 mm，L=15 mm，初始軸向氣隙c0=6 mm。定義永磁環(huán)性能參數(shù)：Br=1.231 T，Hc=917 530 A/m，則磁環(huán)的相對磁導(dǎo)率μr為

(17)

將上述參數(shù)代入到等效磁荷法公式中，編制程序進(jìn)行求解。

3.1 樣本數(shù)量N對磁力值的影響

其他參數(shù)均不變，改變N值得出的永磁軸承磁力曲線如圖4、圖5所示。

圖4 改變N值時(shí)磁力曲線變化示意圖

圖5 不同N值時(shí)z=4 mm處的磁力曲線(重復(fù)運(yùn)行程序6次)

從圖4可以看出，隨著N值的減小，磁力曲線波動(dòng)越來越大，曲線越來越不平滑，計(jì)算數(shù)值準(zhǔn)確性隨之大大降低；并且積分區(qū)域越大，蒙特卡洛法的計(jì)算精度越低，因此，當(dāng)磁環(huán)面積較大時(shí)，應(yīng)該提高抽樣數(shù)量N。

從圖5可以看出，隨著N值的減小，數(shù)據(jù)的可重復(fù)性差；蒙特卡洛積分誤差與抽樣數(shù)的平方根成反比，在滿足工程應(yīng)用(誤差在10%以內(nèi))的前提下，N應(yīng)取一合適值(106)。

3.2 蒙特卡洛法與數(shù)值解法對比

研究表明，在計(jì)算永磁軸承的承載能力時(shí)，有限元數(shù)值解是較為準(zhǔn)確的[2-3，10]。用ANSYS求解時(shí)，若只考慮軸向氣隙的變化，采用軸對稱二維模型，單元類型選擇plane53，磁環(huán)網(wǎng)格尺寸為0.1 mm，空氣網(wǎng)格尺寸為2 mm，給最外層空氣邊界施加通量平行條件(即Az=0)，將動(dòng)磁環(huán)定義為一個(gè)組件，求解即可得到磁力線分布和軸向力的值。

圖6和圖7分別為蒙特卡洛法(N=1 000 000)與有限元數(shù)值解法的軸向力和軸向剛度計(jì)算結(jié)果對比。

圖6 軸向力與軸向位移關(guān)系曲線

由圖可以看出，兩種算法的計(jì)算結(jié)果接近，兩者的磁力值最大偏差對應(yīng)的軸向偏移在5.5 mm附近，誤差約為6.5%。在最常使用的軸向偏移4 mm處，偏差僅為5.3%。軸向剛度與軸向位移關(guān)系曲線的變化趨勢比ANSYS得出的更平緩。

綜合圖6和圖7中的對比結(jié)果可知，蒙特卡洛法完全可以應(yīng)用到永磁軸承的工程設(shè)計(jì)中。

圖7 軸向剛度與軸向位移關(guān)系曲線

4 結(jié)論

(1)通過蒙特卡洛法與有限元數(shù)值解的結(jié)果比較可知，采用蒙特卡洛法結(jié)合等效磁荷法能夠較準(zhǔn)確地計(jì)算出軸承的磁力以及剛度。

(2)在用蒙特卡洛法求解時(shí)，抽樣樣本N值越大，計(jì)算精度越高，但N值過大會(huì)增加計(jì)算周期，N值的選取應(yīng)該使得磁力反復(fù)求解的波動(dòng)不超過1%。