OCPs, Hysteresis, etc

[tommaull]

Hi everyone, I wanted to open up a discussion about how we handle OCPs. It seems like a submodel might be the best way to handle this. Right now OCPs are a function of concentration, and in the sigmoid hys model, current direction. As far as I know, entropic effects aren't included either. I think it'd be quite helpful for users to be able to choose their OCP function based on need, so parameters like decay rate in hysteresis can be included. This might also help with MSMR, as OCP is not longer a function of concentration.

[brosaplanella]

About the entropic term, I think it might be useful to define OCP as a function of stoichiometry and temperature directly (like with electrolyte diffusivity, for example) so it is more general. Then the entropic term can be computed from the OCP. The only issue is that, even though this makes a lot of sense from the modelling point of view, it is not how the parameters are usually measured...

[ejfdickinson]

"can be computed from the OCP" - how, in practice?

If you can do symbolic differentiation, fine. If you are reliant on numerical differentiation, feels dangerous.

[rtimms]

If users provide a function we can do symbolic differentiation, but would need to be numerical for data. My vote would still be to give OCP at fixed T, and then an entropic term -- this is the only way I've ever seen the OCP being given/collected.

[ejfdickinson]

@rtimms I'd agree.

[rtimms]

Just to comment on structure, the OCP is already a separate submodel where the three options are:

1. single_ocp - OCP is given as a function of stoichiometry
2. current_sigmoid_ocp - lithiation and delithiation OCPs are given as separate functions of stoichiometry and switched between using a sigmoid function
3. msmr_ocp - OCP comes from the MSMR formulation (solve for potential in the particles and provide stoichiometry as a function of potential, see e.g. Baker 2018)

In all of these, the entropic term is given separately and the OCP has linear temperature dependence U(x, T) = U(x, T0) + (T - T0)*dU(x)/dT), and in the current sigmoid option, the entropic term is the same in both directions.

It is possible to add new submodels for OCP that have new parameters and solve new equations. In principle, these could be whatever you like (e.g. Plett-style model based on some definition of electrode SOC or one that depends on the local extent of lithiation). The only other place this would impact the code would be the thermal model, which calculates the reversible heat source based on the entropic change. But so long as your OCP had linear temperature dependence with an entropic term, then any new OCP submodels would be fairly self-contained.

[tommaull]

Hahahaha, I'm not proud of this habit I have of requesting things that already exist. Tried searching discussions/docs and poking through submodels before I posted but never opened "interface".

[rtimms]

Haha no worries, it just tells us our docs need improving!

[tommaull]

@rtimms do you agree that the best place to add a more plett-like hysteresis model, running off a single OCP?

[rtimms]

Yes, I would add a new submodel in the open_circuit_potential folder. A rough outline:

• the code would use the current f"{Domain} electrode OCP [V]" parameter as the "equilibrium" OCP
• add new parameters for hysteresis voltage and decay rate
• add a new ODE for a hysteresis variable (details up to the user)
• compute the OCP as "equilibrium" OCP plus f(hysteresis voltage, hysteresis variable) which now depends on time as a result of the ODE for the hysteresis variable

[ejfdickinson]

@rtimms @tommaull

I noticed that the current 0th-order hysteresis model implemented by CurrentSigmoidOpenCircuitPotential references the global scalar "Total current density [A.m-2]" on the cell.

In a spatially-resolved model (DFN), there's a good argument that this variable should probably be evaluated locally, based on the faradaic current density experienced by the specific particle for which open-circuit potential is being computed. Local current density is not necessarily zero, even when cell current is zero.

A local hysteresis model supports better prediction of cell relaxation at rest when different parts of an electrode are in different non-equilibrium states. It's also important if you have a blended electrode in a non-equilibrium state, where one particle type in the blend can lithiate the other, again for cell relaxation at rest. For example with graphite-Si you might want a hysteresis model for the Si but not for the graphite.

• Was this option decided against at original implementation due to numerical stability issues, maybe?
  • A 0th-order hysteresis model with spatial resolution might cause overpotential and current density to change abruptly or in an oscillatory manner, depending on how local mesh size compares to the sigmoid constant k.
• A possible alternative is to leave the existing implementation alone, but implement @tommaull 's suggestion in Decay rate for single state hysteresis models #3375 with spatial resolution. Using a fast decay rate based on local current density would then allow a safe approach to the numerically unstable limit.

[rtimms]

Tagging @mbonkile to comment on why the total current density was chosen for this model. I believe it was simply a choice of a simple model to implement, rather than any numerical reason, but I may be wrong.

I agree it would be good to implement a model that uses the local current density. My suggestion would be to implement this as an extra submodel with its own option (e.g. "interfacial current-sigmoid"). That way the "current-sigmoid" option can be used to reproduce the results of the paper where the implementation came from. As you say, the current 0th-order model does not make good predictions during relaxation. This may be improved for spatially resolved models in the "local" version, but not for averaged models like SPM.

A Plett-style hysteresis model that is a function of local extent lithiation with a decay rate would be a very welcome addition too, and allow for much better predictions of hysteresis during rest and over longer timescales. I don't think it would be too difficult to implement.