基于Matlab的电磁场与波“镜像电荷法”仿真——点电荷位于半球形凸起导体附近

1.目的：

  对镜像法的知识点进行应用，举一反三解决实际问题，并用MATLAB将其形象化。

2.理论分析：

  （1）平面镜像法

      对于一无限大导体平面上方一个点电荷的电场，只要在导体平面的下方点电荷q对称的点处放置一点电荷（-q），并把无限大导体平板撤去，整个空间充满介电常数为的电介质，这两个的电荷在上半平面共同产生的电场即为原电场。

（2）球面镜像法

      理论分析见专栏内第一期分享

  （3）本题模型

      本题为一点电荷位于半球形凸起导体附近

 点电荷q在此模型中会产生三个镜像电荷，其中有：

 导体上方电场由这四个电荷共同组成：

3.源代码与结果

3.1二维图源代码：

k=8.9875e9;
q=1.602e-18;
d=0.2;
a=0.1;
b=a*a/d;
q1=-q;
q2=-a/d*q;
q3=-q2;
theta=linspace(0,pi,50);
th=linspace(0,2*pi,25);
xr=a*cos(theta);
yr=a*sin(theta);
hold on
plot(0,d,'o','MarkerSize',6)
plot(xr,yr,'k','linewidth',2)
plot([-0.4,-0.1],[0,0],'k','linewidth',2)
plot([0.1,0.4],[0,0],'k','linewidth',2)
xlabel('X','fontsize',16)
ylabel('Y','fontsize',16)
title('接地半球突起','fontsize',20);
axis([-0.4,0.4,-0.4,0.4]);


%第一部分
xd1=linspace(-0.4,0.4,100);
yd1=linspace(0,0.4,100);
[Xd1,Yd1]=meshgrid(xd1,yd1);
r1=sqrt(Xd1.^2+(Yd1-d).^2);
r2=sqrt(Xd1.^2+(Yd1+d).^2);
r3=sqrt(Xd1.^2+(Yd1-b).^2);
r4=sqrt(Xd1.^2+(Yd1+b).^2);
U=k*q./r1+k*q1./r2+k*q2./r3+k*q3./r4;
u0=k*q./0.01+k*q1/0.39+k*q2/0.14+k*q3/0.24;
u=linspace(-u0,u0,101);
contour(Xd1,Yd1,U,u,'--');
colormap(jet);
hold on
[Ex,Ey]=gradient(-U);
E=sqrt(Ex.^2+Ey.^2);
Ex=Ex./E;
Ey=Ey./E;
x1=0.001*cos(th);
y1=d+0.001*sin(th);
for i=2:24
    h1=streamline(Xd1,Yd1,Ex,Ey,x1(i),y1(i));
    set(h1,'LineWidth',1); 
    pipi=drawEline(h1,184,185);
    set(pipi,'HeadWidth',4,'HeadLength',4,'Color','b')
end
fill(xr,yr,'w');
hold on


%第二部分
xd1=linspace(-0.4,0.4,100);
yd1=linspace(0,0.4,100);
% for i=1:100
%     for j=1:100
%         if xd1(i)^2+yd1(j)^2 <= a^2
%             xd1(i)=nan;
%             yd1(j)=nan;
%         end
%     end
% end
[Xd1,Yd1]=meshgrid(xd1,yd1);
r1=sqrt(Xd1.^2+(Yd1-d).^2);
r2=sqrt(Xd1.^2+(Yd1+d).^2);
r3=sqrt(Xd1.^2+(Yd1-b).^2);
r4=sqrt(Xd1.^2+(Yd1+b).^2);
U=k*q./r1+k*q1./r2+k*q2./r3+k*q3./r4;
u0=k*q./0.01+k*q1/0.39+k*q2/0.14+k*q3/0.24;
u=linspace(-u0,u0,101);
contour(Xd1,Yd1,U,u,'--');
colormap(jet);
hold on
[Ex,Ey]=gradient(-U);
E=sqrt(Ex.^2+Ey.^2);
Ex=Ex./E;
Ey=Ey./E;
x1=0.001*cos(th);
y1=d+0.001*sin(th);
for i=2:24
    h1=streamline(Xd1,Yd1,Ex,Ey,x1(i),y1(i));
    set(h1,'LineWidth',1,'LineStyle','--'); 
%     pipi=drawEline(h1,164,165);
%     set(pipi,'HeadWidth',4,'HeadLength',4,'Color','b')
end

