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

通電螺線管磁場的雙重數(shù)值積分法和可視化

莫云飛1, 周群益2, 侯兆陽3, 周麗麗4

(1. 長沙學(xué)院 電子信息與電氣工程學(xué)院, 湖南 長沙, 410022; 2. 廣州理工學(xué)院, 廣東 廣州, 510540; 3. 長安大學(xué) 理學(xué)院 應(yīng)用物理系, 陜西 西安, 710064; 4. 贛南醫(yī)學(xué)院 信息工程學(xué)院, 江西 贛州, 341000)

通電螺線管可以當(dāng)作環(huán)電流密繞而成的, 本文利用環(huán)電流的磁場公式推導(dǎo)了通電螺線管磁感應(yīng)強度的積分式, 將公式無量綱化。利用MATLAB的多維矩陣形成多維被積函數(shù), 用二重數(shù)值積分計算場強。利用MATLAB的曲面指令繪制出磁感應(yīng)強度2個分量以及合場強和方向曲面, 利用MATLAB的流線指令畫出管內(nèi)外的磁感應(yīng)線, 充分顯示了通電螺線管場強的分布規(guī)律。

通電螺線管; 磁感應(yīng)強度; 數(shù)值積分; 無量綱作圖

通電螺線管的磁場是電磁場理論中的一個典型問題, 文獻(xiàn)[1–3]只計算軸線上的磁感應(yīng)強度, 繼而得到無限長通電螺線管在軸線上的場強公式; 文獻(xiàn)[4–9]討論了無限長通電螺線管內(nèi)部的場強。文獻(xiàn)[10–13]推導(dǎo)了有限長通電螺線管內(nèi)外的場強的積分公式, 其中有些論文還畫出了曲線。通電螺線管的磁場沒有精確的解析解, 文獻(xiàn)[14]研究了數(shù)值解, 并畫出曲線族。文獻(xiàn)[15–16]用MATLAB和Mathema- tica計算和研究場強、繪制曲線和曲面, 這是2種很好的方法。

1 通電螺線管的磁場

本文建立直角坐標(biāo)系中場強的雙重積分公式, 將公式無量綱化, 通過雙重數(shù)值積分的方法計算了管內(nèi)外磁感應(yīng)強度的大小, 并用曲面和磁感應(yīng)線顯示磁場的分布規(guī)律。

如圖1所示, 設(shè)螺線管的半徑為, 長度為2, 通有電流。在圓環(huán)取一電流元d, 螺線管各匝線圈都是螺線形的, 在密繞的情況下可以當(dāng)作多匝圓形線圈緊密排列而成。設(shè)單位長度上的匝數(shù)為, 在長度處取一線元d, 線圈匝數(shù)為d=d, 電流強度

d=d=d。 (1)

圓環(huán)形電流元d到點的軸坐標(biāo)之間的距離為–, 環(huán)上的點到場點的距離

這里,是坐標(biāo)和以及方位角和長度的四元函數(shù)=(,,,)。根據(jù)環(huán)電流產(chǎn)生的場強的公式, 環(huán)電流元在空間產(chǎn)生的磁感應(yīng)強度的2個分量的微元為[17–19]

場強的2個分量

這是雙重定積分公式,B=B(,)和B=B(,)是坐標(biāo)和的二元函數(shù)。

合場強

方向由角度確定

能量密度

當(dāng)= 0時, 距離

場強

, (12)

這是軸線上的場強。

當(dāng)= 0時, 距離

由于B(, 0)是的奇函數(shù), 所以B(, 0) = 0,B(, 0)由式(6)決定, 這是中垂線上的場強。

2 磁感應(yīng)強度公式的無量綱化

取半徑為長度單位, 則無量綱的坐標(biāo)

*=/,*=/。 (14)

無量綱的距離

其中,*=/是無量綱的長度或相對長度。取0=0為磁感應(yīng)強度單位, 則無量綱的磁感應(yīng)強度的兩個分量

, (16)

其中,*=/是無量綱的半長。

無量綱的合場強

方向角

軸線上無量綱的場強

。 (21)

將磁感應(yīng)強度公式無量綱化是為了便于數(shù)值計算。應(yīng)用MATLAB, 將被積函數(shù)化為多維矩陣, 利用雙重數(shù)值積分的方法求出平面內(nèi)各點的磁感應(yīng)強度。

圖2 通電螺線管軸線和中垂上的磁感應(yīng)強度

3 磁場的可視化

根據(jù)磁感應(yīng)強度的數(shù)值, 利用MATLAB的plot指令可以畫出曲線族, 利用surf指令可以畫出磁感應(yīng)強度的曲面。利用流線指令streamline可以畫出磁感應(yīng)線。

(1) 通電螺線管軸線上的磁場如圖2(a)所示。當(dāng)管的長度一定時, 中間的磁感應(yīng)強度最大; 當(dāng)=時, 最大場強比較小。管越長, 管心處的磁感應(yīng)強度就越接近于0, 軸線上均勻的范圍就越大, 管口的場強約為最大場強的一半。通電螺線管中垂線上的磁場如圖2(b)所示, 當(dāng)=時, 中垂面內(nèi)部的場強比較小, 管心處的場強最小; 管越長, 中垂面內(nèi)部場強越均勻, 越接近于0, 外部場強越接近于0,= ±處兩側(cè)的場強變化越光滑。

(2) 當(dāng)= 5時, 通電螺線管磁場的分量B的曲面如圖3所示。管內(nèi)外的B都很小; 只在在管口邊緣,B有兩對對稱的“峰”和“谷”, “峰”和“谷”的B方向相反。軸線和中垂線上的場強B(, 0)和B(, π/2)為0, 其直線分布在水平面上。

(3) 磁場的分量B的曲面如圖4所示。管內(nèi)B形似一堵“墻”,= ±處是“墻”的邊緣; 在管口,B迅速減小; 在管口的邊緣,B< 0, 場強方向相反; 在管外, 除了管口附近,B很小。軸線和中垂線上的場強B(,0)和B(,π/2)曲線分布在曲面上。

圖3 通電螺線管磁場x分量的分布面(L = 5a)

(4) 合場強的曲面如圖5所示。圖5與圖4類似, 不過,恒大于0; 除了管口邊緣之外, 管外的曲面比較平, 也就是接近于0。軸線和中垂線上的場強(, 0) =B(, 0)和(, π/2) = |B(, π/2)|曲線分布在曲面上。場強的能量密度曲面與圖5類似, 只是管外的曲面更加平, 也就是更加接近于0。

圖4 通電螺線管磁場z分量的分布面(L = 5a)

圖5 通電螺線管合場強的分布面(L = 5a)

(5) 合場強的方向角曲面如圖6所示。在= 0的軸線上,= 0。取軸為極軸, 在管內(nèi),幾乎為0, 表示場強的方向幾乎與軸線方向一致; 在管外,隨極角的增加而增加,在= ± π/2處發(fā)生從-π到π的躍變。

(6) 通電螺線管內(nèi)外的磁感應(yīng)線如圖7所示。在管內(nèi), 磁感應(yīng)線密集而均勻, 可當(dāng)作勻強磁場; 在管處, 磁感應(yīng)線比較稀疏, 磁感應(yīng)強度和能量密度都比較小。因此, 螺線管的能量可以當(dāng)作均勻分布在管內(nèi)處理。

圖6 通電螺線管磁場方向的分布面(L = 5a)

圖7 通電螺線管的磁感應(yīng)線(L = 5a)

4 結(jié)束語

本文解決了通電螺線管磁場的計算和可視化的問題: (1) 建立了直角坐標(biāo)系中通電螺線管的磁感應(yīng)強度的雙重定積分公式; (2) 將公式無量綱化, 以便純數(shù)值計算; (3) 用MATLAB做雙重數(shù)值積分解決了計算問題; (4) 畫出了彩色場強分布圖和磁感應(yīng)線圖, 說明磁場的分布規(guī)律; (5)設(shè)計了2個MATLAB程序(見附錄), 一個專門計算軸線和中垂線上的場強, 畫出曲線族; 一個專門計算空間場強, 畫出曲面和磁感應(yīng)線。2個程序都不到100行。掌握程序設(shè)計方法是提出問題和解決問題的重要手段; (6) 程序采用雙重積分公式。如果利用兩類完全橢圓積分, 可以將雙重積分公式化為單一積分公式, 計算效率更高, 曲面可以畫得更加精致。

[1] 趙凱華, 陳熙謀. 電磁學(xué)[M]. 北京: 高等教育出版社, 2004.

[2] 張三慧. 大學(xué)物理學(xué)(下冊)[M]. 3版. 北京: 清華大學(xué)出版社, 2010.

[3] 周群益. MATLAB可視化大學(xué)物理學(xué)[M]. 北京: 清華大學(xué)出版社, 2011.

[4] 成軍. 無限長直均勻密繞載流螺線管的磁場[J]. 浙江師大學(xué)報(自然科學(xué)版), 2001, 24(4): 325–326.

[5] 王國強. 螺線管磁場的另一計算方法[J]. 高等函授學(xué)報(自然科學(xué)版), 2003, 16(2): 14–16.

[6] 張保軍. 無限長載流螺線管磁場分布特點的證明[J]. 江蘇技術(shù)師范學(xué)院學(xué)報, 2005, 11(6): 17–19.

[7] 丁紅, 藍(lán)海江. 無限長通電直螺線管內(nèi)、外的磁感應(yīng)強度[J]. 欽州學(xué)院學(xué)報, 2010, 25(3): 4–7.

[8] 牛中明, 于欣欣. 無限長直密繞通電螺線管磁場的一種簡單計算方法[J]. 赤峰學(xué)院學(xué)報(自然科學(xué)版), 2018, 34(7): 11–12.

[9] 李永樂. 無限長密繞通電螺線管磁場簡明計算[J]. 物理通報, 2017(11): 24–26.

[10] 程昌林, 王慧, 李業(yè)鳳. 細(xì)導(dǎo)線密繞螺線管的磁場[J]. 物理與工程, 2003, 13(1): 6–8.

[11] 侯昭武. 通電線圈產(chǎn)生的電磁場強度[J]. 欽州學(xué)院學(xué)報, 2007, 22(6): 55–58.

[12] 李春生, 楊中海, 黃桃. 有限長通電螺線管空間磁場分析[J]. 現(xiàn)代電子技術(shù), 2009, 32(11): 28–30.

[13] 梁麗萍. 關(guān)于載流長直螺線管磁場的討論[J]. 教育教學(xué)論, 2014(18): 106–108.

[14] 劉耀康. 螺線管磁場的數(shù)值解[J]. 高師理科學(xué)刊, 2008, 28(1): 60–64.

[15] 惠小強, 陳文學(xué). 有限長通電螺線管空間的磁場分[J]. 物理與工程, 2004, 14(2): 22–25.

[16] 王鍇, 廖斌, 吳先映. 利用 Matlab研究多螺線管磁場分布[J]. 北京師范大學(xué)學(xué)報(自然科學(xué)版), 2013, 49(6): 565– 570.

[17] 向裕民. 圓環(huán)電流磁場的普遍分布[J]. 大學(xué)物理, 1999, 18(1): 14–17.

[18] 李海, 張玉穎. 圓環(huán)線電流的磁感應(yīng)強度[J]. 大學(xué)物理, 1999, 18(6): 20–22.

[19] 朱平. 圓電流空間磁場分布[J]. 大學(xué)物理, 2005, 24(9): 13–17.

ringN1.m%通電螺線管軸線上和中垂線的場強分布(雙重數(shù)值積分)

clear,l=1:2:13;%清除變量,半長度與半徑比的向量

zm=20;z=-zm:0.1:zm;%橫坐標(biāo)的最大范圍,自變量向量

[Z,L]=meshgrid(z,l);%坐標(biāo)和比值矩陣

Bzr=((L-Z)./sqrt(1+(L-Z).^2)+(L+Z)./sqrt(1+(L+Z).^2))/2;%相對磁感應(yīng)強度

fs=16;figure,subplot(2,1,1)%字體大小,創(chuàng)建圖形窗口,選子圖

plot(z,Bzr,'LineWidth',2)%畫B的直角坐標(biāo)曲線族

legend([repmat('itL/a m=',length(l),1),num2str(l')])%加圖例

grid on,hold on,axis([-zm,zm,0,1])%加網(wǎng)格,曲線范圍,保持圖像

plot([-zm;zm],[0.5;0.5],'-.','LineWidth',2)%畫水平零線

xlabel('itz/a','FontSize',fs)%顯示橫坐標(biāo)

ylabel('itB_z m(0,itz m)/itmu m_0itnI','FontSize',fs)%顯示縱坐標(biāo)

title('通電螺線管軸線上的磁場','FontSize',fs)%顯示標(biāo)題

ph=linspace(0,pi);%角度向量

r='(1+l.^2+x.^2-2*x.*cos(ph)).^(3/2)';%距離三次方字符串

ddbzx=inline(['(1-x.*cos(ph))./',r],'x','l','ph');%x軸Bz分量被積內(nèi)線函數(shù)

xm=3;x=linspace(-xm,xm,500);Bzx=[];%橫坐標(biāo)的最大范圍,縱坐標(biāo)向量,矩陣置空

for i=1:length(l)%按長度循環(huán)

li=linspace(-l(i),l(i));%管的長度向量(繞過奇點)

[X,L,PH]=ndgrid(x,li,ph);%設(shè)置數(shù)三維矩陣

dbzx=trapz(ddbzx(X,L,PH),3)/2/pi*ph(2);%x軸Bz分量對角度積分

bzx=trapz(dbzx,2)*(li(2)-li(1));%x軸Bz分量對長度積分

Bzx=[Bzx,bzx];%連接矩陣

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

subplot(2,1,2)%選子圖

plot(x,Bzx,'LineWidth',2)%畫B的直角坐標(biāo)曲線族

legend([repmat('itL/a m=',length(l),1),num2str(l')])%加圖例

grid on%加網(wǎng)格

xlabel('itx/a','FontSize',fs)%顯示橫坐標(biāo)

ylabel('itB_z m(itx m,0)/itmu m_0itnI','FontSize',fs)%顯示縱坐標(biāo)

title('通電螺線管中垂線上的磁場','FontSize',fs)%顯示標(biāo)題

ringN2.m

%通電螺線管的磁場(雙重數(shù)值積分)

clear,lm=5;%清除變量,螺線管半長度與半徑之比

r3='(1+(z-l).^2+x.^2-2*x.*cos(ph)).^(3/2)';%距離3次方字符串

ddBx=inline(['(z-l).*cos(ph)./',r3],'x','z','l','ph');%磁場x分量被積內(nèi)線函數(shù)

ddBz=inline(['(1-x.*cos(ph))./',r3],'x','z','l','ph');%磁場z分量被積內(nèi)線函數(shù)

ph=linspace(0,pi);%角度向量

l=linspace(-lm,lm);dl=l(2)-l(1);%管的長度向量(繞過奇點),管的長度間隔

zm=2*lm;z=linspace(0,zm,30);%z坐標(biāo)范圍,橫坐標(biāo)向量

xm=4/5*zm;x=linspace(0,xm,30);%x坐標(biāo)范圍,縱坐標(biāo)向量

z(1)=sqrt(eps);x(1)=sqrt(eps);%零改為小量

[X,Z,L,PH]=ndgrid(x,z,l,ph);%設(shè)置四維矩陣

dBx=trapz(ddBx(X,Z,L,PH),4)/2/pi*ph(2);%Bx分量對角度積分

dBz=trapz(ddBz(X,Z,L,PH),4)/2/pi*ph(2);%Bz分量對角度積分

Bx=trapz(dBx,3)*dl;Bz=trapz(dBz,3)*dl;%Bx,Bz分量對長度積分

clear PH L%清除兩個矩陣

sx=0.1:0.1:0.9;sz=0*sx+0.001;%磁感應(yīng)線的起點坐標(biāo)

figure,plot([0;0],[-xm;xm])%創(chuàng)建圖形窗口,畫豎直軸線

hold on,plot([-zm;zm],[0;0])%保持圖像,畫水平軸線

ms=6;z0=-lm:0.25:lm;%符號大小,管的橫坐標(biāo)向量

x0=ones(size(z0));%上管的縱坐標(biāo)向量

plot(z0,x0,'ro',z0,x0,'r.','MarkerSize',ms)%畫上面流出屏幕的電流

plot(z0,-x0,'ro',z0,-x0,'rx','MarkerSize',ms)%畫下面流出屏幕的電流

plot([1 1]*lm,[-1+0.1;1-0.1],'r','LineWidth',3)%畫右端的剖面

plot([-1 -1]*lm,[-1+0.1;1-0.1],'r','LineWidth',3)%畫左端的剖面

grid on,axis equal,fs=16;%加網(wǎng)格,使坐標(biāo)刻度相等,字體大小

title(['通電螺線管的磁感應(yīng)線(itL m=',num2str(lm),'ita m)'],...

'FontSize',fs)%顯示標(biāo)題

xlabel('itz/a','FontSize',fs)%顯示橫坐標(biāo)

ylabel('itx/a','FontSize',fs)%顯示縱坐標(biāo)

Z=Z(:,:,1);X=X(:,:,1);%取橫坐標(biāo)矩陣,取縱坐標(biāo)矩陣

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)線

r3='(1+l.^2+x.^2-2*x.*cos(ph)).^(3/2)';%距離3次方字符串

ddbzx=inline(['(1-x.*cos(ph))./',r3],'x','l','ph');%x軸Bz被積內(nèi)線函數(shù)

x=linspace(-xm,xm);%縱坐標(biāo)向量

[XX,L,PH]=ndgrid(x,l,ph); %設(shè)置數(shù)三維矩陣

dbzx=trapz(ddbzx(XX,L,PH),3)/2/pi*ph(2);%x軸Bz分量對角度積分

bzx=trapz(dbzx,2)*dl;%x軸Bz分量對長度積分

Z=[-fliplr(Z),Z;-fliplr(Z),Z];%四象限橫坐標(biāo)矩陣

X=[flipud(X),flipud(X);-X,-X];%四象限縱坐標(biāo)矩陣

Bz=[fliplr(flipud(Bz)),flipud(Bz);fliplr(Bz),Bz];%四象限Bz矩陣

Bx=[-fliplr(flipud(Bx)),flipud(Bx);fliplr(Bx),-Bx];%四象限Bx矩陣

B=sqrt(Bz.^2+Bx.^2);%總場強

A=atan2(Bx,Bz)*180/pi;%場強角度

BC={Bx,Bz,B,A,B.^2};%數(shù)據(jù)元胞

zc={'itB_x/B m_0','itB_z/B m_0',...

'itB/B m_0','italpha m/circ','itomega/omega m_0'};%豎坐標(biāo)元胞

tc={'itx m分量','itz m分量','合場強itB m','方向',... '能量密度itomega m'};%標(biāo)題的一部分

z=linspace(-zm,zm);%橫坐標(biāo)向量

bzr=((z+lm)./sqrt(1+(z+lm).^2)-(z-lm)./sqrt(1+(z-lm).^2))/2;%軸線上磁感應(yīng)強度

for i=1:length(BC)%循環(huán)

figure,mesh(Z,X,BC{i})%創(chuàng)建圖形窗口,畫曲面

box on,axis tight%加框,帖軸

t=['通電螺線管磁場',tc{i},'的分布面(itL m=',...

num2str(lm),'ita m)'];%標(biāo)題

title(t,'FontSize',fs)%顯示標(biāo)題

xlabel('itz/a','FontSize',fs)%顯示橫坐標(biāo)

ylabel('itx/a','FontSize',fs)%顯示縱坐標(biāo)

zlabel(zc{i},'FontSize',fs),hold on%顯示豎坐標(biāo),保持圖像

if i==1,plot3(z,0*z,0*z,'r',0*x,x,0*x,'m','LineWidth',2),end%畫線

if i==2,plot3(z,0*z,bzr,'r',0*x,x,bzx,'m','LineWidth',2),end%畫線

if i==3,plot3(z,0*z,bzr,'r',0*x,x,abs(bzx),'m','LineWidth',2),end%畫線

if i==4,plot3(z,0*z,0*z,'r','LineWidth',2),view(-75,60),end%調(diào)整視角

if i==5,plot3(z,0*z,bzr.^2,'r',0*x,x,bzx.^2,'m','LineWidth',2),end%畫線

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

Calculating the magnetic filed of current solenoid according to a double numerical integral and its visualization

Mo Yunfei1, Zhou Qunyi2, Hou Zhaoyang3, Zhou Lili4

(1. School of Electronic Information and Electrical Engineering, Changsha University, Changsha 410022, China; 2. Guangzhou Institute of Science and Technology, Guangzhou, 510540, China; 3. School of Science, Chang’an University, Xi’an 710064, China; 4. Department of Information Engineering, Gannan Medical University, Ganzhou 341000, China)

The current solenoid can be considered that it is composed of the tightly wound loop current, this work has been obtained the integration formula of the magnetic field for current solenoid depending on the magnetic field formula for ring current, and the formula is also nondimensionalized. After the multi-dimensional integrand has been produced by using the multi-dimensional matrix in MATLAB, then the magnetic field is obtained according to a double numerical integral. The two components of magnetic field, the composite of filed and the direction surface are drawn by using the surface instruction of MATLAB. The magnetic lines inside and outside of the solenoid are also drawn with the streamline instruction of MATLAB, which can fully show the distribution characterization of magnetic field strength in the current solenoid.

current solenoid; magnetic field;numerical integral; dimensionless drawing

O 441

A

1672–6146(2020)04–0020–07

10.3969/j.issn.1672–6146.2020.04.004

周群益,1845901757@qq.com。

2020–04–24

湖南省自科基金(2018JJ3560); 長沙市科技計劃項目(kc1809022); 湖南省教育廳科學(xué)研究項目(19C0176)。

(責(zé)任編校: 劉剛毅)