Are there any problems with the OKane2022 parameter set, or is there something wrong with my simulation experiment settings?

[RobinBluehub]

When I used the parameter set of OKane2022 to reproduce the Craking and SEI experiments in Dr. OKane's paper: Lithium-ion battery degradation: how to model it. If the Negative/Positive electrode cracking rate is increased, for example, 30✖3.9e-20, or greater than 50✖3.9e-20, the cycle will be terminated early because the positive electrode active particle crack is larger than the particle radius, and the cycle can only reach a few hundred times. I tried to find the reason, but found nothing.
Hope you can give some advice, best wishes.

Code show as below：
'import pybamm
import numpy as np
import matplotlib.pyplot as plt
pybamm.set_logging_level("NOTICE")
from pybamm import Parameter, constants, exp

#Used to calculate Total charge throughput (A.h)
def cal_integral(x, y):
import scipy
from scipy.integrate import simps
integrals = []
for i in range(len(y)):
integrals.append(scipy.integrate.trapz(y[:i + 1], x[:i + 1]))
return integrals

option={
"SEI":"solvent-diffusion limited",
"particle mechanics": "swelling and cracking",
"SEI on cracks": "true",
}
model=pybamm.lithium_ion.DFN(options=option)
param1 = pybamm.ParameterValues("OKane2022")

def graphite_cracking_rate_Ai2020(T_dim):
k_cr = 50*3.9e-20
Eac_cr = 0 # to be implemented
arrhenius = exp(Eac_cr / constants.R * (1 / T_dim - 1 / 298.15))
return k_cr * arrhenius

def cracking_rate_Ai2020(T_dim):
k_cr = 50*3.9e-20
Eac_cr = 0 # to be implemented
arrhenius = exp(Eac_cr / constants.R * (1 / T_dim - 1 / 298.15))
return k_cr * arrhenius

param1.update({"Negative electrode cracking rate":graphite_cracking_rate_Ai2020})
param1.update({"Positive electrode cracking rate":cracking_rate_Ai2020})
param1.update({"Outer SEI solvent diffusivity [m2.s-1]":2.5E-22})
param1.update({"Initial outer SEI thickness [m]":5E-9})

experiment_ini = pybamm.Experiment(
[(
"Hold at 4.2V until 50mA ",
"Rest for 2 hours ",
"Discharge at 0.5 A until 2.5V ",
"Charge at 1.5 A until 4.2V ",
"Hold at 4.2V until 50mA ",

)],period="30 seconds"

)

cycle =1000
experiment = pybamm.Experiment(
[(
"Discharge at 5 A until 2.5V",
"Charge at 1.5 A until 4.2V",
"Hold at 4.2V until 50mA"
)] * cycle
)

solver=pybamm.CasadiSolver("safe",dt_max=200,max_step_decrease_count=20)
var_pts={"x_n": 40, "x_s":20, "x_p": 40, "r_n": 40, "r_p":40}
sim1_ini = pybamm.Simulation(model, experiment = experiment_ini,parameter_values = param1,solver = solver,var_pts=var_pts)
solution1_ini = sim1_ini.solve(calc_esoh=False)
N=np.arange(1,1001,50)
sim1=pybamm.Simulation(model, experiment = experiment,parameter_values = param1,solver = solver,var_pts=var_pts)
solution1=sim1.solve(starting_solution=solution1_ini,calc_esoh=False,save_at_cycles=N)
#solution1.save("crack_1.spk")
x1=[]
y1=[]
zong_dis_cap_1=0
for i in range(1,cycle+1):
current1 = abs(solution1.cycles[i].steps[0]["Current [A]"].entries)
time1 = solution1.cycles[i].steps[0]["Time [h]"].entries
biaoding_discharge_capacity_1=cal_integral(time1,current1)[-1]

current_zong_1 = abs(solution1.cycles[i]["Current [A]"].entries)
time_zong_1= solution1.cycles[i]["Time [h]"].entries
zong_dis_cap_dan_1=cal_integral(time_zong_1,current_zong_1)[-1]
zong_dis_cap_1=zong_dis_cap_1+zong_dis_cap_dan_1

x1.append(zong_dis_cap_1)
y1.append(biaoding_discharge_capacity_1)

print(x1)
print(y1)'

The warning is as follows:
2022-08-18 09:51:18.148 - [NOTICE] callbacks.on_cycle_start(174): Cycle 181/1001 (29 minutes, 32 seconds elapsed) -------------------- 2022-08-18 09:51:18.149 - [NOTICE] callbacks.on_step_start(182): Cycle 181/1001, step 1/3: Discharge at 5 A until 2.5V 2022-08-18 09:51:21.575 - [NOTICE] callbacks.on_step_start(182): Cycle 181/1001, step 2/3: Charge at 1.5 A until 4.2V 2022-08-18 09:51:26.532 - [WARNING] callbacks.on_experiment_infeasible(233): Experiment is infeasible: 'event: positive particle crack length larger than particle radius' was triggered during 'Charge at 1.5 A until 4.2V'. The returned solution only contains up to step 2 of cycle 181. 2022-08-18 09:51:26.664 - [NOTICE] callbacks.on_experiment_end(222): Finish experiment simulation, took 29 minutes, 40 seconds

[brosaplanella]

@DrSOKane, what are your thoughts on this?

[DrSOKane]

Hi Robin,

I recommend disabling cracking in the positive electrode, because there aren't actually any degradation mechanisms associated wit it so you're just wasting computational power. You can do this by setting your options to be

option={
"SEI":"solvent-diffusion limited",
"particle mechanics": ("swelling and cracking", "none"),
"SEI on cracks": "true",
}

If you're using the stress-driven loss of active material model, you should set

option={
"SEI":"solvent-diffusion limited",
"particle mechanics": ("swelling and cracking", "swelling only"),
"SEI on cracks": "true",
"loss of active material": "stress-driven",
}

instead. In pull request #2262 I've modified the code so it does both of these automatically and you don't need to set the "particle mechanics" option at all (but you still can if you wish).

I also strongly recommend setting "SEI porosity change" to "true" if you're modelling SEI on cracks.

Sorry about the slow response: I had found a different bug with the SEI on cracks, where the loss of lithium inventory variable was being incremented but the lithium wasn't actually going anywhere. A partial fix for this has already been implemented: the full fix is in pull request #2262.

[RobinBluehub]

Hi, DrSOKane. Thank you very much for your reply, your answer is very helpful to me. In addition, I would like to ask why there is no ”Knees“ in the capacity drop caused by li plating in your paper. Is it because the two parameters are set relatively small, or because the cycle period is relatively short, the ”Knees“ does not appear?

Hope to hear from you, and I wish you all the best in your research work!
Robin

[DrSOKane]

Hi Robin, are you referring to the Yang et al. (2017) paper?

Yang et al.'s simulations did indeed predict that lithium plating occurs after the knee point, but did not predict that lithium plating caused the knee point. What they predict is that pore clogging caused by SEI growth causes the knee point, after which rapid plating occurs, accounting for the sudden drop in capacity.

As you know, one of my simulations predicts a knee point caused by rapid particle cracking leading to pore clogging. However, lithium plating was not enabled for that simulation. None of the simulations where plating was enabled had a large enough cracking rate to cause a knee point.

Many thanks, Simon

[RobinBluehub]

Yes, I refer to the Yang et al. (2017) paper.