二维0V边界平行板电容器电势绘图异常问题排查

问题概述 #

在二维0V边界的箱体中模拟平行板电容器，已知极板电压与间距，采用带超松弛的Gauss-Seidel方法计算电势，初始假设极板电荷密度ρ（代码中为p）为0，绘图结果出现两个异常：

• 电势分布左偏，不符合平行板的对称特征
• 电场矢量指向正极板，与理论上“电场从正极指向负极”的结论相反

原始代码 #

import numpy as np
import matplotlib.pyplot as plt

L=100
a=1
e0=8.85e-12
p=0
w=0.9
iterations=0
tolerance=1e-6
dvmax=100
dec=4
vpot=np.zeros((L,L))
# 设置极板电压：左侧极板（x=20）为1V，右侧极板（x=80）为-1V
vpot[20:80,20:21]=1
vpot[20:80,80:81]=-1
xy = np.linspace(0,L-1,L)
xydec=xy[::dec]

while dvmax>tolerance:
      iterations+=1
      dvmax=0
      for i in range(1,L-1):
          for j in range(1,L-1):
              vnew=(vpot[i-1,j]+vpot[i+1,j]+vpot[i,j-1]+vpot[i,j+1]+(p*(a**2)/e0))/4
              dv=(1+w)*(vnew-vpot[i,j])
              if abs(dv)>dvmax:
                 dvmax=abs(dv)
              vpot[i,j]+=dv

print(f'Number of iterations: {iterations}')

plt.imshow(vpot, origin='lower')
plt.plasma()
plt.colorbar()
plt.title(f'Plot of Positive Potentials and Centered Charge Density with Tolerance {tolerance}')
plt.figure()
cv=plt.contour(vpot,15,colors='k', linestyles='solid')
plt.clabel(cv, fontsize=8, inline=True)
plt.title('Contour of Potentials')
plt.figure()

Ex=np.zeros((L,L))
Ey=np.zeros((L,L))
for j in range(1,L-1):
          for k in range(1,L-1):
              if k==0:
                 Ex[j,k] = -(vpot[j,k+1] - vpot[j,k])/a
              elif k==L-1:
                 Ex[j,k] = -(vpot[j,k] - vpot[j,k-1])/a
              else:
                 Ex[j,k] = -(vpot[j,k+1] - vpot[j,k-1])/(2*a)
              if j==0:
                 Ey[j,k] = -(vpot[j+1,k] - vpot[j,k])/a
              elif j==L-1:
                 Ey[j,k] = -(vpot[j,k] - vpot[j-1,k])/a
              else:
                 Ey[j,k] = -(vpot[j+1,k] - vpot[j-1,k])/(2*a)

Exdec = Ex[::dec,::dec]
Eydec = Ey[::dec,::dec]
plt.quiver(xydec,xydec,Exdec,Eydec)
plt.title('Electric Field Vectors')
plt.show()

问题根源分析 #

1. 极板边界条件被迭代覆盖：
   迭代过程中没有排除极板所在的网格点，导致预先设置的极板电压被Gauss-Seidel迭代修改，破坏了边界条件，这是电势分布左偏的核心原因。

2. 电场计算的索引混淆：
   代码中vpot[i,j]的i对应y轴（行）、j对应x轴（列），但电场计算时循环变量j和k的对应关系错误，导致电场梯度的计算方向反转，最终矢量方向与理论相反。

3. 电荷密度的误解：
   平行板内部区域（真空）电荷密度ρ=0的设定是正确的，极板的电荷是边界条件的结果，不需要预先设定。我们通过固定极板电压（狄利克雷边界条件）已经间接定义了极板电荷，无需额外设置ρ值。

修正后的代码 #

import numpy as np
import matplotlib.pyplot as plt

