Matlab在靜電場中的應用

1、摘 要 本文介紹了用Matlab軟件模擬點電荷電場的方法，詳細介紹了用Matlab模擬單個點電荷、一對點電荷的電勢線和電場線的方法。然后，利用鏡像法求解靜電場，用MATLAB將解出的電場可視畫，并討論了感應電荷是如何影響實際電場的，以此，更加深刻的理解導體表面的感應電荷對求解空間電場的影響。關鍵詞：MATLAB 靜電場 鏡像法引言電磁場是物質世界的重要的組成部分之一，在生產實踐和科學技術領域內，存在著大量和電磁場有關的問題，例如，電力系統(tǒng)，凝聚態(tài)物理，天體物理，粒子加速器等，都涉及不少宏觀電磁場的理論問題。因此掌握電磁場的基本理論對生產實踐和科學實驗都有重大意義。現(xiàn)將電磁場的基本理論應用到最簡

2、單的情況：電荷靜止，相應的電場不隨時間而變化的情況。研究的主要問題是：在給定自由電荷分布以及周圍空間介質和導體分布的情況下，電場的分布問題。由Maxwell方程組我們知道靜電勢滿足的微分方程，并且利用唯一性定理可知在給定自由電荷分布和邊界條件的情況下求解空間內的電場唯一確定。接下來，就是在區(qū)域內對方程進行求解。常用的計算電磁場問題的方法主要有兩大類，一是解析法，二是數(shù)值法。對于那些具有簡單邊界條件和幾何形狀的問題，可用分離變量法，鏡像法和格林函數(shù)法求解電磁場邊值問題的解析解。本文首先介紹了用Matlab軟件模擬點電荷電場的方法，詳細介紹了用Matlab模擬單個點電荷、一對點電荷的電勢線和電場線

3、的方法。然后，利用鏡像法求解靜電場，用MATLAB將解出的電場可視畫，并討論了感應電荷是如何影響實際電場的，以此，更加深刻的理解導體表面的感應電荷對求解空間電場的影響。1 Matlab模擬點電荷電場的分布1.1 原理根據(jù)庫倫定律：在真空中，兩個靜止點電荷之間的作用力與這兩個電荷的電量乘積成正比，與它們之間距離的平方成反比，作用力的方向在兩個電荷的連線上，兩電荷同號為斥力，異號為吸力，它們之間的力F滿足： (1.1) 由電場強度E的定義可知： (1.2)對于點電荷，根據(jù)場論基礎中的定義，有勢場E的勢函數(shù)為 (1.3) (1.4)在Matlab中，由以上公式算出各點的電勢U，電場強度E后，可以用M

4、atlab自帶的庫函數(shù)繪出相應電荷的電場分布情況。1.2 單個點電荷真空中點電荷的場強大小是E=kq /r2，其中k為靜電力恒量，q為電量，r為點電荷到場點P(x,y)的距離。電場呈球對稱分布, 取電量q> 0，電力線是以電荷為起點的射線簇。以無窮遠處為零勢點，點電荷的電勢為U=kq /r，當U取常數(shù)時，此式就是等勢面方程。等勢面是以電荷為中心以r為半徑的球面。1.2.1 單個點電荷平面電場線的畫法在平面上，電力線是等角分布的射線簇，用Matlab畫射線簇很簡單。取射線的半徑為(都取國際制單位) r0=0.18，不同的角度用向量表示(單位為弧度) th=linspace(0,2*pi,1

5、5)。射線簇的終點的直角坐標為：x,y=pol2cart(th,r0)。插入x 的起始坐標x=x; 0.1*x。同樣插入y的起始坐標，y=y; 0.1*y，x和y都是二維數(shù)組，每一列是一條射線的起始和終止坐標。用二維畫線命令plot(x,y)就畫出所有電力線。1.2.2 單個點電荷平面等勢線的畫法在過電荷的截面上，等勢線就是以電荷為中心的圓簇，用MATLAB 畫等勢線更加簡單。靜電力常量為k=9e9，電量可取為q=1e- 9；最大的等勢線的半徑應該比射線的半徑小一點，r0=0.1，其電勢為u0=k8q /r0。如果從外到里取7 條等勢線，最里面的等勢線的電勢是最外面的3倍，那么各條線的電勢用向

