交流异步电机的Modelica建模

编者语：

细数我国工业的短板众多，而工业软件也是其中之一。工业APP是工业软件在工业互联网环境下的新形态，其发展无法回避工业软件本身的核心关键技术问题。2013年4月发布的德国工业4.0认为软件是工业的未来，未来的工业软件应基于模型的理论、方法和工具自动产生，建立全系统模型必须建立跨领域的模型标准，Modelica是工业领域重要的模型表达与互联标准。为此，小编转发由陈冰博士撰写的《交流异步电机的Modelica建模》一文，仅供读者参考。

Modelica是开放式的建模语言，可以用来方便对各种物理对象建模，具体的仿真软件可以采用OpenModelica等开源软件，也可以采用国产MWorks等支持Modelica的仿真软件。

本文用OpenModelica标准库里的基本元件，搭建了一个三相异步电机的模型，用于机电传动控制系统的学习。原理公式参照陈伯时的《电力拖动自动控制系统--运动控制系统》第3版的190页到195页的内容，实际的电机模型参数参照了Novotny和Lipo的《Vector Control and Dynamics of AC Drives》第78页的一个例子参数并稍作修改。这个模型没有使用dq坐标系。

本模型中使用的电机主要参数为：

• 额定电压（相电压）：220 V

• 额定频率：50 Hz

• 极对数：2

• 转动惯量：0.1 kg.m^2

• 定子电阻：0.531 Ohm

• 转子电阻：0.408 Ohm

• 定子漏感：2.52 mH

• 转子漏感：2.52 mH

• 互感：8.47 mH

上述参数可根据实际电机的参数进行修改，负载转矩和负载惯量可根据实际仿真情况加以修改。

