矩形磁铁 磁场计算 matlab源代码

%This code was created by Xiying Wang in 10/01/2011

%calculate magnetic field of a rectangular magnet****%

% **x_dir_len：the width of magnet****** %

% **y_dir_len：the height of magnet***** %

% **z_dir_len：the thickness magnet***** %

% **x,y,z:the magnetic field point****** %

% **the zero point of coordinat is %

% **located in the mass center of magnet %

% ************************************** %

u0 = 4*pi*10^(-7);

M=11.8611e5;%M=m/V m is the magnetic moment and V is the volume of magnet

m=M*x_dir_len*y_dir_len*z_dir_len;

A=x_dir_len*y_dir_len;

I=m/A;

zstep=0.1e-3;

%----------------------------------------

z=z-(-1*z_dir_len/2:zstep:z_dir_len/2);

%----------------------------------------

%B field of y direction wire(left wire)

d=sqrt((x_dir_len/2+x)^2+z.^2);

B_y_dir1=abs(u0*((y_dir_len/2+y)./sqrt((y_dir_len/2+y)^2+d.^2)-((y-y_dir_len/2)./sqrt((y-y_dir_len/2)^2+d.^2)))./(4*pi*d));

B_x1=abs((B_y_dir1.*z)./d);

B_z1=abs(B_y_dir1.*(x_dir_len/2+x)./d);

%B field of y direction wire(right wire)

d=sqrt((x_dir_len/2-x)^2+z.^2);

B_y_dir2=abs(u0*((y_dir_len/2+y)./sqrt((y_dir_len/2+y)^2+d.^2)-((y-y_dir_len/2)./sqrt((y-y_dir_len/2)^2+d.^2)))./(4*pi*d));

B_x2=abs((B_y_dir2.*z)./d);

B_z2=abs(B_y_dir2.*(x_dir_len/2-x)./d);

%B field of x direction wire(bottom wire)

d=sqrt((y_dir_len/2+y)^2+z.^2);

B_x_dir1=abs(u0*((x_dir_len/2+x)./sqrt((x_dir_len/2+x)^2+d.^2)-((x-x_dir_len/2)./sqrt((x-x_dir_len/2)^2+d.^2)))./(4*pi*d));

B_y1=abs((B_x_dir1.*z)./d);

B_z3=abs(B_x_dir1*(y_dir_len/2+y)./d);

%B field of x direction wire(top wire)

d=sqrt((y_dir_len/2-y)^2+z.^2);

B_x_dir2=abs(u0*((x_dir_len/2+x)./sqrt((x_dir_len/2+x)^2+d.^2)-((x-x_dir_len/2)./sqrt((x-x_dir_len/2)^2+d.^2)))./(4*pi*d));

B_y2=abs((B_x_dir2.*z)./d);

B_z4=abs(B_x_dir2*(y_dir_len/2-y)./d);

Bx=B_x2-B_x1;

By=B_y2-B_y1;

if (x<-x_dir_len/2)

B_z1=-1*B_z1;

end

if (x>x_dir_len/2)

B_z2=-1*B_z2;

end

if (y<-y_dir_len/2)

B_z3=-1*B_z3;

end

if (y>y_dir_len/2)

B_z4=-1*B_z4;

end

Bz=B_z1+B_z2+B_z3+B_z4;

Bx=I*sum(Bx)/length(z);By=I*sum(By)/length(z);Bz=I*sum(Bz)/length(z);

end

可以用下面的代码调用以上计算磁场的代码绘制 磁场强度曲线

clc;close all;clear all;

flag_bx_by_bz=0;

flag_bz1_bz2=0;

flag_2magnet_bz=0;

flag_coil_size_optimization=1;

%%

if flag_bx_by_bz==1

z=16e-3;

x=-10e-3:0.1e-3:10e-3;

y=0e-3;

xl=10e-3;yl=10e-3;zl=30e-3;

Bx=zeros(1,length(x));

By=Bx;Bz=Bx;

for i=1:length(x)

[bx by bz]=RectangularMagnet(xl,yl,zl,x(i),y,z);

Bx(i)=bx;

By(i)=by;

Bz(i)=bz;

end

