Zenn
🔋

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

2022/12/05に公開

概要

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

方程式の導出

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

Lx×LyL_x\times L_y の矩形導体の中に半径 rr の円筒形電極2つを 2d2d だけ離して置く. そして, 電極の中心を結ぶ方向に xx 軸, 両電極の中点を通りこれと垂直な方向に yy 軸をとる. また, 両電極の電位差は 2V2V である.

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

2ϕx2+2ϕy2=0(1) \frac{\partial^2\phi}{\partial x^2}+\frac{\partial^2\phi}{\partial y^2}=0 \tag{1}

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

ϕ(xi+Δx,yj)=ϕ(xi,yj)+ϕxΔx+12!2ϕx2Δx2+13!3ϕx3Δx3+(2) \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}

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

ϕ(xiΔx,yj)=ϕ(xi,yj)ϕxΔx+12!2ϕx2Δx213!3ϕx3Δx3+(3) \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)(2),(3)の和をとれば, Δx\Delta x の奇数次の項が消えて

ϕ(xi+Δx,yj)+ϕ(xiΔx,yj)=2ϕ(xi,yj)2ϕx2Δx2+O(Δx2)(4) \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)(4)より

2ϕx2=1Δx2{ϕ(xi+Δx,yj)+ϕ(xiΔx,yj)2ϕ(xi,yj)}(5) \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}

yy 方向にも同様に

2ϕy2=1Δy2{ϕ(xi,yj+Δy)+ϕ(xi,yjΔy)2ϕ(xi,yj)}(6) \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)(5),(6)の和をとれば, 式(1)(1)の左辺を ϕ\phi の値で表すことができる. ϕi jϕ(xi,yj)\phi_{i\ j}\equiv\phi(x_i,y_j), Δx=Δy\Delta x=\Delta y とすれば

2ϕx2+2ϕy2=ϕi+1 j+ϕi1 j+ϕi j+1+ϕi j14ϕi j=0(7) \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}

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

ϕi j=14(ϕi+1 j+ϕi1 j+ϕi j+1+ϕi j1) (i=0,1,...,nx, j=0,1,...,ny)(8) \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}

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

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

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

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

{ϕ(x,y)(xd)2+y2r2=Vϕ(x,y)(x+d)2+y2r2=V \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 (ϕx=0\frac{\partial \phi}{\partial x}=0 または ϕy=0\frac{\partial \phi}{\partial y}=0) になるように領域外の ϕ\phi の値を決めてやればよい. つまり, 境界上の ϕ\phi と領域外の ϕ\phi を等しくすればよい.

i=nxi=n_x のとき, 式(8)(8)

ϕnx j=14(ϕnx j+ϕnx1 j+ϕnx j+1+ϕnx j1)=13(ϕnx1 j+ϕnx j+1+ϕnx j1) \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=nx, j=nyi=n_x,\ j=n_y のときは

ϕnx ny=14(ϕnx ny+ϕnx1 ny+ϕnx ny+ϕnx ny1)=12(ϕnx1 ny+ϕnx ny1) \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)(7)を反復法によって解く. 反復法は適当な初期値から始めて, 漸化式を反復計算していくことで解を得る数値計算の手法である. kk 回目の反復における ϕi j\phi_{i\ j} の値を ϕi j(k)\phi_{i\ j}^{(k)} としたとき, 式(8)(8)

ϕi j(k+1)=14(ϕi+1 j(k)+ϕi1 j(k)+ϕi j+1(k)+ϕi j1(k))(9) \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}

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

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

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

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

ϕi j(k+1)=14(ϕi+1 j(k)+ϕi1 j(k)+ϕi j+1(k)+ϕi j1(k)) \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)})

対して黒点では

ϕi j(k+1)=14(ϕi+1 j(k+1)+ϕi1 j(k+1)+ϕi j+1(k+1)+ϕi j1(k+1)) \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)(9)よりも早く収束すると考えられる. このように k+1k+1 回目の漸化式に k+1k+1 回目の値も用いる方法をガウス=ザイデル法という. さらに収束を早めた手法として SOR 法がある. SOR 法では, 加速パラメータ ω\omegaϕ\phi の修正量に対して次のように導入する.

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

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

μ=12(cosπnx+cosπny)ω=21+1μ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)(10)は桁落ちを考慮して, 次式に変形しておく.

ϕi j(k+1)=(1ω)ϕi j(k)+ωϕi j(k+1) \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

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

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

real    2m28.149s

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

参考文献

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

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

GitHubで編集を提案

Discussion

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