Inability to perform coupled degradation analysis on non OKane2022 parameter sets

[MirAbbasAli2A]

Hello Pybamm Team -

I have been struggling with trying to perform a coupled degradation analysis on parameter sets that are not OKane2022. My motivation for using different parameter sets is to analyze the degradation performance on battery cells that do not have Si-Ox in the anode. I have tried performing the analysis on the Ai2022, Ecker2015 and Mohtat2020 parameter sets all with different types of error messages. With the Ai2022 and Ecker 2015 parameter sets, I kept getting a solver error and/or Q_p doesn't match Q_Li warning message. With the Mohtat2020 dataset, I get the same warning message and my porosity in the negative electrode goes to 0 very quickly.

Here is my code for the Mohtat2020 parameter set:

import pybamm
import matplotlib.pyplot as plt
import numpy as np
from os import path
import pandas as pd
from tec_reduced_model.set_parameters import (
    set_thermal_parameters,
    set_experiment_parameters,
    set_ambient_temperature,
    )
# from tec_reduced_model.process_experimental_data import import_thermal_data, get_idxs
from pybamm import Parameter, constants, exp
# temperature = 25
# Crate = 1
# cells_ignore = ["791"]
#Define the DFN Model
model = pybamm.lithium_ion.DFN(options = {
    "thermal":"x-full",
    "cell geometry": "pouch", 
    "SEI": "solvent-diffusion limited",
    "SEI porosity change": "true",
    "particle mechanics": ("swelling and cracking", "swelling only"),
    "SEI on cracks": "true",
    "loss of active material": "stress-driven",
    "calculate discharge energy": "true",
    "lithium plating": "partially reversible",
    "lithium plating porosity change": "true"    
    },
    name="TDFN",
    )
#Parameter assignment and definition
param = pybamm.ParameterValues("Mohtat2020")
# param = set_thermal_parameters(param, 16, 2.32e6, temperature)
# print(param.search("cut"))
# print(model.variables.search("maximum"))
T = model.variables["Volume-averaged cell temperature [K]"]
c_e = model.variables["Electrolyte concentration [mol.m-3]"]
c_Li = model.variables["Lithium plating concentration [mol.m-3]"]
L_sei = model.variables["Total SEI thickness [m]"]
sto_neg= model.variables["Negative electrode stoichiometry"]
sto_pos = model.variables["Positive electrode stoichiometry"]
cs_max_neg = model.variables["Maximum negative particle concentration"]
cs_max_pos = model.variables["Maximum positive particle concentration [mol.m-3]"]
 
param["Ambient temperature [K]"] = 298.15
param["Initial temperature [K]"] = 298.15
param["Inner SEI reaction proportion"] = 0.0
param["SEI growth activation energy [J.mol-1]"] = 38000.0
param["Ratio of lithium moles to SEI moles"]=1.0
param["Initial inner SEI thickness [m]"]=0.0
param["Initial outer SEI thickness [m]"] = 5e-9
param["SEI growth activation energy [J.mol-1]"]=38000.0

