基于MATLAB的永磁體直線電機的設計

王雪艷

（淄博職業學院，淄博 255314）

0 引言

一臺旋轉運動的電動機沿著半徑方向的剖面圖與直線電動機的相同。其簡單的結構和運動方式來帶了高效率。磁懸浮列車就是用直線電機來驅動的最好的例子。航空母艦上的飛機電磁彈射器也可以使用直線電機（未運用，正處于研究階段）。目前直線電機的應用已經越來越廣泛。本文以永磁電機為例探討直線電機的工作原理和設計思路。隨著計算機技術的發展，電機設計方法的發展日新月異。目前最常用的是用商業的專用軟件來進行設計，但是要找到一種既有針對性又符合自己設計風格的軟件并不容易，自己編程來輔助設計勢在必行。程序需要根據畢奧—薩伐爾定律計算磁場，鑒于計算過程比較復雜，國內的研究成果僅有：蘭州大學茍曉凡等學者的《矩形永磁體磁場分布的解析表達式》（國家自然科學重點基金項目）。經本文研究后發現，借助于MATLAB工具強大的公式推導和數值運算能力，只需用簡單的程序語句便可計算出數據，然后繪制出直觀圖像，此法比傳統的方法高效簡便。MATLAB還提供了各種豐富的工具箱，比如神經網絡工具箱、優化工具箱等等，可以很方便地利用它們去進行更深層次的研究。

1 單塊矩形永磁體磁場計算

永磁直線電機是由通電線圈和若干塊矩形永磁體組成的。要研究直線電機的工作過程，需要先弄清電機磁場的空間分布。一個長、寬、高為a、b、h的矩形永磁體如圖1所示。

圖1 矩形永磁體

其外部空間任一點P（x, y, z）的磁場是由磁體內部分子電流共同疊加而產生的合磁場。當矩形磁體沿一個方向充分磁化并達到飽和后，因磁體內分子電流排列整齊而使磁體內部分子電流全部被抵銷，所以可看成只有磁體表面電流在流動。如圖1所示，矩形永磁體的A’ B’ C’ D’平面處于XY坐標平面，A’點與坐標原點距離為d，根據磁體NS極性，按右手螺旋定則可知其體表電流沿ABCD方向流動，并且只有四個豎直方向（平行Z軸）的平面有電流流動，而上下兩個水平平面（垂直Z 軸）無電流流動。P點磁場的計算只需計算出四個豎直平面電流在P點產生的磁場，然后加起來。

先來計算平面ABB’A’在P點形成的磁場，假設此平面上有一點（xi, yi, zi），可由畢奧—薩伐爾定律來計算（xi, yi, zi）點在P點形成的磁場強度，公式為，式中是沿電流方向的矢量微元。在MATLAB中，用向量來表示它。因為它只有x和y兩個方向的分量，可以記為＝[dxi，dyi，0]，式中r為P (x, y, z)點到(xi, yi,zi)的距離。我們可以將構造成一個MATLAB向量r_matrix，設：

則 r_matrix=[(x-xi)/rr，(y-yi)/rr，(z-zi)/rr]，這里需要引入一個面電流密度J＝I/h，運算式變換為：記為K，它是一個與磁體材料和尺寸等參數有關的常數。通過實測磁體外某一點的磁感應強度，計算出常數K，上面的公式可轉變為：，可編程處理：

clear;

syms x y z xi yi zi dxi dyi dzi K real

rr=((x-xi)^2+(y-yi)^2+(z-zi)^2)^(3/2);

dl=[dxi,dyi,0];

r_matrix=[(x-xi)/rr,(y-yi)/rr,(z-zi)/rr];

dB=K*dzi*cross(dl,r_matrix);

通過如上程序進行符號運算后得到的dB是一個三維向量，dB(1),dB(2),dB(3)三個分量分別代表x，y，z方向的磁感應強度微分元。使用quad2d（）函數對ABB’ A’ 整個平面區間進行數值積分得到B，但dB的向量符號表達式中包含有積分變量微元dxi、dyi、dzi，必須去掉后才能用quad2d（）函數進行數值積分，如何具體處理這三個積分變量微元呢？這里我們可以根據不同的平面來確定。

1）ABB’ A’ 平面，平行于ZY平面，積分元為dyi和dzi，xi=d，dxi=0，yi變量積分區間為[0，b]，zi變量積分區間為[0，h]。處理方法如下：把dB向量中的三個積分元做一個替換，dxi替換為0，dyi和dzi分別替換為1，相當于先得到dB/(dyidzi)的符號表達式，然后用quad2d（）函數進行積分。

2）BCC’ B’ 平面，平行于 ZX平面，積分元為dxi和dzi，yi=b，dyi=0，xi變量積分區間為[d，d+a]，zi變量積分區間為[0，h]。

3）CDD’ C’ 平面，平行于ZY平面，積分元為 dyi和 dzi，xi=d+a，dxi=0，yi變量積分區間為 [b，0]，zi變量積分區間為[0，h]。

4）DAA’ D’ 平面：平行于ZX平面，積分元為 dxi和 dzi，yi=0，dyi=0，xi變量積分區間為 [d+a，d]，zi變量積分區間為[0，h]，編程如下：