6、量表示為：u=linspace(1,3,7)*u0.從- r0 到r0取偶數(shù)個點，例如100個點，使最中心點的坐標繞過0，各點的坐標可用向量表示：x=linspace(- r0,r0,100)，在直角坐標系中可形成網(wǎng)格坐標：X,Y=meshgrid(x)。各點到原點的距離為：r=sqrt(X.2+Y.2)，各點的電勢為U=k8q. /r。用等高線命令即可畫出等勢線contour(X,Y,U,u)。以下是單個點電荷的平面電場線和等勢線的Matlab代碼和仿真結果。q=1e-9;r0=0.18;th=linspace(0,2*pi,15);x,y=pol2cart(th,r0);x=x;0.1*x

7、;y=y;0.1*y;plot(x,y);grid onhold onplot(0,0,'o','MarkerSize',12)xlabel('x','fontsize',16)ylabel('y','fontsize',16)title('單個點電荷的電場線與等勢線','fontsize',16)k=9e9;r0=0.12;u0=k*q/r0;u=linspace(1,3,7)*u0;x=linspace(-r0,r0,100);X,Y=meshgrid(x);r=

8、sqrt(X.2+Y.2);U=k*q./r;hold on;contour(X,Y,U,u)圖1.1單個點電荷的平面電場線和等勢線如上圖1.1所示，可以很清楚的看出，越靠近點電荷的中心，電勢越高，電場強度越大，電力線和等勢線也越密。1.3 一對點電荷 同理畫出一對點電荷的電場線和等勢線，其Matlab程序和仿真結果如下：%同號點電荷對的電場線和等勢線clear q=0.5; xm=2.5; ym=2; x=linspace(-xm,xm); %橫坐標向量y=linspace(-ym,ym); %縱坐標向量X,Y=meshgrid(x,y); R1=sqrt(X+1).2+Y.2); %第一個