# param["Negative electrode LAM constant proportional term [s-1]"]=2.7778e-07
# param["Positive electrode LAM constant proportional term [s-1]"]=2.7778e-07
param["Total heat transfer coefficient [W.m-2.K-1]"]=10.0
# param.update({"Typical plated lithium concentration [mol.m-3]": 1000.0}, check_already_exists=False)
# param.update({"Lithium metal partial molar volume [m3.mol-1]": 1.3e-05}, check_already_exists=False)
param.update({"Lithium plating kinetic rate constant [m.s-1]": 1e-09}, check_already_exists=False)
# param.update({"Initial plated lithium concentration [mol.m-3]": 0.0}, check_already_exists=False)
param.update({"Lithium plating transfer coefficient": 0.65})
param.update({"Dead lithium decay constant [s-1]": 1e-06}, check_already_exists=False)
# param.update({"Negative electrode charge transfer coefficient": 0.5}, check_already_exists=False)
# param.update({"Negative electrode double-layer capacitance [F.m-2]": 0.2}, check_already_exists=False)
param.update({"Negative electrode Poisson's ratio": 0.3}, check_already_exists=False)
param.update({"Negative electrode Young's modulus [Pa]": 15000000000.0}, check_already_exists=False)
param.update({"Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0}, check_already_exists=False)
param.update({"Negative electrode partial molar volume [m3.mol-1]": 3.1e-06}, check_already_exists=False)
param.update({"Negative electrode initial crack length [m]": 2e-08}, check_already_exists=False)
param.update({"Negative electrode initial crack width [m]": 1.5e-08}, check_already_exists=False)
param.update({"Negative electrode number of cracks per unit area [m-2]": 3180000000000000.0}, check_already_exists=False)
param.update({"Negative electrode Paris' law constant b": 1.12}, check_already_exists=False)
param.update({"Negative electrode Paris' law constant m": 2.2}, check_already_exists=False)
param.update({"Negative electrode LAM constant exponential term": 2.0}, check_already_exists=False)
param.update({"Negative electrode critical stress [Pa]": 60000000.0}, check_already_exists=False)
param.update({"Negative electrode LAM constant proportional term [s-1]": 2.7778e-07}, check_already_exists=False)
param.update({"Positive electrode charge transfer coefficient": 0.5}, check_already_exists=False)
param.update({"Positive electrode double-layer capacitance [F.m-2]": 0.2}, check_already_exists=False)
param.update({"Positive electrode Poisson's ratio": 0.2}, check_already_exists=False)
param.update({"Positive electrode Young's modulus [Pa]": 375000000000.0}, check_already_exists=False)
param.update({"Positive electrode reference concentration for free of deformation [mol.m-3]": 0.0}, check_already_exists=False)
param.update({"Positive electrode partial molar volume [m3.mol-1]": 1.25e-05}, check_already_exists=False)
param.update({"Positive electrode initial crack length [m]": 2e-08}, check_already_exists=False)
param.update({"Positive electrode initial crack width [m]": 1.5e-08}, check_already_exists=False)
param.update({"Positive electrode number of cracks per unit area [m-2]": 3180000000000000.0}, check_already_exists=False)
param.update({"Positive electrode Paris' law constant b": 1.12}, check_already_exists=False)
param.update({"Positive electrode Paris' law constant m": 2.2}, check_already_exists=False)
param.update({"Positive electrode LAM constant exponential term": 2.0}, check_already_exists=False)
param.update({"Positive electrode critical stress [Pa]": 375000000.0}, check_already_exists=False)
param.update({"Positive electrode LAM constant proportional term [s-1]": 2.7778e-07}, check_already_exists=False)
param.update({"Cell thermal expansion coefficient [m.K-1]":1.1e-06}, check_already_exists=False)
 
def anode_cracking_rate_Ai2020(T):
    k_cr = 3.9e-20
    T_ref = pybamm.Parameter("Reference temperature [K]")
    Eac_cr = 81044.08
    arrhenius = exp(Eac_cr / constants.R * (1 / T - 1 / T_ref))
    return k_cr * arrhenius
def cathode_cracking_rate_Ai2020(T_dim):
    k_cr = 3.9e-20
    T_ref = pybamm.Parameter("Reference temperature [K]")
    Eac_cr = 81044.08
    arrhenius = exp(Eac_cr / constants.R * (1 / T - 1 / T_ref))
    return k_cr * arrhenius
def solvent_diffusivity_sei(T):
    Ea = 37000
    arrhenius = exp(Ea / constants.R * (1 / 298.15 - 1 / T))
    return 2.50000000000000002e-22 * arrhenius
def plating_exchange_current_density_OKane2020(c_e, c_Li, T):
    k_plating = param["Lithium plating kinetic rate constant [m.s-1]"]
    return pybamm.constants.F * k_plating * c_e
def stripping_exchange_current_density_OKane2020(c_e, c_Li, T):
    k_plating = param["Lithium plating kinetic rate constant [m.s-1]"]
    return pybamm.constants.F * k_plating * c_e
def SEI_limited_dead_lithium_OKane2022(L_sei):
    gamma_0=param["Dead lithium decay constant [s-1]"]
    L_inner_0 = pybamm.Parameter("Initial inner SEI thickness [m]")
    L_outer_0 = pybamm.Parameter("Initial outer SEI thickness [m]")
    L_sei_0 = L_inner_0 + L_outer_0
    gamma = gamma_0 * L_sei_0 / L_sei
    return gamma
def graphite_volume_change_Ai2020(sto_neg, cs_max_neg):
    p1 = 145.907
    p2 = -681.229
    p3 = 1334.442
    p4 = -1415.710
    p5 = 873.906
    p6 = -312.528
    p7 = 60.641
    p8 = -5.706
    p9 = 0.386
    p10 = -4.966e-05
    t_change = (
        p1 * sto_neg**9
        + p2 * sto_neg**8
        + p3 * sto_neg**7
        + p4 * sto_neg**6
        + p5 * sto_neg**5
        + p6 * sto_neg**4
        + p7 * sto_neg**3
        + p8 * sto_neg**2
        + p9 * sto_neg
        + p10
    )
    return t_change
 
def volume_change_Ai2020(sto_pos, cs_max_pos):
    omega = param["Positive electrode partial molar volume [m3.mol-1]"]
    t_change = omega * cs_max_pos * sto_pos
    return t_change

param.update({"Outer SEI solvent diffusivity [m2.s-1]": solvent_diffusivity_sei(T)}, check_already_exists=False)
param.update({"Negative electrode cracking rate":anode_cracking_rate_Ai2020(T)}, check_already_exists=False)
param.update({"Positive electrode cracking rate":cathode_cracking_rate_Ai2020}, check_already_exists=False)
param.update({"Exchange-current density for stripping [A.m-2]": stripping_exchange_current_density_OKane2020}, check_already_exists=False)
param.update({"Exchange-current density for plating [A.m-2]": plating_exchange_current_density_OKane2020})
param.update({"Dead lithium decay rate [s-1]": SEI_limited_dead_lithium_OKane2022}, check_already_exists=False)
param.update({"Negative electrode volume change": graphite_volume_change_Ai2020}, check_already_exists=False)
param.update({"Positive electrode volume change": volume_change_Ai2020}, check_already_exists=False)
#Define experiment
pybamm.set_logging_level("NOTICE")
cycle_number = 1
#Simple cycling experiment
N= 50
cccv_experiment = pybamm.Experiment([
    ("Charge at C/5 until 4.2V",
     "Hold at 4.2V until C/50",
    "Discharge at C/5 until 3.0V",
    )
    ]*N
    )
 

rpt_experiment = pybamm.Experiment([(
    "Rest for 3 hours",
    "Charge at C/20 until 4.2V",
    "Hold at 4.2V until C/50",
    "Rest for 1 hour",
    "Discharge at C/20 until 3.0V",
    "Hold at 3.0V until C/50",
    "Rest for 1 hour",
    )])
#Refine mesh
var_pts = {
    "x_n" : 40,  # x-direction, length, negative electrode
    "x_s" : 40,  # x-direction, length, separator
    "x_p" : 40,  # x-direction, length, positive electrode
    "r_n" : 40,  # number of volumes in the radial direction, negative particle
    "r_p" : 40   # number of volumes in the radial direction, positive particle
}
#Run model
solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-3)
cccv_sols = []
rpt_sols = []
cccv_sol = None
rpt_sol = None
M = 8
for i in range(M):
    sim = pybamm.Simulation(model, experiment=rpt_experiment,parameter_values=param, var_pts=var_pts, solver=solver)
    rpt_sol = sim.solve(starting_solution=cccv_sol)
    sim = pybamm.Simulation(model, experiment=cccv_experiment,parameter_values=param, var_pts=var_pts, solver=solver)
    cccv_sol = sim.solve(starting_solution=rpt_sol, save_at_cycles= 10)
    cccv_sols.append(cccv_sol)
    rpt_sols.append(rpt_sol)
cccv_cycles = []
cccv_capacities = []
rpt_cycles = []
rpt_capacities = []
capacities=[]
chargecapacitites=[]
LLI=[]
LAMne=[]
LAMpe=[]
for i in range(M):
    for j in range(N):
        # print(f"i: {i}, j: {j}, cycle: {i * (N + 1) + j + 2}")
        if cccv_sol.cycles[i * (N + 1) + j + 1] is not None:
            cccv_cycles.append(i * (N + 1) + j + 2)
            start_capacity = cccv_sol.cycles[i * (N + 1) + j + 1].steps[2]["Discharge capacity [A.h]"].entries[0]
            end_capacity = cccv_sol.cycles[i * (N + 1) + j + 1].steps[2]["Discharge capacity [A.h]"].entries[-1]
            cccv_capacities.append(end_capacity - start_capacity)
    rpt_cycles.append(i*(N+1)+1)
    start_capacity = rpt_sol.cycles[(i)*(N+1)].steps[4]["Discharge capacity [A.h]"].entries[0]
    end_capacity = rpt_sol.cycles[(i)*(N+1)].steps[4]["Discharge capacity [A.h]"].entries[-1]
    rpt_capacities.append(end_capacity-start_capacity)
    LLIval=rpt_sol.cycles[(i)*(N+1)]["Loss of lithium inventory [%]"].entries[-1]
    LAMneval=rpt_sol.cycles[(i)*(N+1)]["Loss of active material in negative electrode [%]"].entries[-1]
    LAMpeval=rpt_sol.cycles[(i)*(N+1)]["Loss of active material in positive electrode [%]"].entries[-1]
    LLI.append(LLIval)
    LAMne.append(LAMneval)
    LAMpe.append(LAMpeval)
#Plot Aging vs RPT Cycle Capacities
plt.scatter(cccv_cycles,cccv_capacities,label="Ageing cycles")
plt.scatter(rpt_cycles,rpt_capacities,label="RPT cycles")
plt.xlabel("Cycle number")
plt.ylabel("Discharge capacity [A.h]")
plt.legend()
plt.show()
#SOH Eval
SOH = []
s = pd.Series(rpt_capacities)
divisior = rpt_capacities[0]
result = s.apply(lambda x: x/divisior)
SOH = result.to_list()
plt.figure()
plt.plot(rpt_cycles, SOH)
# # #Plot LLI contributions
Qt = cccv_sol["Throughput capacity [A.h]"].entries
Q_SEI = cccv_sol["Loss of capacity to SEI [A.h]"].entries
Q_SEI_cr = cccv_sol["Loss of capacity to SEI on cracks [A.h]"].entries
Q_plating = cccv_sol["Loss of capacity to lithium plating [A.h]"].entries
Q_side = cccv_sol["Total capacity lost to side reactions [A.h]"].entries
Q_LLI = cccv_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()
# # #Plot LAM contributions and compare to LLI
Qt = cccv_sol["Throughput capacity [A.h]"].entries
LLI = cccv_sol["Loss of lithium inventory [%]"].entries
LAM_neg = cccv_sol["Loss of active material in negative electrode [%]"].entries
LAM_pos = cccv_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()
# # #Evaluate porosity changes due to pore clogging
eps_neg_avg = cccv_sol["X-averaged negative electrode porosity"].entries
eps_neg_sep = cccv_sol["Negative electrode porosity"].entries[-1,:]
eps_neg_CC = cccv_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()

Can someone please tell me what is wrong with my approach and how I can correctly/successfully apply this coupled degradation methodology to other parameter sets? @DrSOKane , I would really appreciate your inputs on this.

[MirAbbasAli2A]

@DrSOKane , do you have any advice on my approach of applying this to the mohtat dataset/how i can fix this zero porosity cutoff issue?

[brosaplanella]

It seems that the issue is related to the parameter values. When coupling many degradation effects, the parameters need to be carefully calibrated. This is extremely case dependent, so unfortunately we cannot provide a way to fix that.

A way to debug this would be to add one mechanism at a time, and adjust the parameters as you go.