Another method for OCP with hysteresis

[chtamar]

Description

This method has been proposed by Axen et al. https://doi.org/10.1016/j.est.2022.103985, and like the Wycisk approach, it employs an ODE system to track the evolution of the ocp between the empirical lithiation and delithiation branches of the hysteresis. The rate at which the ocp moves back and forth between the hysteresis branches is scaled by the applied current and a decay rate.

Motivation

Although the Wycisk and Axen methods perform similarly, the proposed implementation does not need to compute the differentiation self.phase_param.U(sto_surf, T_bulk).diff(sto_surf), and it also allows to choose separate decay rates for the lithiation and delithiation branches of the hysteresis. This option can be implemented for the Wycisk approach as well.

Together with the Wycisk approach, this method could be another alternative to choose from in pybamm.

A figure is attached comparing Wycisk and Axen for the silicon's ocp, sweeping over the decay rate (K) and temperature (T), using the Chen2020_composite parametrization and an OCP entropic change for silicon of 2e-3 V/K, for testing purposes. Here the lithiation and delithiation decay rates for Axen take the same value.

Possible Implementation

The method can be readily implemented using the module wycisk_ocp as starting point.

The method would need the following input parameters:

• The hysteresis branches for lithiation and delithition, ocp_lith and ocp_delith, respectively.
• The hysteresis decay rate for lithiation K_lith
• The hysteresis decay rate for delithiation K_delith
• The initial value for the state variable h_init

Main features of the method:

• As in the Wycisk approach, the state variable would be h, taking values between 0 and 1
• The rhs expression of the state variable would be scaled by the volumetric interfacial current density i_vol and the decay rates K_lith and K_delith
• The magnitude of the OCP hysteresis would be tracked with the variable H = lith_ref - delith_ref, where lith_ref = ocp_lith(sto, T) and delith_ref = ocp_delith(sto, T)
• The resulting OCP is calculated as ocp_surf = delith_ref + H * h

The method would be implemented in the new file pybamm\models\submodels\interface\open_circuit_potential\axen_ocp.py

And the following files would be modified:

• pybamm\models\submodels\interface\open_circuit_potential\__init__.py: to make the axen_ocp module available in the namespace and import the class AxenOpenCircuitPotential
• pybamm\models\full_battery_models\base_battery_model.py and pybamm\models\full_battery_models\lithium_ion\base_lithium_ion_model.py: to included the model's ocp option "Axen"
• pybamm\parameters\lithium_ion_parameters.py: to add the new function U_hysteresis_branch, which includes the effect of temperature in the hysteresis branches, and to define the parameters K_lith and K_delith
• pybamm\CITATIONS.bib: to add the relevant citation

I could contribute with the implementation if you think this is worth including in pybamm.

Additional context

No response

[rtimms]

This would be a very welcome addition, please let me know if you need any help with this!

Also see #4332 for some discussion around refactoring the current implementation of the Wycisk model.

[chtamar]

That's great @rtimms . I'd be glad to contribute to the project.
I'll test the implementation on the model in #4332, see how it goes and then get the pull request ready.

[chtamar]

Hi @rtimms, very interesting discussion in #4332.

I'm not sure if this is still relevant but the problem in #4332 seems to be related to the choice of the solver tolerance. The figure below shows the hysteresis state for the Wycisk approach using the Casadi solver with two different tolerances and initial state of h_init = 0 (pybamm v25.1.1). For a tolerance of 1e-3 (used in the model of #4332), the hysteresis state h can take values larger than 1 (the bug reported in #4332). But when the tolerance is reduced to 1e-6, this problem is not observed.

The Axen method would be a special case of the Wycisk implementation discussed in #4332 when $\gamma$ is a Scalar. However, in contrast to Wycisk, the Axen method skips the average OCP, only the hysteresis branches are required.

I have tested the implementation and it has been pushed to my fork. Should I create a pull request?

[rtimms]

Great, thanks for the update! Yes, please make a pull request.

[ejfdickinson]

@rtimms Let me know if you'd like to pick up on last week's in-person conversation on aligning all the hysteresis models.

I think it would probably make sense to package these under a single 'base hysteresis model' or else bring the base treatment of hysteresis within the base OCP model. That would make it much easier to fix #3867 in the general case. I think you or another of the core team would be the best implementer, but I'd be happy to contribute some suggestions on base mathematical formulation.

[rtimms]

Yes, sounds good. There is an open PR for a new model #4816, so let’s merge that and then consolidate everything.

[ejfdickinson]

@rtimms Yes, agreed that that's the right order.