基于Monte Carlo法和積分法的被動(dòng)永磁軸承磁力計(jì)算

李維，程文杰，肖玲，鐘斌，李明

(西安科技大學(xué) a.理學(xué)院；b.機(jī)械學(xué)院，西安 710054)

被動(dòng)永磁軸承利用永磁體環(huán)產(chǎn)生磁力和剛度，實(shí)現(xiàn)轉(zhuǎn)子懸浮，不需要外界輸入功率進(jìn)行控制，相對(duì)于主動(dòng)電磁軸承，被動(dòng)永磁軸承結(jié)構(gòu)簡(jiǎn)單，具有很高的魯棒性[1]。根據(jù)Earnshaw規(guī)則，被動(dòng)永磁軸承不可能自發(fā)地獲得一個(gè)穩(wěn)定的平衡點(diǎn)，即對(duì)于徑向被動(dòng)永磁軸承，需要對(duì)軸向位置進(jìn)行控制,而對(duì)于軸向被動(dòng)永磁軸承，則需對(duì)其徑向位置進(jìn)行控制[2]。另外，被動(dòng)永磁軸承的阻尼非常小，這種低阻尼、負(fù)剛度的特性使其受到的關(guān)注較少。然而，被動(dòng)永磁軸承的承載力和剛度可以通過采用堆疊結(jié)構(gòu)獲得較大提升，如今隨著高磁能積永磁材料的發(fā)展，其在工業(yè)領(lǐng)域獲得了越來越廣泛的應(yīng)用，比如高速壓縮機(jī)[3]、風(fēng)力發(fā)電機(jī)[4]、儲(chǔ)能飛輪[5]。文獻(xiàn)[6]采用徑向箔片軸承及軸向被動(dòng)永磁軸承結(jié)構(gòu)實(shí)現(xiàn)了轉(zhuǎn)子在40 000 r/min轉(zhuǎn)速下的完全穩(wěn)定懸浮。文獻(xiàn)[7]對(duì)被動(dòng)永磁軸承并聯(lián)油潤(rùn)滑滑動(dòng)軸承的穩(wěn)態(tài)承載力分配進(jìn)行了研究。最簡(jiǎn)單、典型的被動(dòng)永磁軸承是由2個(gè)永磁環(huán)構(gòu)成，其中一個(gè)固定在定子上，另一個(gè)隨轉(zhuǎn)子旋轉(zhuǎn)。這種結(jié)構(gòu)的磁力一般有2種解析法，即偶極子法和等效磁荷法(庫(kù)侖模型)，偶極子法適用于永磁體尺寸不大于氣隙尺寸的情形；等效磁荷法適用于永磁體尺寸大于氣隙尺寸的情形[8]。此外，可以利用磁路法求解這種典型被動(dòng)永磁軸承的磁場(chǎng)分布及剛度[9-10]。文獻(xiàn)[11]研究表明，對(duì)于軸向充磁的徑向永磁軸承，當(dāng)定子為180°的半環(huán)時(shí)承載力最大。采用等效磁荷法獲得的軸承磁力表達(dá)式涉及到四重積分，且難以獲得解析表達(dá)式,文獻(xiàn)[12]采用MATLAB軟件中的梯形積分法進(jìn)行求解,文獻(xiàn)[13]采用Monte Carlo法計(jì)算永磁軸承磁力,文獻(xiàn)[14]專門針對(duì)軸向充磁的徑向永磁軸承，給出了一種基于系數(shù)擬合的軸承磁力解析解，使設(shè)計(jì)更為簡(jiǎn)便，但對(duì)于其他類型的軸承不適用。文獻(xiàn)[12]研究了軸向磁化、徑向磁化和垂直磁化下的堆疊型被動(dòng)永磁軸承的剛度數(shù)學(xué)模型，涉及到大量積分、求和運(yùn)算;文獻(xiàn)[15]提出了堆疊型徑向永磁軸承的剛度、峰值磁力解析式，涉及到指數(shù)、求和運(yùn)算。

為了更簡(jiǎn)便、直觀地獲得四重積分的數(shù)值解，可以離散軸承計(jì)算區(qū)域，將積分轉(zhuǎn)換為求和運(yùn)算,劃分的單元數(shù)越多，計(jì)算結(jié)果越準(zhǔn)確，但計(jì)算速度將變慢。當(dāng)采用Monte Carlo法求解時(shí)，計(jì)算速度也會(huì)隨樣本點(diǎn)數(shù)目增多而下降，但計(jì)算結(jié)果更準(zhǔn)確。現(xiàn)以一種半定子環(huán)軸向充磁徑向永磁軸承為例，分別采用數(shù)值積分法(采用MATLAB軟件中的四重for循環(huán)實(shí)施)和Monte Carlo法進(jìn)行求解，對(duì)比2種方法的計(jì)算結(jié)果和時(shí)間。另外，現(xiàn)有文獻(xiàn)對(duì)堆疊型永磁軸承磁力的求解略顯復(fù)雜，不利于工程中應(yīng)用,為此，將軸承模塊化，提出一種方便、準(zhǔn)確的堆疊型永磁軸承磁力計(jì)算方法。