syms a b d h real

%ABB'A'平面 xi=d, dxi=0, dyi=1, dzi=1

dB_AB_div_dyidzi=subs(dB,{xi,dxi,dyi,dzi},{d,0,1,1});

%BCC'B'平面 yi=b, dxi=1, dyi=0, dzi=1

dB_BC_div_dxidzi=subs(dB,{yi,dxi,dyi,dzi},{b,1,0,1});

%CDD'C'平面 xi=d+a, dxi=0, dyi=1, dzi=1

dB_CD_div_dyidzi=subs(dB,{xi,dxi,dyi,dzi},{d+a,0,1,1});

%DAA'D'平面 yi=0, dxi=1, dyi=0, dzi=1

dB_DA_div_dxidzi=subs(dB,{yi,dxi,dyi,dzi},{0,1,0,1});

此后對四個面的表達式進行合并處理，YZ平面的總磁感應強度dByz由ABB'A'平面和CDD'C'平面相加得到，考慮兩平面電流方向相反，為了把yi的積分區間統一為0到b, 相加時把CDD’ C’段加負號，同樣XZ平面的總磁感應強度dBxz由BCC’ B’ 平面和DAA’ D’ 平面兩部分相加得到，把DAA’ D’ 段也反一下符號，使xi的積分區間統一為d到d+a。對dByz和dBxz分別進行數值積分后再相加，就得到P（x, y, z）點磁感應強度B。為了先求出常數K，測量出磁體表面或者外面一點的磁感應強度是必要的。比如測量出ABCD平面中心點外0.5mm處Z方向的磁感應強度Bz0為0.1T，然后令K=1，計算出Bz，則K=Bz0/Bz, 假設磁體的參數為：a=50mm，b=500mm，h=10mm。下面給出計算程序：

dByz=dB_AB_div_dyidzi-dB_CD_div_dyidzi

dBxz=dB_BC_div_dxidzi-dB_DA_div_dxidzi

a=50;b=500;h=10;d=0;Bz0=0.1;x=a/2;y=b/2;z=1 0.5;K=1;

dBz1=vectorize(char(subs(dByz(3))));%YZ平面在中心處的磁感應強度

dBz2=vectorize(char(subs(dBxz(3))));%XZ平面在中心處的磁感應強度

下面要對dBz1和dBz2進行積分。先用sprintf()函數構造出符號運算表達式，然后再放到eval()中去運算得到結果。

Bz1=sprintf('quad2d(@(yi,zi)%s,0,b,0,h)',dBz1);

Bz2=sprintf('quad2d(@(xi,zi)%s,0,b,0,h)',dBz2);

f=sprintf('%s+%s',Bz1,Bz2); %f為XZ和YZ兩個平面Z方向磁場總和

K=abs(Bz0/eval(f));

接下來，可以利用前面所推導出的符號表達式，進行諸多分析設計。比如要計算分析ABCD平面上在y=b/2,z=10.5處X方向和Z方向的磁場沿坐標X的分布圖，可用如下程序，繪制圖像如圖2所示。

Syms x real;

dBxzz1=vectorize(char(subs(dBxz(3))));

Bxzz1=sprintf('quad2d(@(xi, zi)%s,0,a,0,h)',dBxzz1);

dByzz1=vectorize(char(subs(dByz(3))));

Byzz1=sprintf('quad2d(@(yi,zi)%s,0,b,0,h)',dByzz1);

f=sprintf('%s+%s',Bxzz1,Byzz1);

fun=eval(['@(x)',f]);

xx=-20:70;

Bz1=arrayfun(@(x) fun(x),xx);

dBx1=vectorize(char(subs(dByz(1))));

Bx1=sprintf('quad2d(@(yi,zi)%s,0,b,0,h)',dBx1);

fun=eval(['@(x)',Bx1]);

Bx1=arrayfun(@(x) fun(x),xx);

plot(xx,Bx1,xx,Bz1)

圖2 Bx1、Bz1磁場分布曲線圖

通過計算得出磁場的空間分布圖。如要求出ABCD平面上方0.5mm的Bz分布，編程如下：

syms x y real;

dBxzz1=vectorize(char(subs(dBxz(3))));

Bxzz1=sprintf('quad2d(@(xi,zi)%s,0,a,0,h)',dBxzz1);

dByzz1=vectorize(char(subs(dByz(3))));

Byzz1=sprintf('quad2d(@(yi,zi)%s,0,b,0,h)',dByzz1);

f=sprintf('%s+%s',Bxzz1,Byzz1);

fun=eval(['@(x,y)',f]);

xx=0:50;yy=0:500;

[X,Y]=meshgrid(xx,yy);

Bz=arrayfun(@(x,y) fun(x,y),X,Y)

plot3(X,Y,Bz)

上面程序繪出的圖像如圖3所示。

圖3 磁體表面Bz三維分布圖

如圖3所示，在磁體表面附近的Bz較強，而遠離磁體表面的Bz迅速減小。因此，在制作電機時，應讓通電線圈盡量靠近磁體表面位置以獲得強大的磁場力。一般放在表面0.5mm上方。

