周群益，莫云飛，周麗麗，侯兆陽

（1.廣州理工學院通識教育學院，廣東 廣州 510540；2.長沙學院電子信息與電氣工程學院，湖南 長沙 410022；3.贛南醫(yī)學院醫(yī)學信息工程學院，江西 贛州 341000；4.長安大學理學院，陜西 西安 710064）

亥姆霍茲線圈是由兩個半徑相同、電流大小和方向相同的共軸圓環(huán)組成，不少文獻研究了亥姆霍茲線圈均勻磁場區(qū)域的問題。環(huán)電流磁場的計算涉及兩類完全橢圓積分，早期的學者通過橢圓積分表手工計算兩類橢圓積分，效率低、精度差［1-4］。后來有學者利用電子計算機計算，效率和精度有所提高，還能畫出曲線圖［5-6］。近年來，有學者利用Mathematica計算兩類完全橢圓積分，并繪制曲線和曲面，效率和精度大大提高［7-11］。也有一些學者利用MATLAB計算兩類完全橢圓積分，并繪制曲線和曲面，效果也不錯［12-14］。有的學者設(shè)計了Mathematica和MATLAB的程序，繪制了均勻區(qū)域圖［15-18］。不過，這些學者沒有將公式無量綱化，畫出的圖形不夠規(guī)范。莫云飛等用MATLAB精確地繪制了亥姆霍茲線圈的均勻區(qū)域，但是還不夠系統(tǒng)和全面［19］。文章畫出了亥姆霍茲線圈在臨界條件下磁感應(yīng)強度分量和合磁場以及方向曲面，增加磁感應(yīng)線說明了兩環(huán)之間的距離對均勻區(qū)域大小和形狀的影響。文后附有兩個主程序和函數(shù)文件，其中主程序用參數(shù)調(diào)用函數(shù)文件，第一個函數(shù)文件說明了磁感應(yīng)強度分量與合磁場以及方向曲面的畫法，第二個函數(shù)文件說明了磁感應(yīng)線和均勻區(qū)域的畫法。這兩個函數(shù)文件使用了一定的編程技巧，有一定的參考價值。

1 亥姆霍茲線圈的磁場

如圖1所示，取x軸為縱軸，z軸為橫軸，兩環(huán)共z軸。設(shè)兩個圓環(huán)的半徑為a，通有同方向的電流I，兩環(huán)之間的距離為2L。

圖1 亥姆霍茲線圈

如果將左環(huán)或右環(huán)移到原點，環(huán)電流產(chǎn)生的磁感應(yīng)強度的兩個分量分別為［20-22］：

其中k稱為模數(shù)，其平方為：

K（k）、E（k）分別為第一類完全橢圓積分和第二類完全橢圓積分：

當x= 0時，必 有k= 0。而K（0）= E（0）= π/2。對式（1a）的第一個等式求x→0時的極限，可得：

對式（1b）的第一個等式求x→0時極限，可得：

這是軸線上的磁感應(yīng)強度，與電磁學教材中的公式相同。當z= 0時，可得：

這是環(huán)心處的磁感應(yīng)強度，用B0表示。

將式（1a）（1b）中的z換為z-L，即可得右環(huán)O1產(chǎn)生的磁場的兩個分量B1x（x，z）和B1z（x，z）；將式（1a）（1b）中的z換為z+L，即可得左環(huán)O2產(chǎn)生的磁場的兩個分量B2x（x，z）和B2z（x，z）。兩個線圈產(chǎn)生的合磁場的兩個分量在左右兩環(huán)產(chǎn)生的分量之和分別為Bx（x，z）和Bz（x，z），由此可以計算合磁場的大小和方向。

當z= 0時，可得原點的磁感應(yīng)強度：

這是一個鞍點，其磁感應(yīng)強度用BS表示。

在z= 0處，令，可得a2= 4L2，所以2L=a，這是臨界距離。z= 0處的臨界磁場為：

2 磁感應(yīng)強度的可視化

取半徑a為長度單位，則無量綱的坐標和長度為：