1 被動(dòng)永磁軸承磁力數(shù)學(xué)模型

半定子環(huán)軸向充磁徑向永磁軸承實(shí)際上是軸向充磁徑向永磁軸承的一種特殊情況，其力學(xué)模型如圖1所示。圖中，O1和O2分別為定子和轉(zhuǎn)子坐標(biāo)系原點(diǎn)，且yz平面均建立在定子左端面2上；R1，R2分別為永磁體環(huán)的內(nèi)、外半徑；R3，R4分別為定子環(huán)的內(nèi)、外半徑；l為磁環(huán)寬度；e為O2相對(duì)于O1的偏心距；x0為內(nèi)磁環(huán)相對(duì)于外磁環(huán)的軸向位移；r2為P相對(duì)于O1的位矢長(zhǎng)度；r3為Q相對(duì)于O2的位矢長(zhǎng)度；r0為Q相對(duì)于P的位矢長(zhǎng)度；α為Q在O2yz坐標(biāo)系中的角度坐標(biāo)；β為P在O1yz坐標(biāo)系中的角度坐標(biāo)。根據(jù)其磁化方向，磁荷分布于外磁環(huán)環(huán)面1，2和內(nèi)磁環(huán)環(huán)面3，4上，定子環(huán)面上磁荷微元P與轉(zhuǎn)子環(huán)面上的磁荷微元Q之間相互作用形成磁力。

圖1 半定子環(huán)軸向充磁徑向永磁軸承力學(xué)模型

根據(jù)等效磁荷理論，內(nèi)、外磁環(huán)間磁力分別在3個(gè)坐標(biāo)軸方向上的分力為

r2r3dr2dr3dαdβ，

(1)

(2)

(r3sinα-r2sinβ)2]1/2，

r24=[(l+x0)2+(e+r3cosα-r2cosβ)2+

(r3sinα-r2sinβ)2]1/2，

r13=[(l-x0)2+(e+r3cosα-r2cosβ)2+

(r3sinα-r2sinβ)2]1/2，

式中：Br為永磁體剩余磁通密度；r23,r24,r13為磁荷微元P與磁荷微元Q的相對(duì)位置長(zhǎng)度；σ1，σ2為磁荷的面密度；μ0為真空磁導(dǎo)率。

軸承徑向磁力為

(4)

2 四重積分?jǐn)?shù)值求解方法

積分?jǐn)?shù)值求解是先將被積函數(shù)和積分微元離散化，再求和。(1)～(3)式均含有r23,r24,r13和r2r3dr2dr3dαdβ等變量，因此只需將上述變量進(jìn)行離散。將定子、轉(zhuǎn)子的單邊徑向方向分別分成n1，n2份，則

(5)

(6)

將轉(zhuǎn)子圓環(huán)分成C2份，將定子半圓環(huán)分成C1份，則

(7)

將(5)～(7)式代入(1)～(3)式可得

(8)

3 高維數(shù)值積分的Monte Carlo法

考慮如下形式的n重積分

(9)

式中：x=(x1,x2,…,xn)為n維空間中的點(diǎn)；D為n維積分區(qū)域。Monte Carlo法計(jì)算積分的步驟如下：

1)在D上構(gòu)造分布密度函數(shù)。取D上任一概率密度函數(shù)ρ(x)，滿足條件ρ(x)≠0，當(dāng)x∈D，f(x)≠0時(shí)，令該分布密度函數(shù)為

則(9)式可改寫為

(10)

(10)式表示I是隨機(jī)變量F(x)的期望。

(11)

為便于計(jì)算，概率密度函數(shù)ρ(x)通常選取D上的均勻分布，即

式中：VD為D的體積。此時(shí)

(12)

將(12)式應(yīng)用到軸承磁力計(jì)算所涉及到的四重積分中，有

(13)

式中：ai，bi分別為積分下限和上限。

4 堆疊型永磁軸承模塊化計(jì)算法