Modelica模型如下：model SACIM 'A Simple AC Induction Motor Model' type Voltage=Real(unit='V'); type Current=Real(unit='A'); type Resistance=Real(unit='Ohm'); type Inductance=Real(unit='H'); type Speed=Real(unit='r/min'); type Torque=Real(unit='N.m'); type Inertia=Real(unit='kg.m^2'); type Frequency=Real(unit='Hz'); type Flux=Real(unit='Wb'); type Angle=Real(unit='rad'); type AngularVelocity=Real(unit='rad/s'); constant Real Pi = 3.1415926;
Current i_A'A Phase Current of Stator'; Current i_B'B Phase Current of Stator'; Current i_C'C Phase Current of Stator'; Voltage u_A'A Phase Voltage of Stator'; Voltage u_B'B Phase Voltage of Stator'; Voltage u_C'C Phase Voltage of Stator'; Current i_a'A Phase Current of Rotor'; Current i_b'B Phase Current of Rotor'; Current i_c'C Phase Current of Rotor'; Frequency f_s'Frequency of Stator'; Torque Tm'Torque of the Motor'; Speed n'Speed of the Motor';
Flux Psi_A'A Phase Flux-Linkage of Stator'; Flux Psi_B'B Phase Flux-Linkage of Stator'; Flux Psi_C'C Phase Flux-Linkage of Stator'; Flux Psi_a'a Phase Flux-Linkage of Rotor'; Flux Psi_b'b Phase Flux-Linkage of Rotor'; Flux Psi_c'c Phase Flux-Linkage of Rotor';
Angle phi'Electrical Angle of Rotor'; Angle phi_m'Mechnical Angle of Rotor'; AngularVelocity w'Angular Velocity of Rotor';
Torque Tl'Load Torque'; parameter Resistance Rs = 0.531'Stator Resistance'; parameter Resistance Rr = 0.408'Rotor Resistance'; parameter Inductance Ls = 0.00252'Stator Leakage Inductance'; parameter Inductance Lr = 0.00252'Rotor Leakage Inductance'; parameter Inductance Lm = 0.00847'Mutual Inductance'; parameter Frequency f_N = 50'Rated Frequency of Stator'; parameter Voltage u_N = 220'Rated Phase Voltage of Stator'; parameter Real p =2'number of pole pairs'; parameter Inertia Jm = 0.1'Motor Inertia'; parameter Inertia Jl = 0.1'Load Inertia';
initial equation
Psi_A = 0; Psi_B = 0; Psi_C = 0; Psi_a = 0; Psi_b = 0; Psi_c = 0; phi = 0; w = 0;
equation u_A = Rs * i_A + 1000 * der(Psi_A); u_B = Rs * i_B + 1000 * der(Psi_B); u_C = Rs * i_C + 1000 * der(Psi_C);
0 = Rr * i_a + 1000 * der(Psi_a); 0 = Rr * i_b + 1000 * der(Psi_b); 0 = Rr * i_c + 1000 * der(Psi_c);
Psi_A = (Lm+Ls)*i_A + (-0.5*Lm)*i_B + (-0.5*Lm)*i_C + (Lm*cos(phi))*i_a + (Lm*cos(phi+2*Pi/3))*i_b + (Lm*cos(phi-2*Pi/3))*i_c; Psi_B = (-0.5*Lm)*i_A + (Lm+Ls)*i_B + (-0.5*Lm)*i_C + (Lm*cos(phi-2*Pi/3))*i_a + (Lm*cos(phi))*i_b + (Lm*cos(phi+2*Pi/3))*i_c; Psi_C = (-0.5*Lm)*i_A + (-0.5*Lm)*i_B + (Lm+Ls)*i_C + (Lm*cos(phi+2*Pi/3))*i_a + (Lm*cos(phi-2*Pi/3))*i_b + (Lm*cos(phi))*i_c;
Psi_a = (Lm*cos(phi))*i_A + (Lm*cos(phi-2*Pi/3))*i_B + (Lm*cos(phi+2*Pi/3))*i_C + (Lm+Lr)*i_a + (-0.5*Lm)*i_b + (-0.5*Lm)*i_c; Psi_b = (Lm*cos(phi+2*Pi/3))*i_A + (Lm*cos(phi))*i_B + (Lm*cos(phi-2*Pi/3))*i_C + (-0.5*Lm)*i_a + (Lm+Lr)*i_b + (-0.5*Lm)*i_c; Psi_c = (Lm*cos(phi-2*Pi/3))*i_A + (Lm*cos(phi+2*Pi/3))*i_B + (Lm*cos(phi))*i_C + (-0.5*Lm)*i_a + (-0.5*Lm)*i_b + (Lm+Lr)*i_c; Tm =-p*Lm*((i_A*i_a+i_B*i_b+i_C*i_c)*sin(phi)+(i_A*i_b+i_B*i_c+i_C*i_a)*sin(phi+2*Pi/3)+(i_A*i_c+i_B*i_a+i_C*i_b)*sin(phi-2*Pi/3));
w = 1000 * der(phi_m); phi_m = phi/p; n= w*60/(2*Pi);
Tm-Tl = (Jm+Jl) * 1000 * der(w);

if time <= 100 then u_A = 0; u_B = 0; u_C = 0; f_s = 0; Tl = 0; else f_s = f_N; u_A = u_N * 1.414 * sin(2*Pi*f_s*time/1000); u_B = u_N * 1.414 * sin(2*Pi*f_s*time/1000-2*Pi/3); u_C = u_N * 1.414 * sin(2*Pi*f_s*time/1000-4*Pi/3); Tl = 10; end if;
end SACIM;

在模型中，我们定义了电机的定子电压方程（由于仿真软件的时间单位是毫秒，所以所有求导操作前都乘以1000，下同）：

  u_A = Rs * i_A + 1000 * der(Psi_A);  u_B = Rs * i_B + 1000 * der(Psi_B);  u_C = Rs * i_C + 1000 * der(Psi_C);

转子电压方程（针对鼠笼转子，如果采用绕线式转子，可以把外接的电阻和电抗参数引入方程）：

