1. はじめに
二次元ポアッソン方程式による電場解析例(1)はQiitaにおいても紹介されているが，電場解析に必要な三次元の例はなかったので，トライしてみた．
本稿では，下記を実施する．
(1)三次元ポアソン方程式を離散化
(2)繰り返し法によるプログラムを作成．
①メッシュ幅は1.0と0.5を選択できるようにする．
②電位を計算
③得られた電位から電場を計算
(3)計算結果を理論値と比較することにより解析の妥当性を検証
なお，今後帯電した円盤，球について報告したいと考えています．
2. 例題の問題設定
各軸のメッシュ数は50とし，中心点(25,25,25)に点電荷を配置する．
図1 問題設定-点電荷の配置
3. ポアソン方程式の離散化
下記に手順を記す．
ここで注意が必要な点はρは電荷密度であることである．
4. プログラム
import pprint
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.cm as cm
delta_L=0.5 # delta_L=1.0
LN=50
HLN=int(LN/2)
nx = LN
ny = LN
nz = LN
xmin = 0
xmax = LN*delta_L
ymin = 0
ymax = LN*delta_L
zmin = 0
zmax = LN*delta_L
p = np.zeros((nz, ny, nx))
pd = np.zeros((nz, ny, nx))
charge = np.zeros((nz, ny, nx))
x = np.linspace(xmin, xmax, nx)
y = np.linspace(ymin, ymax, ny)
z = np.linspace(zmin, zmax, nz)
eps0=1 # とりあえず
y00=1
Q=1
q=Q/delta_L**3 # 電荷密度に変換
charge[HLN][HLN][HLN]=q
for I in range(500): #2000 # 大きくすると精度は上がるが時間がかかる
pd = p.copy()
p[1:-1,1:-1,1:-1] = (pd[2:,1:-1,1:-1] + pd[:-2,1:-1,1:-1] +
pd[1:-1,2:,1:-1] + pd[1:-1,:-2,1:-1] +
pd[1:-1,1:-1,2:] + pd[1:-1,1:-1,:-2] + (delta_L**2)*charge[1:-1, 1:-1, 1:-1]/eps0) /6
p[0][:][:]=0
p[nx-1][:][:]=0
p[:][0][:]=0
p[:][ny-1][:]=0
p[:][:][0]=0
p[:][:][nz-1]=0
def plot2D(x, y, p):
fig = plt.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(x, y)
surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)
ax.view_init(30, 225)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('p')
plot2D(x,y,p[:,:,HLN])
def plot_surface(pp):
nx = LN
ny = LN
fig = plt.figure(figsize=(11,7), dpi=100)
fig
xmin = 0
xmax = LN-1
ymin = 0
ymax = LN-1
x = np.linspace(xmin, xmax, nx)
y = np.linspace(ymin, ymax, ny)
X, Y = np.meshgrid(x, y)
plt.contourf(X, Y, pp, alpha=0.1, cmap=cm.viridis)
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')
plot_surface(p[:,:,HLN])
plt.style.use('ggplot')
plt.rcParams["axes.facecolor"] = 'white'
fig = plt.figure()
ax = fig.gca(projection='3d')
L=5
NLS=HLN-2#25
NLY=5
X, Y, Z= np.meshgrid(np.arange(0, L,1), np.arange(0, L,1), np.arange(0, NLY,1))
Ex = np.zeros([L,L,abs(NLY)])
Ey = np.zeros([L,L,abs(NLY)])
Ez = np.zeros([L,L,abs(NLY)])
for i in range(L):
for j in range(L):
for k in range(0,NLY):
ii=i+HLN-2
jj=j+HLN-2
kk=k+HLN-2
Ex[i,j,k]=-(p[ii,jj+1,kk]-p[ii,jj-1,kk])/2/delta_L
Ey[i,j,k]=-(p[ii+1,jj,kk]-p[ii-1,jj,kk])/2/delta_L
Ez[i,j,k]=-(p[ii,jj,kk+1]-p[ii,jj,kk-1])/2/delta_L
# 電位の理論値計算，表示
print("Riron Keisan V=")
for y00 in range(1,5):
x2=1/(4*3.14)/(y00*delta_L)
print('z=',y00*delta_L+25,' V=',x2)
# 電位の計算結果表示
print("V=")
for i in range(26,33):
print('z=',25+(i-25)*delta_L,' V=',p[HLN,HLN,HLN+(i-25)])
# 電場の理論値計算，表示
print("Riron Keisan E=")
for y00 in range(1,6):
x2=1/(4*3.14)/(y00*delta_L)**2
print('z=',y00*delta_L+25,' E=',x2)
# 電場の計算結果表示
i=25
j=25
for k in range(26,31):
print("E=",25+(k-25)*delta_L,"-> ",-(p[i,j,k+1]-p[i,j,k-1])/2/delta_L)
ax.quiver(X,Y,Z, Ex, Ey, Ez, color='red',length=6.0*delta_L**2, normalize=False)
X1,Y1, Z1=2,2,2
ax.scatter3D(X1,Y1,Z1,"o", color='blue')
X1,Y1, Z1=3,2,2
ax.scatter3D(X1,Y1,Z1,"o", color='blue')
X1,Y1, Z1=4,2,2
ax.scatter3D(X1,Y1,Z1,"o", color='blue')
#
plt.grid()
plt.draw()
plt.show()
5. 計算結果
5.1 delta_T=1.0の場合
図2(a)にz=25面の電位分布の2Dﾌﾟﾛｯﾄ図を，図2(b)に電荷点(25,25,25)を中心とした電場図を示す．
図2(a)電位分布
図2(b)電場
5.2 delta_T=0.5の場合
同様に，図3(a)にz=25面の電位分布の2Dﾌﾟﾛｯﾄ図を，図3(b)に電荷点(25,25,25)を中心とした電場図を示す．
電場図は全く同一に表されている．詳細は次章で説明するが本来なら４倍になるはずであるが，描画プログラムの中で電場ベクトルの長さを1/4にしているためである．
ax.quiver(X,Y,Z, Ex, Ey, Ez, color='red',length=6.0*delta_L**2, normalize=False)
図3(a)電位分布
図3(b)電場
6. 解析の妥当性検証
点電荷周りの電位，電場は理論的に下式で計算できる．
表1にdelta_L=1.0の場合の電位，電場の理論値と計算値の比較を示す．
表2にdelta_L=0.5の場合の電位，電場の理論値と計算値の比較を示す．
いずれの場合も，計算値は理論値にほぼ一致している．
また，delta_L=0.5の場合，電荷からの距離は1/2となるので，上記計算式で分かるように，電位は2倍，電場は4倍となっている．
表1 delta_L=1.0の場合の電位，電場の理論値と計算値の比較
表2 delta_L=0.5の場合の電位，電場の理論値と計算値の比較
参考文献
(1)[Pythonによる科学・技術計算] 静電位に対する2次元ラプラス・ポアソン方程式のヤコビ法による数値解法，楕円型偏微分方程式，境界値問題