Tip

An interactive online version of this notebook is available, which can be accessed via

---

Alternatively, you may download this notebook and run it offline.

Hysteresis State models

%pip install "pybamm[plot,cite]" -q    # install PyBaMM if it is not installed
import pybamm

[notice] A new release of pip is available: 25.1.1 -> 25.2
[notice] To update, run: pip install --upgrade pip
Note: you may need to restart the kernel to use updated packages.

Model Equations

Herein, we outline the equations for the “one-state hysteresis” model (Axen, 2022), the “one-state differential capacity hysteresis” model (Wycisk, 2022), and the “current Sigmoid” model (Ai et al., 2022). In this notebook we collectively paraphrase all the original literature models into a common notation and definition. For the original model expressions please refer to the appropriate references.

In all of the models, the open-circuit potential is given by

\[U = \frac{1+h}{2} U_{\text{delith}} + \frac{1-h}{2} U_{\text{lith}},\]

where \(h\) is a hysteresis state variable, \(U_{\text{delith}}\) is the delithiation open-circuit potential, and \(U_{\text{lith}}\) is the lithiation open-circuit potential. The hysteresis state variable \(h\) is defined between -1 and 1, where 1 corresponds to the delithiation branch and -1 corresponds to the lithiation branch.

One-state hysteresis (Axen)

In this model, the state variable \(h(z,t)\) is both stoichiometry and time-dependent, and is governed by the following ODE:

\[\frac{dh(z,t)}{dt} = \gamma \left( \frac{i_{\text{vol}}}{F \, c_{\text{max}} \, \epsilon}\right) \frac{\left(1 - \text{sgn}(i_{\text{vol}}) \, h(z,t)\right)}{2},\]

where \(i_{\text{vol}}\) is the volumetric interfacial current density, \(F\) is the Faraday constant, \(c_{\text{max}}\) is the maximum lithium concentration in the particles and \(\epsilon\) is the active material volume fraction. In this model, the rate parameter \(\gamma\) is given by

\[\gamma = \left( \frac{1-\text{sgn}(i_{\text{vol}})}{2} \gamma_{\text{lith}} + \frac{1+\text{sgn}(i_{\text{vol}})}{2} \gamma_{\text{delith}}\right),\]

which allows for a different decay rate during delithiation (\(\gamma_{\text{delith}}\)) and lithiation (\(\gamma_{\text{lith}}\)). Note that this implementation is a reformulation of the governing equations in the original paper. Also note that \(\gamma\) is dimensionless.

One-state differential capacity hysteresis (Wycisk)

As in the Axen model, the state variable \(h(z,t)\) is both stoichiometry and time-dependent, and is governed by the same ODE

\[\frac{dh(z,t)}{dt} = \gamma \left( \frac{i_{\text{vol}}}{F \, c_{\text{max}} \, \epsilon}\right) \frac{\left(1 - \text{sgn}(i_{\text{vol}}) \, h(z,t)\right)}{2},\]

where now \(\gamma(z)\) is given by

\[\gamma(z) = \Gamma(z) \cdot \left(\frac{C_{\text{diff}}\left(z\right)}{C_{\text{ref,vol}}}\right)^{-x}.\]

Here \(C_{\text{diff}}(z)\) is the differential capacity with respect to potential, expressed as

\[C_{\text{diff}}(z) = \frac{dQ}{dU_{\text{eq}}(z)},\]

where \(Q\) is the capacity of the phase of active material experiencing the voltage hysteresis and \(U_{\text{eq}}\) is the average open-circuit potential of the phase of active material experiencing the voltage hysteresis. The remaining parameters are \(\Gamma\) and \(x\) which are both fitting parameters that affect the response of the hysteresis state decay when passing charge in either direction. The parameter \(C_{\text{ref,vol}}\) is the reference differential capacity, which in this implementation is hard-coded to 1 \(\text{F}\).

Current Sigmoid (Ai)

This approach uses a sigmoid function to switch between the delithiation and lithiation branches of the open-circuit potential, depending on the sign of the current. The original paper gives the following expression for the open-circuit potential:

\[U = \text{sigmoid}\left(-\frac{KI(t)}{Q_{\text{cell}}}\right) U_{\text{delith}} + \text{sigmoid}\left(\frac{KI(t)}{Q_{\text{cell}}}\right) U_{\text{lith}}\]

Where \(K\) is a fitting parameter (in the current implementation this is hard-coded to 100, as in the original publication), \(I\) is the cell current, and \(Q_{\text{cell}}\) is the cell capacity. To simplify the comparison with the other models, we can rewrite this expression in terms of the hysteresis state variable \(h\) given by

\[h(t) = 1 - 2 \, \text{sigmoid}\left(-\frac{KI(t)}{Q_{\text{cell}}}\right).\]

Comparing the different hysteresis state models

We now compare the behavior of the different hysteresis state models. First we define the models, which are all based on the DFN model with a composite negative electrode comprised of two particles phases. One phase is modelled using a single open-circuit potential, while the other is modelled using a hysteresis state model. The positive electrode is modelled using a single open-circuit potential.

