promgram for linear case.m

永磁同步电机输出特性模拟程序。此处电机材料为线性。

%N=input(' input coil turns = ');
 %Np=input(' Pairs_of_Poles_Number= ');
 %D=input(' Motor_Depth= ');
 %Dso=input('Stator_outer_Diameter=');
 %Dsi=input('Inner_Stator_Diameter=');
 %Dro=input('Out_Rotor_Diameter=');
 %Dri=input('Inner_Rotor_Diameter=');
 %Lm=input('Magnet_Thickness=');
 %Lt=input('Length_of_Teeth=');
 %Tsp=input('Tooth_Width/Tooth_Pitch=');               %propably 1/2
 %Msp=input('Magnet_span=');                           %probably 2/3                 
 %Uo=4*3.14*10^-7;                                     % permiability of free space
 
 N=100;
 Np=1;
 D=0.1;
 Dso=0.8;
 Dsi=0.44;
 Dro=0.4;
 Dri=0.12;
 Lm=0.017;
 Lt=0.13;
 Tsp=0.5;
 Msp=2/3;
 Uo=4*3.14*10^-7; 
 
 Wt=(Dsi+Dsi+2*Lt)*pi*Tsp/2/2/2/Np;
 Wair=(Dsi+Dro+2*Lm)*pi/2/2/2/Np;
 Lair=(Dsi-Dro-2*Lm)/2;
 Wsbi=(Dso-Dsi-2*Lt)/2;
 Lsbi=(Dso+Dsi+2*Lt)*pi/2/2/Np;
 Wrbi=(Dro-Dri)/2;
 Lrbi=(Dro+Dri)*pi/2/2/Np;
 Wm=(Dro+Dro+2*Lm)*pi*Msp/2/2/2/Np;
 
 b=[0.00:0.02:0.4];
    h=1000*b;                     %for linear material
 plot(b,h,'b.')
    for i=1:21                     % bi are airgap flux densities
    h1=1000*b(i)*Wair*D/Wt/D;     % teeth      
    h2=1000*b(i)*Wair*D/Wsbi/D;   % stator back iron 
    h3=1000*b(i)*Wair*D/Wrbi/D;   % rotor back
    h4=1000*b(i)*Wair*D/Wm/D;     % Magnet
    I(i)=1/N*((h1*2*Lt+h2*Lsbi+h3*Lrbi)+2*b(i)*(Lair+Lm)/Uo);
    end
    for i=1:21                         % calculate various H around the circuit for given airgap flux density
    h1=1000*b(i)*Wair*D/Wt/D;     % teeth      
    h2=1000*b(i)*Wair*D/Wsbi/D;   % stator back iron 
    h3=1000*b(i)*Wair*D/Wrbi/D;   % rotor back
    h4=1000*b(i)*Wair*D/Wm/D;     % Magnet
    hc(i)=1/2/Lm*((h1*2*Lt+h2*Lsbi+h3*Lrbi)+2*b(i)*Wair*D/Wm/D*Lm/Uo+2*b(i)*Lair/Uo);  %calculate 20 different hc=[.........]
    end;                    


    n=21
    I1 = I(1); In = I(n);
    q=2*N*Wair*D*b;        % calculation of flux linkage
    for i=1:21                                        
    LL(i)=q(i)/I(i);     %calculate 18 diffent inductance=[.........]
    end;                 %LLi is inductance

   figure
   Ii=[I1:0.1:In];             %create a look up table to find a certain LLi for a given Ii, the length of the step is 0.1
   LLi = interp1(I,LL,Ii);   
plot(I,LL,'o',Ii,LLi);      %plot the curve

   hc1=hc(1);hcn=hc(n);
figure
   hci=[hc1:100:hcn];             %create a look up table to find a certain bi for a given hc, the length of the step is 0.1
   bf = interp1(hc,b,hci);   
plot(hc,b,'o',hci,bf);      %plot the curve


I0=input(' input current= ');     
i=1;
while I0>Ii(i)              % to choose Ii(i) from the look up table, make sure it is very closed to the input current I0. The errer is less than 0.1(the length of the step)
     i=i+1;
end 
LL0=LLi(i);

hc0=input(' input hc= ');     
i=1;
while hc0>hci(i)              % to choose hc(i) from the look up table, make sure it is very closed to the input 
     i=i+1;
end 
bf0=bf(i);



bf=bf0;
fprintf (1,' The inductance LL0 = %3.4f \n',LL0)
w=input(' input Angular velocity（rpm）= ');
E=w*3.1415926*(Dsi+Lt)/60*D*bf*N*cos(30/180*3.1415926);               
Force=D*1.5*I0*bf*N*2;
Torque=Force*(Dro/2+Lm);

fprintf (1,' The voltage E = %3.4f \n',E)
fprintf (1,' The Force = %3.4f \n',Force)
fprintf (1,' The Torque = %3.4f \n',Torque)