MATLAB語(yǔ)言課程論文 基于MATLAB的電磁場(chǎng)數(shù)值圖像分析

1、基于matlab的電磁場(chǎng)數(shù)值分析應(yīng)用 摘要 matlab使用計(jì)算機(jī)進(jìn)行電磁場(chǎng)數(shù)值分析已成為電磁場(chǎng)的工程開發(fā)、科研和教學(xué)的重要手段。編程實(shí)現(xiàn)從電磁場(chǎng)微分方程到有限元求解全過(guò)程需要很好的理論基礎(chǔ)和編程技巧，難度較高。該文介紹了電磁場(chǎng)數(shù)值分析的基本理論并通過(guò)幾個(gè)實(shí)例介紹了使用matlab 實(shí)現(xiàn)電磁場(chǎng)偏微分方程的有限元解法。 實(shí)驗(yàn)結(jié)果表明這一方法具有操作簡(jiǎn)單明了!運(yùn)算速度快，計(jì)算誤差可控制等優(yōu)點(diǎn)關(guān)鍵詞 電磁場(chǎng)數(shù)值分析 matlab 麥克斯韋方程1、 問(wèn)題的提出電磁學(xué)是物理學(xué)的一個(gè)分支，是研究電場(chǎng)和電磁的相互作用現(xiàn)象。電磁學(xué)從原來(lái)互相獨(dú)立的兩門科學(xué)(電學(xué)、磁學(xué))發(fā)展成為物理學(xué)中一個(gè)完整的分支學(xué)科，主要是

2、基于電流的磁效應(yīng)和變化的磁場(chǎng)的電效應(yīng)的發(fā)現(xiàn)。這兩個(gè)實(shí)驗(yàn)現(xiàn)象，加上麥克斯韋關(guān)于變化電場(chǎng)產(chǎn)生磁場(chǎng)的假設(shè)，奠定了電磁學(xué)的整個(gè)理論體系，發(fā)展了對(duì)現(xiàn)代文明起重大影響的電工和電子技術(shù)。針對(duì)電磁場(chǎng)學(xué)習(xí)理論性強(qiáng)、概念抽象等特點(diǎn)，利用matlab強(qiáng)大的數(shù)值計(jì)算和圖形技術(shù)，通過(guò)具體實(shí)例進(jìn)行仿真，繪制相應(yīng)的圖形，使其形象化，便于對(duì)其的理解和掌握。將matlab引入電磁學(xué)中，利用其可視化功能對(duì)電磁學(xué)實(shí)驗(yàn)現(xiàn)象進(jìn)行計(jì)算機(jī)模擬，可以提高學(xué)習(xí)效率于學(xué)習(xí)積極性，使學(xué)習(xí)效果明顯。通過(guò)matlab軟件工具，對(duì)點(diǎn)電荷電場(chǎng)、線電荷產(chǎn)生的電位、平面上n個(gè)電荷之間的庫(kù)侖引力、仿真電荷在變化磁場(chǎng)中的運(yùn)動(dòng)等問(wèn)題分別給出了直觀形象的的仿真圖和數(shù)

3、值分析，形實(shí)現(xiàn)了可視化學(xué)習(xí)，豐富了學(xué)習(xí)內(nèi)容，提高了對(duì)電磁場(chǎng)理論知識(shí)的興趣。從而更好地解決電磁場(chǎng)中數(shù)值分析的問(wèn)題。二、電磁場(chǎng)數(shù)值解法 麥克斯韋方程組是電磁場(chǎng)理論的基礎(chǔ)，也是電磁場(chǎng)數(shù)值分析的出發(fā)點(diǎn)。它的微分形式方程： （1）式中磁場(chǎng)強(qiáng)度電通密度電場(chǎng)強(qiáng)度磁感應(yīng)強(qiáng)度。電磁場(chǎng)中各種場(chǎng)量之間的關(guān)系由媒質(zhì)的特性確定。在各向同性媒質(zhì)中，由下列結(jié)構(gòu)方程組確定 (2)為獲得電磁場(chǎng)問(wèn)題的唯一解!除上述方程組之外尚需給出定解條件，對(duì)靜態(tài)場(chǎng)和穩(wěn)態(tài)場(chǎng)只需加邊界條件，對(duì)時(shí)變場(chǎng)還需另加初始條件。邊界條件包括：（1）第一類邊界條件是給定邊界上的值，其中是邊界點(diǎn)的函數(shù)或常數(shù)；（2）第二類邊界條件給定邊界上法向?qū)?shù)的值；（3）第三

