Issues with Aging Simulation of LGM50 Cell: Current Behavior, Voltage Rise, and Feature Extraction

Hi,

I am simulating an aging experiment for an LGM50 cell using PyBaMM, involving constant current discharge and CCCV charging cycles. My goal is to analyze the degradation behavior over 500 cycles, including feature extraction during CCCV cycles. However, I am encountering several issues that I would appreciate your insights on.

Simulation Setup
I am using the DFN model with side reaction options for SEI growth, lithium plating, particle mechanics, and thermal effects. The parameter set combines ORegan2022 and OKane2022, with additional updates for aging parameters. I am using IDAKLUSolver instead of CasadiSolver as it has been suggested for aging simulations. Here’s the initialization and model setup:

# ------------- Initialization --------------------------
parameter_values_ORegan = pybamm.ParameterValues("ORegan2022")
parameter_values_okane = pybamm.ParameterValues("OKane2022")
combined_dict = {**parameter_values_okane, **parameter_values_ORegan}
combined_parameters = pybamm.ParameterValues(combined_dict)
parameter_values_298 = combined_parameters

# Aging parameters (from https://arxiv.org/pdf/2311.05482)
parameter_values_298.update({"Ratio of lithium moles to SEI moles": 2})
parameter_values_298.update({"Negative electrode cracking rate": 5.29e-25})
parameter_values_298.update({"Positive electrode cracking rate": 5.29e-25})

#Considering irreversible Li-plating from 1e-10 -> 5e-11
parameter_values_298.update({"Lithium plating kinetic rate constant [m.s-1]": 5e-11})

# Model 
model = pybamm.lithium_ion.DFN(options={
    "SEI": "interstitial-diffusion limited",
    "SEI porosity change": "true",
    "lithium plating": "irreversible",
    "lithium plating porosity change": "true",
    "particle mechanics": ("swelling and cracking", "swelling only"),
    "SEI on cracks": "true",
    "loss of active material": "stress-driven",
    "thermal": "lumped",
})

# Spatial discretization
var_pts = {
    "x_n": 5,  # negative electrode
    "x_s": 5,  # separator
    "x_p": 30,  # positive electrode
    "r_n": 30,  # negative particle
    "r_p": 30,  # positive particle
}

# Experiment
cycle_number = 500
experiment_298 = pybamm.Experiment(
    [
        "Hold at 4.4 V until C/100 (5 minute period)",
        "Rest for 4 hours (5 minute period)",
        "Discharge at 0.1C until 2.5 V (5 minute period)",  # initial capacity check
        "Charge at 0.3C until 4.4 V (5 minute period)",
        "Hold at 4.4 V until C/100 (5 minute period)",
    ] + [
        (
            "Discharge at 1C until 2.5 V",  # ageing cycles 1C discharge
            "Charge at 0.3C until 4.4 V (5 minute period)",  # Using 4.4V since Upper voltage cut-off = 4.4V
            "Hold at 4.4 V until C/100 (5 minute period)",
        )
    ] * cycle_number + ["Discharge at 0.1C until 2.5 V (5 minute period)"],  # final capacity check
)

solver = pybamm.IDAKLUSolver()
sim_298 = pb.Simulation(model, parameter_values=parameter_values_298, experiment=experiment_298, solver=solver, var_pts=var_pts)

Here are some plots :

Early cycles:

Abrupt change of voltage / current profile:

Current increasing during cv phase :

Issue 1: Current Increasing During CV Phase
In a typical CCCV charging protocol, I expect the current to decrease during the constant voltage phase as the battery approaches full charge. However, in my simulation, in the last cycles, the current magnitude during the CV phase increases over time, as shown in the attached plots. This behavior is unexpected and seems unphysical. Could this be related to the model options (e.g., lithium plating or SEI growth), the solver I am using or the experiment definition (CV phase at 4.4V, which is specified as the upper voltage cut-off in O’Regan 2022.)?

