🔋

数値計算による等電位線のシミュレーション

2022/12/05に公開

Rust

数値計算

tech

概要

一様な厚さの導体に電極を取り付けて定常電流を流し, 電位の等しい点を結ぶことで等電位線を求める実験がある. 本記事では, この等電位線の実験を数値計算によりシミュレーションする. なお, 実装は Rust, 描画は gnuplot を用いて行う.

方程式の導出

電位を求める方程式を導出する.

$L_x\times L_y$ の矩形導体の中に半径 $r$ の円筒形電極2つを $2d$ だけ離して置く. そして, 電極の中心を結ぶ方向に $x$ 軸, 両電極の中点を通りこれと垂直な方向に $y$ 軸をとる. また, 両電極の電位差は $2V$ である.

スカラーポテンシャル (= 電位) を $\phi$ とすれば次のラプラス方程式が成り立つ.

\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}=0 \tag{1}

ここで, $\phi$ は $x,y$ の関数で, $\phi(x,y)$ と書ける. 数値計算のために上式を離散化する. すなわち, 格子点 $(x_i,y_j)$ における $\phi$ の値を求める. $n_x+1\times n_y+1$ の格子点を考え, 格子点 $(x_{i+1},y_j)$ における $\phi(x_i+\Delta x,y_j)$ を $x$ 方向にテイラー展開すれば

\phi(x_i+\Delta x,y_j) = \phi(x_i,y_j) + \frac{\partial\phi}{\partial x}\Delta x + \frac{1}{2!}\frac{\partial^2\phi}{\partial x^2}\Delta x^2 + \frac{1}{3!}\frac{\partial^3\phi}{\partial x^3}\Delta x^3 + \cdots \tag{2}

格子点 $(x_{i-1},y_j)$ においても同様に

\phi(x_i-\Delta x,y_j) = \phi(x_i,y_j) - \frac{\partial\phi}{\partial x}\Delta x + \frac{1}{2!}\frac{\partial^2\phi}{\partial x^2}\Delta x^2 - \frac{1}{3!}\frac{\partial^3\phi}{\partial x^3}\Delta x^3 + \cdots \tag{3}

式 $(2),(3)$ の和をとれば, $\Delta x$ の奇数次の項が消えて

\phi(x_i+\Delta x,y_j)+\phi(x_i-\Delta x,y_j) = 2\phi(x_i,y_j) - \frac{\partial^2\phi}{\partial x^2}\Delta x^2 + O(\Delta x^2) \tag{4}

式 $(4)$ より

\frac{\partial^2\phi}{\partial x^2} = \frac{1}{\Delta x^2}\{\phi(x_i+\Delta x,y_j)+\phi(x_i-\Delta x,y_j)-2\phi(x_i,y_j)\} \tag{5}

$y$ 方向にも同様に

\frac{\partial^2\phi}{\partial y^2} = \frac{1}{\Delta y^2}\{\phi(x_i,y_j+\Delta y)+\phi(x_i,y_j-\Delta y)-2\phi(x_i,y_j)\} \tag{6}

式 $(5),(6)$ の和をとれば, 式 $(1)$ の左辺を $\phi$ の値で表すことができる. $\phi_{i\ j}\equiv\phi(x_i,y_j)$ , $\Delta x=\Delta y$ とすれば

\frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}=\phi_{i+1\ j}+\phi_{i-1\ j}+\phi_{i\ j+1}+\phi_{i\ j-1}-4\phi_{i\ j}=0 \tag{7}

$\phi_{i\ j}$ について解くと

\phi_{i\ j}=\frac{1}{4}(\phi_{i+1\ j}+\phi_{i-1\ j}+\phi_{i\ j+1}+\phi_{i\ j-1})\ (i=0,1,...,n_x,\ j=0,1,...,n_y) \tag{8}

よって, 格子点 $(x_i,y_j)$ での $\phi$ の値は, その周囲4格子点の平均値であると示された.

また, 式 $(8)$ は差分方程式であり, 右辺には境界値 $\phi_{0\ j},\phi_{n_x\ j},\phi_{i\ 0},\phi_{i\ n_y}$ が含まれる. よって, 解くためにはこれらの値を定める境界条件が必要である. 今回は次の2条件を考える.

電極内の電位.
金属箔の端面.

電極については, 両電極の電位差が $2V$ であるから, 左の電極を正とすれば