plot(x,Bx,x,By,x,Bz);legend('x

','y','z');

hold on;

end

%%

if flag_bz1_bz2==1

z=6e-3;

x=-10e-3:0.1e-3:10e-3;

y=0e-3;

xl=10e-3;yl=10e-3;zl=10e-3;

Bz1=zeros(1,length(x));

Bz2=zeros(1,length(x));

for i=1:length(x)

[bx by bz]=RectangularMagnet(xl,yl,zl,x(i),y,z);

Bz1(i)=bz;

[bx by bz]=RectangularMagnet(xl,yl,zl,x(i),2e-3,z);

Bz2(i)=bz;

end

plot(x,Bz1,x,Bz2);legend('z1','z2');

end

%%

%两平行磁铁气隙磁场

if flag_2magnet_bz==1

magnet_w=1*0.0254;magnet_h=0.75*0.0254;magnet_t=0.5*0.0254;

gap=10e-3;

D=0.5*(magnet_t+magnet_t)+gap;%the center distance of two face to face parallel magnets;

x=-30e-3:0.1e-3:30e-3;

y=0;

z1=0.5*magnet_t+0.1*gap;

z2=D-z1;

len=length(x);

Bz1=zeros(1,len);Bz2=Bz1;

for i=1:len

[bx by bz1]=RectangularMagnet(magnet_w,magnet_h,magnet_t,x(i),y,z1);

[bx by bz2]=RectangularMagnet(magnet_w,magnet_h,magnet_t,x(i),y,z1);

Bz1(i)=abs(bz1);

Bz2(i)=abs(bz2);

end

plot(x,Bz1+Bz2,x,Bz1,x,Bz2);

legend('1+2','1','2');

end

%%

%Optimization of coil size. To find out the best thickness of coil for

%given size rectangular magnets

if flag_coil_size_optimization==1

%******************Parameters 1****************************************

magnet_w=25.4e-3;magnet_h=19.05e-3;magnet_t=12.7e-3;%the size of magnet

coil_w=magnet_w;coil_h=13e-3;%the size of coil

useless_coil_length=43.76e-3;%the wires that are not in the magnetic field

coil_velocity=0.6;%m/s f=50hz Amplitude=2mm peak(v)=2*pi*f*Amplitude=6.28*50*2e-3=0.6

resistance_ratio=0.0172;%ohm/m

%******************Parameters 2****************************************

wire_diameter=0.127e-3;%wire diameter

wire_space_ratio=0.3*wire_diameter;%the space between two wires in a coil

layer_num=100;%the number of wire layer

x_dir_grid=100;%devide x direction into 100 elements

magnet_coil_gap=2e-3;%the gap between magnet surface and coil

%**********************************************************************

coil_t=1*0.5*(wire_diameter+wire_space_ratio):...

(wire_diameter+wire_space_ratio):...

(0.5*(wire_diameter+wire_space_ratio)...

+(layer_num-1)*(wire_diameter...

+wire_space_ratio));%the thickness of coil(from inner face to outer face)

x=-coil_w*0.5:(coil_w/100):coil_w*0.5;

lenx=length(x);

y=-coil_h*0.5:(wire_diameter+wire_space_ratio):coil_h*0.5;%the increment of y direction is the diameter of wire plus its space

leny=length(y);

lent=length(coil_t);

wire_number=0;%total wire number

voltage_out=0;%total output voltage

wire_number_per_layer=leny;

voltage_out_total=zeros(1,lent);

power_out_total=zeros(1,lent);

total_resistance=zeros(1,lent);

for i=1:lent

clc

display(sp

rintf('Calculation finished %g%% .....' ,i*100/lent));

z=coil_t(i)+magnet_t*0.5+magnet_coil_gap;

wire_number=wire_number+wire_number_per_layer;

wire_length=wire_number*coil_w;

wire_resistance=(2*wire_number*coil_w+4*i*wire_diameter+useless_coil_length*2)*...

resistance_ratio*(1e-6)*4/(pi*wire_diameter^2);%total resistance of coil

for j=1:leny

for k=1:lenx

[bx by bz]=RectangularMagnet(magnet_w,magnet_h,magnet_t,x(k),y(j),z);

voltage_out=voltage_out+bz*(coil_w/100)*coil_velocity;

end

end

voltage_out_total(i)=4*voltage_out;

power_out_total(i)=(0.5*voltage_out)^2/wire_resistance;

total_resistance(i)=wire_resistance;

end

plot(coil_t.*1000,voltage_out_total./sqrt(2));legend('voltage');

xlabel('thickness (mm)');ylabel('voltage (V)');

figure;

plot(coil_t.*1000,power_out_total./sqrt(2));%RMS power of 4 coils

legend('power');xlabel('thickness (mm)');ylabel('power (w)');

[maxpower t]=max(power_out_total./sqrt(2));

hold on;

plot(coil_t(t)*1000,maxpower,'r*');

text(coil_t(t)*1000,maxpower*1.04,['thickness:' num2str(coil_t(t)*1000) ' mm' ', power:' num2str(maxpower) ' W']);

display(sprintf('Best wire loops: %g' ,round(wire_number_per_layer*coil_t(t)/(wire_diameter+wire_space_ratio))));

figure;

plot(coil_t.*1000,total_resistance);%total resistance

end