4種典型的被動(dòng)永磁軸承結(jié)構(gòu)如圖2所示。其中，斥力型徑向軸承徑向剛度為正，軸向剛度為負(fù)；吸力型軸向軸承徑向剛度為正，軸向剛度為負(fù)[10]。前2種軸承結(jié)構(gòu)采用軸向充磁；后2種軸承采用徑向充磁。堆疊型永磁軸承可以視為由這4種基本模塊組合而成，由于磁力是矢量，滿足線性疊加原理，因此堆疊型軸承的磁力是各模塊軸承磁力之和。

圖2 4種常見永磁軸承結(jié)構(gòu)

5 算例分析

考慮定子環(huán)為全環(huán)和半環(huán)的2種情形，其中全環(huán)模型尺寸參考文獻(xiàn)[7]446，2種模型的主要參數(shù)見表1。

表1 軸向充磁徑向永磁軸承主要參數(shù)

5.1 for循環(huán)計(jì)算結(jié)果分析

使用for循環(huán)方法進(jìn)行高維數(shù)值積分計(jì)算，剖分參數(shù)見表2。對(duì)于本算例，當(dāng)C1=C2=120，n1=n2=14時(shí)，全定子環(huán)模型豎直方向磁力Fy的計(jì)算結(jié)果如圖3所示。由圖可知，計(jì)算結(jié)果與文獻(xiàn)[7]的結(jié)果非常吻合，證明了此方法的正確性。C1，C2，n1，n2取值越大，計(jì)算結(jié)果越精確，但結(jié)果之間的差別非常小。

表2 for循環(huán)方法進(jìn)行計(jì)算的剖分參數(shù)

對(duì)于半定子環(huán)模型，豎直方向偏心距e=3 mm時(shí)，計(jì)算結(jié)果如圖4所示。由圖可知，當(dāng)剖分單元數(shù)目小于(28×180)2時(shí)，剖分單元數(shù)對(duì)軸承徑向磁力影響顯著，且隨著剖分?jǐn)?shù)目的增加，軸承徑向磁力逐漸增大；當(dāng)剖分單元數(shù)目大于(28×180)2時(shí)，軸承徑向磁力雖仍隨剖分?jǐn)?shù)目的增加而增大，但增幅變小，軸承徑向磁力趨于穩(wěn)定。另外，剖分單元數(shù)的增加直接影響程序運(yùn)行效率，當(dāng)剖分單元數(shù)大于(28×180)2時(shí)，剖分單元數(shù)對(duì)軸承徑向磁力的影響不顯著，但程序運(yùn)行耗時(shí)t(由MATLAB軟件自動(dòng)記錄)隨單元數(shù)的增加而急劇增長(zhǎng)。由此可見，應(yīng)選擇合適的剖分單元數(shù)來提高計(jì)算效率。

圖4 e=3 mm處for循環(huán)計(jì)算結(jié)果

5.2 Monte Carlo計(jì)算結(jié)果分析

Monte Carlo法是一種近似方法，在使用過程中涉及隨機(jī)變量的產(chǎn)生，且對(duì)于不同隨機(jī)變量個(gè)數(shù)N，其計(jì)算的積分值一般不同。為確保Monte Carlo法的效果，降低偶然性的干擾，在不同隨機(jī)變量個(gè)數(shù)N(1.0×106，5.0×106，1.0×107)下分別運(yùn)行7次，然后取均值，獲得永磁軸承磁力。當(dāng)e=3 mm時(shí)，半定子環(huán)模型的計(jì)算結(jié)果如圖5所示。

圖5 e=3 mm處Monte Carlo法計(jì)算結(jié)果

由圖5a可知，N越大，軸承徑向磁力值波動(dòng)幅度越小，且軸承徑向磁力均值隨N的增大而增加。由圖5b可知，程序耗時(shí)隨N的增大而增加，當(dāng)N分別取1.0×106，1.0×107時(shí)，單次Monte Carlo法所需時(shí)間分別約為6，60 s。

不同樣本數(shù)下，軸承徑向磁力的均值隨偏心距的變化曲線如圖6所示。由圖可知，軸承徑向磁力隨偏心距的增加而增大，不同樣本數(shù)下的曲線非常接近。說明當(dāng)樣本數(shù)量達(dá)到一定值后(如N=1.0×106)，再增加樣本數(shù)不會(huì)導(dǎo)致計(jì)算精度發(fā)生大幅度變化，因此應(yīng)當(dāng)選擇合適的樣本數(shù)來提高計(jì)算效率。

圖6 偏心距對(duì)軸承徑向磁力的影響

5.3 for循環(huán)方法和Monte Carlo法的比較

