Solver error when increasing degradation of graphite parameters

alex.0249 · 25 July 2025 09:54

Hi,
I keep getting the solver error:
2025-07-25 10:44:02.905 - [ERROR] callbacks.on_experiment_error(233): Simulation error: FAILURE IDA_CONV_FAIL: Convergence test failures occurred too many times during one internal time step or minimum step size was reached.

[IDAS ERROR] IDASolve
At t = 254878 and h = 1e-06, the corrector convergence failed repeatedly or with |h| = hmin. This occurs after a certain number of cycles. I am trying to emulate the parameters from this paper: Lithium-ion battery degradation: how to model it - Physical Chemistry Chemical Physics (RSC Publishing)

import pybamm 
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pouch import get_parameter_values_pouch
import matplotlib.ticker as mticker
from scipy.signal import find_peaks
# Load the CSV file
graphite_deg = np.loadtxt("graphite_deg.csv", delimiter=",", skiprows=1)
cycle = graphite_deg[:, 0]  # Cycle number
E_dis = graphite_deg[:, 1]  # Discharged energy
cycle = cycle
E_dis = E_dis # Discharged energy

model1 = pybamm.lithium_ion.SPMe(
    {

        "SEI": "solvent-diffusion limited",
        "SEI porosity change": "true",
        "lithium plating": "partially reversible",
        "lithium plating porosity change": "true",  # alias for "SEI porosity change"
        "particle mechanics": "swelling and cracking",
        "SEI on cracks": "true",
        "loss of active material": "stress-driven",
        "thermal": "lumped",
    }
)



param = pybamm.ParameterValues("OKane2022")

cycle_number = 350
def graphite_cracking_rate_Ai2020(T_dim):

    k_cr = 3.9e-20 * 25
    Eac_cr = 0  # to be implemented
    arrhenius = np.exp(Eac_cr / pybamm.constants.R * (1 / T_dim - 1 / 298.15))
    return k_cr * arrhenius


param.update({
    "SEI solvent diffusivity [m2.s-1]":5e-22,
    "Dead lithium decay constant [s-1]": 1e-5,
    "Negative electrode cracking rate": graphite_cracking_rate_Ai2020,
    "Negative electrode LAM constant proportional term [s-1]": 2.5e-06 ,
    "Ambient temperature [K]": 273.15+37,  # accelerate degradation temperature
})
Q_nom = param["Nominal cell capacity [A.h]"]

param.update({"Total heat transfer coefficient [W.m-2.K-1]" : 150, })

param.set_initial_stoichiometries(1)
exp = pybamm.Experiment(
    # [
    #     "Hold at 4.2 V until C/100",
    #     "Rest for 4 hours",
    #     "Discharge at 0.1C until 2.5 V",  # initial capacity check
    #     "Charge at 0.3C until 4.2 V",
    #     "Hold at 4.2 V until C/100",
    # ]
    [
        (
            "Discharge at 1.5C until 2.5 V",  # ageing cycles
            # "Charge at 3C until 4.2 V",
            "Charge at 1.5C until 4.2 V",
            "Hold at 4.2 V until C/3",
        )
    ]
    * cycle_number
    # + ["Discharge at 0.1C until 2.5 V"],  # final capacity check
)

solver = pybamm.IDAKLUSolver(
    rtol=3e-3,
    atol=3e-3,
    options={
        "dt_max": 10,
        "dt_min": 1e-6,
        "max_num_steps": 5000,
        "print_stats": False,
    },
)

var_pts = {"x_n": 6, "x_s": 6, "x_p": 6, "r_n": 6, "r_p": 6}



sim1 = pybamm.Simulation(
    model1, parameter_values=param, experiment=exp, solver=solver, var_pts=var_pts
)

If anyone has any ideas about why this happens, please let me know!
Thanks

Marc · 25 July 2025 15:08

Hi @alex.0249, are you on the latest version of pybamm? The example you posted runs without error for me, using both your solver settings and the default settings.

alex.0249 · 25 July 2025 15:55

I’ve reinstalled Pybamm but get the same error just before 500 cycle. All my solutions spike like this :

. I’ve tried increasing the mesh and changing solvers too.

Marc · 28 July 2025 12:29

If you are seeing the same result with all solvers, I’d suggest looking at the model to see if something unphysical is occurring at high degradation. Some common places of model failure are the active material fractions or porosities going below 0, for example.

alex.0249 · 28 July 2025 18:21

The porosity does indeed spike to 0 suddenly at this point, but I am unsure this is the root cause of the issue. I have followed the process from the notebook: Modelling coupled degradation mechanisms in PyBaMM — PyBaMM v25.6.0 Manual and changed to an SPMe model and discharge/charging at 1.5C. When cycling to 2000 cycles I get the same error: 2025-07-28 19:15:59.030 - [ERROR] callbacks.on_experiment_error(233): Simulation error: FAILURE IDA_ERR_FAIL: Error test failures occurred too many times during one internal time step or minimum step size was reached. The porosity and LLI spike just before the failure but how can I be sure this is the root cause of the issue?

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