4、類邊界條件給定邊界值與法向?qū)?shù)的線性組合 (3)根據(jù)麥克斯韋方程組和結(jié)構(gòu)方程組!在靜電場(chǎng)中，以電位 為求解對(duì)象，在各向同性介質(zhì)中，電位滿足泊松方程： (4)在恒定磁場(chǎng)中，取矢量磁位為求解對(duì)象，令有 (5)考慮電磁波動(dòng)方程： (6)在正弦穩(wěn)態(tài)條件下,上式可分別導(dǎo)出亥姆霍茲方程: (7)如上述幾個(gè)例子,對(duì)于不同的電磁場(chǎng)實(shí)際應(yīng)用問(wèn)題求解可以得到對(duì)應(yīng)的電磁場(chǎng)偏微分方程,直接用解析法求解這些方程組往往會(huì)遇到很多困難甚至無(wú)法求解,電磁場(chǎng)數(shù)值分析方法已成為求解電磁場(chǎng)問(wèn)題的重要方法。數(shù)值分析方法將原來(lái)連續(xù)的場(chǎng)域離散化求解， 再用離散點(diǎn)的結(jié)果近似逼近連續(xù)場(chǎng)域的解&常用的電磁場(chǎng)數(shù)值分析方法包括有限差分法、邊界元法

5、和有限元法。 有限差分法是以差分原理為基礎(chǔ)的一種數(shù)值計(jì)算方法，即用差分方程代替偏微分方程，把要求解的邊值問(wèn)題轉(zhuǎn)化為一組相應(yīng)的差分問(wèn)題，將求解區(qū)域劃分，求解差分方程組從而得出各網(wǎng)格單元的場(chǎng)值。這種方法的特點(diǎn)是方法簡(jiǎn)單，網(wǎng)格劃分容易，但對(duì)不規(guī)則邊界處理不便，網(wǎng)格劃分缺乏靈活性。邊界元法以麥克斯韋積分方程為基礎(chǔ)，它采用分步積分如格林定理等在一定條件下把該積分方程轉(zhuǎn)化為關(guān)于邊界的積分方程，并據(jù)此進(jìn)行離散獲得對(duì)應(yīng)的代數(shù)方程!求出場(chǎng)域中變量的數(shù)值。 它的特點(diǎn)是將數(shù)值法與解析法相結(jié)合， 在數(shù)學(xué)上起到降維的作用，減少了計(jì)算量，但對(duì)非線性情況失去了高精度特點(diǎn)，有局限性。有限元法由于單元定義靈活，處理邊界條件容易

6、，具有正定對(duì)稱系數(shù)矩陣而占據(jù)主導(dǎo)地位&有限元法是根據(jù)變分原理和離散化而取得近似解的方法。是先從偏微分方程邊值問(wèn)題出發(fā)，找出一個(gè)能量泛函的積分式，并令其在滿足第一類邊界條件的前提下取極值，即構(gòu)成條件變分問(wèn)題& 例如與上式對(duì)應(yīng)的二維泊松場(chǎng)第三類邊界條件下的泛函極值問(wèn)題為： (8)與上式對(duì)應(yīng)的亥姆霍茲方程泛函為： (9) 這個(gè)條件變分問(wèn)題是和偏微分方程邊值問(wèn)題等價(jià)的。有限元法便是以條件變分問(wèn)題為對(duì)象來(lái)求解電磁場(chǎng)問(wèn)題。在求解過(guò)程中，將場(chǎng)的求解區(qū)域剖分成有限個(gè)單元，因此在單元中構(gòu)造出插值函數(shù)，將插值函數(shù)代入能量泛函的積分式，再把泛函離散化成多元函數(shù)。通過(guò)多元函數(shù)極值求極值方法得到一個(gè)代數(shù)方程組。最后由此