0 = Rr * i_a + 1000 * der(Psi_a); 0 = Rr * i_b + 1000 * der(Psi_b); 0 = Rr * i_c + 1000 * der(Psi_c);

定子磁链方程：

  Psi_A =   (Lm+Ls)*i_A + (-0.5*Lm)*i_B + (-0.5*Lm)*i_C +        (Lm*cos(phi))*i_a + (Lm*cos(phi+2*Pi/3))*i_b + (Lm*cos(phi-2*Pi/3))*i_c;  Psi_B = (-0.5*Lm)*i_A +   (Lm+Ls)*i_B + (-0.5*Lm)*i_C + (Lm*cos(phi-2*Pi/3))*i_a +        (Lm*cos(phi))*i_b + (Lm*cos(phi+2*Pi/3))*i_c;  Psi_C = (-0.5*Lm)*i_A + (-0.5*Lm)*i_B +   (Lm+Ls)*i_C + (Lm*cos(phi+2*Pi/3))*i_a + (Lm*cos(phi-2*Pi/3))*i_b +        (Lm*cos(phi))*i_c;

转子磁链方程：

Psi_a = (Lm*cos(phi))*i_A + (Lm*cos(phi-2*Pi/3))*i_B + (Lm*cos(phi+2*Pi/3))*i_C + (Lm+Lr)*i_a + (-0.5*Lm)*i_b + (-0.5*Lm)*i_c;Psi_b = (Lm*cos(phi+2*Pi/3))*i_A + (Lm*cos(phi))*i_B + (Lm*cos(phi-2*Pi/3))*i_C + (-0.5*Lm)*i_a + (Lm+Lr)*i_b + (-0.5*Lm)*i_c;Psi_c = (Lm*cos(phi-2*Pi/3))*i_A + (Lm*cos(phi+2*Pi/3))*i_B + (Lm*cos(phi))*i_C + (-0.5*Lm)*i_a + (-0.5*Lm)*i_b + (Lm+Lr)*i_c;

电磁转矩方程：

Tm =-p*Lm*((i_A*i_a+i_B*i_b+i_C*i_c)*sin(phi)+(i_A*i_b+i_B*i_c+i_C*i_a)*sin(phi+2*Pi/3)+(i_A*i_c+i_B*i_a+i_C*i_b)*sin(phi-2*Pi/3));

转子的角速度和机械角位移存在导数关系：

w = 1000 * der(phi_m);

下面两个公式实现转子电角度与机械角度，角速度和电机转速之间的单位转换。

phi_m = phi/p; n= w*60/(2*Pi);

电机的实际速度由电机的电磁转矩、负载转矩以及电机和负载的共同负载惯量决定(采用最简化的刚体动力学模型，可逐渐扩展):

Tm-Tl = (Jm+Jl) * 1000 * der(w);

在此基础上，我们就可以通过设定外部条件的变化来仿真电机的运行，如改变定子电压和频率，定子串电阻，转子串电阻，改变负载大小等。

这里只给出了最简单的额定电压直接启动的例子：

if time <= 100 then u_A = 0; u_B = 0; u_C = 0; f_s = 0; Tl = 0; else f_s = f_N; u_A = u_N * 1.414 * sin(2*Pi*f_s*time/1000); u_B = u_N * 1.414 * sin(2*Pi*f_s*time/1000-2*Pi/3); u_C = u_N * 1.414 * sin(2*Pi*f_s*time/1000-4*Pi/3); Tl = 10; end if;

在OpenModelica运行如下命令可进行模型的仿真：

simulate(SACIM,startTime=0,stopTime=2000)

作者简介

陈冰，男，工学博士，华中科技大学机械科学与工程学院副教授，中国自动化学会边缘计算专业委员会委员。主要从事机电装备运动控制、网络化控制和工业边缘计算等领域的研发工作，承担“新工科”机电教学改革和实践教学系统设计相关工作。