9、正電荷到場點的距離R2=sqrt(X-1).2+Y.2); %第二個正電荷到場點的距離U=1./R1+q./R2; %計算電勢u=1:0.5:4; %等勢線的電勢向量figure contour(X,Y,U,u) %畫等勢線grid on hold on plot(-xm;xm,0;0) %畫水平線plot(0;0,-ym;ym) %畫豎直線plot(-1,0,'o','MarkerSize',12) plot(1,0,'o','MarkerSize',12) Ex,Ey=gradient(-U,x(2)-x(1),y(2)-y(

10、1); %用電勢梯度求場強的兩個分量dth1=20; th1=(dth1:dth1:180-dth1)*pi/180; %電場線的起始角度r0=0.1; x1=r0*cos(th1)-1; %電場線的起點橫坐標y1=r0*sin(th1); %電場線的起點縱坐標streamline(X,Y,Ex,Ey,x1,y1) %畫左上電場線streamline(X,-Y,Ex,-Ey,x1,-y1) %畫左下電場線dth2=dth1/q; %右邊電場線角度間隔th2=(180-dth2:-dth2:dth2)*pi/180; x2=r0*cos(th2)+1; %電場線的起點橫坐標y2=r0*sin(t

11、h2); %電場線的起點縱坐標streamline(X,Y,Ex,Ey,x2,y2) %畫右上電場線streamline(X,-Y,Ex,-Ey,x2,-y2) %畫右下電場線axis equal tight title('同號點電荷的電場線和等勢線','fontsize',16)xlabel('x','fontsize',16) ylabel('y','fontsize',16) txt='電荷比:itQrm_2/itQrm_1=' num2str(q);text(-xm,-ym+

12、0.5,txt,'fontsize',16) 圖1.2圖1.3(1)如圖1.2所示，電場線從正電荷出發(fā)，終止在無窮遠處。電場線與等勢線垂直，任何兩條電場線都不相交。(2)電勢較高的等勢線分別包圍著電荷，電勢較低等勢線包圍著兩個電荷。電場強度大的地方，電場線較密，等勢線也較密。(3)當兩個電荷的電量相等時，電場線和等勢線對中垂線是對稱的。當兩個電荷的電量不相等時，電場線和等勢線對中垂線都不對稱，如圖1.3所示，但是電場線與等勢線仍然垂直。2 鏡像法點電荷Q附近有一導體，在點電荷的電場作用下，導體面上出現(xiàn)感應電荷。而真實的電場是自由電荷和感應電荷所激發(fā)的電場的疊加。設想感應電荷對空

13、間中電場的影響能否用導體內的一個或幾個假想電荷代替，由唯一性定理可知只要我們能利用假象電荷構造出相同的邊界條件，在求解空間中的電場是相同的。假設接地的導體平面上有一個半徑為a的半球突起，球心在導體平面上，一自由電荷在上半空間任意位置，求電場。為此，我門專門做了一個函數(shù)lilei,如下：function =lilei(a,r0,phi0,q)%在接地的導體平面上有一半徑為a的半球凸部，半球的球心在導體平面%上,電電荷4*pi*p*q(p為真空電介常數(shù))在導體上半空間。%本函數(shù)以過電荷和球心并垂直于導體平面的平面上%畫出電勢線和電場%lilei(a,r0,phi0,q)%a為導體球的半徑%(r0,

14、phi0)為電荷坐標r0>a,0<phi0<pi%q為電荷量4*pi*p*q if nargin=4 disp('請輸入a,ro,phi0,q')elseif a<=0 disp('a>0');elseif r0<=a disp('r0>a')elseif or(phi0>=pi,phi0<=0) disp('phi0>=pi,phi0<=0')elsex,y=meshgrid(-(2*a):0.01:2*a,0:0.01:4*a);Q,R=cart2pol(x,y

15、);R(R<=a)=NaN;ar=a/r0; figure(1) subplot(221)hold onu1=q./sqrt(r02-2*r0*R.*cos(abs(Q-phi0)+R.2);contour(x,y,u1,-1:20,300);ex,ey=gradient(-u1);t=0:pi/10:2*pi;sx=r0*cos(phi0)+0.1*cos(t);sy=r0*sin(phi0)+0.1*sin(t);streamline(x,y,ex,ey,sx,sy)axis equal;tt=0:pi/30:pi;plot(a*exp(i*tt),'r')plot(

16、-3*a:0.1:-a,zeros(size(-3*a:0.1:-a),'r')plot(a:0.1:3*a,zeros(size(a:0.1:3*a),'r')title('孤立電荷產生的電場線和電勢線')hold offsubplot(222)u2=-ar*q./sqrt(a*ar)2-2*a*ar.*R.*cos(Q-phi0)+R.2);contour(x,y,u2,20);hold onex,ey=gradient(-u2);axis equal;tt=0:pi/30:pi;plot(a*exp(i*tt),'r')pl

17、ot(-3*a:0.1:-a,zeros(size(-3*a:0.1:-a),'r')plot(a:0.1:3*a,zeros(size(a:0.1:3*a),'r')k,l=pol2cart(phi0,a*ar);plot(k,l,'or')text(k+0.1,l,'象電荷1')title('象電荷1產生的電場線和電勢線')axis(-3*a,3*a,-0.3,-0.3+6*a)hold off subplot(223)hold onu3=-q./sqrt(r02-2*r0*R.*cos(abs(Q+phi0)

18、+R.2);contour(x,y,u3,20);ex,ey=gradient(-u3);t=0:pi/10:2*pi;sx=r0*cos(-phi0)+0.1*cos(t);sy=r0*sin(-phi0)+0.1*sin(t);streamline(x,y,ex,ey,sx,sy)axis equal;tt=0:pi/30:pi;plot(a*exp(i*tt),'r')plot(-3*a:0.1:-a,zeros(size(-3*a:0.1:-a),'r')plot(a:0.1:3*a,zeros(size(a:0.1:3*a),'r')k

19、,l=pol2cart(-phi0,r0);plot(k,l,'or')text(k+0.1,l,'象電荷2');title('象電荷2產生的電場線和電勢線')axis(-3*a,3*a,-r0*sin(phi0)-0.3,-r0*sin(phi0)-0.3+6*a)hold offsubplot(224)u4=ar*q./sqrt(a*ar)2-2*a*ar.*R.*cos(Q+phi0)+R.2);contour(x,y,u4,20);hold onex,ey=gradient(-u4);axis equal;tt=0:pi/30:pi;pl

20、ot(a*exp(i*tt),'r')plot(-3*a:0.1:-a,zeros(size(-3*a:0.1:-a),'r')plot(a:0.1:3*a,zeros(size(a:0.1:3*a),'r')k,l=pol2cart(-phi0,a*ar);plot(k,l,'or')text(k+0.1,l,'象電荷3');axis(-3*a,3*a,-a*ar*sin(phi0)-0.3,-a*ar*sin(phi0)-0.3+6*a)title('象電荷3產生的電場線和電勢線')hold off figure(2)u=u1+u2+u3+u4;c,h=contour(x,y,u, 0.01 0.1 0.4 1 2 5 9 ,300);hold onset(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)ex,ey=gradient(-u);t=0:pi/10:2*pi;sx=r0*cos(phi0)+0.1*cos(t);sy=r0*sin(phi0)+0.1*sin(t);streamline(x,y,ex,e