7、方程組求解得到數(shù)值解。下面舉幾例論證： 例證1設(shè)圓線圈的中心為o，半徑為r，放置于y-z平面，線圈通過(guò)的電流為i0，如右圖所示。用畢奧薩伐爾定律計(jì)算載流圓線圈在z=0處x-y平面上的磁場(chǎng)分布。解題分析根據(jù)畢奧薩伐爾定律， (10)線圈上任一點(diǎn)處的電流元在x-y平面上一點(diǎn)p產(chǎn)生的元磁場(chǎng)為db。在編制程序時(shí)，將電流環(huán)分為n段，每一小段視為一電流元，然后求出每一電流元在觀察點(diǎn)處的磁場(chǎng)分量，求出總磁場(chǎng)，最后疊加。 matlab程序如下： clear all %清除 r=input(請(qǐng)輸入圓環(huán)半徑,r=); %定義輸入變量 i0=input(請(qǐng)輸入電流,i0=); %定義輸出變量mu0=4*pi*1e-

8、7; c0=mu0/(4*pi); %歸并常數(shù)n=20; %電流環(huán)分段數(shù)x=linspace(-3,3,n); y=x; %確定觀測(cè)點(diǎn)范圍theta0=linspace(0,2*pi,n+1); %環(huán)的圓周角分段 theta1=theta0(1:n); y1=r*cos(theta1); z1=r*sin(theta1); %環(huán)各段矢量的起始坐標(biāo)y1,z1theta2=theta0(2:n+1); y2=r*cos(theta2); z2=r*sin(theta2); %環(huán)各段矢量的終點(diǎn)坐標(biāo)y2,z2xc=0; yc=(y2+y1)./2; zc=(z2+z1)./2; %計(jì)算環(huán)各段矢量中點(diǎn)的

9、三個(gè)坐標(biāo)分量xc,yc,zcdlx=0;dly=y2-y1;dlz=z2-z1; %計(jì)算環(huán)各段矢量dl的三個(gè)長(zhǎng)度分量，其中x1=x2=0。 ngx=n; ngy=ngx; %網(wǎng)格線數(shù)for i=1:ngy %循環(huán)計(jì)算各網(wǎng)點(diǎn)上的b（x,y）值 for j=1:ngxrx=x(j)-xc; ry=y(i)-yc; rz=0-zc; %計(jì)算徑矢r的3個(gè)長(zhǎng)度分量，r在z=0的面上。 r3=sqrt(rx.2+ry.2+rz.2).3; %計(jì)算r3 dlxr_x=dly.*rz-dlz.*ry; %計(jì)算叉乘dlr的x和y分量，z分量為0 dlxr_y=dlz.*rx-dlx.*rz; bx(i,j)=s

10、um(c0*i0.*dlxr_x./r3); %把環(huán)各段產(chǎn)生的磁場(chǎng)分量累加by(i,j)=sum(c0*i0.*dlxr_y./r3);b=(bx.2+by.2).0.5; %計(jì)算b的大小endendsubplot(1,2,1), quiver(x,y,bx,by), %畫矢量場(chǎng)圖hold on plot(0,1,ro,0,-1,bo),xlabel(x),ylabel(y), %修飾圖形，標(biāo)注坐標(biāo)軸axis(-3,3,-3,3), subplot(1,2,2)mesh(x,y,b);axis(-3,3,-3,3,0,1e-4) %畫磁場(chǎng)大小分布圖xlabel(x),ylabel(y),zla

11、bel(b)運(yùn)行該程序，例如，在命令窗中的r和i0的提示后分別鍵入1和100，運(yùn)行結(jié)果如圖1所示圖1 載流圓線圈的磁場(chǎng)分布度及位置曲線例證2一對(duì)相同的圓形線圈，彼此平行而共軸。設(shè)兩線圈內(nèi)的電流都是i，且回繞方向一致，線圈的半徑為r，二者的間距為a（當(dāng)a=r時(shí)，稱為亥姆霍茲線圈），求軸線附近的磁場(chǎng)分布。解題分析 本題是把觀測(cè)區(qū)域取在兩線圈之間的小范圍內(nèi)。線圈b生成的左邊的磁場(chǎng)等于線圈a的左邊磁場(chǎng)。因 此，a、b兩線圈在中間部分的合成磁場(chǎng)等于a線圈的右磁場(chǎng)與其左磁場(chǎng)平移r后的和。matlab程序如下：clear all %清除 mu0=4*pi*1e-7;c0=mu0/(4*pi); %歸并常數(shù)i

