Pythonを使ったDCモータのシミュレーション

Pythonを使ってブラシ付きDCモータのシミュレーションをやってみたいと思います。

(1) まずは、ブラシ付きDCモータの電気特性について整理します。

モータの内部に踏み入れずに、モータの入力電圧と出力回転速度のみについてみると、
　　 $v_a - R_a \cdot i_a - V_{B} = K_E \cdot \omega$ ...........①
が成り立つ。

ただし、

$v_a$ : モータ端子電圧　(V)
$R_a$ : 電機子抵抗 $(\Omega)$
$i_a$ : モータ電流 (A)
$V_B$ : ブラシ接触抵抗による電圧降下 (V)
$K_E$ : 逆起電力定数 (Vs/rad)
$\omega$ : 角速度 (rad/s)

また、モータの入力電流と出力トルク間についてみると、
　　 $T = K_T \cdot i_a - T_F$ ..........................②
が成り立つ。

ただし、
$T$ : 出力トルク　(Nm)
$K_T$ : トルク定数 (Vs/rad)
$T_F$ : 摩擦などによるトルク損失 (Nm)

さらにモータの入出力は

モータ入力 $P_{in} = v_a \cdot i_a$ ..........................③
モータ出力 $P_{out} = T \cdot \omega$ ..........................④
モータ効率 $\eta=\frac{P_{out}}{P_{in}}$ ..........................⑤

式①と②を式④に代入すると
$P_{out} = -\frac {R_a \cdot K_T} {K_E} i^2_a + \frac {K_T v_a - K_T V_B + T_F R_a} {K_E} \cdot i_a + \frac {T_F V_B - T_F v_a} {K_E}$ .....⑥
によって表現されます。式⑥が2次関数なので、出力が最大値をとるための入力値 $i_a$ を求めることができます。

#!/usr/bin/env python

import sympy as sym

Ra, Kt, Ke, ia, va, Tf, Vb, = sym.symbols('Ra Kt Ke ia va Tf Vb')
expr = -Ra * Kt * ia**2/Ke  + (Kt*va - Kt*Vb + Tf*Ra)*ia/Ke +(Tf*Vb - Tf*va)/Ke
d = sym.Derivative(expr, ia) 
f = d.doit()
print(sym.solve(sym.Eq(f, 0), ia))

$i_a = \frac {K_T \cdot (-V_B + v_a) + R_a \cdot T_F)} {2 \cdot K_T \cdot R_a}$

のときに、最大値をとることがわかります。

(2) 続いてモータにつながる機械系の特性をも考慮に入れます。

トルクは慣性モーメントJ、回転角 $\theta$ および角速度 $\omega$ 間の関係は、
　　 $T=J \cdot \frac{\mathrm{d}^2 \theta}{\mathrm{d}t^2} = J \cdot \frac{\mathrm{d} \omega}{\mathrm{d}t}$ ...................⑦
式②から損失部分を取り除き、式③に代入すると、モータ電流と出力トルクの関係が得られる。
　　 $T = K_{T} \cdot i_a = J \cdot \frac{\mathrm{d} \omega}{\mathrm{d}t}$ ..................⑧
モータ逆起電力と回転速度間の関係は、
　　 $e = K_{E} \cdot \omega$ ...............................⑨
によって与えられます。

また、モータの中身に踏み入れてコイルに存在するL成分コイルインダクタンス[ $L_a$ の影響をモデリングすると、モータ中電圧の平衡状態は
　　 $v_a(t) = L_a \cdot \frac{\mathrm{d} i_a(t)}{\mathrm{d}t} + R_a \cdot i_a(t) + V_{B} + e(t)$ ............⑩
式⑥のブラシ接触抵抗による電圧降下を無視し、式④⑤と連立してラプラス変換すれば、
　　 $\{\array{lcl$V(s) = (s \cdot L_a + R_a) \cdot I_a(s) + E(s)\\T(s) = K_T \cdot I_a(s) \\ \Omega(s) = T(s) \cdot \frac {1} {J \cdot s} \\ E(s) = K_E \cdot \Omega(s) }$ .......................⑪
となります。

式⑦の連列方程式をブロック線図で表現すると下図となり、モータ特性を表すものとなります。

式⑦の時間ドメイン方程式として表現すると
$\frac {d} {dt} \left(\large\begin{array}{GC+40}\omega(t)\\i_a(t)\end{array}\right){\Large=} \left( \begin{array}{CC} 0 & \frac {K_T} {J} \\ -\frac {K_E} {L_a} & -\frac {R_a} {L_a} \end{array}\right) \cdot \left(\large\begin{array}\omega(t)\\i_a(t) \end{array}\right) + \left(\begin{array}0 \\ \frac {v_a(t)} {L_a} \end{array}\right)$ ............⑫
となります。

逆起電力定数 $K_E$ とトルク定数 $K_T$ を0.175、慣性モーメントJを0.00005、コイルインダクタンス $L_a$ を3.1mH、電機子抵抗 $R_a$ を2.5ohmとおいて、3.0vのステップ電圧を与えた場合のモータ回転速度を計算してみます。

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

Kt = 0.175     # Vs/rad
Ke = 0.175    # Nm/A
J  = 0.00005  # kg m2
La = 3.1e-3   # H
Ra = 2.5      # ohm

def model(z, t, u):
    w = z[0]
    ia = z[1]
    dwdt = Kt*ia/J
    didt = -Ke*w/La-Ra*ia/La + u/La
    dzdt = [dwdt, didt]
    return dzdt

# initial condition
z0 = [0,0]

# number of time points
n = 50001

# time points
t = np.linspace(0,2.5,n)

# step input
u = np.zeros(n)
# change to 3.0 at time = 1.0s
u[20001:] = 3.0

# store solution
w = np.empty_like(t)
ia = np.empty_like(t)
# record initial conditions
w[0] = z0[0]
ia[0] = z0[1]

# solve ODE
for i in range(1,n):
    # span for next time step
    tspan = [t[i-1],t[i]]
    # solve for next step
    z = odeint(model,z0,tspan,args=(u[i],))
    # store solution for plotting
    w[i] = z[1][0]
    ia[i] = z[1][1]
    # next initial condition
    z0 = z[1]

# plot results
plt.plot(t,u,'g:',label='Input v(t)')
plt.plot(t,w,'b-',label='Output w(t)')
plt.plot(t,ia,'r--',label='ia(t)')
plt.ylabel('values')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()