Issue with numerical drift with two Neumann BCs

[muntasirms]

I've been working on a 1D transient PDE and have been running into some solver issues. I was wondering if I could ask for some advice. The problem describes the discharge of a biphasic electrode which is generally described as a fixed flux entering through one phase, and the same flux leaving at the other end of the domain out of the other phase. Here is the problem:

$\frac{dC_s}{dt} = \underbrace{k_m(C^\infty-C_a)}_{convective \space transport}-\underbrace{\frac{{nFk^0\left(C_ae^{\alpha_anF(\eta - E_{o})/RT} - (C^\infty-C_a) e^{\alpha_cnF(\eta - E_{o})/RT}\right)}}{nF}}_{source/sink}$

$\underbrace{(\frac{1}{\sigma}+\frac{1}{\kappa})^{-1}\frac{\partial^2 \eta}{\partial x^2} = -aC_{dl} \frac{\partial\eta}{\partial t}}_{electrochemical \space transport} + \underbrace{{anFk^0\left(C_ae^{\alpha_anF(\eta - E_{o})/RT} - (C^\infty-C_a) e^{\alpha_cnF(\eta - E_{o})/RT}\right)}}_{source/sink}$

The first equation is the mass balance around a particle and the second is a charge balance across the porous electrode.

with BCs

1. At (x = 0):
   $\left. \frac{\partial \eta}{\partial x} \right|_{x=0} = -\frac{I}{\sigma}$
   $\left. \frac{\partial {C_s}}{\partial x} \right|_{x=0} = 0$

2. At (x = L) :
   $\left. \frac{\partial \eta}{\partial x} \right|_{x=L} = \frac{I}{\kappa}$
   $\left. \frac{\partial {C_s}}{\partial x} \right|_{x=L} = 0$

and ICs:
$C_s = 0.5$
$\eta = E_o$ (generally set to zero)

I understand that invariant problems with two neumann boundary conditions are generally susceptible to this numerical drift in their solutions. My main question is how to address this: I've found suggestions for "pinning" a value to create a gauge - but this seems to change the physics by creating a Dirichlet condition at one end. This inherently changes the problem physics, which is an issue since I am looking for information about the charge balance alongside the final solution. Do you have any suggestions for other ways to address this drift? I've also attached a functional version of my code that shows the drift when run.

Thanks so much for your help and take care!
fipyHybridCapacitorBattery.txt

[guyer]

The internal pinning is not implemented correctly:

• The explicit term mask * largeValue * value must be combined with its implicit counterpart - ImplicitSourceTerm(coeff=mask * largeValue, var=etaBar).
• Even with the implicit part, the mask doesn't actually include any cells, so this condition doesn't do anything. Changing the mask to mask = (mesh.x >= 1 - dx / 2) & (mesh.x <= 1 + dx / 2) or mask = (mesh.x >= 0.5 - dx/2) & (mesh.x <= 0.5 + dx/2).
• Because of the way it works, pinning an internal value to zero doesn't really work. Try setting something very small, like value = 1e-5, combined with largeValue = 1e20 (1e12 evidently isn't big enough for this problem).

All that being said, you don't need this term. Because of the TransientTerm, these equations are well posed and don't have.

Regardless, I see no evidence of a Dirichlet condition at one end. I'm also not sure what drift you are experiencing. The solutions seems stable to me. Can you describe what you expect to see?