model = pybamm.lithium_ion.SPMe(
    {
        "SEI": "solvent-diffusion limited",
        "SEI porosity change": "true",
        "lithium plating": "partially reversible",
        "lithium plating porosity change": "true",  # alias for "SEI porosity change"
        "particle mechanics": ("swelling and cracking", "swelling only"),
        "SEI on cracks": "true",
        "loss of active material": "stress-driven",
    }
)
param = pybamm.ParameterValues("OKane2022")
var_pts = {
    "x_n": 5,  # negative electrode
    "x_s": 5,  # separator
    "x_p": 5,  # positive electrode
    "r_n": 30,  # negative particle
    "r_p": 30,  # positive particle
}

cycle_number = 2000
exp = pybamm.Experiment(
    [
        "Hold at 4.2 V until C/100",
        "Rest for 4 hours",
        "Discharge at 0.1C until 2.5 V",  # initial capacity check
        "Charge at 0.3C until 4.2 V",
        "Hold at 4.2 V until C/100",
    ]
    + [
        (
            "Discharge at 1.5C until 2.5 V",  # ageing cycles
            "Charge at 1.5C until 4.2 V",
            "Hold at 4.2 V until C/100",
        )
    ]
    * cycle_number
    + ["Discharge at 0.1C until 2.5 V"],  # final capacity check
)
solver = pybamm.IDAKLUSolver()
sim = pybamm.Simulation(
    model, parameter_values=param, experiment=exp, solver=solver, var_pts=var_pts
)
sol = sim.solve()

Qt = sol["Throughput capacity [A.h]"].entries
Q_SEI = sol["Loss of capacity to negative SEI [A.h]"].entries
Q_SEI_cr = sol["Loss of capacity to negative SEI on cracks [A.h]"].entries
Q_plating = sol["Loss of capacity to negative lithium plating [A.h]"].entries
Q_side = sol["Total capacity lost to side reactions [A.h]"].entries
Q_LLI = (
    sol["Total lithium lost [mol]"].entries * 96485.3 / 3600
)  # convert from mol to A.h
plt.figure()
plt.plot(Qt, Q_SEI, label="SEI", linestyle="dashed")
plt.plot(Qt, Q_SEI_cr, label="SEI on cracks", linestyle="dashdot")
plt.plot(Qt, Q_plating, label="Li plating", linestyle="dotted")
plt.plot(Qt, Q_side, label="All side reactions", linestyle=(0, (6, 1)))
plt.plot(Qt, Q_LLI, label="All LLI")
plt.xlabel("Throughput capacity [A.h]")
plt.ylabel("Capacity loss [A.h]")
plt.legend()
plt.show()

Qt = sol["Throughput capacity [A.h]"].entries
LLI = sol["Loss of lithium inventory [%]"].entries
LAM_neg = sol["Loss of active material in negative electrode [%]"].entries
LAM_pos = sol["Loss of active material in positive electrode [%]"].entries
plt.figure()
plt.plot(Qt, LLI, label="LLI")
plt.plot(Qt, LAM_neg, label="LAM (negative)")
plt.plot(Qt, LAM_pos, label="LAM (positive)")
plt.xlabel("Throughput capacity [A.h]")
plt.ylabel("Degradation modes [%]")
plt.legend()
plt.show()

eps_neg_avg = sol["X-averaged negative electrode porosity"].entries
eps_neg_sep = sol["Negative electrode porosity"].entries[-1, :]
eps_neg_CC = sol["Negative electrode porosity"].entries[0, :]
plt.figure()
plt.plot(Qt, eps_neg_avg, label="Average")
plt.plot(Qt, eps_neg_sep, label="Separator", linestyle="dotted")
plt.plot(Qt, eps_neg_CC, label="Current collector", linestyle="dashed")
plt.xlabel("Throughput capacity [A.h]")
plt.ylabel("Negative electrode porosity")
plt.legend()
plt.show()

alex.0249 · 28 July 2025 18:22

Marc · 28 July 2025 21:07

The porosity approaching zero is a common root cause of degradation model failure. The model breaks at zero porosity because of the electrolyte concentration PDE, which has a factor of 1/porosity on the right-hand side. The solver fails because there is no solution for porosity = 0. Physically, porosity < 0 is not meaningful.

I’d recommend trying other degradation mechanisms or adjusting the degradation prefactors to match experimental measurements.

alex.0249 · 29 July 2025 11:13

Thank you for your responses! I have turned off the lithium plating & SEI porosity change, and kept the temperature constant but I still see a large spike in the LLi and side reactions with the base parameters of the OKane set:

I am trying to fit the degradation parameters so the capacity reaches 80% after 1500-2000 cycles from experimental data but I just can’t work out how to do this and avoid the solver error/ reach physical limits. I have found if I remove the SEI cracking then it runs. What parameters could I change for the model to still run with this submodel? I have tried reducing “Initial SEI on cracks thickness” but this still caused a failure

Topic		Replies	Views
Cycle Ageing with PyBaMM EIS PyBaMM Q&A	1	268	26 January 2025
PyBaMM EIS with coupled degradation models PyBaMM EIS	2	268	13 February 2025
Degradation of Lg 18650 MH1 cell General	2	189	14 January 2025
Using the Idaklu Solver causing a memory leak Solvers	8	184	28 January 2025
Error plotting the different contributions to LLI against throughput capacity in coupled degradation model Models	1	116	9 January 2025

Solver error when increasing degradation of graphite parameters

Related topics