L=100
a=1
e0=8.85e-12
p=0
w=0.9
iterations=0
tolerance=1e-6
dvmax=100
dec=4
vpot=np.zeros((L,L))
# 定义极板区域：左侧极板（x=20，y从20到79），右侧极板（x=80，y从20到79）
left_plate = (slice(20,80), slice(20,21))
right_plate = (slice(20,80), slice(80,81))
vpot[left_plate] = 1
vpot[right_plate] = -1

xy = np.linspace(0,L-1,L)
xydec=xy[::dec]

while dvmax>tolerance:
      iterations+=1
      dvmax=0
      for i in range(1,L-1):
          for j in range(1,L-1):
              # 跳过极板区域，保持固定电压
              if (20<=i<80 and j==20) or (20<=i<80 and j==80):
                  continue
              vnew=(vpot[i-1,j]+vpot[i+1,j]+vpot[i,j-1]+vpot[i,j+1]+(p*(a**2)/e0))/4
              dv=(1+w)*(vnew-vpot[i,j])
              if abs(dv)>dvmax:
                 dvmax=abs(dv)
              vpot[i,j]+=dv

print(f'Number of iterations: {iterations}')

# 绘制电势图，指定extent确保x/y轴对应正确
plt.imshow(vpot, origin='lower', extent=[0, L-1, 0, L-1])
plt.plasma()
plt.colorbar(label='Potential (V)')
plt.title('Parallel Plate Capacitor Potential Distribution')
plt.xlabel('X')
plt.ylabel('Y')
plt.figure()

# 绘制等势线
cv=plt.contour(vpot, 15, colors='k', linestyles='solid', extent=[0, L-1, 0, L-1])
plt.clabel(cv, fontsize=8, inline=True)
plt.title('Equipotential Lines')
plt.xlabel('X')
plt.ylabel('Y')
plt.figure()

# 计算电场：Ex = -dV/dx，Ey = -dV/dy
Ex=np.zeros((L,L))
Ey=np.zeros((L,L))
for i in range(1, L-1):  # i对应y轴（行）
    for j in range(1, L-1):  # j对应x轴（列）
        # 计算x方向电场Ex
        if j == 0:
            Ex[i,j] = -(vpot[i,j+1] - vpot[i,j])/a
        elif j == L-1:
            Ex[i,j] = -(vpot[i,j] - vpot[i,j-1])/a
        else:
            Ex[i,j] = -(vpot[i,j+1] - vpot[i,j-1])/(2*a)
        # 计算y方向电场Ey
        if i == 0:
            Ey[i,j] = -(vpot[i+1,j] - vpot[i,j])/a
        elif i == L-1:
            Ey[i,j] = -(vpot[i,j] - vpot[i-1,j])/a
        else:
            Ey[i,j] = -(vpot[i+1,j] - vpot[i-1,j])/(2*a)

# 抽取稀疏的电场矢量用于绘图
Exdec = Ex[::dec,::dec]
Eydec = Ey[::dec,::dec]
plt.quiver(xydec, xydec, Exdec, Eydec, scale=50)
plt.title('Electric Field Vectors')
plt.xlabel('X')
plt.ylabel('Y')
plt.xlim(0, L-1)
plt.ylim(0, L-1)
plt.show()

修正说明 #

1. 保护极板边界条件：在迭代循环中加入判断，跳过极板所在的网格点，确保极板电压始终保持预设值，保证电势分布符合平行板的对称特征。
2. 修正电场索引对应关系：明确i对应y轴、j对应x轴，正确计算电势的梯度方向，使得电场矢量从正极板（1V）指向负极板（-1V），符合理论预期。
3. 优化绘图坐标对应：在imshow和contour中添加extent参数，确保x/y轴与网格坐标一致，避免视觉上的混淆。
4. 电荷密度的正确处理：内部区域ρ=0的设定保持不变，极板电荷由边界条件自动隐含，无需额外设置，因为我们求解的是拉普拉斯方程（ρ=0），极板为狄利克雷边界。

内容的提问来源于stack exchange，提问作者Giau Diep