2 多塊永磁體磁場中線圈受力計算

直線電機是由多塊尺寸相同的永磁體按一定間距并排，并且相鄰的兩塊極性相反而構成的。如圖4所示。

圖4 多塊磁體排列示意圖

假設各磁體之間距離為c，在磁體上方放置一個平行于磁體頂端面的通電線圈EFGH，線圈EF邊長b，FG邊長＝a+c，線圈距離磁體表面0.5mm，假設電流I沿線圈FGHE方向流動，現在來分析線圈沿X軸方向所受的磁場力Fx，磁場力可以根據進行計算。由于磁體的對稱性，沿X軸向位于磁體中心線兩側的B是對稱的，而GF邊和EH邊電流相反，所以兩側邊導線受力抵銷，只需計算導線GH和EF上的受力，并且只要計算出了EF段的磁場力，可通過坐標平移變換得出GH導線段的受力。下面就來對EF的受力進行計算：電流沿Y方向流動，＝I×[0，dy, 0]×[Bx, By, Bz]，可以推導出dFx= -I×Bz×dy，這里的Bz是多塊磁體在Z方向的磁場總和。根據此式，通過數值積分就能得到Fx，并畫出其在X方向的變化曲線圖。dy積分區間應該沿電流方向，在[0，b]區間。

現在以8塊磁鐵為例來編程說明，對于各磁體的磁場先分別計算，再累加起來。各磁體主要差別在X軸上放置的位置不一樣，就是前面求出的dB向量表達式中的d參數不同，另外還有偶數位置的磁體極性反向，也就是向量正負號要反向。第1塊磁鐵放在坐標原點處d=0,第2塊d=a+c，第3塊d=2 (a+c),第i塊磁鐵d=(i-1) (a+c)。寫出程序如下：

Syms c d x y real;

p=-1;dByzz=0;dBxzz=0;

for i=1:8;

p=-p;

dd=(i-1)*(a+c);

dByzz=dByzz+p*subs(dByz(3),d,dd);

dBxzz=dBxzz+p*subs(dBxz(3),xi,xi+dd)；%

此處對xi坐標變換，以便將積分區間由[d,d+a]變為[0, a]

end

上面這段程序得出的dBxzz和dByzz分別是8塊磁鐵在XZ平面和YZ平面Z方向的總磁感應強度的微分符號表達式。下面進一步計算Fx，程序仍然使用前面的磁體參數，假設I=2A。

a=50;b=500;h=10;d=0;z=10.5;c=5;

dBxzz=vectorize(char(subs(dBxzz)));

Bxzz=sprintf('quadl(@(yy) arrayfun(@(y)quad2d (@(xi,zi)%s,0,a,0,h),yy),0,b,0.1)',dBxzz);

dByzz=vectorize(char(subs(dByzz)));

Byzz=sprintf('quadl(@(yy) arrayfun(@(y)quad2d(@(yi,zi)%s,0,b,0,h),yy),0,b,0.1)',dByzz);

f=sprintf('%s+%s',Bxzz,Byzz)

fun=eval(['@(x)',f]);

xx=0:300;I=2;

Fx=-I*(arrayfun(@(x) fun(x),xx));

這里得出的Fx是一維向量，坐標變化單位為1mm，如圖5所示。

圖5 線圈EF段導線受力曲線圖

GH段受力計算通過坐標變換x=x+a+c,由于電流反向，-Fx(x+a+c)，所以線圈總的受力為F=Fx (x) -Fx(x+a+c)，繪制圖像如圖6所示。

3 結論

1） 如圖6所示，線圈的受力是不斷正負變化的。為保證線圈沿X軸方向持續運動，需在F變為負值時及時將線圈中的電流反向一次，這樣就能使線圈的受力始終為正值，大小等于| F |。現代電子技術的發展使得對電流的控制非常精確，制造直線電機也變得比以往更加容易。電機由旋轉運動回歸到直線運動，從而帶來更高的效率。

2）多個線圈的共同應用使受力變得平滑。只需要根據位置做坐標變換就可以對多個線圈進行分析，不用重復計算。另外，由于磁場的對稱性，只需要計算前四塊磁鐵的情況。

3） 本文通過改變磁體的參數分析磁場力與各參數之間的關系曲線。MATLAB編程的使用使得最優化設計方案的尋求更為方便。

圖6 線圈受力圖

4） MATLAB是一種解釋型語言。雖然運行效率不高，但可以通過優化、使用并行運算等方法來提高速度，還可以通過改變分析計算的計算精度、取樣點間距等來提高運行速度。本文程序在MATLAB R2009b上運行通過。

5）為了驗證磁場分布計算的正確性，將本文計算結果與美國希捷公司高級顧問工程師侯春洪博士設計軟件計算繪制的圖表進行對比，其一致性證明了本文研究方法和程序運算的真實性和可信度。

[1] 茍曉凡, 楊勇, 鄭曉靜. 矩形永磁體磁場分布的解析表達式[J]. 應用數學和力學, 2004,25(3).