#1. はじめに
Python三次元ポアソン方程式による静電場解析に関しては先の投稿において報告した(1)(2)(3).
本稿では,二次元ポアソン方程式による静電場解析として十分長い直線状電荷分布を対象に下記を実施する.
(1)二次元ポアソン方程式を差分化
(2)繰り返し法によるプログラムを作成.
①電位を計算
②得られた電位から電場を計算
(3)計算結果を理論値と比較することにより解析の妥当性を検証
#2. 例題の問題設定
各軸のメッシュ数は100とし,中心に一様に帯電した十分長い直線を配置する.
図1 問題設定-十分長い直線状電荷の配置
#3. ポアソン方程式の離散化
下記に手順を記す.
ここで注意が必要な点はρは線密度であることである.
#4. プログラム
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
import math
from numpy import log as ln
delta_L=1.0
nx = 100
ny = 100
nt = 10000
xmin = 0
xmax = 100*delta_L
ymin = 0
ymax = 100*delta_L
dx = (xmax - xmin) / (nx - 1)
dy = (ymax - ymin) / (ny - 1)
p = np.zeros((ny, nx))
pd = np.zeros((ny, nx))
b = np.zeros((ny, nx))
x = np.linspace(xmin, xmax, nx)
y = np.linspace(xmin, xmax, ny)
eps0=1
charge= np.zeros((ny, nx))
Q1=1
Q2=Q1/delta_L**2
charge[50,50] = Q2
for it in range(10000):
pd = p.copy()
p[1:-1,1:-1] = (pd[1:-1, 2:] + pd[1:-1, :-2] +
pd[2:, 1:-1] + pd[:-2, 1:-1] + (delta_L**2/(1*eps0))*charge[1:-1, 1:-1] ) / 4
p[0, :] = 0
p[ny-1, :] = 0
p[:, 0] = 0
p[:, nx-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)
fig = plt.figure(figsize=(11,7), dpi=100)
fig
xmin = 0
xmax = 100
ymin = 0
ymax = 100
x = np.linspace(xmin, xmax, nx)
y = np.linspace(xmin, xmax, ny)
X, Y = np.meshgrid(x, y)
plt.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Y')
plt.ylim([35,65])
plt.xlim([35,65])
fig = plt.figure(figsize=(11,7), dpi=100)
fig
x1 = np.linspace(0, 100, nx)
y1 = np.linspace(xmin, xmax, ny)
for i in range(nx):
y1[i]=p[i,50]
plt.plot(x1,y1)
plt.title("A simple line graph")
plt.xlim([44,56])
L= 99
Ex = np.zeros([L,L])
Ey = np.zeros([L,L])
for i in range(L):
for j in range(L):
Ey[i,j] = -(p[i+1,j]-p[i-1,j])/2/delta_L
Ex[i,j] = -(p[i,j+1]-p[i,j-1])/2/delta_L
print("Riron Keisan E")
for i in range(1,6):
print("x=",i*delta_L,"E=",Q1/(2*np.pi)/(i*delta_L))
print("Riron Keisan V")
for i in range(1,6):
print("x=",i*delta_L,"V=",-Q1/(2*np.pi)*ln(i*delta_L))
print("Keisan E")
for i in range(1,6):
print("x=",i*delta_L,"E=",Ey[50+i,50])
print("Keisan V")
for i in range(1,6):
print("x=",i*delta_L,"V=",p[50+i,50]-p[50+int(1/delta_L),50])
fig = plt.figure(figsize=(11,7), dpi=100)
fig
X, Y= np.meshgrid(np.arange(0, L,1), np.arange(0, L,1)) # メッシュ生成
plt.quiver(X,Y,Ex,Ey,color='red',angles='xy',scale_units='xy', scale=0.25/delta_L) # ベクトル場をプロット
plt.xlim([45,55])
plt.ylim([45,55])
plt.grid()
plt.draw()
plt.show()
#5. 計算結果-delta_T=1.0の場合
図2(a)にz=50面の電位分布の2Dプロット図を,図2(b)に電場図を示す.
参考に,直線状の電荷が二本の場合のベクトル図を図3に示す.
#6. 解析の妥当性検証
##6.1十分長い直線状電荷周りの電位,電場の理論的計算
図4に理論計算のためのモデルを示す.
##6.2 電位,電場の理論値と計算値の比較
表1にdelta_L=1.0の場合の電位,電場の理論値と計算値の比較を示す.
計算値は理論値にほぼ一致している.
表1 delta_L=1.0の場合の電位,電場の理論値と計算値の比較
参考文献
(1)Python三次元ポアソン方程式による静電場解析 (1)点電荷の電場解析 Qiita
(2)Python三次元ポアソン方程式による静電場解析 (2)円板状電荷の電場解析 Qiita
(3)Python三次元ポアソン方程式による静電場解析 (3)一様に帯電した球の電場解析 Qiita