单匝线圈的磁感应强度模拟

系统：Windows 10 专业版 1903
matlab版本：Windows版 2018a

使用说明

详细的推导过程见这篇文章，其中magnetic.m用来计算磁场分布，Single_Coil_3D.m 用来画三维的磁力线， Single_Coil_2D.m 用来画二维的磁力线。使用时需要保证magnetic.m和 Single_Coil_3D.m或 Single_Coil_2D.m 在同一个文件夹下。在近环处，可以利用极限求出磁感应强度，程序中没有涉及。

程序

magnetic.m

%c0=u0/2pi,R:线圈半径（m），I：电流（A）
%Bx,By,Bz：磁感应强度在三个坐标上的分量
function [Bx,By,Bz]=magnetic(c0,R,I,x,y,z)
rho=sqrt(x^2+y^2);
theta=atan2(y,x);      %tan(y/x)=theta,theta为幅度值   
k2=4*R*rho/((R+rho)^2+z^2);	%椭圆函数中的参数k方
k2=k2-eps;                  
[K,E]=ellipke(k2);
dd=(R-rho)^2+z^2;
dd=dd+eps; %防止分母为零
B_rho=c0*I*z/((rho+1e-3)*sqrt((R+rho)^2+z^2))*((R^2+rho^2+z^2)/dd*E-K);
B_z=c0*I/sqrt((R+rho)^2+z^2)*((R^2-rho^2-z^2)/dd*E+K);
B_theta=0;
Bx=B_rho*cos(theta)-B_theta*sin(theta);
By=B_rho*sin(theta)+B_theta*cos(theta);
Bz=B_z;
end

Single_Coil_3D.m

tic;clear;clc;
%I:电流(A),c0=u0/2pi,R:线圈半径(m),dl:步长
I=1e6;c0=2e-7;R=0.5;dl=0.001;x=[];y=[];z=[];
r=[0.1;0.15;0.2;0.25;0.3;0.35;0.4;0.45];
%用这么多点画磁力线的一半
n=[510;525;545;600;800;880;500;220];
for theta=0:30:360      %在不同theta角画磁力线
for j=1:length(r)   %在某个角度下，磁力线条数
x(1)=r(j)*cos(theta*pi/180);y(1)=r(j)*sin(theta*pi/180);z(1)=0; %磁力线的起点
for pm=[-1,1] %磁力线是对称的，先画上半部分，在画下半部分 
for i=1:n(j)   
[Bx,By,Bz]=magnetic(c0,R,I,x(i),y(i),z(i));
B=sqrt(Bx^2+By^2+Bz^2);
x(i+1)=x(i)+dl*Bx/B*pm;
y(i+1)=y(i)+dl*By/B*pm;
z(i+1)=z(i)+dl*Bz/B*pm;
end
plot3(x,y,z);hold on;box on;
end
x=[];y=[];z=[]; 
end
end
% 画线圈
h=0; t=0:0.1:(2*pi);t=[t,0];
plot3(0+0.5*sin(t),0+0.5*cos(t),h*ones(size(t)),'LineWidth',4)
xlabel(['\fontname{Arial}x']);
ylabel(['\fontname{Arial}y']);
zlabel(['\fontname{Arial}z']);
set(gca,'FontSize',20);
disp(['计算时间' num2str(toc) 's']);
%写入到文件
time=datestr(now,30);
print(gcf,['Single_Coil_3D_',time,'.eps'],'-depsc','-r600')
savefig(['Single_Coil_3D',time,'.fig'])
print(gcf,['Single_Coil_3D',time,'.png'],'-dpng')
disp(['耗时' num2str(toc) 's']);

得到的结果

Single_Coil_2D.m

tic;clc;clear;
%x-z平面磁力线分布，当y=0,\rho=x>0;Bx=B\rho
%c0=u0/2pi,I:线圈电流(A),R:线圈半径(m)
c0=2e-7;I=1e6;R=0.5;
[x,z]=meshgrid(0:0.01:1.0);  %网格划分
k2=4*R*x./((R+x).^2+z.^2);k2=k2-eps*(k2==1);
[K,E]=ellipke(k2);  
dd=(R-x).^2+z.^2;dd=dd+eps*(dd==1); %防止分母为零
Bx=c0*I*z./(x+eps*(x==0))./sqrt((R+x).^2+z.^2).*((R^2+x.^2+z.^2)./dd.*E-K);
Bz=c0*I*1./sqrt((R+x).^2+z.^2).*((R^2-x.^2-z.^2)./dd.*E+K);
px=[0:R/50:R];
pz=zeros(1,length(px)); %磁力线起点
streamline(x,z,Bx,Bz,px,pz);hold on;
streamline(-x,z,-Bx,Bz,-px,pz)
streamline(x,-z,Bx,-Bz,px,-pz)
streamline(-x,-z,-Bx,-Bz,-px,-pz)
hold on;
xlabel(['\fontname{Arial}x']);
ylabel(['\fontname{Arial}z']);
set(gca,'FontSize',20);
plot([-R R],[0 0],'r.','markersize',20);
disp(['计算时间' num2str(toc) 's']);
%写入到文件
time=datestr(now,30);
print(gcf,['Single_Coil_2D',time,'.eps'],'-depsc','-r600')
savefig(['Single_Coil_2D',time,'.fig'])
print(gcf,['Single_Coil_2D',time,'.png'],'-dpng')
disp(['耗时' num2str(toc) 's']);

得到的结果

源码

欢迎转载，不需注明出处，就说是你写的