Python三次元ポアソン方程式による静電場解析 (1)点電荷の電場解析

#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次元ラプラス・ポアソン方程式のヤコビ法による数値解法，楕円型偏微分方程式，境界値問題