How do I deal with nonlinear Neumann BC and Laplace equation?

Hi everyone,

I’m trying to simulate something similar to what has been done in this paper (it’s a rather short read):

Briefly summarized, the problem looks like this:

My aim is to calculate an electric potential in an electrolyte which is adjacent an anode (index a) and a cathode (index c) in a rectangular domain (the bottom boundary of the domain is split equally into a cathode and an anode boundary). The electric potential in the electrolyte domain is calculated by the Laplace equation \nabla^{2} \varphi=0. The two Neumann boundary conditions (for anode and cathode) are represented by \frac{\partial \varphi}{\partial n}=\frac{i_{a,c}}{\sigma}, where i_{a,c} is the current density at the anode and cathode, respectively, \sigma is the conductivity of the electrolyte. The current densities of anode and cathode are calculated by the following equations
i_{a,c}=i_{0,a,c}(\mathrm{exp}(\frac{\eta}{\alpha})-\mathrm{exp}(\frac{\eta}{-\alpha})). The remaining boundary conditions are no flux. There are no Dirichlet boundaries for now.

This is where my problem arises. \eta is defined as \varphi_{a,c}-\varphi. Here, \varphi_{a,c} are material constants which makes my Neumann boundary condition depend on the solution of the Laplace equation. I have trouble implementing this nonlinearity and get my simulation to converge.

Things I have tried

  • using a constant \eta, which works fine but this would just be a simplification of the physical process.
  • using the solution of the constant \eta as an initial guess. Does not work.

I’m using the latest fenics debian package and I also used the latest stable docker image. I’m new to fenics and programming, so any advice would be helpful and is greatly appreciated. Thank you for reading.

A .geo file for the geometry can be found at the end of the post. I have provided two sets of Neumann BC; one set for a constant \eta, and one set where \eta is a function of u. Only use one set in the simulations and comment out the other set. Here is my code so far:

# import necessary libraries
from fenics import *
import matplotlib.pyplot as plt

# defining constants
sigma_1 = 10  # 1/(ohm*m), conductivity of electrolyte
alpha = 0.05
i0_anode = 1.0  # A/m², exchange current density
i0_cathode = 1.0  # A/m², exchange current density
pot_cathode = 0.5  # V, potential of anode
pot_anode = -0.5  # V, potential of cathode
eta = 0.28  # V, constant overpotential

# import mesh and boundaries with tags
mesh = Mesh()

with XDMFFile("mesh.xdmf") as infile:
	infile.read(mesh)

mvc = MeshValueCollection('size_t', mesh, 1) 
with XDMFFile('mf.xdmf') as infile:
	infile.read(mvc, 'Boundary')

boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc) 

# define boundaries in the mesh
ds_anode = Measure("ds", domain=mesh, subdomain_data=boundaries, subdomain_id=4)  # call boundary markers
ds_cathode = Measure("ds", domain=mesh, subdomain_data=boundaries, subdomain_id=5)  # call boundary markers

# set up function space
V = FunctionSpace(mesh, 'P', 1)
v = TestFunction(V)
u = Function(V)
f = Constant(0)

# Neumann BC
#g_anode = Expression('i0_anode*exp((pot_anode - u)/alpha)-exp((pot_anode - u)/-alpha)/sigma_1', u=u, i0_anode=i0_anode, pot_anode=pot_anode, alpha=alpha, sigma_1=sigma_1, degree=2)
#g_cathode = Expression('i0_cathode*exp((pot_cathode - u)/-alpha)-exp((pot_cathode - u)/alpha)/sigma_1', u=u, i0_cathode=i0_cathode, pot_cathode=pot_cathode, alpha=alpha, sigma_1=sigma_1, degree=2)

g_anode = Expression('(i0_anode*(exp((eta)/alpha)-exp((eta)/-alpha)))/sigma_1',  i0_anode=i0_anode, pot_anode=pot_anode, alpha=alpha, eta=eta, sigma_1=sigma_1, degree=2)
g_cathode = Expression('(i0_cathode*(exp((eta)/-alpha)-exp((eta)/alpha)))/sigma_1', i0_cathode=i0_cathode, pot_cathode=pot_cathode, alpha=alpha, eta=eta, sigma_1=sigma_1, degree=2)

# solve
bc = []
F = dot(grad(u), grad(v)) * dx - f * v * dx + g_anode * v * ds_anode + g_cathode * v * ds_cathode
solve(F == 0, u, bc)

# some plotting
electric_field = project(-grad(u))
plt.figure()
p = plot(project(u, V))
plot(electric_field)
plt.colorbar(p)
#plt.savefig('output/potential.eps')
plt.show()

.geo file:

// Gmsh project created on Tue Sep  8 13:53:25 2020
//+
SetFactory("OpenCASCADE");
//+
Point(1) = {-0.0001, 0, 0, 1.0};
//+
Point(2) = {-0.0001, 0.0001, -0, 1.0};
//+
Point(3) = {0.0001, 0.0001, -0, 1.0};
//+
Point(4) = {0.0001, 0, -0, 1.0};
//+
Point(5) = {0, 0, -0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 5};
//+
Line(5) = {5, 1};
//+
Curve Loop(1) = {2, 3, 4, 5, 1};
//+
Plane Surface(4) = {1};
//+
Physical Line(1) = {1};
//+
Physical Line(2) = {2};
//+
Physical Line(3) = {3};
//+
Physical Line(4) = {4};
//+
Physical Line(5) = {5};
//+
Physical Surface(0) = {4};
//+
Transfinite Curve {1, 3, 4, 5} = 20 Using Progression 1;
//+
Transfinite Curve {2} = 39 Using Progression 1;