取B0=μ0I/2a為磁感應(yīng)強度單位，即可將磁感應(yīng)強度的公式無量綱化。注意，L*是可調(diào)節(jié)的參數(shù)，當L*= 1/2時，它就是臨界無量綱距離。

根據(jù)磁感應(yīng)強度的無量綱公式，利用ellipke函數(shù)可以計算橢圓積分，利用surf指令可以繪制曲面（見附錄1）［23］49-52。

（1）亥姆霍茲線圈的臨界磁感應(yīng)強度的z分量Bz的曲面如圖2所示，點（±a，±L）是奇點，表示環(huán)所在的位置，奇點處出現(xiàn)4個“背靠背的峰”和4個“面對面的谷”。

圖2 臨界磁感應(yīng)強度z分量的分布面

（2）磁感應(yīng)強度的x分量Bx的曲面如圖3所示，奇點處出現(xiàn)4個“面朝背的峰”和4個“面朝背的谷”。

圖3 臨界磁感應(yīng)強度x分量的分布面

（3）合磁場B的曲面如圖4所示，中心處比較平坦，說明此處接近于勻強磁場，圓環(huán)處是4個很高的峰，外面磁場比較弱。

圖4 臨界合磁場的分布面

（4）合磁場B的方向角α的曲面如圖5所示，在x= 0的z軸上，α= 0；在兩環(huán)之間，α隨著極角θ的增加變化較小，說明磁場比較均勻；在兩環(huán)之外，α隨著θ的增加而增加；當θ取某些值的時候，α發(fā)生從π到-π的躍變。

圖5 臨界合磁場方向的分布面

3 磁場均勻區(qū)域的可視化

利用三維等高線指令contour3可畫等值線，利用流線指令streamline可以繪制磁感應(yīng)線（見附錄2）［23］。

取|B–BS| ≤ 0.01BS為線圈的均勻范圍，即0.99BS≤B≤ 1.01BS。

（1）當L=LC= 0.5a時，亥姆霍茲線圈臨界截面和磁場的均勻區(qū)域以及磁感應(yīng)線如圖6所示。兩個電流從下面流入，從上面流出。磁感應(yīng)線在原點附近是平行的，因而是勻強磁場。磁場BS=BC= 1.431 1B0的等值線恰好經(jīng)過原點并包圍4個電流；B= 1.01BC的4條等值線分別包圍了4個電流，4條曲線并不連接；B= 0.99BC的等值線是一條曲線，同時包圍了4個電流，曲線如同“X”形狀。B= 1.01BC的等值線與B= 0.99BC的等值線之間的區(qū)域就是線圈臨界磁場的均勻范圍。均勻區(qū)域呈現(xiàn)“八爪”形狀，一對“爪”包圍一個電流。

圖6 臨界磁場的均勻范圍和磁感應(yīng)線（L=0.5a）

（2）當L=0.45a時，兩環(huán)距離比較近，線圈磁場的均勻范圍和磁感應(yīng)線如圖7所示。鞍點磁場增加為BS=1.516 7B0，磁場的均勻范圍左右收縮而上下拉長。B=BS的左右曲線經(jīng)過原點，B=1.01BS的等值線與B=0.99BS的等值線分別向上下移動而不經(jīng)過原點，均勻區(qū)域仍然呈現(xiàn)“八爪”形狀，范圍稍微縮小。

圖7 磁場的均勻范圍和磁感應(yīng)線（L=0.45a）

（3）當L=0.4a時，兩環(huán)距離更加小，磁場的均勻范圍和磁感應(yīng)線如圖8所示。鞍點磁場增加為BS=1.600 8B0，B=1.01BS的曲線合并成上下兩條曲線，分別包圍上面兩個電流和下面兩個電流，上下各有一個小小的均勻區(qū)域，中間均勻區(qū)域呈“四爪”形狀，上下一對“爪”分別包圍上下兩個電流。雖然鞍點的磁場增強了，但是均勻區(qū)域變得更小了。如果距離繼續(xù)變小，磁場均勻區(qū)域會變得更小而形狀不變（圖略）。

