同軸圓環永磁體磁力計算與實驗對比研究

程書晗，高 琳，蘇宇鋒，王天航

（鄭州大學機械與動力工程學院，河南 鄭州 450001）

1 引言

隨著稀土永磁材料的發展，得益于磁力機械結構簡單、環保、非接觸、損耗小等優點，磁力在機械系統中有了越來越廣泛的應用。磁力齒輪［1］、磁力泵［2］和磁力彈簧［3-4］等磁力機械都是通過永磁鐵或電磁鐵的磁場耦合實現力或力矩的傳遞，通過無接觸傳動的方式可以有效解決傳統機械結構存在的問題。電磁式能量采集裝置［4-10］也得到了廣泛的研究，這種裝置一般都通過不同形狀、大小的永磁體布置使磁場相互耦合，在外界激勵的作用下產生相應的響應，利用法拉第電磁感應原理進行能量的采集。為了對磁力機械進行結構設計和響應分析需要引入磁力的計算，磁力的計算一般抽象為兩個永磁體之間力的計算。文獻［11］基于磁路分析，使用虛位移方法得出了軸向混合磁軸承的徑向承載力解析數學模型，得到了永磁體徑向承載力與徑向位移的關系。文獻［12-13］采用等效磁荷法模型對同軸或偏心時兩圓環磁體之間的磁力進行了建模，得到含有四重積分的磁力計算公式。文獻［14］通過磁偶極子之間的力模型計算兩磁體之間的磁力。文獻［15］對文獻［12］的計算公式進行簡化并增加校正項，得到相對精確的同軸圓柱磁體之間磁力計算結果。文獻［16-17］通過等效電流法對同軸圓柱磁體之間磁力進行求解，得到含有四重積分的解析表達式。

引入了文獻［18］的磁場求解方法，對其進行補充分析，使用等效電流法得到圓環永磁體含有廣義完全橢圓積分的磁場求解公式，分別使用自適應辛普森法則和文獻［19］中的數值解法進行Matlab 編程計算，得到圓環永磁體在空間中任意位置產生的磁場。進一步對同軸布置的兩個圓環磁體進行建模求解磁力，通過與COMSOL中仿真結果進行對比驗證模型的正確性，最后通過該模型在一種抗磁懸浮能量采集器［8］中的應用，在實驗中驗證這種磁力計算方法的準確性。

2 等效電流法進行磁場計算

2.1 圓環永磁體磁場計算

根據等效電流模型，均勻磁化的圓環磁體內部等效圓電流相互抵消，只剩余圓環邊界上的電流。圓環厚度為H，外圓半徑為R，內圓半徑為R’，磁場強度矢量為M，磁化方向為軸向，如圖1所示。

圖1 磁場求解示意圖Fig.1 Schematic Diagram of Magnetic Field Solution

把坐標系Oxyz建立在圓環磁體中心，圓環磁體等效為一個個邊界電流環，每個邊界電流環的寬度為dz'，電流環L的參數方程可以表示為：

