I have scoured the documentation and internet for examples, but cannot find one that explains how to impose a boundary condition that obeys the Tafel or Butler-Volmer equations. This occurs wherever an electrolyte is in contact with an electrode and passes current, so it must be pretty commonplace.

The potential is proportional to the log of current, and the current is proportional to the slope of the potential. Therefore, I would want to solve for u where the boundary obeys:

u = a + b log(du/dn) (Tafel equation)

where u is electrical potential, n is the vector normal to the boundary, and a and b are constants. The best I have come up with is this:

mesh Th = square(10, 10);

// The boundaries are labelled 1, 2, 3, 4 from bottom anti-clockwise

fespace Vh(Th, P1);

Vh u, v;

func f = 0.0;

func a = 1.0;

func b = 1.0;

// Copied from FreeFEM documentation Page 197

problem Pw (u, v)

= int2d(Th)(dx(u) * dx(v) + dy(u) * dy(v))

- int1d(Th, 1)(exp(a * u) * v)

- int2d(Th)(f * v)
- int1d(Th, 4)(exp(b) * v)

- on(2, 3, u = 0);

Pw;

// Plot the result

plot(u);

But it raises an error on line 14:

14 : + int1d(Th, 1)(exp(a * u) error operator ( <10LinearCombI7MGauche4C_F0E>

which seems to imply that you can’t have a non-linear function like exp() in a boundary condition. And yet I know of a couple of academic papers where they claim to have solved this problem using FreeFEM. However, they didn’t publish their source code!