models = {
    "current sigmoid": pybamm.lithium_ion.DFN(
        {
            "open-circuit potential": (("single", "current sigmoid"), "single"),
            "particle phases": ("2", "1"),
            "thermal": "lumped",
        }
    ),
    "one-state differential capacity hysteresis": pybamm.lithium_ion.DFN(
        {
            "open-circuit potential": (
                ("single", "one-state differential capacity hysteresis"),
                "single",
            ),
            "particle phases": ("2", "1"),
            "thermal": "lumped",
        }
    ),
    "one-state hysteresis": pybamm.lithium_ion.DFN(
        {
            "open-circuit potential": (("single", "one-state hysteresis"), "single"),
            "particle phases": ("2", "1"),
            "thermal": "lumped",
        }
    ),
}

Next we define the parameter set, which is for a graphite and silicon negative electrode with a NMC positive electrode. We add some thermal parameters.

parameters = pybamm.ParameterValues("Chen2020_composite")
parameters.update(
    {
        "Negative current collector density [kg.m-3]": 8933.0,
        "Negative current collector specific heat capacity [J.kg-1.K-1]": 385.0,
        "Negative current collector thermal conductivity [W.m-1.K-1]": 398.0,
        "Negative electrode density [kg.m-3]": 1555.0,
        "Negative electrode specific heat capacity [J.kg-1.K-1]": 1437.0,
        "Negative electrode thermal conductivity [W.m-1.K-1]": 1.58,
        "Separator density [kg.m-3]": 1017.0,
        "Separator specific heat capacity [J.kg-1.K-1]": 1978.0,
        "Separator thermal conductivity [W.m-1.K-1]": 0.34,
        "Positive electrode density [kg.m-3]": 3699.0,
        "Positive electrode specific heat capacity [J.kg-1.K-1]": 1270.0,
        "Positive electrode thermal conductivity [W.m-1.K-1]": 1.04,
        "Positive current collector density [kg.m-3]": 2702.0,
        "Positive current collector specific heat capacity [J.kg-1.K-1]": 903.0,
        "Positive current collector thermal conductivity [W.m-1.K-1]": 238.0,
        "Total heat transfer coefficient [W.m-2.K-1]": 5.0,
    }
)

Next, we need to add the additional parameters required by the model. In this example, \(K_{\text{lith}} = K_{\text{delith}} = 100\) for the HS model, but it is recommended to test different values for the decay rate aiming for good agreement to experiment if the information is available, or a smooth transition between the lithiation and delithiation branches of the hysteresis. For example, depending on the case study, an acceptable value for \(K_{\text{lith}}\) and \(K_{\text{delith}}\) could be \(10\) or lower.

# Wycisk parameters
parameters_wycisk = parameters.copy()
parameters_wycisk.update(
    {
        "Secondary: Negative particle hysteresis decay rate": 0.1,
        "Secondary: Negative particle hysteresis switching factor": 10,
        "Secondary: Initial hysteresis state in negative electrode": 1,
    }
)

# Axen parameters
parameters_axen = parameters.copy()
parameters_axen.update(
    {
        "Secondary: Negative particle lithiation hysteresis decay rate": 100,
        "Secondary: Negative particle delithiation hysteresis decay rate": 100,
        "Secondary: Initial hysteresis state in negative electrode": 1,
    }
)

parameter_values = {
    "current sigmoid": parameters,
    "one-state differential capacity hysteresis": parameters_wycisk,
    "one-state hysteresis": parameters_axen,
}

Then we define the experiment, which is a discharge and charge cycle at 1C.

experiment = pybamm.Experiment(
    [
        (
            "Discharge at 1C for 1 hour or until 2.5 V",
            "Rest for 15 minutes",
            "Charge at 1C until 4.2 V",
            "Hold at 4.2 V until 0.05 C",
            "Rest for 15 minutes",
        ),
    ]
)

We then loop over the models and solve them.

solutions = {}
for name, model in models.items():
    sim = pybamm.Simulation(
        model, experiment=experiment, parameter_values=parameter_values[name]
    )
    solutions[name] = sim.solve(calc_esoh=False)

Finally, we plot the results.

output_variables = [
    "Current [A]",
    "Negative electrode secondary stoichiometry",
    "Voltage [V]",
    "Volume-averaged cell temperature [K]",
    "Total heating [W.m-3]",
    "Hysteresis electrochemical heating [W.m-3]",
    "X-averaged negative electrode secondary hysteresis state",
    "X-averaged negative electrode secondary open-circuit potential [V]",
]

pybamm.dynamic_plot(
    list(solutions.values()),
    labels=list(solutions.keys()),
    colors=["black", "blue", "red"],
    linestyles=["-", "--", "-."],
    output_variables=output_variables,
)

<pybamm.plotting.quick_plot.QuickPlot at 0x32a3cf2f0>

We see that both one-state hysteresis models include a decay time for the hysteresis state variable \(h\), and that they “stick” the (de)lithiation branches during rest. This is in contrast to the current sigmoid model, which does not stick to the (de)lithiation branches during rest; it instantly switches back the the equilibrium open-circuit potential. The effect of this can clearly be seen in the voltage and open-circuit potential plots.