設單位寬度上電流匝數為n，外圈每匝電流為I，所以外圈每個電流環的電流大小為nIdz'。內圈每匝電流為I'，有I'=-I，內圈每個電流環的電流大小為nI'dz'。設圓環頂面外圓面積為S，內圓面積為S'，則有：MH(S-S')=ISnH+I'S'nH即：

因此外圈、內圈電流環電流分別可以寫成Mdz'、-Mdz'。

對圓環外側面進行求解，P(x'，y'，z')為電流環上任意一點，則可用柱坐標表示向量OP為(Rcosφ'，Rsinφ'，z')。設待求解磁場點Q(x，y，z)，則根據畢奧-薩伐爾定律，首先對單個電流環每個位置進行積分，再對所有軸向上所有電流環累加磁場進行求解，得到如下等式：

式中：μ0—真空磁導率。

在電流環L上的P點位置，dL∕dφ'可以表示為：

為方便求解，調整坐標系的x軸位置，使待求點Q在xOz平面上，則向量OQ可用柱坐標表示為(ρ，0，z)。

將dL、OP和OQ代入方程（2）中，且將線積分轉換為對φ'的積分后可以得到：

由于電流環的對稱性，且點Q在xOz平面上，對于電流恒定的電流環來說不存在軸向的電流和位移電流，所以Q點在y軸方向的磁感應強度為0，即對于柱坐標來說，切向磁感應強度大小為0。因而可以舍掉等式（3）中的y方向磁場，可將等式（3）分解為徑向磁場Bρ(Q)和軸向磁場Bz(Q)。

將等式（4）和等式（5）對z'以及φ'積分，擴展到空間中任意一點可以得到：

廣義完全橢圓積分：

同理推出圓環內側面產生的磁場。得到圓環產生的總磁場Bρtot和Bztot為：

2.2 廣義完全橢圓積分求解

對于圓環永磁體所產生的磁場為式（14）和式（15），從解析式中可以看出，其中包含廣義完全橢圓積分項，想要求解該方程，需要對積分項進行求解。對于廣義完全橢圓積分，一般采用數值法進行求解，這里分別采用文獻［19］中的算法與自適應辛普森法則［20］進行求解。

自適應辛普森法則是一種對常規積分求解的數值解法，也可應用于橢圓積分。與之相比，文獻［19］提出的算法是一種基于拉登變換，利用高斯和拉格朗日等人提出的算術幾何平均的快速收斂尺度，只適用于求解橢圓積分的快速收斂算法。通過對兩種方法在Matlab中編程對比發現，對于同樣的求解精度，快速收斂算法求解速度遠快于自適應辛普森法，為提升后續磁力求解速度與求解精度，采用快速收斂算法求解磁場。

3 磁力計算推導

根據圓環磁體的磁場計算式（14）與式（15），進行兩同軸圓環磁體之間磁力的推導。

如圖2所示，兩個圓環磁體同軸，磁體1外圈半徑為R1，內圈半徑為R1'，厚度為H1，磁場強度矢量為M1；磁體2 外圈半徑為R2，內圈半徑為R2'，厚度為H2，磁場強度矢量為M2；兩磁體磁化方向均為軸向，且磁體間距為d。待求磁體2對磁體1的磁力。

圖2 圓環磁體示意圖Fig.2 Schematic Diagram of Ring Magnet

由于對于永磁體的參數已知的經常是磁體的剩余磁通密度Br，為簡化計算，當不考慮兩磁體之間磁場相互作用的對磁化的影響時，可以通過M=Br∕u0推導出磁體的磁化強度M。

由于磁場具有疊加性，磁體1所受的力為：

式中：Fij—磁體i所產生的磁場對磁體j的力。易得F11=0，所以磁體1所受的力為磁體2產生的磁場對其上等效電流的力。

將坐標軸建立在磁體2中心處將磁體1軸向劃分為一個個邊界電流環，對應的外環與內環分別為L1、L'1，根據安培力公式，單個電流環受磁力為：

式中：n1—磁體1單位長度電流環匝數；dz1—電流環寬度；I1—每匝外圈電流大小；I2—每匝內圈電流大小；B2—磁體2在磁體1處磁感應強度。

由等式（6）和等式（7）推導出磁體2在磁體1側面上各位置磁感應強度：

電流環在軸向和徑向磁場中的位置示意圖，如圖3所示。

圖3 電流環受力示意圖Fig.3 Schematic Diagram of Current Loop Force

由圖3（a）可知，徑向磁場對電流環產生的安培力合力方向垂直電流環平面。由圖3（b）可知，軸向磁場所產生的安培力方向都為徑向，且由于電流環的中心對稱性，所有軸向磁場所產生的安培力合力為0。所以將等式（18）帶入等式（17），取z軸方向為正方向，并將軸向磁場舍去，對圓環積分后可以得到：

將M1=n1I1=-n1帶入等式（19），對z1進行積分后可以獲得磁體2對磁體1的磁力F：

等式（20）是一個含有二重積分的解析式，需要使用數值法進行求解，內部的廣義橢圓積分可以使用快速收斂法進行求解，外部對z1的積分可以使用自適應辛普森法則進行求解，這樣只相當于一重積分的求解，與含有四重積分的磁力求解公式相比，速度大有提高。采用Matlab數值計算編程語言利用等式（20）可以對同軸圓環磁體間的相互作用力計算求解。

4 COMSOL磁力仿真結果對比研究

本節分別在有限元軟件COMSOL中進行磁力仿真分析，與Matlab求解出來的磁力結果進行比較。

4.1 COMSOL中磁場及磁力計算方法

在COMSOL中進行永磁體之間磁力的仿真需要使用磁場、無電流接口。在COMSOL軟件中內置的磁場求解方法為等效磁荷模型與有限元相結合。通過

式中：Vm—磁標勢；Br—剩余磁化強度；μrec—回復磁導率。

COMSOL中對磁力的求解采用Maxwell應力張量法。

式中：n1—表面的法向量；H—磁場強度；T2—空氣的應力張量。通過磁體所有邊界表面應力張量的積分可以求得磁力為：

式中：S—磁體邊界表面。

表1 圓柱磁體組參數Tab.1 Cylindrical Magnet Group Parameters

表2 圓環磁體組參數Tab.2 Ring Magnet Group Parameters

使用COMSOL進行仿真時，設置各組圓環磁體的尺寸，為保證求解精度，空氣域半徑統一設置為55mm，無限元域半徑設置為25mm，為保證網格質量，對圓環磁體進行結構化網格劃分，根據各組具體尺寸調節結構化網格最大單元大小，保證收斂精度以及平衡運算時間，其他部分采用自由四面體網格，網格大小采用超細化，網格劃分后模型，如圖4所示。穩態求解器的相對容差設置為1e-6。

圖4 COMSOL仿真模型Fig.4 COMSOL Simulation Model

由于使用COMSOL進行磁力仿真時，會有迭代誤差和空氣域和無限元域網格劃分所產生的誤差，使用文獻［21］中的研究方法進行校正，通過計算含誤差力與誤差力之差，得到兩磁體之間真實作用力，增加仿真結果的準確度。

4.2 兩種計算方法的對比

經過仿真與計算后結果對比，如圖5所示。可以看出由等效電流法求解出的兩圓環永磁體之間的磁力與使用COMSOL仿真得到的磁力結果非常相近。使用等式（20）擬合COMSOL中求解的離散點，六組數據得到的R2值分別為：0.9958，0.9999，0.9998，0.99997，0.9991，0.9996。可以看出擬合效果非常好，關于圓環磁體之間的磁力求解，使用等式（20）完全可以代替COMSOL軟件建模仿真。在求解速度上，同樣精度的情況下，使用COMSOL求解一個位置磁力的平均時間約2min，使用Matlab 模型進行求解的平均時間小于1ms，極大地提高了計算效率。

圖5 各組磁體不同間距下磁力對比圖Fig.5 Comparison of the Magnetic Force of Each Group of Magnets at Different Distances

5 實驗驗證

為了對磁力求解公式進行進一步驗證，對一種抗磁懸浮振動能量采集結構［22］進行分析，該結構主要由一個圓環提升磁鐵、一個圓柱懸浮磁鐵以及一對抗磁體構成，通過提升磁體對懸浮磁鐵的吸引力平衡懸浮磁鐵的重力判斷懸浮間距，將該懸浮間距求解問題抽象為磁力求解模型，分別進行模型求解和實驗求解的對比。實驗平臺，如圖6所示。實驗采用15組數據進行計算，為保證實驗數據的完備性，實驗中包含提升磁體和懸浮磁體分別為圓柱或圓環的四種情況。實驗使用的磁體具體尺寸，如表3所示。

表3 磁力驗證實驗中磁體參數Tab.3 Magnet Parameters in the Magnetic Verification Experiment

圖6 實驗平臺Fig.6 Experiment Platform

磁鐵的材料為不同牌號的燒結釹鐵硼，密度為7.5g∕cm2，提升磁體與懸浮磁體磁化方向相同且為軸向。尋找間距采用二分法，結合前面磁力計算模型自動尋找懸浮間距，最終懸浮間距的實驗測量結果和模型計算結果及誤差，如圖7所示。

圖7 磁體懸浮間距對比Fig.7 Comparison of Magnet Suspension Spacing

通過對比圖可以看出，實驗結果和模型預測結果之間有一定誤差，經計算十五組數據使用公式求解出的懸浮間距與實驗得到的結果誤差都小于5%，由于實際磁體存在磁化強度不準確，磁化不均勻，實驗磁體布置不對中等問題，該計算結果在可接受范圍之內。

6 結論

首先通過等效電流法求解同軸圓環永磁體之間的磁力，得到含有廣義完全橢圓積分的解析式，通過Matlab編程對橢圓積分部分采用快速數值解法進行求解。其次通過解析方程對不同參數下通過有限元仿真軟件獲得的仿真結果進行擬合，得到解析方程的R2值均大于0.995，驗證了解析式的準確性，并對比兩種求解方式的求解速度，分析結果表明此磁力計算模型速度遠優于有限元軟件仿真計算速度。最后將磁力求解模型應用于一種抗磁懸浮能量采集裝置中，用于求解提升磁體與懸浮磁體之間的懸浮間隙，通過建立15組實驗的實驗測量結果與模型預測結果進行對比，得到的模型計算誤差均小于5%，進一步驗證了同軸磁體間相互作用力的有效便捷計算方法。這種同軸圓環磁體磁力計算方法，由于求解速度快，精度較高的特點，可用于基于磁力作用機理的驅動器和傳感器分析與設計。