%第三部分
xd2=linspace(-0.4,0.4,100);
yd2=linspace(-0.2,0.2,100);
[Xd2,Yd2]=meshgrid(xd2,yd2);
r1=sqrt(Xd2.^2+(Yd2-d).^2);
r2=sqrt(Xd2.^2+(Yd2+d).^2);
r3=sqrt(Xd2.^2+(Yd2-b).^2);
r4=sqrt(Xd2.^2+(Yd2+b).^2);
U=k*q./r1+k*q1./r2+k*q2./r3+k*q3./r4;
u0=k*q./0.01+k*q1/0.39+k*q2/0.14+k*q3/0.24;
u=linspace(-u0,u0,101);
% contour(Xd2,Yd2,U,u,'--');
colormap(jet);
hold on
[Ex,Ey]=gradient(-U);
E=sqrt(Ex.^2+Ey.^2);
Ex=Ex./E;
Ey=Ey./E;
x2=0.001*cos(th);
y2=-d+0.001*sin(th);
for i=2:24
    h1=streamline(Xd2,Yd2,-Ex,-Ey,x2(i),y2(i));
    set(h1,'LineWidth',1,'LineStyle','--'); 
%     pipi=drawEline(h1,164,163);
%     set(pipi,'HeadWidth',4,'HeadLength',4,'Color','b')
end
x4=0.001*cos(th);
y4=-b+0.001*sin(th);
for i=2:11
    h1=streamline(Xd2,Yd2,Ex,Ey,x4(i),y4(i));
    set(h1,'LineWidth',1,'LineStyle','--'); 
%     pipi=drawEline(h1,80,81);
%     set(pipi,'HeadWidth',4,'HeadLength',4,'Color','b')
end

fill([-0.4,0.4,0.4,-0.4],[0,0,-0.4,-0.4],'w');

%第四部分
xd2=linspace(-0.4,0.4,100);
yd2=linspace(-0.4,0,100);
[Xd2,Yd2]=meshgrid(xd2,yd2);
r1=sqrt(Xd2.^2+(Yd2-d).^2);
r2=sqrt(Xd2.^2+(Yd2+d).^2);
r3=sqrt(Xd2.^2+(Yd2-b).^2);
r4=sqrt(Xd2.^2+(Yd2+b).^2);
U=k*q./r1+k*q1./r2+k*q2./r3+k*q3./r4;
u0=k*q./0.01+k*q1/0.39+k*q2/0.14+k*q3/0.24;
u=linspace(-u0,u0,101);
contour(Xd2,Yd2,U,u,'--');
colormap(jet);
hold on
[Ex,Ey]=gradient(-U);
E=sqrt(Ex.^2+Ey.^2);
Ex=Ex./E;
Ey=Ey./E;
x2=0.001*cos(th);
y2=-d+0.001*sin(th);
for i=2:24
    h1=streamline(Xd2,Yd2,-Ex,-Ey,x2(i),y2(i));
    set(h1,'LineWidth',1,'LineStyle','--'); 
    pipi=drawEline(h1,184,183);
    set(pipi,'HeadWidth',4,'HeadLength',4,'Color','b')
end
x4=0.001*cos(th);
y4=-b+0.001*sin(th);
for i=2:11
    h1=streamline(Xd2,Yd2,Ex,Ey,x4(i),y4(i));
    set(h1,'LineWidth',1,'LineStyle','--'); 
    pipi=drawEline(h1,120,121);
    set(pipi,'HeadWidth',4,'HeadLength',4,'Color','b')
end

text(0.03,d,'q=+1.602e-18C')

3.2二维图结果：

修改前

 修改后

 3.3三维图源代码：

k=8.9875e9;
q=1.602e-18;
d=0.2;
a=0.1;
b=a*a/d;
q1=-q;
q2=-a/d*q;
q3=-q2;
theta=linspace(0,pi,50);
xr=a*cos(theta);
yr=a*sin(theta);
x=linspace(-0.2,0.2,100);
y=linspace(0.1,0.5,100);
[X,Y]=meshgrid(x,y);
r1=sqrt(X.^2+(Y-d).^2);
r2=sqrt(X.^2+(Y+d).^2);
r3=sqrt(X.^2+(Y-b).^2);
r4=sqrt(X.^2+(Y+b).^2);
U=k*q./r1+k*q1./r2+k*q2./r3+k*q3./r4;
u0=k*q./0.01+k*q1/0.39+k*q2/0.14+k*q3/0.24;

u=linspace(-4*u0,4*u0,400);
[Ex,Ey]=gradient(-U);
E=sqrt(Ex.^2+Ey.^2);
Ex=Ex./E;
Ey=Ey./E;
surf(U);
hold on
shading interp;

subplot(2,1,1);
[c,h]=contour3(U,u);
colormap(jet);
h.EdgeColor='c';

subplot(2,1,2);
[c,h]=contour3(U,u);
colormap(white);
h.EdgeColor='w';
% % h.LineStyle='--';
h1=streamslice(Ex,Ey);
set(h1,'Color',[0.1 0.9 0.8]);
for j=1:length(h1)
    ui = interp2(U,h1(j).XData,h1(j).YData);
    h1(j).ZData = ui;
end

3.4三维图结果

修改前

 修改后