msp1.m

基于MATLAB电机磁场程序

% 1:air-gap, 2:culasse stator, 3: fer rotor,  
% 4-7: stator conductor,  28 - 29: rotor conductor

tic; clear;clc; close all;



%     ------------- Draw the geometry of stator, rotor et air-gap --------------

entref=1e-3;    % air-gap depth
Rs=39.385e-3;	% radius of (rotor+airgap) 
Rr=Rs-entref;	% radius of rotor
Rc=71.75e-3; 	% radius of stator

Nb=140;         % number of turns per phase
Long=0.125;     % longitude of machine in rotation axis direction
p=1;            % ??
f=50;           % ??

beta=35*pi/180; beta=beta/p;    % culasse rotor
beta1=50*pi/180; beta1=beta1/p; % epanouissement polaire

dessinP1;                        % a sub-programme to draw the geometry of the machine section with the parameters defined above
%     --------------


%     ------------- Rotor positioning and initial mesh --------------

aa=-30*pi/180;         % initial angle of rotor
gr=groto(gr,(aa));     % rotating the rotor to the angle specified above using the sub-programme 'groto' 
g=[gsf gr];            % stator + rotor given the whole machine geometry
limites;               % a sub-programme defining the bondary condition
[p,e,t]=initmesh(g);   % initial mesh for the partial differential equation resolution

 %   ------------- mesh refine if necessary ------------
% [p,e,t]=refinemesh(g,p,e,t,[1 30]);
% pdegplot(g); axis equal, hold on ; pdemesh(p,e,t), hold off
% pause
%    -------------------------------------



%     -------------  specification of parameters required for calculation -------
    kexc=0;  
    sigmaexc=kexc*5e6;	               % conductivity of rotor coil conductor
    nuo=1/(4e-7*pi);                   % inverse of permeability of vacuum 
    murs=500; murr=500;                % relative permeability of stator and rotor
                                    
    ks=1; kr=0;                        % ks=1 or 0: stator excited or not; kr=1 or 0: rotor excited or not.
    
    Jex1=0;    Jex11=0;
    Jex2 = 10e6*ks;     Jex21 = -Jex2;
    Jex3 = 0;     Jex31 = 0;
%     Jex3=10e6*ks; Jex31=-Jex3;         % current density of phase a (stator)
%     Jex1=-Jex3/2; Jex11=-Jex1;         % current density of phase b (stator)
% 	Jex2=-Jex3/2; Jex21=-Jex2;         % current density of phase c (stator)
    Jexr1=10e6*kr; Jexr2=-Jexr1;       % current density of rotor excitation

    garnissageP1;	                   % a sub-programme for the giving the current density and permeability at each mesh element
%     ------------- 



%     -------------   resolution of the partial differential equation and figure ploting ---------

u=assempde(bond,p,e,t,nu,0,J);         % resolution of the partial differential equation
% U=ASSEMPDE(B,P,E,T,C,A,F) assembles and solves the PDE problem -div(c*grad(u))+a*u=f

figure (1), pdemesh(p,e,t), axis equal, title ('Mesh of the machine section')
figure (2), pdegplot(g); axis equal; hold on; pdecont(p,t,full(real(u)),30), hold off, title('Magnetic field distribution'), % plots using 30 levels.
%    ------------------------------



% -------------- calculation of the energy stocked in the machine (energy sum in the magnetic field) --------------------

% [ux,uy]=pdegrad(p,t,u);
% bx=uy;
% by=-ux;
% for i=1:length(bx);
%     b(i)=(bx(i)^2+by(i)^2)^.5;
% end
% nu_e=1/(4e-7*pi);
% nu_c=nu_e;
% nu_s=1/500*nu_e;
% nu_r=1/500*nu_e;
% h=sparse(1,ntrg);
% ind_e=find(t(4,:)==2);
% ind_s=find(t(4,:)==1);
% ind_r=find(t(4,:)==3);
% ind_c=find((t(4,:)>=4) & (t(4,:)<=27));
% h(ind_e)=nu_e*b(ind_e)';
% h(ind_s)=nu_s*b(ind_s);
% h(ind_r)=nu_r*b(ind_r);
% h(ind_c)=nu_c*b(ind_c)';
% aire=pdetrg(p,t);
% vol_trg=Long*aire;
% energ_elem_e=1/2*(h(ind_e).*b(ind_e)').*vol_trg(ind_e);
% energ_tot_e=sum(energ_elem_e);
% energ_elem_s=1/2*(h(ind_s).*b(ind_s)').*vol_trg(ind_s);
% energ_tot_s=sum(energ_elem_s);
% energ_elem_r=1/2*(h(ind_r).*b(ind_r)').*vol_trg(ind_r);
% energ_tot_r=sum(energ_elem_r);
% energ_elem_c=1/2*(h(ind_c).*b(ind_c)').*vol_trg(ind_c);
% energ_tot_c=sum(energ_elem_c);
% energ_total=energ_tot_e+energ_tot_s+energ_tot_r+energ_tot_c;
% enr=energ_total;

toc;