12、0=5.0;r=1; ngx=21;ngy=21; %設(shè)定網(wǎng)格線數(shù)x=linspace(-1,1,ngx); %確定觀測(cè)點(diǎn)范圍y=linspace(-1,1,ngy); n=20; %電流環(huán)分段數(shù) theta0=linspace(0,2*pi,n+1); %環(huán)的圓周角分段 theta1=theta0(1:n);y1=r*cos(theta1); z1=r*sin(theta1); %環(huán)各段矢量的起始坐標(biāo)y1,z1theta2=theta0(2:n+1);y2=r*cos(theta2); z2=r*sin(theta2); %環(huán)各段矢量的終點(diǎn)坐標(biāo)y2,z2dlx=0;dly=y2-y1;dlz

13、=z2-z1; %計(jì)算環(huán)各段矢量dl的三個(gè)長(zhǎng)度分量，其中x1=x2=0。xc=0; yc=(y2+y1)/2; zc=(z2+z1)/2;for i=1:ngy for j=1:ngx rx=x(j)-xc; ry=y(i)-yc; rz=0-zc; %計(jì)算徑矢r的3個(gè)長(zhǎng)度分量,r在z=0的面上。 r3=sqrt(rx.2+ry.2+rz.2).3; dlxr_x=dly.*rz-dlz.*ry; %計(jì)算叉乘dlr的x和y分量，z分量為0 dlxr_y=dly.*rx-dlx.*rz; bx(i,j)=sum(c0*i0*dlxr_x./r3); %把環(huán)各段產(chǎn)生的磁場(chǎng)分量累加 by(i,j)=

14、sum(c0*i0*dlxr_y./r3);endendbax=bx(:,11:21)+bx(:,1:11); %把x0區(qū)域bay=by(:,11:21)+by(:,1:11);subplot(1,2,1),mesh(x(11:21),y,bax);xlabel(x);ylabel(y);zlabel(b);subplot(1,2,2),quiver(x,y,bx,by,1.5), axis(square),axis(-1,1,-1,1),xlabel(x);ylabel(y); 運(yùn)行結(jié)果如圖2所示。可以看出，在軸線附近磁場(chǎng)大小均勻且沿x方向。圖2 亥姆霍茲線圈軸線附近bx在x-y平面上的分布

