🔋

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

2022/12/04に公開約7,700字

概要

一様な厚さの導体に電極を取り付けて定常電流を流し, 電位の等しい点を結ぶことで等電位線を求める実験がある. 本記事では, この等電位線の実験を数値計算によりシミュレーションする. なお, 実装は 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}

ここで, \phix,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*}

実装

ソースコードは次のリポジトリに置いておく.
https://github.com/k-kuroguro/rust-electric-potential-solver

環境

  • Windows 11
  • Rust 1.64.0
  • gnuplot 5.4.5

格子間隔による計算精度

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

https://github.com/k-kuroguro/rust-electric-potential-solver/blob/554bc96d7e75ecf0d374fc4cc199efea16685217/src/solver.rs#L44-L94

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)}

https://github.com/k-kuroguro/rust-electric-potential-solver/blob/554bc96d7e75ecf0d374fc4cc199efea16685217/src/solver.rs#L168-L195

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

https://github.com/k-kuroguro/rust-electric-potential-solver/blob/554bc96d7e75ecf0d374fc4cc199efea16685217/src/solver.rs#L130-L150

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

https://github.com/k-kuroguro/rust-electric-potential-solver/blob/554bc96d7e75ecf0d374fc4cc199efea16685217/src/solver.rs#L201-L215

出力

得られた計算結果から等電位線を出力する. 計算結果は .dat ファイルとして書き出す.
https://github.com/k-kuroguro/rust-electric-potential-solver/blob/554bc96d7e75ecf0d374fc4cc199efea16685217/src/solver.rs#L106-L128

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

https://github.com/k-kuroguro/rust-electric-potential-solver/blob/554bc96d7e75ecf0d374fc4cc199efea16685217/src/main.rs#L17-L55

実際に L_x=30,\ L_y=25,\ r=0.4,\ d=3,\ V=20,\ \epsilon=0.001 で, \Delta0.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} となっており, 比較的高速に計算できていることが分かる.

参考文献

脚注
  1. 幸谷智紀. ソフトウェアとしての数値計算, 第9章 連立一次方程式の解法2 — 反復法. 2007. https://na-inet.jp/nasoft/chap09.pdf, (参照 2022-11-17). ↩︎

  2. 山崎郭滋. 偏微分方程式の数値解法入門. 森北出版株式会社, 2015. ↩︎

GitHubで編集を提案

Discussion

ログインするとコメントできます