[Bug]: Hysteresis state changes under no current + OCP artifacts

[Daniel-Nicolae23]

PyBaMM Version

25.4.2

Python Version

3.12.7

Describe the bug

I'm seeing a weird behaviour when applying current pulses to an SPMe model with anode hysteresis. There is no multiphase, just 2 anode branches I made up:

I'm charging the cell at 1C in steps of 10%, then resting it for 10 minutes.

Using Wycisk and the algebraic surface form

model = pybamm.lithium_ion.SPMe(
    {
        "working electrode": "both",
        "surface form": "algebraic",
        "contact resistance": "false",
        "open-circuit potential": ("Wycisk", "single")
    }
)

First of all, the h state (which is initialised at -1) changes when the current is 0, which should not happen under the Wycisk ODE. Also, there are some bumps in the anode OCP (sol["Battery negative electrode bulk open-circuit potential [V]"].entries):

This was using decay rate 0.1 and switching factor 0. Increasing the decay rate affects the way the h-state changes in those regions with 0 current:

Surface form false

This way the h state no longer changes, but the artifacts in the OCP are still there:

Using current sigmoid

No more OCP artifacts, here the curves match well during the relaxations and we also get the expected jumps due to the current.

Steps to Reproduce

import numpy as np
import pybamm
import matplotlib.pyplot as plt

def graphite_ocp(sto):
    return (graphite_ocp_lit(sto) + graphite_ocp_delit(sto)) / 2

def graphite_ocp_delit(sto):
    u_eq = (
        1.9793 * np.exp(-15.3631 * (sto))
        + 0.2482
        - 0.0709 * np.tanh(29.8538 * (sto - 0.234))
        - 0.04478 * np.tanh(14.9159 * (sto - 0.3769))
        - 0.0205 * np.tanh(30.4444 * (sto - 0.7103))
        - 0.0205 * np.tanh(30.4444 * (sto - 0.9103))
    )

    return u_eq

def graphite_ocp_lit(sto):
    u_eq = (
        1.9793 * np.exp(-50.3631 * sto)
        + 0.2482
        - 0.0909 * np.tanh(35.8538 * (sto - 0.0834))
        - 0.04478 * np.tanh(14.9159 * (sto - 0.1769))
        - 0.0105 * np.tanh(30.4444 * (sto - 0.6103))
        - 0.0105 * np.tanh(30.4444 * (sto - 0.4103))
    )

    return u_eq


parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.update({
    "Initial hysteresis state in negative electrode": -1,
    "Negative particle hysteresis decay rate": 1,
    "Negative particle hysteresis switching factor": 0,
    "Negative electrode OCP [V]": graphite_ocp,
    "Negative electrode lithiation OCP [V]": graphite_ocp_lit,
    "Negative electrode delithiation OCP [V]": graphite_ocp_delit,
}, check_already_exists=False)


model = pybamm.lithium_ion.SPMe(
    {
        "working electrode": "both",
        "surface form": "false",
        "contact resistance": "false",
        "open-circuit potential": ("Wycisk", "single")
    }
)
model.insert_reference_electrode()


loop = [
    f"Charge at 1C for 6 minutes",
    "Rest for 10 minutes"
]

exp = pybamm.Experiment(["Rest for 1 second"] + loop * 9, 
period = '1 second')
sim = pybamm.Simulation(model, parameter_values=parameter_values, experiment=exp, solver=pybamm.IDAKLUSolver())
sol = sim.solve(initial_soc=0)


fig, axs = plt.subplots(1, 2, figsize=(18, 6), width_ratios=[1, 2])

axs[0].plot(sol["Time [s]"].entries, sol["X-averaged negative electrode hysteresis state"].entries)

ax = axs[0].twinx()
ax.plot(sol["Time [s]"].entries, sol["Current [A]"].entries)

axs[1].plot(sol["Time [s]"].entries,
            sol["Negative electrode 3E potential [V]"].entries)  
axs[1].plot(sol["Time [s]"].entries,
            sol["Battery negative electrode bulk open-circuit potential [V]"].entries)

Relevant log output

[Daniel-Nicolae23]

An extra note on this. Maybe it should be a different issue, but I'll start it here.

In the Wycisk documentation, the rate of the ODE is proportional to I/Q, which I guess is a way of normalising the current (the behavior should depend on the C-rate, but not on the cell dimensions). However, looking at the WyciskOpenCircuitPotential implementation in PyBaMM:

i_surf = variables[
            f"{Domain} electrode {phase_name}interfacial current density [A.m-2]"
        ]
h = variables[f"{Domain} electrode {phase_name}hysteresis state"]

dhdt = (
    K * (i_surf / (Q_cell * (dQdU**K_x))) * (1 - pybamm.sign(i_surf) * h)
)  #! current is backwards for a halfcell
self.rhs[h] = dhdt

Using i_surf is already a way of normalising for the cell dimensions, but we are still dividing by the capacity. This means that for a certain parameter set, if we change the ["Electrode height [m]", "Electrode width [m]", "Nominal cell capacity [A.h]"] to simulate a larger cell, that cell will have a lower h decay rate. For me, it would make the most sense to say:

$$\frac{dh}{dt} = (K * \text{C-rate}) * (1 - sgn ...)$$

In a very thin strip of electrode of capacity dQ, the reaction current is dI:

$$\text{C-rate}(x) = \frac{dI(x)}{dQ} = \frac{j dV \frac{3 \varepsilon}{R}} {dV \varepsilon c_{max} F} = \frac{3j}{RFc_{max}}$$

So we could have:

dhdt = (
    K * (i_surf / (R * c_max * pybamm.constants.F) * (dQdU**K_x)) * (1 - pybamm.sign(i_surf) * h)
) 

This should work for the DFN, since the $j$ is calculated at every x. And it could also work for the MPM if it's done separately for each $R$ in the radius domain discretisation, I don't know exactly how that works in the backend.

[rtimms]

Thanks for pointing this out.

1. I think this may be a bug in the way SPMe with capacitance plays with hysteresis. The variable "Negative electrode hysteresis state" in general depends on x so that each particle can have its own hysteresis state. When surface form is false in SPMe, then "Negative electrode hysteresis state" is equal everywhere, but not when it is true, see:
   
   So there is some kind of equilibration of the hysteresis state during rest. I need to look into this in more detail to provide a better answer.
2. The bumps appear to be a bug that is fixed in this Issue 4884 unify hysteresis #4893, I haven't dug into what the issue is, but here is the same plot made on that branch (with surface form false):
   
3. See Issue 4884 unify hysteresis #4893 and the linked issues for a discussion of the unusual units in the Wycisk model

[rtimms]

Zooming in on that plot reveals this isn't a bug. In SPMe, the interfacial current density can be a function of x, and therefore the hysteresis state can be too. Note that PyBaMM implements a "particle-wise" implementation of the Wycisk model, so the variable in the rhs is the interfacial current density (a function of x and t, which is changing even during rest).

[Daniel-Nicolae23]

Thanks for your answers @rtimms! Indeed, given these profiles for the interfacial current density, the model is behaving correctly (except for those OCP bumps you're already working on).

But I'm still a bit confused about this variation of the current density (this is turning from an issue into a discussion, I'm happy to take it elsewhere). During the charging, j varies quite a lot along x:

And this doesn't happen without the capacitance (surface form false). Wouldn't the double layer capacitors be fully charged everywhere after such a long time? I would expect that including them would only introduce some transient effects, but the j profile from above stays the same with time (even for longer charges).

[rtimms]

In the below example, you can see the capacitance timescale right at the start of the discharge. You end up with a nonuniform (in x) interfacial current density at steady state because, for SPMe, j+j_dl=i_app/a/L and j_dl is found by integrating phi_s-phi_e in time. In the below, "Negative electrode surface potential difference [V]" == phi_s-phi_e` is non-uniform in space and increasing with time.

surface_form = ["false", "differential", "algebraic"]
models = [
    pybamm.lithium_ion.SPMe(
        {
            "surface form": form,
        }
    )
    for form in surface_form
]
parameter_values = pybamm.ParameterValues("Chen2020")
parameter_values.update({"Negative electrode double-layer capacity [F.m-2]": 1e-8})
sols = []
for model in models:
    sim = pybamm.Simulation(model, parameter_values=parameter_values)
    sol = sim.solve([0, 3600])
    sols.append(sol)

pybamm.dynamic_plot(
    sols,
    [
        "Negative electrode volumetric interfacial current density [A.m-3]",
        "X-averaged negative electrode interfacial current density [A.m-2]",
        "Negative electrode surface potential difference [V]",
        "X-averaged negative electrode surface potential difference [V]",
    ],
    labels=surface_form,
)

[rtimms]

For the x-averaged models, the hysteresis eqn should use the x-averaged interfacial current density. I will fix this in #4893

[Daniel-Nicolae23]

Thanks @rtimms! Just to make sure I got the nomenclature right: the interfacial current density variable from PyBaMM is just the j from j+j_dl=i_app/a/L, right?

Also, I think you meant differentiate

and j_dl is found by integrating phi_s-phi_e in time.

I was indeed wrong about the role of the double layer capacitor: I thought it only charges up to the voltage $\eta_r$ from the electrokinetics, i.e. $\phi_e - \phi_s - U(x_s)$, and that it goes back to zero when no current is applied.

But its voltage does actually include the surface OCV component (in the PyBaMM convention, V_anode(x) = 0 + ohmic_anode(x) and the surface OCV is included in the electrolyte voltage, but I often think about it in the 3E convention of V_anode containing the OCV and the electrolyte being 0 + $c_e$ effects + ohmic, there is a surface jump to account for either way).

So it makes sense that the "X-averaged negative electrode surface potential difference [V]" changes, it's because the particles are being discharged of course. And the dl capacitors have to charge so they can produce the new potential difference due to OCV. And it also makes sense that the x-dependent difference is not constant, since that's the effect of concentration overpotential in the electrolyte (and ohmic).

But I'm still not convinced this explains it. If phi_s-phi_e has four components:

• Electrolyte concentration overpotential
• Electrolyte ohmic
• Solid ohmic (very small for graphite anyway)
• Anode OCV

Then the first 3 have zero derivative at steady-state (electrolyte concentration is steady and the ohmic is pre-determined for the SPMe). Plotting the negative particle surface stoichiometry during the charge shows a constant with x, meaning that the OCV should be constant with x at all times, and same its d/dt. That would mean j_dl is constant with x since all the dl capacitors charge equally.
So j too should remain constant, and adding the dl should only slightly lower it. Does that make sense or am I getting something wrong about the model?