CODEWITHSUNDEEP
pythonnumerical-methodspdechemistryfipy
Nadav
Nadav

Reputation: 51

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

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:

Equations

And the boundary conditions:

Boundary Conditions

Problem

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

The solution I look for should be that 𝜙 is more like a sigmoid, and Cp and Cn 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

Upvotes: 5

Views: 320

Answers (0)

Related Questions