MATLAB在静电场中的应用

MATLAB在静电场中的应用何凯摘要：利用数值和解析的方法求解静电场，并利用MATLAB将解出的电场可视画。数值法主要应用有限元法（偏微分方程工具箱）和差分法，解析法主要应用镜像法，并讨论了想电荷是如何影响实际电场的，这其实也就是导体表面的感应电荷对求解空间电场的影响。最后讨论了导体上半空间和球外空间的格林函数。关键词：MATLAB静电场一引言电磁场是物质世界的重要的组成部分之一，在生产实践和科学技术领域内，存在着大量和电磁场有关的问题，例如电力系统，凝聚态物理，天体物理，粒子加速器等，都涉及不少宏观电磁场的理论问题。因此掌握电磁场的基本理论对生产实践和科学实验都有重大意义。1864年Maxwell把电磁规律总结为Maxwell方程组。这是电磁理论中最基本的关系。现在我们将电磁场的基本理论应用到最简单的情况：电荷静止，相应的电场不随时间而变化的情况。我们研究的主要问题是：在给定自由电荷分布以及周围空间介质和导体分布的情况下，电场的分布问题。由Maxwell方程组我们知道静电势满足的微分方程为：并且利用唯一性定理可知在给定自由电荷分布和边界条件的情况下求解空间内的电场唯一确定。下面的问题就是在区域内对方程进行求解。常用的计算电磁场问题的方法主要有两大类，一是解析法，二是数值法。对于那些具有简单边界条件和几何形状的问题，可用分离变量法，镜像法和格林函数法求解电磁场边值问题的解析解。但是在许多实际问题中往往由于边界条件过于复杂而无法求解，在这种情况下，一般用数值法求解。二解析法—镜像法设电点荷Q附近有一导体。在电电荷的电场作用下，导体面上出现感应电荷。而真实的电场是自由电荷和感应电荷所激发的电场的叠加。我们设想感应电荷对空间中电场的影响能否用导体内的一个或几个假想电荷代替。由唯一性定理可知只要我们能利用假象电荷构造出相同的边界条件，在求解空间中的电场是相同的。例1接地的导体平面上有一个半径为a的半球突起，球心在导体平面上，一自由电荷在上半空间任意位置，求电场。为此，我门专门做了一个函数hekai,如下：function[]=hekai(a,r0,phi0,q)%在接地的导体平面上有一半径为a的半球凸部，半球的球心在导体平面%上,电电荷4*pi*p*q(p为真空电介常数)在导体上半空间。%本函数以过电荷和球心并垂直于导体平面的平面上%画出电势线和电场%hekai(a,r0,phi0,q)%a为导体球的半径%(r0,phi0)为电荷坐标r0a,0phi0pi%q为电荷量4*pi*p*qifnargin~=4disp('请输入a,ro,phi0,q')elseifa=0disp('a0');elseifr0=adisp('r0a')elseifor(phi0=pi,phi0=0)disp('phi0=pi,phi0=0')else[x,y]=meshgrid(-(2*a):0.01:2*a,0:0.01:4*a);[Q,R]=cart2pol(x,y);R(R=a)=NaN;ar=a/r0;figure(1)subplot(221)holdonu1=q./sqrt(r0^2-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)axisequal;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')title('孤立电荷产生的电场线和电势线')holdoffsubplot(222)u2=-ar*q./sqrt((a*ar)^2-2*a*ar.*R.*cos(Q-phi0)+R.^2);contour(x,y,u2,20);holdon[ex,ey]=gradient(-u2);axisequal;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,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])holdoffsubplot(223)holdonu3=-q./sqrt(r0^2-2*r0*R.*cos(abs(Q+phi0))+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)axisequal;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,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])holdoffsubplot(224)u4=ar*q./sqrt((a*ar)^2-2*a*ar.*R.*cos(Q+phi0)+R.^2);contour(x,y,u4,20);holdon[ex,ey]=gradient(-u4);axisequal;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,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产生的电场线和电势线')holdofffigure(2)u=u1+u2+u3+u4;[c,h]=contour(x,y,u,[[0.010.10.41259],300]);holdonset(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,ey,sx,sy)axisequal;tt=0:pi/30:pi;plot(a*exp(i*tt),'b')plot(-3*a:0.1:-a,zeros(size(-3*a:0.1:-a)),'b')plot(a:0.1:3*a,zeros(size(a:0.1:3*a)),'b')colormapflagtitle('真实电场')holdoffend以hekai(1,2,pi/2,1)的运行结果为：由结果可以明显的看出真实场可有几个假想电荷和自由电荷叠加而成。以上是一个电场的刨面图，下面求体的格林函数。例2：真空中有一半径为R0的接地导体球，距球心为a处有一点电荷Q，求空间电场。作程序如下：clearclcR0=1;a=1.5;q=50;[x,y,z]=meshgrid(-2:0.03:2,-2:0.03:2,-1:0.03:1.6);[t,p,R]=cart2sph(z,y,x);R(R=1)=NaN;ar=R0/a;clearxyzu2=-ar*q./sqrt(R.^2+(R0*ar)^2-2.*R.*cos(p).*(R0*ar).*cos(t));u1=q./sqrt(a^2+R.^2-2*R.*cos(p).*a.*cos(t));u=u1+u2;clearu1u2cleartpR[x,y,z]=meshgrid(-2:0.03:2,-2:0.03:2,-1:0.03:1.6);isosurface(x,y,z,u,5)holdongridonisosurface(x,y,z,u,15)isosurface(x,y,z,u,25)isosurface(x,y,z,u,35)isosurface(x,y,z,u,45)daspect([111])view(3);axistightcamlightlightinggouraudsphere(16);holdoff图中的彩色面是等势面数值从下向上为：5v,15v,25v,35v,45v。底下带格子的球体就是导体球，半径为1，球心坐标为（0，0，0），自由电荷坐标为（0，0，1.5）。等势面之间的间距代表场强的大小，可知越靠近自由电荷的地方场强越大。起初，在做这个电场时是希望用streamline流线函数计算电场的电场线，但由于计算空间电场线的数据量太大，所以只好用等势面表征电场的特性。程序中的clear…也是为了释放内存空间。三数值计算在电磁场边值问题的分析计算中，场域边界形状及介质分布比较复杂的问题，用解析方法无法分析求解，而用图解法作场图又相当困难，且一般情况下精度较底，难以满足工程要求，而必须使用数值计算方法来求解。有限差分法和有限单元法，原则上都是将一个连续域离散化成有限个分区，利用计算机求数值解，理论上可达到任意要求的精度。1有限差分法利用有限差分法计算如下电场，横截面为矩形的无限长槽由三块接地导体构成，槽盖接电压50V。程序如下：clearclchx=17;hy=11;v1=ones(hy,hx);v1(hy,:)=ones(1,hx)*50;v1(1,:)=zeros(1,hx);fori=1:hyv1(i,1)=0;v1(i,hx)=0;endv2=v1;maxt=1;t=0;k=0;while(maxt1e-6)k=k+1;maxt=0;fori=2:hy-1forj=2:hx-1v2(i,j)=(v1(i,j+1)+v1(i+1,j)+v2(i-1,j)+v2(i,j-1))/4;t=abs(v2(i,j)-v1(i,j));if(tmaxt)maxt=t;endendendv1=v2;endsubplot(121)mesh(v2)axis([0,17,0,11,0,50])subplot(122)contour(v2,15)holdonx=1:1:hx;y=1:1:hy;[xx,yy]=meshgrid(x,y);[gx,gy]=gradient(v2,1,1);axis([-1.5,hx+2.5,-2,13])plot([1,1,hx,hx,1],[1,hy,hy,1,1],'k');2有限元法在MATLAB中提供一个图形界面的偏微分方程工具箱来进行数值求解，它使用的就是有限元法。我们利用这个工具来考察两个问题。一是考察镜像法在实际中的近似程度。因为我们知道镜像法求得的是无限大导体平面上半空间的精确解，而实际中并没有无限大的导体。首先看下模拟的由自由电荷和像电荷激发的电场f