\left\{ \begin{array}{l} \phi(x,y)|_{(x-d)^2+y^2\leq r^2}=-V \\ \phi(x,y)|_{(x+d)^2+y^2\leq r^2}=V \end{array} \right.

端面では, 境界上と領域外の電位差が0 ( $\frac{\partial \phi}{\partial x}=0$ または $\frac{\partial \phi}{\partial y}=0$ ) になるように領域外の $\phi$ の値を決めてやればよい. つまり, 境界上の $\phi$ と領域外の $\phi$ を等しくすればよい.

$i=n_x$ のとき, 式 $(8)$ は

\begin{align*} \phi_{n_x\ j}&=\frac{1}{4}(\phi_{n_x\ j}+\phi_{n_x-1\ j}+\phi_{n_x\ j+1}+\phi_{n_x\ j-1})\\ &=\frac{1}{3}(\phi_{n_x-1\ j}+\phi_{n_x\ j+1}+\phi_{n_x\ j-1})\\ \end{align*}

また, $i=n_x,\ j=n_y$ のときは

\begin{align*} \phi_{n_x\ n_y}&=\frac{1}{4}(\phi_{n_x\ n_y}+\phi_{n_x-1\ n_y}+\phi_{n_x\ n_y}+\phi_{n_x\ n_y-1})\\ &=\frac{1}{2}(\phi_{n_x-1\ n_y}+\phi_{n_x\ n_y-1})\\ \end{align*}

となる. 他の境界上の点でも同様となるため, 端面では領域外の点を除いた2または3点の平均をとればよいと分かる.

数値計算手法

式 $(7)$ を反復法によって解く. 反復法は適当な初期値から始めて, 漸化式を反復計算していくことで解を得る数値計算の手法である. $k$ 回目の反復における $\phi_{i\ j}$ の値を $\phi_{i\ j}^{(k)}$ としたとき, 式 $(8)$ は

\phi_{i\ j}^{(k+1)}=\frac{1}{4}(\phi_{i+1\ j}^{(k)}+\phi_{i-1\ j}^{(k)}+\phi_{i\ j+1}^{(k)}+\phi_{i\ j-1}^{(k)}) \tag{9}

$\phi_{i\ j}^{(0)}$ に適当な値を与え, $\phi_{i\ j}^{(1)}, \phi_{i\ j}^{(2)}, ..., \phi_{i\ j}^{(k)}$ と計算していく. そうすれば徐々に解は収束していく. そして誤差が許容誤差 $\epsilon$ 以下になったときの $\phi_{i\ j}^{(k)}$ を解とする. なお, 誤差は次式で与える.

\mathrm{max}\ |\phi_{i\ j}^{(k+1)}-\phi_{i\ j}^{(k)}|

式 $(9)$ のように, $k+1$ 回目の値を $k$ 回目の値で求める反復法をヤコビ法という. 式 $(9)$ より, $\phi_{i\ j}$ は周囲4格子点だけに依存しているため, 下図のような格子点を考えれば

黒点は赤点のみに, 赤点は黒点のみに依存することになる. したがって, 赤点においては

\phi_{i\ j}^{(k+1)}=\frac{1}{4}(\phi_{i+1\ j}^{(k)}+\phi_{i-1\ j}^{(k)}+\phi_{i\ j+1}^{(k)}+\phi_{i\ j-1}^{(k)})

対して黒点では

\phi_{i\ j}^{(k+1)}=\frac{1}{4}(\phi_{i+1\ j}^{(k+1)}+\phi_{i-1\ j}^{(k+1)}+\phi_{i\ j+1}^{(k+1)}+\phi_{i\ j-1}^{(k+1)})

と計算できる. これは式 $(9)$ よりも早く収束すると考えられる. このように $k+1$ 回目の漸化式に $k+1$ 回目の値も用いる方法をガウス＝ザイデル法という. さらに収束を早めた手法として SOR 法がある. SOR 法では, 加速パラメータ $\omega$ を $\phi$ の修正量に対して次のように導入する.

\phi_{i\ j}^{(k+1)}=\phi_{i\ j}^{(k)}+\omega(\phi_{i\ j}^{(k+1)}-\phi_{i\ j}^{(k)}) \tag{10}

上式の収束性は少なくとも $0<\omega<2$ でなければ保証されず^[1], $\omega=1$ のときガウス＝ザイデル法と同じになる. よって, $1<\omega<2$ の範囲で最適な $\omega$ を選ぶ必要がある. 矩形領域のラプラス方程式に対しては最適な $\omega$ が分かっており, 次式で与えられる^[2].

\begin{align*} \mu&=\frac{1}{2}\left(\cos\frac{\pi}{n_x}+\cos\frac{\pi}{n_y}\right)\\ \omega&=\frac{2}{1+\sqrt{1-\mu^2}} \end{align*}

実装

ソースコードは次のリポジトリに置いておく.

環境

Windows 11
Rust 1.64.0
gnuplot 5.4.5

格子間隔による計算精度

格子間隔 $\Delta$ を小さくするほど計算結果は正確になる. しかし, その分計算時間は長くなる. これを避けるために最初は大きな $\Delta$ から始め, その条件で得られた値を初期値として小さな $\Delta$ で計算することにする.

main.rs

fn main() {
   let stop: usize = 1;
   let mut solver = Solver::new(0.1, 1e-3, 30., 25., 0.4, 3., 20.);
   for i in 0..=stop {
      solver.solve();
      if i != stop {
         solver = Solver::create_double(&solver);
      }
   }
}

計算部分

実装は, 数値計算手法節の数式をそのままループで計算させればよい. 但し, 式 $(10)$ は桁落ちを考慮して, 次式に変形しておく.

\phi_{i\ j}^{(k+1)}=(1-\omega)\phi_{i\ j}^{(k)}+\omega\phi_{i\ j}^{(k+1)}

$\phi$ が変化しない電極内では初期化時にマスクしておき, 計算時に更新されないようにしている.

境界上の点では, インデックスの加算・減算時に格子数以上および 0 以下にならないようにすることで対応する.

出力

得られた計算結果から等電位線を出力する. 計算結果は .dat ファイルとして書き出す.

描画には gnuplot の等高線プロット機能を使う.

実際に $L_x=30,\ L_y=25,\ r=0.4,\ d=3,\ V=20,\ \epsilon=0.001$ で, $\Delta$ を $0.1$ から $0.125$ に変化させたときの計算結果を出力すると, 以下の図が得られた.

$ time cargo run
error = 6.235707296588218
error = 4.928699908037741
...
error = 0.0010016977552442796
error = 0.0009947205673150883

real    2m28.149s

最終誤差は $0.000995$ , 計算時間は $148.149\ \mathrm{s}$ となっており, 比較的高速に計算できていることが分かる.

参考文献

脚注

幸谷智紀. ソフトウェアとしての数値計算, 第9章連立一次方程式の解法2 — 反復法. 2007. https://na-inet.jp/nasoft/chap09.pdf, (参照 2022-11-17). ↩︎
山崎郭滋. 偏微分方程式の数値解法入門. 森北出版株式会社, 2015. ↩︎

GitHubで編集を提案