内置式电机解析计算

yunyunfeng

我利用子域法和磁路法对内置式电机进行了计算，但是做出的磁密波形不对称，求解过程中用到的公式和矩阵是从一篇文献中拿来的，计算了其中的2维部分，没有考虑齿槽效应，最后用matlab进行了计算。以前做过考虑齿槽的模型，出的图也不对称。
以下为编的程序：        最后附带了参考的文献
%-----------------主要参数
Rr=0.0791;
Rs=0.0798;
p=4;
n=5100;
bsa=0.05;
boa=0.011;
Lz=0.13;
hM=0.006;
LM=0.017;
miu0=4*pi*10^(-7);
miur=1.059;
BM=1.282;
Bs=2;
alpha_p=0.558;
alpha_b=0.035;
hb=0.0021;
hb_1=0.0028;

%-----------------------磁路计算
FM=BM*hM/miu0/miur;
RM=hM/miu0/miur/LM/Lz;
phi_a=Bs*hb*Lz;
phi_b=Bs*hb_1*Lz/2;

%-----------------------主程序
K=100;
for k=1:K
   
    Ik_1(k)=1;
    Ik=diag(Ik_1);
    K11=Ik;
    K23=Ik;
    K32=-1*Ik;
    K44=-1*Ik;
   
    Gk_1(k)=(Rr/Rs)^k;
    Gk=diag(Gk_1);
    K12=-1*Gk;
    K24=-1*Gk;
    K31=Gk;
    K43=Gk;
   
    tau(k)=8*p*sin(k*(alpha_p+alpha_b)/2)*sin(k*alpha_b)/pi/alpha_b/k/k;
    w=2*pi*340*pi/180;            % 角速度;  频率：f=340
    alpha_0=7.5*pi/180;           %初始角度：360/Ns/2
    t=1;                          %选择时间为1 的时刻
    alpha_t=w*t+alpha_0;
   
    K35_1(k)=-1*miu0*(tau(k)*sin(k*alpha_t));
    K35=K35_1';
    K45_1(k)=-1*miu0*(tau(k)*cos(k*alpha_t));
    K45=K45_1';
   
    Gak(k)=-2*Lz*Gk(k)*sin(k*alpha_t);
    Gbk(k)=-2*Lz*sin(k*alpha_t);
    Gck(k)=2*Lz*Gk(k)*cos(k*alpha_t);
    Gdk(k)=2*Lz*sin(k*alpha_t);
    alpha_w=(alpha_p+2*alpha_b)/2;
    K51(k)=Gak(k)*sin(k*alpha_w);
    K52(k)=Gbk(k)*sin(k*alpha_w);
    K53(k)=Gck(k)*sin(k*alpha_w);
    K54(k)=Gdk(k)*sin(k*alpha_w);
    K55=1/RM;
    Y5=FM/RM-phi_a-phi_b;
   
    K13=zeros(K,K);
    K14=zeros(K,K);
    K21=zeros(K,K);
    K22=zeros(K,K);
    K33=zeros(K,K);
    K34=zeros(K,K);
    K41=zeros(K,K);
    K42=zeros(K,K);
    K15=zeros(K,1);
    K25=zeros(K,1);
   
    Y1=zeros(K,1);
    Y2=zeros(K,1);
    Y3=zeros(K,1);
    Y4=zeros(K,1);
end
x=[K11,K12,K13,K14,K15
   K21,K22,K23,K24,K25
   K31,K32,K33,K34,K35
   K41,K42,K43,K44,K45
   K51,K52,K53,K54,K55];
y=[Y1
   Y2
   Y3
   Y4
   Y5];
b=x\y;
Ak=b(1:K,1);
Bk=b((K+1):(2*K),1);
Ck=b((2*K+1):(3*K),1);
Dk=b((3*K+1):(4*K),1);
ao=b((4*K+1),1);

r=(Rs+Rr)/2;
alpha=0:pi/2/251:pi/2;
Br=0;
for k=1:K
    Br_1=-1*(k/r)*(Ak(k)*(r/Rs)^k+Bk(k)*(r/Rr)^(-1*k))*sin(k*alpha)+(k/r)*(Ck(k)*(r/Rs)^k+Dk(k)*(r/Rr)^(-1*k))*cos(k*alpha);
    Br=Br+Br_1;
end
B=Br';
%--------------------作图

figure
plot(alpha,B,'b-');
xlabel('Rotor position angle (Mech.Deg)');
ylabel('Air gap flux density (T)');
grid;

sakymuni16

楼主问题解决了吗；虽然我没找到问题所在。。。

哥的一把刀

楼主解决了么，有联系方式可以一起讨论下么