电磁场仿真——绘制电场线和等势线

空间中有三个电荷量分别为+1C，-1C，+1C的点电荷，三个点电荷在边长为10mm的正三角形顶点处，绘制电场线的空间二维分布。

一：基本原理

在空间二维平面中建立二维坐标系，如图1所示在（5，0），（-5，0）和（0，5）处分别有1C的正电荷，1C的正电荷和-1C的负电荷。

（一）电势原理

根据空间中某一点距离点电荷的电势知道：

U=kqr

根据空间某点处电势的叠加原理，在场点P（x，y）处产生的电势U为：

 

其中，

知道了空间某一点P（x，y）处三个点电荷产生的电势的叠加值，我们就可以找出某一电势所对应的多个空间点，将这些点连接成线即可以得到等势线。

（二）电场原理

根据库仑定律知道，空间某一点相对某一位置的电荷的电场强度为 :

根据点电荷电场强度的叠加原理，在场点P(x,y)处产生的电场强度E为：

通过计算可以知道：

。

大小分别为1C,1C,-1C。

知道了空间某一点处的各分量的电场强度的大小，我们就推导出下一段电场线上的点的各分量电场强度大小。即满足 在模型中我们可以将各个电荷看成一定半径的球体，按照一定的间距将每一个球体表面圆周分成若干个点，这些点就是电场线的起始点。通过从各个起始点取一定间隔逐次推导出下一个坐标点的位置，知道满足一定的条件停止推导，将这些点连线即可得到电场线。

二：仿真数据设置

将x的范围设置成[-2.5,2.5]/cm。其中取样间隔为0.01cm，将y的范围同样设置成[-2.5,2.5]/cm，取样间隔同样为0.01。这样就形成了一个501*501的网格。我们取k介质为真空中，化简得到。在等势线绘制的时候我们把网格内的电势线分成了100份。我们设置各个电荷的半径0.05*2^(1/2)m,各个电荷表面有50个起始点。起始点到下一个推导电场点的距离间隔为0.001。最终电场线终止要满足的条件有：

1.下一个电场点靠近其他电荷。

2.下一个电场点的位置超出了网格点的边界。

三．结果分析及总结

最终仿真图像结果如下：

四：MATLAB代码运行

clear,number=3;        %设定空间电荷数，并用charge矩阵来保存各个电荷的电荷量和位置 
charge=[1,0.5,0;1,-0.5,0;-1,0,0.86]
k=9*10^9;%这里取介质为无限大真空，化简得到k=9*10^9
x=-2.5:0.01:2.5;
y=-2.5:0.01:2.5;
[X,Y]=meshgrid(x,y); %建立X_Y坐标系
U=zeros(501,501);
clf  %clf指令是清除之前的图像
for j=1:number
    U=U+k*charge(j,1)./sqrt((X-charge(j,2)).^2+(Y-charge(j,3)).^2);
end%得到电势的值 注意U与XY的对应关系是U[Y,X]
for a=1:number
for b=-10:10 %横坐标变化
for c=-10:10 %纵坐标变化
if abs(U(charge(a,3)*100+251+b,charge(a,2)*100+251+c))>abs(U(charge(a,3)*100+251+5,charge(a,2)*100+251+5))
U(charge(a,3)*100+251+b,charge(a,2)*100+251+c)=U(charge(a,3)*100+251+5,charge(a,2)*100+251+5);
end
end
end
end%由于在电荷处，电势理论上无穷大，所以通过设置阈值对其进行“削顶”
contour(X,Y,U,100)%画等势面
hold on
%下面是对电场线的绘制
total=sum(abs(charge(:,1)));%电荷模值求和
px=[];py=[];%存放电场线数据
for a=1:number%画由第a个电荷发出的电场线
    for b=1:round(abs(charge(a,1))/total*50*number)%电场线数目由归一化电荷决定 其中round函数是四舍五入的意思
        j=round(abs(charge(a,1))/total*50*number);i=1;%j是单个电荷表面电场线起始点的数目
        px(a,b,i)=charge(a,2)+0.05*2^(1/2)*cos(b*2*pi/j);
        py(a,b,i)=charge(a,3)+0.05*2^(1/2)*sin(b*2*pi/j);%分配电场线起点坐标 其中ab表示的是第a个电荷的第b个起始点的位置
        t=1;
        while t==1
            Ex=0;Ey=0;
            for c=1:number
                Rx=px(a,b,i)-charge(c,2);
                Ry=py(a,b,i)-charge(c,3);
                E=k*charge(c,1)./(Rx^2+Ry^2);
                Ex=Ex+E*Rx/sqrt((Rx^2+Ry^2));
                Ey=Ey+E*Ry/sqrt((Rx^2+Ry^2));%得到线端点电出场的X,Y方向上的大小
            end
            rx=0.001*Ex/sqrt(Ex^2+Ey^2);
            ry=0.001*Ey/sqrt(Ex^2+Ey^2);%得到下一步要连接的点跨度
            if charge(a,1)>0
            else
                rx=-rx;
                ry=-ry;
            end              %判断切线的方向
            px(a,b,i+1)=px(a,b,i)+rx;
            py(a,b,i+1)=py(a,b,i)+ry;
            i=i+1;
            for c=1:number
                if ((((px(a,b,i)+rx-charge(c,2))^2+(py(a,b,i)+ry-charge(c,3))^2)>0.005)&&(px(a,b,i)+rx<2.5)&&(py(a,b,i)+ry<2.5)&&(px(a,b,i)+rx>-2.5)&&(py(a,b,i)+ry>-2.5))
                    %上面这一行是判断下一个连接点即没有出坐标范围也没连接到另一个电荷
                else 
                    t=0;
                end
            end
        end
    end
end%我们终于得到所有电场线的信息了
hold on
for a=1:number
    for b=1:round(abs(charge(a,1))/total*50*number)
        t=1;
        hold on
        pxp=squeeze(px(a,b,:));
        pyp=squeeze(py(a,b,:));
        pxp(find(pxp==0))=[];
        pyp(find(pyp==0))=[];
       for c=1:number
        if ((((pxp(end)-charge(c,2))^2+(pyp(end)-charge(c,3))^2)<=0.01)&&charge(a,1)<0)
            t=0;
        else
            hold on
        end
        end
        if t==1
        plot(pxp,pyp)
        end
    end
end%一根根的画出来
xlabel('X'),ylabel('Y'),title('电场线等势面绘制图'),text(1.8,2.3,'电荷量单位 C')
text(2.1,-2.35,'坐标单位cm')
text(1.9,-2.6,date)
for a=1:number
    text(charge(a,2)-0.04,charge(a,3),num2str(charge(a,1)*1))
end