はじめに
RLC回路の過渡応答解析は一般的には微分方程式を用いて考察する。しかし、その微分方程式を解析的に解くのは難易度が高く、ましては応答をグラフ化するのは手計算では限界がある。そこで、微分方程式をプログラムを用いて解くことでRLC回路の過渡応答の様子を伺うことを本記事では試みる。具体的には、微分方程式を状態方程式に変換した上で、簡易な差分法を用いて解析する。以下は解析結果のグラフである。
問題設定
以下のようなRLC直列回路に正弦波電圧を加えたときのコンデンサに蓄積された電荷の応答を調べる。ただし、R=2[Ω],L=1[mH],C=10[μF]として、電源周波数を50[Hz]、電源電圧を1[V]とした。
微分方程式と状態方程式
RLC直列回路の微分方程式は以下のようになる。ただし、電源電圧とコンデンサに蓄積された電荷の瞬時値をそれぞれ$e(t),q(t)$とした。
e(t)=R\frac{dq(t)}{dt}+L\frac{d^2 q(t)}{dt^2}+\frac{q(t)}{C}
ゆえに、
\frac{d^2 q(t)}{dt^2}=-\frac{q(t)}{LC}-\frac{R \frac{dq(t)}{dt}}{L}+\frac{e(t)}{L}
これを状態方程式にすると以下のようになる。
\begin{pmatrix}
\frac{dq(t)}{dt} \\
\frac{d^2 q(t)}{dt^2} \\
\end{pmatrix}
=
\begin{pmatrix}
0 & 1 \\
-\frac{R}{L} & -\frac{1}{CL} \\
\end{pmatrix}
\begin{pmatrix}
q(t) \\
\frac{d q(t)}{dt} \\
\end{pmatrix}
+
\begin{pmatrix}
0 \\
\frac{1}{L} \\
\end{pmatrix}
e(t)
したがって、以下のように表す。
\dot{\textbf{x}}(t)
=
\begin{pmatrix}
0 & 1 \\
-\frac{R}{L} & -\frac{1}{CL} \\
\end{pmatrix}
\textbf{x(t)}
+
\begin{pmatrix}
0 \\
\frac{1}{L} \\
\end{pmatrix}
e(t)
ただし、
\textbf{x}(t)=\begin{pmatrix}
q(t) \\
\frac{d q(t)}{dt} \\
\end{pmatrix}
\dot{\textbf{x}}(t)=\textbf{v}(t)=
\begin{pmatrix}
\frac{dq(t)}{dt} \\
\frac{d^2 q(t)}{dt^2} \\
\end{pmatrix}
とする。
ここで、位置ベクトル、速度ベクトルをそれぞれ$\textbf{x}(t),\dot{\textbf{x}}(t)$として、差分法で微分方程式を解くことを考える。
アルゴリズム
微分方程式(状態方程式)を解くアルゴリズムは以下の通りである。
ここでは、ある物体がx軸に対して直線運動していると考える。
変位と速度(差分法の原理)
微小区間($t=t$から$t=t+\Delta t$)に速度$v(t)$（変位の微分）と変位$x(t)$が以下のよう与えられたとする。
$x(t+\Delta t)$は微分の定義式から以下のように考えることができる。
これは、$t=t$のときの情報のみを頼りに$t=t+\Delta t$のときの情報を推定することができることを意味している。
\textbf{x}(t+\Delta t)=\textbf{x}(t)+\textbf{v}(t)\Delta t
これを$t=0$から$t=t$まで加算することで$t=t$での位置$\textbf{x}(t)$を求めることができる。
プログラム
さて、以上の考察をもとにして以下のようなプログラムを作成した。
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
# 微分方程式:V =R Q'*L Q''+Q/C
#パラメータ
R = 2 # 抵抗値
L = 1e-3 # インダクタンス
C = 10e-6 # 静電容量
#初期条件 （電荷と電荷の微分量について）
x=np.array([0,0])
#初期時間
t=0
#時間と状態量（ここでは電荷と電荷の微分量）
time_ary=[]
X=[]
#微小時間の時間間隔(小さくするほど精度は上がるが計算時間が比例して増大する)
delta_t=0.00001
while(t<0.02):
E=np.sin(2*math.pi*50*t)
v=np.array([[0, 1], [-1/(L*C),-R/L]])@x+np.array([0,1/L])*E
#位置の微小変化
x=x+v*delta_t
#print(v)
#演算結果を記録
time_ary.append(t)
X.append(x[1])
#時間の更新
t=t+delta_t
plt.xlabel("応答時間[s]")
plt.ylabel("コンデンサに蓄えられた電荷量[C]")
plt.plot(time_ary,X)
#グラフを保存する
plt.savefig("RLC_状態方程式_電荷.png")
#グラフの表示
plt.show()
これを実行すると以下のようなグラフが出力される。
まず、初期電荷量は0なので、電荷は急速に電源から供給され蓄積されるが、共振周波数付近で振動した。次に、十分に時間が経過すると電源周波数で緩やかな振動をした。
まとめ
今回は、RLC直列回路の過渡応答について、状態方程式を用いることでプログラミングによりシミュレーションした。結果、振動解を実用計算時間内において十分に再現することが可能であるということが分かった。また電流値の応答自体は、初期は共振周波数付近の振動をするが、時間が経つに従って電源周波数での緩やかな振動に落ち着くということも確認することができた。