【MATLAB】差分法による電位分布計算

ここでは球ー球ギャップの電位分布をMATLABで計算して可視化してみた

高圧電極の電位で規格化してある

差分法による電界分布計算

数式は数理工学社の『高電圧工学』を，プログラムは大学の「電磁界応用工学」という講義で配布されたC言語のプログラムを参考にしました

高電圧工学 (新・電気システム工学)

作者:日高邦彦
出版社/メーカー: 数理工学社
発売日: 2009/01/01
メディア: 単行本

差分法とは
プログラム
改良したコード

差分法とは

電位分布を数値計算するための手法の1つ

電位分布を知りたい領域を細かく分け、ある点での電位はその隣の電位との関係から計算できることを利用する

まずは二次元平面の場合を考える

点$(x_0,y_0)$の電位を$\phi_0=\phi(x_0,y_0)$、その周りの電位を $\phi_1=\phi(x_0,y_0+h), \phi_2=\phi(x_0-h,y_0), \phi_3=\phi(x_0,y_0-h), \phi_4=\phi(x_0+h, y_0)$とおく

点$(x_0, y_0)$のまわりのテイラー展開やラプラスの方程式 $\phi_{xx}(x_0, y_0)+\phi_{yy}(x_0, y_0)=0$ を用いると式(1)の関係が導かれる

（詳しい導出過程はページ下の参考文献を参照されたい）

\phi_0=\frac{1}{4}(\phi_1+\phi_2+\phi_3+\phi_4) \tag{1}

この場合z方向に隣り合った電位は考えていないし、実際の電極は回転体で、円筒座標系の方が表しやすい

上の直交座標系でのx, y をr, zに変えるが、円筒座標系のラプラスの方程式が $\frac{\partial^2\phi}{\partial r^2}+\frac{1}{r}\frac{\partial\phi}{\partial r}+\frac{\partial^2\phi}{\partial z^2}=0$ となることから次式のように項が1つ増えることに注意する

\begin{aligned} \phi_0 &= \frac{1}{4}\left(\phi_1+\phi_2+\phi_3+\phi_4-\frac{h^2}{r_0}\phi(r_0)'\right) \\ &= \frac{1}{4}\left\{\phi_1+\phi_2+\phi_3+\phi_4-\frac{h}{2r_0}(\phi_4-\phi_2)\right\}\,\, \left(\because \phi(r_0)'=\frac{\phi_4-\phi_2}{2h}\right) \end{aligned}

以上で領域の内部の点の電位が隣の点から近似計算できることが分かったが、これだけでは値が定まらないので境界部分に条件が必要である

境界条件には次の2種類がある

ディリクレ条件

電位の値を決める

ループ内でディリクレ条件が決まっている場所の電位は計算しなくてよい

例: 高圧電極、低圧電極、接地面、r方向遠方

ノイマン条件

電界が0で電位の傾き$\phi'=0$

$r$方向のノイマン条件なら、$\phi_r=0$より$\phi(r_0+h,z_0)=\phi(r_0-h,z_0)$となるので, 点 $(r_0,z_0)$の電位は

\begin{aligned} \phi_0 = \frac{1}{4}\left\{2\phi(r_0+h,z_0)+\phi(r_0,z_0+h)+\phi(r_0,z_0-h)\right\} \end{aligned}

z方向では$\phi_z=0$より$\phi(r_0,z_0+h)=\phi(r_0,z_0-h)$となるので

\begin{aligned} \phi_0 = \frac{1}{4}\left\{\phi(r_0+h,z_0)+\phi(r_0-h,z_0)+2\phi(r_0,z_0-h)-\frac{h}{2r_0}(\phi(r_0+h,z_0)-\phi(r_0-h,z_0) )\right\} \end{aligned}

例: 回転体の中心軸上など