Coupling of reaction limited model

General
1

Hi,

I am trying to couple the reaction limited SEI growth submodel with LAM and cracking for graphite and another parameter set to fit the degradation using the SEI reaction exchange current. For the OKane paramete set I get an error when combining with other submodels:

2025-08-08 11:15:39.885 - [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.

import pybamm
import numpy as np
import matplotlib.pyplot as plt
# 1) Base model, geometry & experiment ---------------------------------------
import matplotlib.ticker as mticker
model = pybamm.lithium_ion.SPMe(
    {
        "SEI": "reaction limited",
        "SEI porosity change": "true",
        "particle mechanics": "swelling and cracking",
        "SEI on cracks": "true",
        "loss of active material": "stress-driven",
        "thermal": "lumped",
    }
)

cycle_number = 150
var_pts = {"x_n": 10, "x_s": 10, "x_p": 10, "r_n": 10, "r_p": 10}

exp = pybamm.Experiment(
    # [
    #     "Hold at 3 V until C/100",
    #     "Rest for 4 hours",
    #     "Discharge at 0.001C until 1.5 V",  # initial capacity check
    #     "Charge at 0.3C until 3 V",
    #     "Hold at 3 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 0.3C",
        )
    ]
    * cycle_number
    # + ["Discharge at 0.1C until 1.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": 10000},
# )

cycles = np.arange(1, cycle_number + 1)


# 2) Base‐update 
base_updates = {
    "Ambient temperature [K]": 273.15 + 45,
}


# 3) Helper to run a single sim and return capacity array -------------------

def run_capacity(param_override):
    pv = pybamm.ParameterValues("OKane2022")
    pv.update(base_updates)
    pv.update(param_override)
    pv.set_initial_stoichiometries(0.99)
    sim = pybamm.Simulation(
        model,
        parameter_values=pv,
        experiment=exp,
        solver=solver,
        var_pts=var_pts,
    )
    sol = sim.solve()

    cap = sol.summary_variables["Capacity [A.h]"]
    cap = np.array(cap)
    return cap


# 5) Sweep 2: SEI j_sei---------------------------------------

j_sei = [1.5e-7]
fig, ax = plt.subplots(figsize=(10, 6), dpi=500)
for J in j_sei:
    cap = run_capacity({ "SEI reaction exchange current density [A.m-2]": J })
    plt.plot(cycles, cap, label=f"j_SEI = {J:.1e}")

plt.xlabel("Cycle")
plt.ylabel("Capacity [A·h]")
plt.title("Effect of SEI exchange current density on Capacity Fade")
plt.legend()
plt.tight_layout()
ax.yaxis.set_major_formatter(mticker.FormatStrFormatter("%.4f"))
ax.set_box_aspect(1)
plt.show()

For my parameter set with a Wadsley roth anode I see very odd behaviour where j_SEI makes no difference to the capacity fade and the Loss of capacity to negative SEI [A.h] becomes negative. I was wondering if anyone had any ideas how this might occur and what parameter might be at fault?

image
image578×448 16 KB

Any help will be much appreciated :slight_smile: