pybamm-team
Unify hysteresis models and add heat generation terms
Description
We now have three options for modelling open-circuit potential hysteresis in PyBaMM: 1) "current sigmoid"; 2) "Wycisk"; and 3) "Axen", all of which have slightly different implementations.
Both "Wycisk" and "Axen" are single-state hysteresis models which contain a single ODE for a hysteresis state variable, h, of the form dh/dt = K(1-sgn(i))h, but again with slightly different implementations. The current sigmoid model is a zero-state model (the OCP just jumps between branches depending on the sign of the applied current), and doesn't track the variable h
We propose we unify these models by (see https://github.com/pybamm-team/PyBaMM/issues/4332#issuecomment-2296943415):
-
Making sure all hysteresis models use the state variable
-1<h<1(in the current sigmoid modelhwill just jump from -1 to 1) -
Add a base one-state model that solves a general one-state equation:
$\frac{\partial h}{\partial t} = \gamma\cdot\left(\frac{i_\mathrm{surf}a_\mathrm{vol}}{\phi_\mathrm{act}Fc_\mathrm{max}}\right)\cdot \left(1 - \mathrm{sgn}(i_\mathrm{surf})h\right)$
for $i_\mathrm{surf}$ the interfacial phase-specific current density (anodic current density positive).
- Allow a scalar input or arbitrary functional input for $\gamma$, or allow $\gamma$ to be computed using the specific Wycisk formula with parameters $K$, $n$ and the OCP:
$\gamma = K \cdot \left(\frac{\frac{\partial q_\mathrm{vol}}{\partial U}}{\left.\frac{\partial q_\mathrm{vol}}{\partial U}\right|_\mathrm{ref}}\right)^{-n}$
where the evaluation at a reference point (implicitly at $x_\mathrm{Li}=x_\mathrm{Li,ref}$ which in turn defines $K$ as the decay rate at the reference stoichiometry) is to ensure that we take an arbitrary power only of a unitless quantity.
- Allow general (particle-wise) initialisation of the hysteresis state, not necessarily = 0.
Unifying the models will also allow us to more reasonably include hysteresis heating (see https://github.com/About-Energy-OpenSource/AEPyBaMM/blob/main/src/aepybamm/pybamm_tools.py#L91).
Motivation
No response
Possible Implementation
No response
Additional context
No response
@rtimms Looks good.
I think we should try to find a formulation here that enables us to also solve #3867 in a general way.
Specifically, in terms of $h$, how are we relating this to the OCP of the branches and the "average" OCP? By my understanding, we would have:
$U_\mathrm{eff} = \frac{1+h}{2}U_\mathrm{delit}+ \frac{1-h}{2}U_\mathrm{lit}$
But $h = 0$ need not correspond to $U_\mathrm{eqm}$, the actual equilibrium condition from a thermodynamic perspective. The general hysteresis heat source would be:
$Q_\mathrm{hys,vol} = i_\mathrm{far,vol}\left(U_\mathrm{eff} - U_\mathrm{eqm}\right)$
How should we parameterise this, and do you think the back-end variable implementation of of Wycisk with a variable "Negative electrode OCP hysteresis [V]", which is then scaled by $h$, is the best way?
How should we parameterise this, and do you think the back-end variable implementation of of Wycisk with a variable "Negative electrode OCP hysteresis [V]", which is then scaled by h, is the best way?
One approach would be to have the user pass all three of "Negative electrode OCP [V]", "Negative electrode lithation OCP [V]" and "Negative electrode delithation OCP [V]" and treat the first as the actual equilibrium condition.
That makes sense to me. Could we make the former parameter optional and compute it as the mean of the two branches if unspecified, so that the breaking change arises only when it was actively specified?