【オイラー法】非線形素子を含んだLR直列回路の過渡解析

【問】以下の電気回路に流れる電流\(I\)の時間変化を、オイラー法を用いて求めよ。ただし、\(E=10[V]\)、\(R_{1}=1[Ω]\)、\(L=1[H]\)、\(R_{2}=(\frac{I}{5})^{10}[Ω])\)とする。計算には、Excelを用いて良く、時間刻み幅\(h=0.1[s]\)とせよ。

はじめに
オイラー法とは
オイラー法を用いた回路方程式の変形
発展：電気と熱の連成解析
最後に

はじめに

電気回路の過渡解析を、オイラー法を用いて行います。

非線形素子\(R_{2}\)の無い、一般的なLR直列回路は

回路方程式

\begin{eqnarray}E=L\dfrac{dI}{dt}+R_{1}I\end{eqnarray}が成立し、ラプラス変換を使用すると

\begin{eqnarray}I\left( t\right) =\dfrac{E}{R_{1}}\left( 1-e^{-\frac{R_{1}}{L}t}\right)\end{eqnarray}

を手計算で難なく導けると思います。

しかし、今回のような、非線形素子\(R_{2}=(\frac{I}{5})^{10}\)が入った場合はどうなるのでしょうか。

イメージとしては、ダイオードです。ある閾値電流までは殆ど電圧は立ちませんが、閾値を超えると非常に大きな電圧が立ち、これ以上回路に電流が流れないような振る舞いになります。

このような、変数を何乗もした式を手計算でラプラス変換することは非常に厳しく、厳密解を求めることは困難だと考えられます。

このようなとき、計算機を使用して微分方程式をt領域のまま近似解を求める方法があり、オイラー法が最もオーソドックスな解法になります。

本記事では、オイラー法の使い方を、電気回路の過渡現象を題材に考えたいと思います。

そして、発展版として、電気回路で発生したジュール熱を次回の時間ステップで抵抗値にフィードバックする「連成解析」の手法についても説明します。

オイラー法とは

下記の１～３一連の流れを言います。

求める変数に対し、t=0の初期値\(I(0)\)を与えて、その時の微分係数\(\frac{dI}{dt}\)を求める
求めた微分係数を時間の刻み幅分掛け算し、次の時間ステップ\(I(h)\)における電流値の近似解を出す
次の時間ステップにおける微分係数を同様に算出し、その次の時間ステップで同様の操作を繰り返す

イメージ論ですが、ある微分方程式の解(真値)は下記のグラフ(桃色)になるとします。

オイラー法で求める近似解は青色のようになります。t=0で求めた傾きを時間幅分掛け算し、その次の時間でも同様のことを行うと、真値に近い離散的な値を得られることを示しています。(線形近似とも言われることがあります。)

オイラー法を用いた回路方程式の変形

概要を説明しても、実際に計算してみないことにはイメージを掴めないと思います。

問題文で与えられた回路方程式を、オイラー法を用いて変形してみましょう。

まず、キルヒホッフの電圧則により

\begin{eqnarray}E=L\dfrac{dI(t)}{dt}+R_{1}I+\left( \dfrac{I(t)}{5}\right) ^{10}\cdot I(t)\qquad・・・(1)\end{eqnarray}

を立式することができます。問題文で与えられたパラメータを代入し

\begin{eqnarray}10=\dfrac{dI(t)}{dt}+1+\left( \dfrac{I(t)}{5}\right) ^{10}+I(t)\qquad・・・(2)\end{eqnarray}

\begin{eqnarray}\dfrac{dI(t)}{dt}=10-I(t)-\left( \dfrac{I(t)}{5}\right) ^{10}\cdot I(t)\qquad・・・(3)\end{eqnarray}

となります。ここで、オイラー法の出番です。

まず、適当な初期値を代入します。ここでは、t=0[s]のとき、回路に電流が流れていないことを利用します。

(i)t=0のとき、\(I(t)=0\)だから、(3)式に値を代入することにより

\begin{eqnarray}\dfrac{dI(t)}{dt}=10\end{eqnarray}

