用Python+SymPy实战推导电机电磁转矩：从磁共能到动态方程

电机控制工程师常被教科书上复杂的电磁转矩公式困扰——那些看似魔术般的推导过程，往往让人知其然而不知其所以然。今天我们将用Python的SymPy库，通过 符号计算 重新演绎经典电磁转矩公式的诞生过程。这不是简单的代码翻译，而是一次 理论重建 的思维实验。

1. 电磁系统建模基础

1.1 建立符号系统

我们首先定义核心物理量的符号变量。在Jupyter Notebook中运行以下代码：

from sympy import symbols, Function, cos, sin, diff

# 定义基本符号
i_A, i_B = symbols('i_A i_B', real=True)  # 绕组电流
θ_r = symbols('theta_r', real=True)       # 转子位置角
L_A, L_B = symbols('L_A L_B', real=True)  # 自感系数
M_AB = symbols('M_AB', real=True)         # 最大互感系数

关键点在于理解这些符号的物理意义：

• i_A , i_B 代表定子和转子绕组电流
• θ_r 是转子机械角度（电角度需要乘以极对数）
• 电感参数 L_A , L_B , M_AB 由电机结构决定

1.2 互感函数建模

实际电机中，互感随转子位置变化。我们建立周期性变化的互感模型：

# 定义位置相关的互感函数
L_AB = M_AB * cos(θ_r)  # 定转子绕组间互感
L_BA = L_AB             # 互易关系

这个简化的余弦模型捕捉了交流电机的基本特征。在真实电机中，可能需要考虑谐波分量：

# 考虑三次谐波的扩展模型（可选）
L_AB_ext = M_AB * cos(θ_r) + 0.1*M_AB * cos(3*θ_r)

2. 磁能与磁共能计算

2.1 磁链方程构建

根据电磁感应定律，建立完整的磁链表达式：

psi_A = L_A * i_A + L_AB * i_B  # 定子磁链
psi_B = L_B * i_B + L_BA * i_A  # 转子磁链

2.2 能量计算实现

磁共能是推导电磁转矩的关键。我们实现磁共能的双重积分计算：

from sympy import integrate

# 计算磁共能（忽略漏感）
W_m_prime = integrate(psi_A, (i_A, 0, i_A)) + integrate(psi_B, (i_B, 0, i_B))
W_m_prime = W_m_prime.simplify()

运行结果应显示：

L_A*i_A**2/2 + L_B*i_B**2/2 + M_AB*i_A*i_B*cos(theta_r)

这个结果与教科书中的经典表达式完全一致，验证了我们的建模正确性。

3. 电磁转矩推导

3.1 符号微分计算

根据虚位移法，电磁转矩是磁共能对位置角的偏导数：

T_e = diff(W_m_prime, θ_r)  # 电磁转矩基本表达式

得到的解析式为：

-M_AB*i_A*i_B*sin(theta_r)

3.2 物理意义解读

这个负号有明确的物理意义：

• 当θ_r=90°时，转矩达到最大值
• 负号表示转矩方向趋向于减小θ_r（磁路最小 reluctance 原理）

用Python可视化转矩特性：

import numpy as np
import matplotlib.pyplot as plt

theta = np.linspace(0, 2*np.pi, 100)
torque = -np.sin(theta)

plt.figure(figsize=(8,4))
plt.plot(theta, torque, label='Electromagnetic Torque')
plt.xlabel('Rotor Position [rad]')
plt.ylabel('Torque [p.u.]')
plt.grid(True)
plt.legend()
plt.show()

4. 动态系统仿真

4.1 运动方程构建

结合机械运动方程，建立完整系统模型：

from sympy import Eq, Derivative

# 定义时间变量和函数
t = symbols('t')
i_A_t = Function('i_A')(t)
i_B_t = Function('i_B')(t)
theta_t = Function('theta')(t)

# 机械运动方程
J, B = symbols('J B', real=True)  # 转动惯量和阻尼系数
T_L = symbols('T_L', real=True)   # 负载转矩

mech_eq = Eq(J*Derivative(theta_t, t, t) + B*Derivative(theta_t, t), 
             -M_AB*i_A_t*i_B_t*sin(theta_t) - T_L)

4.2 状态空间模型

将系统转换为状态空间形式便于仿真：

# 定义状态变量
omega = Derivative(theta_t, t)
states = [theta_t, omega]

# 状态方程
state_eq = [
    omega,
    (-B*omega + (-M_AB*i_A_t*i_B_t*sin(theta_t) - T_L))/J
]

5. 实际应用案例

5.1 永磁同步电机简化模型

对于PMSM，可固定转子电流为等效永磁电流：

I_m = symbols('I_m', real=True)  # 等效永磁电流
W_pm = L_A*i_A**2/2 + M_AB*i_A*I_m*cos(θ_r)  # 简化磁共能
T_pm = diff(W_pm, θ_r)  # 永磁转矩

得到经典的永磁转矩公式：

-M_AB*I_m*i_A*sin(theta_r)

5.2 转矩优化控制

通过电流控制实现最大转矩输出：

# 定义最优电流角
gamma = symbols('gamma', real=True)

# 采用id=0控制策略
i_d, i_q = symbols('i_d i_q', real=True)
i_A_opt = i_q * cos(θ_r + gamma)
i_B_opt = I_m  # 永磁等效

T_opt = (-M_AB*i_A_opt*i_B_opt*sin(θ_r)).simplify()

6. 验证与扩展

6.1 能量守恒验证

检查机电能量转换的完整性：

# 电能输入
dWe = i_A * diff(psi_A, t) + i_B * diff(psi_B, t)

# 机械能输出
dWmech = T_e * diff(θ_r, t)

# 磁能变化
dWm = diff(W_m_prime, t)

# 验证能量守恒
energy_balance = (dWe - dWm - dWmech).simplify()
assert energy_balance == 0  # 应返回True

6.2 非线性效应扩展

考虑磁饱和效应时，电感变为电流的函数：

L_A_sat = L_A / (1 + 0.1*i_A**2)  # 简单饱和模型
psi_A_sat = L_A_sat * i_A + L_AB * i_B
W_prime_sat = integrate(psi_A_sat, i_A) + integrate(psi_B, i_B)

这种建模方式可以更精确地预测实际电机的非线性特性。

7. 工程实践建议

1. 参数获取 ：实际电机参数可通过：

   • 空载实验测定互感曲线
   • 锁转子测试获取自感
   • 有限元分析验证理论模型

2. 控制实现 ：

   # 典型FOC电流控制伪代码
   def current_controller(theta, i_q_ref):
       i_a = i_q_ref * cos(theta)
       i_b = i_q_ref * sin(theta)
       return convert_to_phase_currents(i_a, i_b)
   

3. 实时仿真 ：将符号模型导出为C代码：

   from sympy.utilities.codegen import codegen
   codegen(('torque_eq', T_e), language='C')
   

通过这种符号计算驱动的建模方法，我们不仅复现了教科书公式，更建立了一个可扩展的分析框架。下次当你面对电磁转矩公式时，不妨打开Python环境，亲手重建这些方程——这或许是最深刻的学习方式。