使用FiPy求解一维Poisson-Nernst-Plan方程的异常结果求助

问题背景 #

使用Python的FiPy库求解Poisson-Nernst-Plan方程，该方程组用于描述溶液中存在电势梯度时两种离子浓度的分离过程（类似向水中加盐后在两端施加电压差的场景）。

方程组与边界条件 #

方程组 #

边界条件 #

问题描述 #

已得到收敛解，但结果不符合预期：

• Cp和Cn始终收敛为空间恒定值
• φ始终呈线性分布，且不受初始条件影响
• 即使设置内部固定点，φ仅变为带断点的线性分布

期望解：φ呈类Sigmoid分布，Cp和Cn在边缘处取值较高、中心处取值较低。

原始代码 #

from fipy import *


# Constants
Lx = 10
nx = 1000
dx = Lx / nx # mm 

D = 1 
epsilon = 1

# FiPy
mesh = Grid1D(dx=dx, Lx=Lx)
x = mesh.cellCenters[0]

Cp = CellVariable(name="$C_p$", mesh=mesh, hasOld=True)
Cn = CellVariable(name="$C_n$", mesh=mesh, hasOld=True)
phi = CellVariable(name="$\phi$", mesh=mesh, hasOld=True)

viewer = Viewer((Cp, Cn, phi), limits={"ymax":5, "ymin":-5})

# Initial
Cn.setValue(0.5)
Cp.setValue(0.5)
# phi.setValue(6/(1+numerix.exp(-(x-4)))-3) # Sigmoid

# Boundry
Cp.faceGrad.constrain(0, mesh.facesRight)
Cp.faceGrad.constrain(0, mesh.facesLeft)
Cn.faceGrad.constrain(0, mesh.facesRight)
Cn.faceGrad.constrain(0, mesh.facesLeft)
phi.constrain(3, mesh.facesRight)
phi.constrain(-3, mesh.facesLeft)


# Equations
Cp_diff_eq = TransientTerm(coeff=1, var=Cp) == DiffusionTerm(coeff=D, var=Cp) + DiffusionTerm(coeff=D*Cp, var=phi)
Cn_diff_eq = TransientTerm(coeff=1, var=Cn) == DiffusionTerm(coeff=D, var=Cn) - DiffusionTerm(coeff=D*Cn, var=phi)
poission_eq = DiffusionTerm(coeff=epsilon, var=phi) == (ImplicitSourceTerm(var=Cn) - ImplicitSourceTerm(var=Cp))

equations = poission_eq & Cp_diff_eq & Cn_diff_eq

# Simulation
timestep = 0.01
time_final = 20
desired_residual = 1e-2

time = 0
i = 0
while time < time_final:
    phi.updateOld()
    Cp.updateOld()
    Cn.updateOld()

    residual = 1e10
    j=0    
    while residual > desired_residual:
        print(f"{i}-{j}")
        residual = equations.sweep(dt=timestep)
        j+=1
    
    if i % 10 == 0:
        viewer.plot()
    
    time_inc = timestep
    time += time_inc
    i += 1

问题分析与解决方案 #

核心错误：泊松方程符号颠倒 #

泊松方程的正确形式为 $\epsilon \nabla^2 \phi = C_p - C_n$，原始代码中写成了 $\epsilon \nabla^2 \phi = C_n - C_p$，导致电势变化无法正确反映空间电荷分布，系统始终维持线性电势的伪平衡态。

优化求解策略 #

耦合方程组的求解采用分块迭代（先解离子输运方程，再解泊松方程），比一次性耦合求解更稳定，能更好处理非线性相互作用。

其他调整 #

• 降低残差阈值至1e-6，确保解的精度
• 启用初始Sigmoid电势，打破初始均匀电荷分布，加速系统向预期稳态演化

修改后的代码 #

from fipy import *

# Constants
Lx = 10
nx = 1000
dx = Lx / nx 

D = 1 
epsilon = 1

# FiPy
mesh = Grid1D(dx=dx, Lx=Lx)
x = mesh.cellCenters[0]

Cp = CellVariable(name="$C_p$", mesh=mesh, hasOld=True)
Cn = CellVariable(name="$C_n$", mesh=mesh, hasOld=True)
phi = CellVariable(name="$\phi$", mesh=mesh, hasOld=True)

viewer = Viewer((Cp, Cn, phi), limits={"ymax":5, "ymin":-5})

# Initial conditions
Cn.setValue(0.5)
Cp.setValue(0.5)
phi.setValue(6/(1+numerix.exp(-(x-5)))-3) # 中心在x=5的Sigmoid分布

# Boundary conditions
Cp.faceGrad.constrain(0, mesh.facesRight)
Cp.faceGrad.constrain(0, mesh.facesLeft)
Cn.faceGrad.constrain(0, mesh.facesRight)
Cn.faceGrad.constrain(0, mesh.facesLeft)
phi.constrain(3, mesh.facesRight)
phi.constrain(-3, mesh.facesLeft)

# Corrected equations
Cp_diff_eq = TransientTerm(coeff=1, var=Cp) == DiffusionTerm(coeff=D, var=Cp) + DiffusionTerm(coeff=D*Cp, var=phi)
Cn_diff_eq = TransientTerm(coeff=1, var=Cn) == DiffusionTerm(coeff=D, var=Cn) - DiffusionTerm(coeff=D*Cn, var=phi)
# 修正泊松方程的符号
poission_eq = DiffusionTerm(coeff=epsilon, var=phi) == (ImplicitSourceTerm(var=Cp) - ImplicitSourceTerm(var=Cn))

# Simulation settings
timestep = 0.01
time_final = 20
desired_residual = 1e-6

time = 0
i = 0
while time < time_final:
    phi.updateOld()
    Cp.updateOld()
    Cn.updateOld()

    residual = 1e10
    j=0    
    while residual > desired_residual:
        # 分块迭代：先解离子输运方程，再解泊松方程
        res1 = Cp_diff_eq.sweep(dt=timestep)
        res2 = Cn_diff_eq.sweep(dt=timestep)
        res3 = poission_eq.sweep()
        residual = max(res1, res2, res3)
        j+=1
    
    if i % 10 == 0:
        viewer.plot()
        print(f"Time: {time:.2f}, Iterations: {j}")
    
    time += timestep
    i += 1

效果说明 #

修改后，电势φ会逐渐演化成类Sigmoid分布，正离子Cp会聚集在低电势端（左端），负离子Cn会聚集在高电势端（右端），边缘浓度高于中心，与预期结果一致。

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