Problem while solving a minimal 1D Poisson-Nernst-Planck equation

[nadav7679]

I am trying to solve the Poisson-Nernst-Planck equation in Python with the library FiPy. It is basically a set of equations that describes - in my example - the separation of two ion concentrations in a solution with a potential gradient. Like mixing salt in water and then applying a voltage-difference between both ends of the solution.
The equations are:

$$\frac{\partial C_p}{\partial t} = \nabla^2 C_p + \nabla(C_p\nabla \phi)$$

$$\frac{\partial C_n}{\partial t} = \nabla^2 C_n - \nabla(C_n\nabla \phi)$$

$$\nabla^2 \phi = C_n - C_p$$

And the boundary conditions:
$$\nabla C_p (x=L_x) = \nabla C_p (x=0) = \nabla C_n (x=L_x) = \nabla C_n (x=0)=0$$

$$\phi(L_x) = 1, \phi(0) = -1$$

Problem

I manage to get a converging solution, but it is not the solution i'd like. Specifically, I get that C_p and C_n always converge to be spatially constant, and $\phi$ is always linear. It seems like $\phi$ is indifferent to initial conditions. Even setting internal fixed point (based on this answer) turns $\phi$ to be linear with a break point, which doesn't help.

The solution I look for should be that $\phi$ is more like a sigmoid, and C_p and C_n get high values around the edges and low values in the center.

I really would appreciate help!

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

[guyer]

I assume you intended this?
$$\frac{\partial C_p}{\partial t} = \nabla^2 C_p + \nabla(C_p\nabla \phi)$$
$$\frac{\partial C_n}{\partial t} = \nabla^2 C_n - \nabla(C_n\nabla \phi)$$

[guyer]

A few things:

• Your initial condition of $C_n = C_p = 0.5$ means that net charge is everywhere zero, so no driving force for curvature in $\phi$ or migration of $C_n$ or $C_p$.

  Initializing to, e.g.,

  Cn.setValue(0.5, where=x < Lx / 2)
  Cp.setValue(0.5, where=x >= Lx / 2)

  will produce some non-linear solution for $\phi$.

• The default coefficient of ImplicitSourceTerm is 0. This is non-obvious and non-helpful. Hopefully we'll change the default sooner rather than later. Set those coefficients to 1.

• Do you know the Debye length for your system? You may need to experiment with epsilon to get the sigmoid you're looking for.

[nadav7679]

I assume you intended this? ∂Cp∂t=∇2Cp+∇(Cp∇ϕ) ∂Cn∂t=∇2Cn−∇(Cn∇ϕ)

Thanks! I Changed it accordingly.

[nadav7679]

A few things:

• Your initial condition of Cn=Cp=0.5 means that net charge is everywhere zero, so no driving force for curvature in ϕ or migration of Cn or Cp.
  Initializing to, e.g.,

  Cn.setValue(0.5, where=x < Lx / 2)
  Cp.setValue(0.5, where=x >= Lx / 2)
  

  will produce some non-linear solution for ϕ.

• The default coefficient of ImplicitSourceTerm is 0. This is non-obvious and non-helpful. Hopefully we'll change the default sooner rather than later. Set those coefficients to 1.

• Do you know the Debye length for your system? You may need to experiment with epsilon to get the sigmoid you're looking for.

Thanks a lot!
The ImplicitSourceTerm coefficient really did the trick.

I still didnt achive my goal though. The goal is to to start with uniform concentrations and coverge to a solution where the positive and negative charges seperate to the sides.

I expirimented with the value of epsilon, the magnitude of the potential drop, and the size of the system. As of now nothing seems to work. Any further ideas would be welcome.

[guyer]

I still didnt achive my goal though. The goal is to to start with uniform concentrations and coverge to a solution where the positive and negative charges seperate to the sides.

I understand why you would expect that to happen, but I think the math doesn't support this. In your initial condition, $\nabla^2 c_n = \nabla^2 c_p = 0$, so there's no diffusion. $\nabla\phi = \text{constant}$, as are $c_n$ and $c_p$, so there's no divergence in electromigration flux, hence no electromigration. Therefore $\partial c_n / \partial t = \partial c_p / \partial t = 0$. Your initial condition is a steady-state solution.

There should be a way to induce charge separation, but I'm going to have to think about it. I've always dealt with electrode-electrolyte interfaces or pn junctions, where there was an initial heterogeneity.