Issue 2: Voltage Rise During Constant Current Discharge
During the constant current discharge at 1C, I observe an unusual voltage rise at the beginning of each discharge cycle (e.g., see the voltage subplot in the 4th figure, where the voltage briefly increases before dropping). This is not typical for a discharge process, where I would expect a monotonic voltage decrease. Is this behavior due to solver issues or side reactions considered?

Issue 3: CCCV Detection Algorithm and Feature Extraction
I have implemented a feature extraction algorithm to detect CCCV cycle boundaries, CC index at 3.8 V, CV start (at 4.4 V), and CV end. Below is the algorithm:

# Feature extraction algorithm
eps = 1e-2
v_low_cutoff = parameter_values_298["Lower voltage cut-off [V]"]
v_cc_start_counting = 3.8
voltage_cv_threshold = 4.4 - 0.01
c_100_threshold = parameter_values_298["Nominal cell capacity [A.h]"] / 100
slope_window = 5

cccv_starts = []
cc_starts38 = []
cv_starts = []
cv_ends = []
cc_durations = []
cv_durations = []
cc_end_slopes = []
soh_at_cycle_ends = []

for i in range(1, len(voltage)):
    if voltage[i - 1] < v_low_cutoff + eps and voltage[i] > v_low_cutoff + eps and current[i] < 0:
        cccv_starts.append(i)

for start_idx in cccv_starts:
    cc_start38_idx = None
    for i in range(start_idx, len(voltage)):
        if voltage[i] >= v_cc_start_counting and current[i] < 0:
            # Check if current remains negative until CV voltage is reached
            valid_cc_start = True
            for j in range(i, len(voltage)):
                if voltage[j] >= voltage_cv_threshold - eps:
                    break  # Stop checking once CV voltage is reached
                if current[j] >= 0:  # Current becomes non-negative before CV
                    valid_cc_start = False
                    break
            if valid_cc_start:
                cc_start38_idx = i
                break

    cv_start_idx = None
    for i in range(start_idx, len(voltage) - 1):
        if (voltage[i] >= voltage_cv_threshold - eps and
                abs(voltage[i + 1] - voltage[i]) < 0.001 and
                current[i] < 0):
            cv_start_idx = i
            break

    cv_end_idx = None
    for i in range(cv_start_idx or start_idx, len(voltage) - 1):
        if (abs(current[i]) <= c_100_threshold or
                (current[i] > 0 and voltage[i] < 4.0)):
            cv_end_idx = i
            break

    if (cc_start38_idx and cv_start_idx and cv_end_idx and
            cc_start38_idx < cv_start_idx < cv_end_idx):
        cc_starts38.append(cc_start38_idx)
        cv_starts.append(cv_start_idx)
        cv_ends.append(cv_end_idx)
        cc_duration = time[cv_start_idx] - time[cc_start38_idx]
        cv_duration = time[cv_end_idx] - time[cv_start_idx]
        cc_durations.append(cc_duration)
        cv_durations.append(cv_duration)

        window_start = max(cc_start38_idx, cv_start_idx - slope_window)
        window_end = cv_start_idx + 1
        time_window = time[window_start:window_end]
        voltage_window = voltage[window_start:window_end]
        slope = np.gradient(voltage_window, time_window)[-1]
        cc_end_slopes.append(slope)

I’ve noticed that the algorithm detects the CC start at a voltage of approximately 4.1 V instead of the intended 3.8 V in the last cycles (see the third figure, where the green dotted lines for “CC Start (3.8 V)” appear closer to 4.1 V). I suspect this might be related to the solver’s sampling time being too coarse, causing the algorithm to miss the exact point where the voltage crosses 3.8 V. Could this be the case? If so, how can I adjust the solver’s time step or sampling frequency to improve the resolution? Alternatively, is there a better way to detect these indices more accurately?

I would be very grateful for any kind of help related to any of these problems.