Python三次元ポアソン方程式による静電場解析 (3)一様に帯電した球の電場解析

#1. はじめに
　Python三次元ポアソン方程式による静電場解析により，基本的課題として空間上の点電荷の電場解析が可能なことを先の投稿において報告した(1)．
　本稿では，応用課題として一様に帯電した球を対象に下記を実施する．
(1)三次元ポアソン方程式を差分化（再掲）
(2)繰り返し法によるプログラムを作成．
①メッシュ幅は1.0と0.5を選択できるようにする．
②電位を計算
③得られた電位から電場を計算
(3)計算結果を理論値と比較することにより解析の妥当性を検証
#2. 例題の問題設定
各軸のメッシュ数は50とし，中心に一様に帯電した球を配置する．
![設定.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/eff958cf-cedb-c527-06e8-9e7b836401f1.jpeg)
　　　　　　図1　問題設定-一様に帯電した球の配置

#3. ポアソン方程式の離散化(再掲)
　下記に手順を記す．
　ここで注意が必要な点はρは電荷密度であることである．
![poisson.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/1144eef4-3011-63e4-fa5c-fb224351f62f.jpeg)
 
#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
Q1=1/delta_L**3
          
for i in range(0, LN):
    for j in range(0,LN):
        for k in range(0,LN):
            rr=np.sqrt((i-HLN)**2+(j-HLN)**2+(k-HLN)**2)
            if rr <= 3.5 :
                charge[i,j,k] = Q1
            elif rr <= 4.5:
                charge[i,j,k] = Q1*(4.5-rr)
#
print("Riron")
for i in range(1,5):
    r=i*delta_L
    x2=-Q1*r**2/6+Q1*(4*delta_L)**2/2
    print('z=',r,' V=',x2)
for i in range(4,10):
    r=i*delta_L
    x2=Q1*(4*delta_L)**3/(3*r)
    print('z=',r,' V=',x2)    
for i in range(1,5):
    r=i*delta_L
    x2=Q1*r/3
    print('z=',r,'Ez=',x2)
for i in range(4,10):
    r=i*delta_L
    x2=Q1*(4*delta_L)**3/3/r**2
    print('z=',r,'Ez=',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
cs=p[HLN-4:HLN+8,HLN-4:HLN+8,HLN-4:HLN+8]

X, Y, Z= np.meshgrid(np.arange(0, L,1), np.arange(0, L,1), np.arange(0, L,1))  

Ex = np.zeros([L,L,L])
Ey = np.zeros([L,L,L])
Ez = np.zeros([L,L,L])
EE = np.zeros([L,L,L])

for i in range(L):
  for j in range(L):
    for k in range(L):
      Ex[i,j,k]=-(cs[i,j+1,k]-cs[i,j-1,k])
      Ey[i,j,k]=-(cs[i+1,j,k]-cs[i-1,j,k])
      Ez[i,j,k]=-(cs[i,j,k+1]-cs[i,j,k-1]) 

for i in range(L):
  for j in range(L):
    for k in range(L):
      ii=i*2+HLN-8
      jj=j*2+HLN-8
      kk=k*2+HLN-8     
               
      Ex[i,j,k]=-(p[ii,jj+1,kk]-p[ii,jj-1,kk])/2/delta_L#*2#*0.004#*10
      Ey[i,j,k]=-(p[ii+1,jj,kk]-p[ii-1,jj,kk])/2/delta_L#*2#*0.004#*10
      Ez[i,j,k]=-(p[ii,jj,kk+1]-p[ii,jj,kk-1])/2/delta_L#*2#*0.004#*10

print("V=")
for i in range(1,12):
  print('z=',i*delta_L,'Vz=',p[HLN,HLN,HLN+i])
  
i=25
k=25
for j in range(26,35):
    print("x=",(j-25)*delta_L,"Ex=",-(p[i,j+1,k]-p[i,j-1,k])/2/delta_L)

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

for i in range(-2,2):
    for j in range(-2,2):
        for k in range(-2,2):
            X1,Y1, Z1=int(L/2)+i,int(L/2)+j,int(L/2)+k
            ax.scatter3D(X1,Y1,Z1,"o", color='blue')

for i in range(5):
    X1,Y1, Z1=int(L/2)-i,int(L/2),int(L/2)
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2)+i,int(L/2),int(L/2)
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2),int(L/2)-i,int(L/2)
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2),int(L/2)+i,int(L/2)
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2),int(L/2),int(L/2)-i
    ax.scatter3D(X1,Y1,Z1,"o", color='blue')
    X1,Y1, Z1=int(L/2),int(L/2),int(L/2)+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)に電場図を示す．
![plot_2D.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/c97c2bac-7a4c-1cd3-9710-f47a879fa141.jpeg)

　　　　　　　　　　　　　　　図2(a)電位分布　　　　　　　　　　　　　　　

![E.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/0d392277-ecd1-6224-4558-fb8808564334.jpeg)

　　　　　　　　　　　　　　　図2(b)電場
　

#6. 解析の妥当性検証
6.1円板状電荷周りの電位，電場の理論的計算
図3(a)(b)に理論計算のためのモデルを示す．

![r_big_a.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/e62d30cf-b79d-3b1f-c521-001971a3ec00.jpeg)

　　図3(a) 理論計算のためのモデル　r>=aの場合
理論値は下記で計算できる．

![r_big_a.jpg_equation.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/b44f70f2-54c6-b7db-9777-374c577043f4.jpeg)


![r_small_a.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/0a1d351e-e121-446f-7d69-309a56dc82e2.jpeg)

　　図3(b) 理論計算のためのモデル　r<aの場合
理論値は下記で計算できる．
![r_small_a_equation.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/c2ced63a-cb6d-ba0f-8afb-d93d6abd1329.jpeg)

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

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

![table.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010670/a18e99d0-a402-3ba1-88f1-a04ecba1a7c1.jpeg)


#参考文献
(1)Python三次元ポアソン方程式による静電場解析 (1)点電荷の電場解析 Qiita
(2)Python三次元ポアソン方程式による静電場解析 (2)円板状電荷の電場解析 Qiita