|

磁镜的磁场位形模拟

目录
    本文上次更新于 1591 天前,其内容可能已经过时,如果文章内容或图片资源失效,请留言反馈,我会及时处理,谢谢!

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

    使用说明

    详细的推导过程见这篇文章,程序需要用到这篇文章中的 magnetic.m 程序,请将它与下面的程序放在一个文件夹下。 Simple_Mirror_3D.m 用来计算模拟磁镜的三维磁力线分布。Simple_Mirror_2D.m用来模拟磁场的二维磁力线分布,在模拟二维磁力线分布时,线圈半径取1.0m。

    程序

    Simple_Mirror_3D.m

    tic;clc;clear;
    %c0:u0/2pi,I:线圈电流(A),R:线圈半径(m),L:线圈间距(m),dl:步长
    c0=2e-7;I=1e6;L=2.0;R=0.2;dl=0.001;
    x=[];y=[];z=[];
    r=[0.01;0.02;0.05;0.1;0.2;0.3;0.4;0.5];r=[-r;r];
    for theta=0:30:180
        for j=1:length(r)
            for pm=-1:1
                for i=1:600    
                    x(1)=r(j)*cos(theta*pi/180);y(1)=r(j)*sin(theta*pi/180);z(1)=0;%磁力线起点
                    [Bx1,By1,Bz1]=magnetic(c0,R,I,x(i),y(i),z(i)+L/2);%线圈1产生的磁场
                    [Bx2,By2,Bz2]=magnetic(c0,R,I,x(i),y(i),z(i)-L/2);%线圈2产生的磁场
                    Bx=Bx1+Bx2;By=By1+By2;Bz=Bz1+Bz2;
                    B=Bx^2+By^2+Bz^2;
                    x(i+1)=x(i)+Bx/B*dl*pm;
                    y(i+1)=y(i)+By/B*dl*pm;
                    z(i+1)=z(i)+Bz/B*dl*pm;
                end
                subplot(121);plot3(x,y,z,'b-');hold on;box on;
                if theta==0;
                    subplot(122);plot(x,z,'b-');hold on;
                end
            end
        end
    end
    %画线圈1
    h=-1;t=0:0.1:(2*pi);t=[t,0];
    plot3(0+0.2*sin(t),0+0.2*cos(t),h*ones(size(t)),'LineWidth',4)
    %画线圈2
    h=1;t=0:0.1:(2*pi);t=[t,0];
    plot3(0+0.2*sin(t),0+0.2*cos(t),h*ones(size(t)),'LineWidth',4)
    xlabel('x');ylabel('y');zlabel('z');set(gca,'FontSize',20);
    subplot(122);xlabel('x');ylabel('z');set(gca,'FontSize',20);
    disp(['计算时间' num2str(toc) 's']);
    %写入到文件
    time=datestr(now,30);
    print(gcf,['Simple_Mirror_3D_',time,'.eps'],'-depsc','-r600')
    savefig(['Simple_Mirror_3D_',time,'.fig'])
    print(gcf,['Simple_Mirror_3D_',time,'.png'],'-dpng')
    disp(['耗时' num2str(toc) 's']);

    得到的结果

    Simple_Mirror_2D.m

    tic;clc;clear;
    %c0:u0/2pi,I:线圈电流(A),I:线圈距离(m),R:线圈半径(m)
    c0=2e-7;I=1e6;L=2;R=1.0;
    n=100; x1=0;x2=6*R;z1=-1.0; z2=1.0;
    dx = (x2-x1)/n; dz = 2*(z2-z1)/n;
    [x,z]=meshgrid(x1:dx:x2,3*z1:dz:3*z2);
    k2_1=4*R*x./((R+x).^2+(z+L/2).^2);
    k2_1=k2_1-eps*(k2_1==1);
    k2_2=4*R*x./((R+x).^2+(z-L/2).^2);
    k2_2=k2_2-eps*(k2_2==1);
    [K1,E1]= ellipke(k2_1);
    [K2,E2]=ellipke(k2_2);
    dd1=(R-x).^2+(z+L/2).^2;
    dd1=dd1+eps*(dd1==0);
    dd2=(R-x).^2+(z-L/2).^2;
    dd2=dd2+eps*(dd2==0);
    %线圈1产生的磁场
    Bx1=c0*I*(z+L/2)./(x+eps*(x==0))./sqrt((R+x).^2+(z+L/2).^2).*((R^2+x.^2+(z+L/2).^2)./dd1.*E1-K1);
    Bz1=c0*I./sqrt((R+x).^2+(z+L/2).^2).*((R^2-x.^2-(z+L/2).^2)./dd1.*E1+K1);
    %线圈2产生的磁场
    Bx2=c0*I*(z-L/2)./(x+eps*(x==0))./sqrt((R+x).^2+(z-L/2).^2).*((R^2+x.^2+(z-L/2).^2)./dd2.*E2-K2);
    Bz2=c0*I./sqrt((R+x).^2+(z-L/2).^2).*((R^2-x.^2-(z-L/2).^2)./dd2.*E2+K2);
    Bx=Bx1+Bx2;
    Bz=Bz1+Bz2;
    vz =[0,1]; vx=[0:0.2:4*R];
    [px,pz]=meshgrid(vx,vz);
    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);
    plot([-R -R R R],[z1 z2 z1 z2],'r.','markersize',40);
    plot([0 0],[3*z1 3*z2],'b-','linewidth',1);
    xlabel('x');ylabel('z');set(gca,'FontSize',20);
    disp(['计算时间' num2str(toc) 's']);
    %写入到文件
    time=datestr(now,30);
    print(gcf,['Simple_Mirror_2D_',time,'.eps'],'-depsc','-r600')
    savefig(['Simple_Mirror_2D_',time,'.fig'])
    print(gcf,['Simple_Mirror_2D_',time,'.png'],'-dpng')
    disp(['耗时' num2str(toc) 's']);

    得到的结果,这里线圈的半径不是0.2m,而是1.0m了。半径为0.2时画出来的图形比较难看。

    源码

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

    如果在这个过程中遇到了其它问题,欢迎在评论区留言,如果你已解决,也欢迎把具体的解决方法留在评论区,以供后来者参考
    ×

    感谢您的支持,请扫码打赏

    微信打赏 支付宝打赏
    guest

    这个站点使用 Akismet 来减少垃圾评论。了解你的评论数据如何被处理

    0 评论
    内联反馈
    查看所有评论