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.