A possible bug related to the treatment of boundary storage and boundary species

[JanSuntajs]

Hello to all,

after a while, I got back to the simulations of the diffusion phenomena in batteries using VoronoiFVM (see also issue #131). I am happy to report that the simulation of the galvanostatic mode which I was dealing with in the said issue now works just fine, however, I am having issues trying to implement terms containing storage terms defined for the boundary species.

Problem formulation

The problem I am describing is a standard diffusion equation with the Butler-Volmer flux at the boundary (see the attached MWE for more information and a longer description) and included double layer effects. The latter states
the relation between the voltage on the electrode and the actual flux entering the electrode particle as:

$$ \frac{\mathrm{d} U}{\mathrm{d} t} = \frac{1}{C} \left(j_\mathrm{app} - j_\mathrm{BV}(c_\mathrm{s}, U)\right). $$

Here, $U$ is the voltage, $C$ the double layer capacitance, $j_\mathrm{app}$ the driving flux (in the MWE, it is sinusoidal in time), $j_\mathrm{BV}(c_s, U)$ is the Butler-Volmer flux, which is a function of the surface concentration $c_s$ and voltage $U$.
In the absence of the double layer effects, the driving flux and the Butler-Volmer flux are the same; in this setting, however, the Butler-Volmer flux becomes the Neumann boundary condition for the concentration variable.

Attempt at a solution

I have tried to solve this problem by defining the voltage as surface species. Then, I defined the corresponding boundary storage operators and boundary conditions (see the attached MWEs in fvm_05_double_layer_debug.jl or fvm_05_double_layer_debug.ipynb). It seems, however, that the boundary storage terms are neglected altogether, as the solutions seem to be the same as in the case where no double layer effects are considered. Note that a clear indicator of the double layer effects is a mismatch between the applied and Butler-Volmer fluxes at higher frequencies. I'd assume that the culprit could be in the assembly procedure, perhaps the boundary storage is not considered for surface species? In the absence of the double layer effects, the approach works just fine, which I've benchmarked by performing some EIS calculations and verified them against known data.

Attempt at a workaround

Surprisingly, I get more sensible results if I define the voltage as a global species on all the control volumes (which is a bit redundant, since we assume the voltage is homogenous across the particle). In this formulation, I also dropped the specification of the boundary storage. In this case, I clearly get the desired mismatch between the applied and Butler-Volmer currents, however, a solution where this would also work for surface species would be desired.

The attached MWEs (fvm_05_double_layer_debug_semi_working.ipynb and fvm_05_double_layer_semi_working.jl) will also produce animations and relevant plots. I am using VoronoiFVM 2.3.0, the solver of choice in both cases is FBDF() from OrdinaryDiffEq. For convenience, MWE's are prepared both as Julia scripts or Jupyter notebooks.

@j-fu
VoronoiFVM_MWE.zip

[j-fu]

Hi, thanks for the update.
Just after having a short look: did you try the example with the inbuilt implicit Euler solver ?
Not that it is better suited, but it is better tested. If that works, then there may be an error in the adapter to DifferentialEquations.

Some other hints/questions:

• You seem to use an uniform grid. But the higher the frequency, the thinner the boundary layer at the excited electrode, and IMHO the grid should be finer there. ExtendableGrids has a geomspace method which could help with this.
• Can you provide me with an evironment (Project.toml+ Manifest.toml) ? This makes it easier to run your testst. (see my take on julia worflows: https://j-fu.github.io/marginalia/julia/)
• How did you perform EIS ? Did you use the impedance mode of VoronoiFVM ?

[JanSuntajs]

Hi,

I did not try the built-in implicit Euler solver, as the FBDF() and QNDF() seemed to be the only reliable solvers able to catch the quick oscillations at higher frequencies. But I agree, this is indeed worth exploring, if I recall correctly, my original problem was rather similar - surface species worked just fine with the built-in ImplicitEuler solver, but they failed once I tried to use DifferentialEquations. For now, it looks as if the storage term for the boundary species is overlooked for some reason.

Regarding other hints:

• in an actual simulation setting, I have the possibility of generating grids either equidistant in $r$ or of equal volume. In a simulation, I typically use the latter which are finer at the boundary and thus more suitable for the finite penetration depth. However, my simulation workflow is a part of a larger DrWatson-managed project and I wanted to avoid all the complexities of sharing the whole project in preparing a MWE. The problem persists regardless of the choice of the mesh.
• since I was unsure how to incorporate the Butler-Volmer fluxes (with their possible surface concentration dependence) into the impedance mode of VoronoiFVM, I did it my own way - namely by numerically solving the diffusion equation at different driving frequencies and then performing the Fourier analysis of the obtained time signals (with some windowing applied to account for initial transient effects). I also wanted to keep the workflow rather general so that one would just have to change the type of boundary driving if one wanted to try out a different type of simulation, for instance DC (dis)charging.
• Again, since my actual simulation code is a part of a larger project, I am including just the bits needed to get the MWEs up and running in the Project.toml and Manifest.toml files:

Archive.zip

Thanks for your prompt response and I hope this helps.
Best,
Jan

[j-fu]

Ok there were indeed bugs with mass matrix calculation. If the CI is ok, this can be released ASAP as v2.9.1.

Concerning the equivalence with the "semi working example" one should be careful - the geometric factors multiplying the time derivative of the voltage are different.

[JanSuntajs]

Ok there were indeed bugs with mass matrix calculation. If the CI is ok, this can be released ASAP as v2.9.1.

Concerning the equivalence with the "semi working example" one should be careful - the geometric factors multiplying the time derivative of the voltage are different.

Thank you for prompt treatment of the issue. Since the results of the two approaches I've sent do not match even after updating to the latest version of VoronoiFVM, I was wondering about the last point you made - could you be more specific about that? What should I multiply the voltage term with?

[j-fu]

A simplified version of what you seem to intend to do involves diffusion a nonlinear Robin boundary condition in a spherically symmetric domain $\Omega= (0,R)$:

• $\partial_t u + j = 0$ in $\Omega$
• $ j = -D \nabla u $ in $\Omega$
• $ j \cdot n + b(u) = 0 $ on $\partial \Omega$ (equivalent to $r=R$)

Now the way this is implemented in VoronoiFVM takes a subdivision $0=r_1 \dots r_n=R$ and assumes that the control volumes are the following:

• $\omega_1$ is a ball around 0 of radius $r_1/2$
• $\omega_2 \dots \omega_{n-1}$ are spherical shells $\omega_i$ with inner radius $\frac {r_{i-1} +r_i}{2}$ and outer radius $\frac {r_{i+1} +r_i}{2}$
• $\omega_n$ is a spherical shell of inner radius $\frac {r_{n-1} +r_n}{2}$ and outer radius $r_n$.
  The volumes of $\omega_i$ are calculated as the volumes of spherical shells.

The boundary condition at $r=R$ is integrated over the full sphere of radius $R$, so it is multiplied by $4\pi R^2$

With your volume trick, you just placed the boundary reaction in the control volume $\omega_n$, so it is multiplied with the spherical shell volume $|\omega_n|$.

If, in your trick example, you multiply the reaction term by $4\pi R^2 / |\omega_n|$, the results should be the same as with the boundary reaction.

[JanSuntajs]

Thank you, after correcting for the geometrical factors, the proposed approach with the boundary reaction now indeed works!

Thank you once again.

[j-fu]

Ok @JanSuntajs I close this - feel free to open other issues if needed.