圖8 磁場的均勻范圍和磁感應(yīng)線（L=0.4a）

（4）當L=0.55a時，兩環(huán)距離比較大，線圈磁場的均勻范圍和磁感應(yīng)線如圖9所示。鞍點磁場減小為BS=1.3454B0，磁場的均勻范圍上下收縮而左右拉長。B=BS的上下曲線經(jīng)過原點，但是左右部分曲線分別向左右移動而不經(jīng)過原點，均勻區(qū)域仍然呈現(xiàn)“八爪”形狀，范圍稍微縮小。

圖9 磁場的均勻范圍和磁感應(yīng)線（L=0.55a）

（5）當L=0.6a時，兩環(huán)距離更大，磁場的均勻范圍和磁感應(yīng)線如圖10所示。鞍點磁場減小為BS=1.261B0，B=1.01BS的曲線合并成右邊兩條曲線，分別包圍左邊兩個電流和右邊兩個電流，左右各有一個小小的均勻區(qū)域，中間均勻區(qū)域呈“四爪”形狀，左右一對“爪”分別包圍左右兩個電流。鞍點的磁場減弱了，均勻區(qū)域也變得更小了。如果距離繼續(xù)變大，磁場均勻區(qū)域會變得更小而形狀不變（圖略）。

圖10 磁場的均勻范圍和磁感應(yīng)線（L= 0.6a）

4 結(jié)束語

文章說明了亥姆霍茲線圈在直角坐標系中磁感應(yīng)強度的解析式，將公式無量綱化，繪制出了亥姆霍茲線圈磁場的分量和合磁場的分布曲面。文章精確地繪制了亥姆霍茲線圈均勻磁場的分布范圍和磁感應(yīng)強度分布線，說明了兩環(huán)之間的距離對均勻范圍的影響。為了提高讀者應(yīng)用MATLAB解決問題的能力，文章附錄提供兩個完整的主程序和函數(shù)文件供讀者參考。

附錄1

ring2B1main.m

%亥姆霍茲線圈磁場的主程序

clear，l=0.5；%清除變量，距離

ring2B1fun（l）%調(diào)用函數(shù)文件

ring2B1fun（0.25）%調(diào)用函數(shù)文件

ring2B1fun（1）%調(diào)用函數(shù)文件

ring2B1fun.m

%亥姆霍茲線圈的磁場的函數(shù)文件

function fun（l）

rm=2；r=linspace（-rm，rm，40）；%最大坐標，坐標向量（繞過奇點）

［X，Z］=meshgrid（r）；R1=sqrt（（Z-l）.^2+（abs（X）+1）.^2）；%坐標矩陣，最大距離

KK1=4*abs（X）./R1.^2；［K1，E1］=ellipke（KK1）；%模數(shù)，兩個完全橢圓積分

Bx1=（（1+X.^2+（Z-l）.^2）./（（abs（X）-1）.^2+（Z-l）.^2）.*E1-K1）.*（Z-l）./X/pi./R1；%B的x分量

Bz1=（（1-X.^2-（Z-l）.^2）./（（abs（X）-1）.^2+（Z-l）.^2）.*E1+K1）/pi./R1；%B的z分量

R2=sqrt（（Z+l）.^2+（abs（X）+1）.^2）；%最大距離

KK2=4*abs（X）./R2.^2；［K2，E2］=ellipke（KK2）；%模數(shù)，兩個完全橢圓積分

Bx2=（（1+X.^2+（Z+l）.^2）./（（abs（X）-1）.^2+（Z+l）.^2）.*E2-K2）.*（Z+l）./X/pi./R2；%B的x分量

Bz2=（（1-X.^2-（Z+l）.^2）./（（abs（X）-1）.^2+（Z+l）.^2）.*E2+K2）/pi./R2；%B的z分量

Bx=Bx1+Bx2；Bz=Bz1+Bz2；%合磁場的x分量，z分量

B=sqrt（Bz.^2+Bx.^2）；A=atan2（Bx，Bz）；%合磁場，磁場角度

BC={Bz，Bx，B，A}；%數(shù)據(jù)元胞