e=3 mm處，通過Monte Carlo法計(jì)算不同隨機(jī)變量個(gè)數(shù)下運(yùn)行多次的平均磁力分別為32.20，32.36，32.40 kN,與for循環(huán)法計(jì)算的磁力(圖4)比較可知，2種方法的計(jì)算結(jié)果吻合良好，隨著樣本數(shù)的增大，Monte Carlo法計(jì)算結(jié)果越來越逼近for循環(huán)方法的極限值。由圖4可知，剖分單元數(shù)為(28×180)2時(shí)，for循環(huán)的計(jì)算精度可滿足要求，但其運(yùn)行耗時(shí)在170 s左右，不僅遠(yuǎn)高于圖5b中N=1.0×106的Monte Carlo法所需時(shí)間，甚至是N=1.0×107的Monte Carlo法所需時(shí)間的2.5倍，因此，與Monte Carlo法相比，for循環(huán)精度的提升需要耗費(fèi)更多的時(shí)間。對(duì)于本算例，當(dāng)N=5.0×106時(shí)，Monte Carlo法計(jì)算結(jié)果與for循環(huán)計(jì)算結(jié)果相差在0.185%之內(nèi)，但其耗時(shí)是后者的10%，說明Monte Carlo法具有更高的計(jì)算效率。

需要指出的是，采用四重for循環(huán)計(jì)算時(shí)，為了獲得更精確結(jié)果，需要將網(wǎng)格劃分更密，從而使得計(jì)算時(shí)間呈幾何級(jí)數(shù)增長(zhǎng)；而Monte Carlo法采用打點(diǎn)法，當(dāng)樣本點(diǎn)數(shù)增大時(shí)，計(jì)算時(shí)間呈線性增大,這使得Monte Carlo法在初次試探計(jì)算時(shí)具有更大的時(shí)間優(yōu)勢(shì)。另外，選擇的樣本點(diǎn)數(shù)N=5.0×106可以作為其他模型初次計(jì)算的參考，并不具有普適性，但樣本點(diǎn)數(shù)的選擇應(yīng)以多次單次計(jì)算結(jié)果不出現(xiàn)大幅波動(dòng)為宜。

5.4 堆疊型永磁軸承磁力計(jì)算

某堆疊型徑向永磁軸承如圖7所示，其中R1=12 mm，R2=23.5 mm，R3=29.5 mm，R4=41 mm；定轉(zhuǎn)子永磁環(huán)的厚度、寬度均相等，分別為d=11.5 mm，λ=2d；最大偏心距e=2 mm。定子環(huán)編號(hào)分別為A～D，轉(zhuǎn)子環(huán)編號(hào)分別為a～d。以轉(zhuǎn)子環(huán)a為例，aA和aC組成圖2a所示徑向軸承(基本模塊)，aB和aD組成圖2b所示的軸向軸承(基本模塊)，其他定轉(zhuǎn)子環(huán)的組合與此類似，所有組合情況見表3。根據(jù)模塊理論，軸承磁力為表3中所有模塊軸承磁力之和。實(shí)際上，只需計(jì)算4種組合的磁力，即aA，aB，aC，aD，根據(jù)結(jié)構(gòu)的對(duì)稱性，其他組合都屬于以上4種情形之一。這4種組合的磁力可以通過Monte Carlo法快速計(jì)算出來。

圖7 堆疊型徑向永磁軸承

表3 模塊軸承匹配表及磁力值

根據(jù)表3，圖7中軸承最大徑向磁力Fymax=32.63×4+12.01×6+4.47×4+0.55×2=221 N。根據(jù)文獻(xiàn)[15]737中的圖表數(shù)據(jù)估算圖7所示的軸承最大徑向磁力為207.2 N，兩者相對(duì)誤差僅6.2%，證明了所提出方法的正確性。

6 結(jié)論

1)將被動(dòng)永磁軸承磁力的四重積分離散化，利用MATLAB軟件中的四重for循環(huán)進(jìn)行求和運(yùn)算，求解結(jié)果與Monte Carlo法計(jì)算結(jié)果吻合。

2)與Monte Carlo法相比，for循環(huán)精度的提升需要耗費(fèi)更多的時(shí)間。對(duì)于文中算例，當(dāng)N=5.0×106時(shí)，Monte Carlo法計(jì)算結(jié)果與for循環(huán)計(jì)算結(jié)果(穩(wěn)定值)相差在0.185%之內(nèi)，但Monte Carlo法具有更高的計(jì)算效率。

3)根據(jù)疊加原理，將堆疊型軸承的磁力等效為各組合模塊軸承磁力之和，通過算例驗(yàn)證了所提方法的正確性，給出了一種適用于工程用的堆疊型軸承磁力的計(jì)算方法。