1. はじめに

Python三次元ポアソン方程式による静電場解析により，基本的課題として空間上の点電荷の電場解析が可能なことを先の投稿において報告した(1)．
本稿では，応用課題として円板状電荷を対象に下記を実施する．
(1)三次元ポアソン方程式を差分化（再掲）
(2)繰り返し法によるプログラムを作成．
①メッシュ幅は1.0と0.5を選択できるようにする．
②電位を計算
③得られた電位から電場を計算
(3)計算結果を理論値と比較することにより解析の妥当性を検証
2. 例題の問題設定
各軸のメッシュ数は50とし，中心に帯電した円板状電荷を配置する．

　　　　　図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=1.0
# delta_L=0.5
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
Q1=1
q=Q1/(delta_L**3)
             
for i in range(0, LN):
    for j in range(0,LN):
        rr=np.sqrt((i-HLN)**2+(j-HLN)**2)
        if rr <= 3:
            charge[i,j,HLN] = q#/(delta_L**3)#*(delta_L**2)

print("Riron")
q1=Q1/(delta_L**2)
for i in range(0,8):
    y00=i*delta_L
    x2=q1*(np.sqrt((3*delta_L)**2+y00**2)-y00)/2
    print('z=',y00+25,' V=',x2)
    
for i in range(0,8):
    y00=i*delta_L
    x2=q1*(1-y00/np.sqrt((3*delta_L)**2+y00**2))/2
    print('Ez=',y00+25,'->',x2)

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]) /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.xlim([HLN-5,HLN+5])
plt.ylim([HLN-5,HLN+5])

plt.style.use('ggplot')
plt.rcParams["axes.facecolor"] = 'white'
fig = plt.figure()
ax = fig.gca(projection='3d')

L=9
NLS=HLN#25
NLY=2
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-4#21
      jj=j+HLN-4#21
      kk=k+NLS
      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
      Ez[i,j,k]=-(p[ii,jj,kk+1]-p[ii,jj,kk])/delta_L

print("Keisan")

for i in range(20,33):
  print('Vz=',25+(i-25)*delta_L,' ',p[HLN,HLN,i])
  
i=HLN
j=HLN
  
for k in range(HLN,HLN+6):
    print("Ez=",k,"-> ",-(p[i,j,k+1]-p[i,j,k-1])/2)
     print("Ez=",25+(k-25)*delta_L,"-> ",-(p[i,j,k+1]-p[i,j,k])/delta_L)

ax.quiver(X,Y,Z, Ex, Ey, Ez, color='red',length=1.0*delta_L**2, normalize=False)

for i in range(NLY):
    X1,Y1, Z1=4,4,i
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    

plt.grid()
plt.draw()
plt.show()

5. 計算結果-delta_T=1.0の場合

図2(a)にz=25面の電位分布の2Dﾌﾟﾛｯﾄ図を，図2(b)にz=25（電荷を配置）とz=26の電場図を示す．

　　　　　　　　　　　　　　図2(a)　電位分布　　　　　　　　　　　　　　　

　　　　　　　　　　　　　　図2(b)　電場

　また，図3にz=23,24,25,26層の電場を示す．上下対象になっていることが分かる．

　　　　　　　　　図3 z=23,24,25,26層の電場

6. 解析の妥当性検証

6.1円板状電荷周りの電位，電場の理論的計算

図4に理論計算のためのモデルを示す．

　　図4 理論計算のためのモデル

　理論値は下記で計算できる．

6.2 電位，電場の理論値と計算値の比較

表1にdelta_L=1.0の場合の電位，電場の理論値と計算値の比較を示す．
計算値は理論値にほぼ一致している．

　　　表1 delta_L=1.0の場合の電位，電場の理論値と計算値の比較

参考文献

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