zc={’itB_z/B m_0’，’itB_x/B m_0’，...

’itB/B m_0’，’italpha m’}；%豎坐標元胞

tc={’itz m分量’，’itx m分量’，’總量itB m’，’方向’}；%標題的一部分

txt=’itB m_0=itmu m_0itI m/2ita’；%磁感應(yīng)強度單位文本

t=’亥姆霍茲線圈’；%標題

if l==0.5，t=［t，’臨界’］；end，fs=16；%如果距離為0.5，修改標題，結(jié)束循環(huán)，字體大小

for i=1：4，figure%循環(huán)，創(chuàng)建圖形窗口

surf（Z，X，BC{i}），box on%畫曲面，加框

title（［t，’磁感應(yīng)強度的’，tc{i}，’的分布面（itL m=’，...

num2str（l），’ita m）’］，’FontSize’，fs）%顯示標題

xlabel（’itz/a’，’FontSize’，fs）%顯示x坐標

ylabel（’itx/a’，’FontSize’，fs）%顯示y坐標

zlabel（zc{i}，’FontSize’，fs）%顯示z坐標

end，view（-45，60）%結(jié)束循環(huán)，設(shè)置角度曲面的視角

附錄2

ring2B2main.m

%亥姆霍茲線圈均勻區(qū)域的主程序

clear，l=0.5；%ring2B2fun（l）%清除變量，距離，調(diào)用函數(shù)文件

ring2B2fun（0.45）%調(diào)用函數(shù)文件

ring2B2fun（0.4）%調(diào)用函數(shù)文件

ring2B2fun（0.55）%調(diào)用函數(shù)文件

ring2B2fun（0.6）%調(diào)用函數(shù)文件

ring2B2fun.m

%亥姆霍茲線圈的磁場均勻區(qū)域的范圍的函數(shù)文件

function fun（l）

rm=2；r=linspace（-rm，rm，200）；%最大坐標，坐標向量（繞過奇點）

［Z，X］=meshgrid（r）；R1=sqrt（（Z-l）.^2+（abs（X）+1）.^2）；%坐標矩陣，最大距離

KK1=4*abs（X）./R1.^2；［K1，E1］=ellipke（KK1）；%模數(shù)，兩個完全橢圓積分

Bx1=（（1+X.^2+（Z-l）.^2）./（（abs（X）-1）.^2+（Z-l）.^2）.*E1-K1）.*（Z-l）./X/pi./R1；%B的x分量

Bz1=（（1-X.^2-（Z-l）.^2）./（（abs（X）-1）.^2+（Z-l）.^2）.*E1+K1）/pi./R1；%B的z分量

R2=sqrt（（Z+l）.^2+（abs（X）+1）.^2）；%最大距離

KK2=4*abs（X）./R2.^2；［K2，E2］=ellipke（KK2）；%模數(shù)，兩個完全橢圓積分

Bx2=（（1+X.^2+（Z+l）.^2）./（（abs（X）-1）.^2+（Z+l）.^2）.*E2-K2）.*（Z+l）./X/pi./R2；%B的x分量

Bz2=（（1-X.^2-（Z+l）.^2）./（（abs（X）-1）.^2+（Z+l）.^2）.*E2+K2）/pi./R2；%B的z分量

Bx=Bx1+Bx2；Bz=Bz1+Bz2；B=sqrt（Bz.^2+Bx.^2）；%合磁場的x分量，z分量，合磁場

m=3；B（B>m）=m；fs=16；%截斷磁場，大于截斷磁場者取截斷磁場，字體大小

t=’亥姆霍茲線圈’；%標題

if l==0.5，t=［t，’臨界’］；end %如果距離為0.5，修改標題，結(jié)束循環(huán)

figure，surf（Z，X，B）%開創(chuàng)圖形窗口，畫B解析解的網(wǎng)格曲線

shading interp，alpha（0.8）%染色，半透明

bs=2/（1+l^2）^（3/2）；hold on%原點磁場，保持屬性

contour3（Z，X，B，［1，1］*bs，’b’）%畫中等值線

