基于Monte Carlo法和積分法的被動永磁軸承磁力計算

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

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

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

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

1 被動永磁軸承磁力數(shù)學模型

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

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

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

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的相對位置長度；σ1，σ2為磁荷的面密度；μ0為真空磁導率。

軸承徑向磁力為

(4)

2 四重積分數(shù)值求解方法

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

(5)

(6)

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

(7)

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

(8)

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

考慮如下形式的n重積分

(9)

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

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

則(9)式可改寫為

(10)

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

(11)

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

式中：VD為D的體積。此時

(12)

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

(13)

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

4 堆疊型永磁軸承模塊化計算法

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

圖2 4種常見永磁軸承結構

5 算例分析

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

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

5.1 for循環(huán)計算結果分析

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

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

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

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

5.2 Monte Carlo計算結果分析

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

圖5 e=3 mm處Monte Carlo法計算結果

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

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

圖6 偏心距對軸承徑向磁力的影響

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

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

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

5.4 堆疊型永磁軸承磁力計算

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

圖7 堆疊型徑向永磁軸承

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

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

6 結論

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

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

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