求解最小一维泊松-能斯特-普朗克方程时出现的问题

Question

我正在尝试使用 FiPy 库在 Python 中求解 Poisson-Nernst-Planck 方程。在我的例子中，它基本上是一组方程，描述了具有电位梯度的溶液中两个离子浓度的分离。就像将盐混合在水中，然后在溶液两端施加电压差一样。方程为：

以及边界条件：

问题

我设法得到一个收敛的解决方案，但这不是我想要的解决方案。具体来说，我发现 C _p和 C _n始终收敛为空间常数，并且始终是线性的。似乎对初始条件漠不关心。即使设置内部固定点（基于此答案）也会变得与断点成线性，这没有帮助。

我寻找的解决方案应该更像是一个 sigmoid，并且 C _p和 C _n在边缘周围获得高值，在中心获得低值。

我真的很感激帮助！

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