C1=contour3（Z，X，B，［0.99，0.99］*bs，’r’）；%畫下等值線取坐標

C2=contour3（Z，X，B，［1.01，1.01］*bs，’k’）；%畫上等值線取坐標

title（［t，’的合磁場（itL m=’，num2str（l），’ita m）’］，’FontSize’，fs）%顯示標題

xlabel（’itz/a’，’FontSize’，fs）%顯示x坐標

ylabel（’itx/a’，’FontSize’，fs）%顯示y坐標

zlabel（’itB/B m_0’，’FontSize’，fs）%顯示z坐標

view（-45，45），box on%設(shè)置視角，加框

plot3（0，0，bs，’r.’，’MarkerSize’，10）%畫原點磁場

txt=［’itL m=’，num2str（l），’ita m’］；%距離字符串

text（0，0，m，txt，’FontSize’，fs）%顯示距離

x1=C1（1，2：end）；y1=C1（2，2：end）；%取下線橫，縱坐標

n=C2（2，1）；x2=C2（1，2：1+n）；y2=C2（2，2：1+n）；%上線坐標對數(shù)，取上線坐標

figure，fill（x1，y1，’r’）%開創(chuàng)圖形窗口，畫下線填色圖

alpha（0.5），hold on%平面半透明，保持屬性

fill（x2，y2，’w’），fill（x2，-y2，’w’）%畫左下方上線，左上方上線白色圖

fill（-x2，-y2，’w’），fill（-x2，y2，’w’）%畫右上方上線，右下方上線白色圖

contour（Z，X，B，［1，1］*bs，’b’）%畫等值線

plot（［0；0］，［-rm；rm］），plot（［-rm；rm］，［0；0］）%畫豎直軸線，畫水平軸線

ms=10；%符號大小

plot（l，1，’ro’，l，1，’r.’，’MarkerSize’，ms）%畫右邊流出屏幕的電流

plot（l，-1，’ro’，l，-1，’rx’，’MarkerSize’，ms）%畫右邊流進屏幕的電流

plot（-l，1，’ro’，-l，1，’r.’，’MarkerSize’，ms）%畫左邊流出屏幕的電流

plot（-l，-1，’ro’，-l，-1，’rx’，’MarkerSize’，ms）%畫左邊流進屏幕的電流

plot（［l，l］，［-1+0.05；1-0.05］，’r’，’LineWidth’，4）%畫右環(huán)的剖面

plot（［-l，-l］，［-1+0.05；1-0.05］，’r’，’LineWidth’，4）%畫左環(huán)的剖面

sx=0.1：0.1：1.2；sz=0*sx+0.01；%磁感應(yīng)線的起點縱，橫坐標

Bx（Z<=-0.1）=NaN；%左邊Bx值改為非數(shù)

streamline（Z，X，Bz，Bx，sz，sx）%畫第一象限磁感應(yīng)線

streamline（-Z，X，-Bz，Bx，-sz，sx）%畫第二象限磁感應(yīng)線

streamline（Z，-X，Bz，-Bx，sz，-sx）%畫第三象限磁感應(yīng)線

streamline（-Z，-X，-Bz，-Bx，-sz，-sx）%畫第四象限磁感應(yīng)線

plot（0，0，’r.’，’MarkerSize’，10）%畫原點

title（［t，’合磁場的均勻范圍’］，’FontSize’，fs）%顯示標題

xlabel（’itz/a’，’FontSize’，fs）%顯示x坐標

ylabel（’itx/a’，’FontSize’，fs）%顯示y坐標

grid on，axis equal%加網(wǎng)格，使坐標刻度相等

text（0，1.5，txt，’FontSize’，fs）%顯示距離

if l==0.5%如果距離為0.5

text（0，0，［’itB m_C=’，num2str（bs），’itB m_0’］，’FontSize’，fs）%顯示鞍點磁場值

else%否則

text（0，0，［’itB m_S=’，num2str（bs），’itB m_0’］，’FontSize’，fs）%顯示鞍點磁場值

end%結(jié)束循環(huán)