\begin{eqnarray}I(0.1)&=&I(0)+\dfrac{dI(t)}{dt}*h \\ &=&0+10*0.1 \\&=&1\end{eqnarray}

(ii)t=0.1のとき、(i)より、\(I(t)=1\)だから、(3)式に値を代入することにより

\begin{eqnarray}\dfrac{dI(t)}{dt}=10-1=9\end{eqnarray}

\begin{eqnarray}I(0.2)=I(0.1)+9*0.1=1.9\end{eqnarray}

(iii)t=0.2のとき、\(I(t)=1.9\)だから・・・

以上を繰り返していくと、近似解を求めることができます。

結局、Excelを使用して解いてみると、以下の表(画像)になります。

0.7秒あたりで、回路に流れる電流値が5Aに到達し、非線形素子(第3項)から急激に電圧が立つことで電流がこれ以上流れなくなりました。

電流波形をグラフにすると、以下のようになります。

非線形素子が無い場合は10Aへ漸近しますが、有る場合は、先ほどの表のように5A以上流れないことが視覚的にも確認できました。

以上の反復計算を行うことで、非線形項がある微分方程式の場合でも近似解を求めることができます。

発展：電気と熱の連成解析

研究レベルのより発展的な内容になってくると、電流の過渡解析だけ行えば良いものでは無くなってきます。現実的には、抵抗素子に電流が流れるとジュール熱が発生するはずです。

これにより、抵抗素子の温度が上昇し、抵抗値も上昇することが予想されます。

より詳細な解析を行うならば、各時間ステップにおいて抵抗温度値を求め、回路パラメータにフィードバックすることが求められます。

一見、別の領域である2つの世界(本記事では電気と熱)を統合してシミュレーションを行う手法を連成解析と言います。モータなど、電気機器における解析で広く用いられています。

電気と熱の連成解析の手法

まず、ある時間ステップにおいて、電流が流れるとジュール熱\(Q[W]\)が発生します。これは瞬時値であるため、オイラー法で取っている時間幅\(h\)をかけると、時間ステップ間で\(Qh[J]\)の熱エネルギーが発生したことになります。

これを熱伝導方程式に代入することで、抵抗素子の温度上昇を計算することができます。

\begin{eqnarray}\dfrac{∂T}{∂t}=\dfrac{λ_{x}}{C} \dfrac{∂^{2} T}{∂x^{2 }}+\dfrac{Qh}{C}\end{eqnarray}

※簡単化のため、一次元で考えています。(\lambda_{x}\)は熱伝導率、\(C\)は熱容量を表します。

メッシュ間の距離\(⊿d\)とおく。このとき、中心差分近似により、抵抗\(R_{2}\)での温度上昇\(T_{0}(t+1)-T_{0}(t)\)は、以下の式に近似できます。

\begin{eqnarray}\dfrac{T_{0}\left( t+1\right) -T_{0}\left( t\right) }{h}=\dfrac{\lambda_{ x}}{C}\dfrac{T_{1}\left( t\right) -2T_{0}\left( t\right) +T_{1}(t)}{\left( \Delta d\right) ^{2}}+\dfrac{Qh}{C} \\ T_{0}\left( t+1\right) =T_{0}\left( t\right) +\dfrac{\lambda _{x}h}{C} \dfrac{T_{1}\left( t\right) -2T_{0}\left( t\right) +T_{1}(t)}{\left( \Delta d\right) ^{2}}+\dfrac{Qh^{2}}{C}\end{eqnarray}

で表すことができます。前の時間ステップの値は既知であるため、左辺の未知数を求めることができますね。

ある位置で時間ステップ\(h\)の間に発生する温度上昇量のイメージは以下になります。

「隣の地点から伝搬する熱量+ある位置で発生するジュール熱に熱容量を割った値」

それぞれ、式(13)の右辺第1項、第2項に対応しています。

次回時間ステップ時、まず電流値の値を求めますが、この時に前回時間ステップで算出した温度を温度に依存する抵抗値の式に代入することで連成解析を実現できます。