15、及矢量場(chǎng)例證3設(shè)帶電粒子質(zhì)量為m，帶電量為，電場(chǎng)強(qiáng)度e沿方向，磁感應(yīng)強(qiáng)度b沿qyz方向. 則帶電粒子在均勻電磁場(chǎng)中的運(yùn)動(dòng)微分方程為 (11)令， 則上面微分方程可化作： (12)選擇和為參量，就可以分別研究0e，0=b和，等情況. 編寫matlab程序如下：clear %清除syms w x y z t b e m q; %定義變量e=input(e=);b=input(b=); %輸入e和b值x,y,z=dsolve(d2x=q*b/m*dy,d2y=q*e/m-q*b/m*dx,d2z=0,x(0)=0,y(0)=0,z(0)=0,dx(0)=0.01,dy(0)=6,dz(0)=0.01

16、) ; %初始條件取x(0)=y(0)=z(0)=0,dx(0)=0.01,dy(0)=6,dz(0)=0.01q=1.6e-2; m=0.02; %賦值x=subs(x y z); x=x(1),y=x(2),z=x(3), ezplot3(x(1),x(2),x(3) %繪圖函數(shù)調(diào)用運(yùn)行上述程序，例如，取e=4, b=8可得下列特解并給出圖運(yùn)行結(jié)果如圖3所示。 研究時(shí)可以采用不同的初始條件和不同的參量觀察不同的現(xiàn)象。運(yùn)算特解如下：x =(t*exp(1)/8 - (100*i*exp(1) - 8*i + 4800)/(10240*exp(32*i*t)/5) - (exp(32*i*t)

17、/5)*(8*i - 100*i*exp(1) + 4800)/10240 + 15/16y =(5*exp(1)/256 + (i*(100*i*exp(1) - 8*i + 4800)/(10240*exp(32*i*t)/5) - (i*exp(32*i*t)/5)*(8*i - 100*i*exp(1) + 4800)/10240 - 1/640z =t/100圖3現(xiàn)有e=4, b=8參數(shù)運(yùn)行結(jié)果例證4如右圖所示，求垂直于無(wú)限長(zhǎng)載流直導(dǎo)線的平面內(nèi)磁感應(yīng)強(qiáng)度的分布。解題分析設(shè)場(chǎng)點(diǎn)p的位置為，電流元位置為，電流元矢量為。由此，場(chǎng)點(diǎn)p相對(duì)于電流元的位置矢量為 利用行列式計(jì)算idlr ，可寫為

18、 （13）也可利用matlab中的det 命令函數(shù)來(lái)求該行列式，matlab程序如下：syms dx dy dz x0 x y0 y z0 z; %定義變量dl=dx,dy,dz; %定義行列式r=x0-x,y0-y, z0-z;d1cr=cross(dl,r) %求dlr的積運(yùn)行結(jié)果為d1cr = dy*(z0-z)-dz*(y0-y), dz*(x0-x)-dx*(z0-z), dx*(y0-y)-dy*(x0-x)即又，r的大小為 設(shè)載流導(dǎo)體通過(guò)坐標(biāo)原點(diǎn)垂直于 x-y平面放置，電流元 idl沿z軸正向，場(chǎng)點(diǎn)p位于x-y平面上。對(duì)本題目而言，dx=dy=0,x=y=0, z0=0, 矢量叉

19、乘積為，r的大小為 由畢奧薩伐爾定律 （14）有； ； （15） matlab程序 1、 用符號(hào)運(yùn)算求b的表達(dá)式syms c0 i z x y r r0; %定義變量bx=c0.*i.*int(-y./(x.2+y.2+z.2).(3/2),z,-inf,inf)by=c0.*i.*int(x./(x.2+y.2+z.2).(3/2),z,-inf,inf)b=(bx.2+by.2).0.5 運(yùn)行結(jié)果：bx =-2*c0*i*y/(x2+y2)(3/2)/(1/(x2+y2)(1/2)by =2*c0*i*x/(x2+y2)(3/2)/(1/(x2+y2)(1/2)b =(4*c02*i2*y

20、2/(x2+y2)2+4*c02*i2*x2/(x2+y2)2)(1/2) 即 ； （16）2、繪制磁場(chǎng)大小分布圖和矢量場(chǎng)圖x=-0.5:0.05:0.5;y=x; %賦值i=input(請(qǐng)輸入電流i=); %設(shè)置輸入mu0=4*pi*1e-7; c0=mu0/(4*pi); %設(shè)置函數(shù)x,y=meshgrid(x,y);bx =-2.*c0.*i.*y./(x.2+y.2).(3./2.)./(1./(x.2+y.2).(1./2);by =2.*c0.*i.*x./(x.2+y.2).(3./2)./(1./(x.2+y.2).(1./2);b=(4.*c0.2.*i.2.*y.2./(x

21、.2+y.2).2+4.*c0.2.*i.2.*x.2./(x.2+y.2).2).(1./2); subplot(1,2,1) %選擇12區(qū)域中的1號(hào)區(qū)quiver(x,y,bx,by,2), axis(-0.5,0.5,-0.5,0.5), axis(square),subplot(1,2,2) %選擇12區(qū)域中的2號(hào)區(qū)mesh(x,y,b) %繪圖運(yùn)行該程序，在命令窗中的提示后鍵入i0值 (例如，取i0100a)，便得到圖4所示圖形。圖4長(zhǎng)載流導(dǎo)線的磁場(chǎng)在x-y平面上的分布三、結(jié)論從以上利用matlab語(yǔ)言對(duì)幾種電磁場(chǎng)模型的分析我們不難的出以下結(jié)論：電磁場(chǎng)與電磁波理論作為電子信息類專業(yè)的一門重要