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

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

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

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

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

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

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

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

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

9、三個坐標分量xc,yc,zcdlx=0;dly=y2-y1;dlz=z2-z1; %計算環(huán)各段矢量dl的三個長度分量，其中x1=x2=0。 ngx=n; ngy=ngx; %網(wǎng)格線數(shù)for i=1:ngy %循環(huán)計算各網(wǎng)點上的b（x,y）值 for j=1:ngxrx=x(j)-xc; ry=y(i)-yc; rz=0-zc; %計算徑矢r的3個長度分量，r在z=0的面上。 r3=sqrt(rx.2+ry.2+rz.2).3; %計算r3 dlxr_x=dly.*rz-dlz.*ry; %計算叉乘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)生的磁場分量累加by(i,j)=sum(c0*i0.*dlxr_y./r3);b=(bx.2+by.2).0.5; %計算b的大小endendsubplot(1,2,1), quiver(x,y,bx,by), %畫矢量場圖hold on plot(0,1,ro,0,-1,bo),xlabel(x),ylabel(y), %修飾圖形，標注坐標軸axis(-3,3,-3,3), subplot(1,2,2)mesh(x,y,b);axis(-3,3,-3,3,0,1e-4) %畫磁場大小分布圖xlabel(x),ylabel(y),zla

11、bel(b)運行該程序，例如，在命令窗中的r和i0的提示后分別鍵入1和100，運行結(jié)果如圖1所示圖1 載流圓線圈的磁場分布度及位置曲線例證2一對相同的圓形線圈，彼此平行而共軸。設兩線圈內(nèi)的電流都是i，且回繞方向一致，線圈的半徑為r，二者的間距為a（當a=r時，稱為亥姆霍茲線圈），求軸線附近的磁場分布。解題分析 本題是把觀測區(qū)域取在兩線圈之間的小范圍內(nèi)。線圈b生成的左邊的磁場等于線圈a的左邊磁場。因 此，a、b兩線圈在中間部分的合成磁場等于a線圈的右磁場與其左磁場平移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; %設定網(wǎng)格線數(shù)x=linspace(-1,1,ngx); %確定觀測點范圍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)各段矢量的起始坐標y1,z1theta2=theta0(2:n+1);y2=r*cos(theta2); z2=r*sin(theta2); %環(huán)各段矢量的終點坐標y2,z2dlx=0;dly=y2-y1;dlz

13、=z2-z1; %計算環(huán)各段矢量dl的三個長度分量，其中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; %計算徑矢r的3個長度分量,r在z=0的面上。 r3=sqrt(rx.2+ry.2+rz.2).3; dlxr_x=dly.*rz-dlz.*ry; %計算叉乘dlr的x和y分量，z分量為0 dlxr_y=dly.*rx-dlx.*rz; bx(i,j)=sum(c0*i0*dlxr_x./r3); %把環(huán)各段產(chǎn)生的磁場分量累加 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); 運行結(jié)果如圖2所示?？梢钥闯?，在軸線附近磁場大小均勻且沿x方向。圖2 亥姆霍茲線圈軸線附近bx在x-y平面上的分布

15、及矢量場例證3設帶電粒子質(zhì)量為m，帶電量為，電場強度e沿方向，磁感應強度b沿qyz方向. 則帶電粒子在均勻電磁場中的運動微分方程為 (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)用運行上述程序，例如，取e=4, b=8可得下列特解并給出圖運行結(jié)果如圖3所示。 研究時可以采用不同的初始條件和不同的參量觀察不同的現(xià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ù)運行結(jié)果例證4如右圖所示，求垂直于無限長載流直導線的平面內(nèi)磁感應強度的分布。解題分析設場點p的位置為，電流元位置為，電流元矢量為。由此，場點p相對于電流元的位置矢量為 利用行列式計算idlr ，可寫為

18、 （13）也可利用matlab中的det 命令函數(shù)來求該行列式，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的積運行結(jié)果為d1cr = dy*(z0-z)-dz*(y0-y), dz*(x0-x)-dx*(z0-z), dx*(y0-y)-dy*(x0-x)即又，r的大小為 設載流導體通過坐標原點垂直于 x-y平面放置，電流元 idl沿z軸正向，場點p位于x-y平面上。對本題目而言，dx=dy=0,x=y=0, z0=0, 矢量叉

19、乘積為，r的大小為 由畢奧薩伐爾定律 （14）有； ； （15） matlab程序 1、 用符號運算求b的表達式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 運行結(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、繪制磁場大小分布圖和矢量場圖x=-0.5:0.05:0.5;y=x; %賦值i=input(請輸入電流i=); %設置輸入mu0=4*pi*1e-7; c0=mu0/(4*pi); %設置函數(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號區(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號區(qū)mesh(x,y,b) %繪圖運行該程序，在命令窗中的提示后鍵入i0值 (例如，取i0100a)，便得到圖4所示圖形。圖4長載流導線的磁場在x-y平面上的分布三、結(jié)論從以上利用matlab語言對幾種電磁場模型的分析我們不難的出以下結(jié)論：電磁場與電磁波理論作為電子信息類專業(yè)的一門重要