概要
一様な厚さの導体に電極を取り付けて定常電流を流し, 電位の等しい点を結ぶことで等電位線を求める実験がある. 本記事では, この等電位線の実験を数値計算によりシミュレーションする. なお, 実装は Rust, 描画は gnuplot を用いて行う.
方程式の導出
電位を求める方程式を導出する.
Lx×Ly の矩形導体の中に半径 r の円筒形電極2つを 2d だけ離して置く. そして, 電極の中心を結ぶ方向に x 軸, 両電極の中点を通りこれと垂直な方向に y 軸をとる. また, 両電極の電位差は 2V である.
![](https://res.cloudinary.com/zenn/image/fetch/s--aKJd6huv--/c_limit%2Cf_auto%2Cfl_progressive%2Cq_auto%2Cw_1200/https://storage.googleapis.com/zenn-user-upload/deployed-images/61d49ed71011714f3b381447.png%3Fsha%3D2100c629f712b6d90e6da3730460ae0d47e289b9)
スカラーポテンシャル (= 電位) を ϕ とすれば次のラプラス方程式が成り立つ.
∂x2∂2ϕ+∂y2∂2ϕ=0(1)
ここで, ϕ は x,y の関数で, ϕ(x,y) と書ける. 数値計算のために上式を離散化する. すなわち, 格子点 (xi,yj) における ϕ の値を求める. nx+1×ny+1 の格子点を考え, 格子点 (xi+1,yj) における ϕ(xi+Δx,yj) を x 方向にテイラー展開すれば
ϕ(xi+Δx,yj)=ϕ(xi,yj)+∂x∂ϕΔx+2!1∂x2∂2ϕΔx2+3!1∂x3∂3ϕΔx3+⋯(2)
格子点 (xi−1,yj) においても同様に
ϕ(xi−Δx,yj)=ϕ(xi,yj)−∂x∂ϕΔx+2!1∂x2∂2ϕΔx2−3!1∂x3∂3ϕΔx3+⋯(3)
式(2),(3)の和をとれば, Δx の奇数次の項が消えて
ϕ(xi+Δx,yj)+ϕ(xi−Δx,yj)=2ϕ(xi,yj)−∂x2∂2ϕΔx2+O(Δx2)(4)
式(4)より
∂x2∂2ϕ=Δx21{ϕ(xi+Δx,yj)+ϕ(xi−Δx,yj)−2ϕ(xi,yj)}(5)
y 方向にも同様に
∂y2∂2ϕ=Δy21{ϕ(xi,yj+Δy)+ϕ(xi,yj−Δy)−2ϕ(xi,yj)}(6)
式(5),(6)の和をとれば, 式(1)の左辺を ϕ の値で表すことができる. ϕi j≡ϕ(xi,yj), Δx=Δy とすれば
∂x2∂2ϕ+∂y2∂2ϕ=ϕi+1 j+ϕi−1 j+ϕi j+1+ϕi j−1−4ϕi j=0(7)
ϕi j について解くと
ϕi j=41(ϕi+1 j+ϕi−1 j+ϕi j+1+ϕi j−1) (i=0,1,...,nx, j=0,1,...,ny)(8)
よって, 格子点 (xi,yj) での ϕ の値は, その周囲4格子点の平均値であると示された.
![](https://res.cloudinary.com/zenn/image/fetch/s--rgDbYfHY--/c_limit%2Cf_auto%2Cfl_progressive%2Cq_auto%2Cw_1200/https://storage.googleapis.com/zenn-user-upload/deployed-images/ada8bd3df573ffeca1f62adf.png%3Fsha%3Dd73b8f3130d882d1c09a384ffb4d5ad4a7bf3d22)
また, 式(8)は差分方程式であり, 右辺には境界値 ϕ0 j,ϕnx j,ϕi 0,ϕi ny が含まれる. よって, 解くためにはこれらの値を定める境界条件が必要である. 今回は次の2条件を考える.
電極については, 両電極の電位差が 2V であるから, 左の電極を正とすれば
{ϕ(x,y)∣(x−d)2+y2≤r2=−Vϕ(x,y)∣(x+d)2+y2≤r2=V
端面では, 境界上と領域外の電位差が0 (∂x∂ϕ=0 または ∂y∂ϕ=0) になるように領域外の ϕ の値を決めてやればよい. つまり, 境界上の ϕ と領域外の ϕ を等しくすればよい.
![](https://res.cloudinary.com/zenn/image/fetch/s--K0JOs2ZK--/c_limit%2Cf_auto%2Cfl_progressive%2Cq_auto%2Cw_1200/https://storage.googleapis.com/zenn-user-upload/deployed-images/f2fe6bcba07c5a5bdec0c42f.png%3Fsha%3D55f3e44a16848c1c5a342922442099b80a322ee2)
i=nx のとき, 式(8)は
ϕnx j=41(ϕnx j+ϕnx−1 j+ϕnx j+1+ϕnx j−1)=31(ϕnx−1 j+ϕnx j+1+ϕnx j−1)
また, i=nx, j=ny のときは
ϕnx ny=41(ϕnx ny+ϕnx−1 ny+ϕnx ny+ϕnx ny−1)=21(ϕnx−1 ny+ϕnx ny−1)
となる. 他の境界上の点でも同様となるため, 端面では領域外の点を除いた2または3点の平均をとればよいと分かる.
数値計算手法
式(7)を反復法によって解く. 反復法は適当な初期値から始めて, 漸化式を反復計算していくことで解を得る数値計算の手法である. k 回目の反復における ϕi j の値を ϕi j(k) としたとき, 式(8)は
ϕi j(k+1)=41(ϕi+1 j(k)+ϕi−1 j(k)+ϕi j+1(k)+ϕi j−1(k))(9)
ϕi j(0) に適当な値を与え, ϕi j(1),ϕi j(2),...,ϕi j(k) と計算していく. そうすれば徐々に解は収束していく. そして誤差が許容誤差 ϵ 以下になったときの ϕi j(k) を解とする. なお, 誤差は次式で与える.
max ∣ϕi j(k+1)−ϕi j(k)∣
式(9)のように, k+1 回目の値を k 回目の値で求める反復法をヤコビ法という. 式(9)より, ϕi j は周囲4格子点だけに依存しているため, 下図のような格子点を考えれば
![](https://res.cloudinary.com/zenn/image/fetch/s--XQlCQ8V---/c_limit%2Cf_auto%2Cfl_progressive%2Cq_auto%2Cw_1200/https://storage.googleapis.com/zenn-user-upload/deployed-images/cfaa781a13a44c9fe6503f6a.png%3Fsha%3D896519f7ccac90280d7b7c8c0936be7ed51f41b2)
黒点は赤点のみに, 赤点は黒点のみに依存することになる. したがって, 赤点においては
ϕi j(k+1)=41(ϕi+1 j(k)+ϕi−1 j(k)+ϕi j+1(k)+ϕi j−1(k))
対して黒点では
ϕi j(k+1)=41(ϕi+1 j(k+1)+ϕi−1 j(k+1)+ϕi j+1(k+1)+ϕi j−1(k+1))
と計算できる. これは式(9)よりも早く収束すると考えられる. このように k+1 回目の漸化式に k+1 回目の値も用いる方法をガウス=ザイデル法という. さらに収束を早めた手法として SOR 法がある. SOR 法では, 加速パラメータ ω を ϕ の修正量に対して次のように導入する.
ϕi j(k+1)=ϕi j(k)+ω(ϕi j(k+1)−ϕi j(k))(10)
上式の収束性は 少なくとも 0<ω<2 でなければ保証されず, ω=1 のときガウス=ザイデル法と同じになる. よって, 1<ω<2 の範囲で最適な ω を選ぶ必要がある. 矩形領域のラプラス方程式に対しては最適な ω が分かっており, 次式で与えられる.
μω=21(cosnxπ+cosnyπ)=1+1−μ22
実装
ソースコードは次のリポジトリに置いておく.
https://github.com/k-kuroguro/rust-electric-potential-solver
環境
- Windows 11
- Rust 1.64.0
- gnuplot 5.4.5
格子間隔による計算精度
格子間隔 Δ を小さくするほど計算結果は正確になる. しかし, その分計算時間は長くなる. これを避けるために最初は大きな Δ から始め, その条件で得られた値を初期値として小さな Δ で計算することにする.
https://github.com/k-kuroguro/rust-electric-potential-solver/blob/554bc96d7e75ecf0d374fc4cc199efea16685217/src/solver.rs#L44-L94
計算部分
実装は, 数値計算手法 節の数式をそのままループで計算させればよい. 但し, 式(10)は桁落ちを考慮して, 次式に変形しておく.
ϕi j(k+1)=(1−ω)ϕi j(k)+ωϕi j(k+1)
https://github.com/k-kuroguro/rust-electric-potential-solver/blob/554bc96d7e75ecf0d374fc4cc199efea16685217/src/solver.rs#L168-L195
ϕ が変化しない電極内では初期化時にマスクしておき, 計算時に更新されないようにしている.
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.001 で, Δ を 0.1 から 0.125 に変化させたときの計算結果を出力すると, 以下の図が得られた.
![](https://res.cloudinary.com/zenn/image/fetch/s--6GSmXENy--/c_limit%2Cf_auto%2Cfl_progressive%2Cq_auto%2Cw_1200/https://storage.googleapis.com/zenn-user-upload/deployed-images/5b061b660722867fdfa2d5f1.png%3Fsha%3D2092860c9f2ae7c1c694f6c945d5a6816f2c3f42)
最終誤差は 0.000995, 計算時間は 148.149 s となっており, 比較的高速に計算できていることが分かる.
参考文献
Discussion