Error plotting the different contributions to LLI against throughput capacity in coupled degradation model

[Akila1993]

Hello Everyone,

I used the following code to run the degradation submodels suggested by @DrSOKane, However, when I tried to obtain the plotting as per the example code I got an empty plot. Also, I cannot see the following as model variables.

• Loss of capacity to negative SEI [A.h].entries
• Q_SEI_cr = sol["Loss of capacity to negative SEI on cracks [A.h]
• Q_plating = sol["Loss of capacity to negative lithium plating [A.h]"]

If @DrSOKane or someone can guide me about the mistake I'm making here it will be really appreciated.

Thank you.

import matplotlib.pyplot as plt
import numpy as np
from scipy.fft import fft
import pybamm
import time as timer

# Set random seed
np.random.seed(0)

# Original function for degradation parameters
pp = pybamm.ParameterValues("OKane2022")

# Function to get Chen2020 parameters
def get_chen2020_parameters():
    parameter_values = pybamm.ParameterValues("Chen2020")
    return parameter_values

# Function to get OKane2022 degradation parameters with k_cr included
def get_okane2022_degradation_parameters():
    degradation_parameters = {
        'Dead lithium decay constant [s-1]': 1e-06,
        'Dead lithium decay rate [s-1]': pp["Dead lithium decay rate [s-1]"],
        'Initial plated lithium concentration [mol.m-3]': 0.0,
        'Lithium metal partial molar volume [m3.mol-1]': 1.3e-05,
        'Lithium plating kinetic rate constant [m.s-1]': 1e-09,
        'Lithium plating transfer coefficient': 0.65,
        'Negative electrode LAM constant exponential term': 2.0,
        'Negative electrode LAM constant proportional term [s-1]': 2.7778e-07,
        "Negative electrode Paris' law constant b": 1.12,
        "Negative electrode Paris' law constant m": 2.2,
        "Negative electrode Poisson's ratio": 0.3,
        "Negative electrode Young's modulus [Pa]": 15000000000.0,
        'Negative electrode cracking rate': pp["Negative electrode cracking rate"],  # Function
        'Negative electrode critical stress [Pa]': 60000000.0,
        'Negative electrode initial crack length [m]': 2e-08,
        'Negative electrode initial crack width [m]': 1.5e-08,
        'Negative electrode number of cracks per unit area [m-2]': 3180000000000000.0,
        'Negative electrode partial molar volume [m3.mol-1]': 3.1e-06,
        'Negative electrode reference concentration for free of deformation [mol.m-3]': 0.0,
        'Negative electrode volume change': pp["Negative electrode volume change"],
        'Positive electrode LAM constant exponential term': 2.0,
        'Positive electrode LAM constant proportional term [s-1]': 2.7778e-07,
        'Positive electrode OCP entropic change [V.K-1]': 0.0,
        "Positive electrode Paris' law constant b": 1.12,
        "Positive electrode Paris' law constant m": 2.2,
        "Positive electrode Poisson's ratio": 0.2,
        "Positive electrode Young's modulus [Pa]": 375000000000.0,
        'Positive electrode active material volume fraction': 0.665,
        'Positive electrode cracking rate': pp["Positive electrode cracking rate"],  # Function
        'Positive electrode critical stress [Pa]': 375000000.0,
        'Positive electrode initial crack length [m]': 2e-08,
        'Positive electrode initial crack width [m]': 1.5e-08,
        'Positive electrode number of cracks per unit area [m-2]': 3180000000000000.0,
        'Positive electrode partial molar volume [m3.mol-1]': 1.25e-05,
        'Positive electrode reference concentration for free of deformation [mol.m-3]': 0.0,
        'Positive electrode volume change': pp["Positive electrode volume change"],
        'Typical plated lithium concentration [mol.m-3]': 1000.0,
        'Exchange-current density for stripping [A.m-2]': pp["Exchange-current density for stripping [A.m-2]"],
        'Exchange-current density for plating [A.m-2]': pp["Exchange-current density for plating [A.m-2]"],
        "k_cr": pybamm.Parameter("k_cr"),  # Define as a parameter object
    }
    return degradation_parameters

# Graphite cracking rate function
def graphite_cracking_rate_Ai2020(T_dim):
    """
    Calculates the graphite cracking rate using an Arrhenius-like relationship.
    T_dim: Dimensional temperature (in Kelvin)
    """
    Eac_cr = 0  # Activation energy for cracking [J/mol], update if needed
    arrhenius = pybamm.exp(Eac_cr / pybamm.constants.R * (1 / T_dim - 1 / 298.15))
    return pybamm.Parameter("k_cr") * arrhenius

# Function to update parameters
def update_parameters():
    # Load Chen2020 parameters
    chen2020_params = get_chen2020_parameters()

    # Load OKane2022 degradation parameters
    okane2022_degradation_params = get_okane2022_degradation_parameters()

    # Update Chen2020 parameters with OKane2022 degradation parameters
    chen2020_params.update(okane2022_degradation_params, check_already_exists=False)

    # Add graphite cracking rate function
    chen2020_params["k_cr"] = 50*3.9e-20

    return chen2020_params

# Update the parameters
updated_parameters = update_parameters()

# Define the model with configuration
model = pybamm.lithium_ion.DFN(
    {   
        #"thermal": "lumped",  # Lumped thermal submodel
        "surface form": "differential",  # Adding capacitance to the model
        "SEI": "solvent-diffusion limited",
        "SEI porosity change": "true",
        "particle mechanics": ("swelling and cracking", "swelling only"),
        "SEI on cracks": "true",
    }
)

# Variable points
var_pts = {
    "x_n": 5,
    "x_s": 5,
    "x_p": 5,
    "r_n": 30,
    "r_p": 30,
}


# Simulation parameters
cycle_numbers = [10]


# Store results for Nyquist plot
impedance_data = {}

# Loop through each cycle count
for cycle_number in cycle_numbers:
    exp = pybamm.Experiment(
        [
            (
                "Discharge at 1C until 2.5 V",  # Ageing cycles
                "Charge at 0.6C until 4.2 V (5 minute period)",
                "Hold at 4.2 V until C/100 (5 minute period)",
            )
        ]
        * cycle_number
    )

    sim = pybamm.Simulation(
        model,
        parameter_values=updated_parameters,
        experiment=exp,
        solver=pybamm.CasadiSolver(mode='safe without grid'),
        var_pts=var_pts,
    )

    sol = sim.solve()


Qt = sol["Throughput capacity [A.h]"].entries
Q_SEI = sol["Loss of capacity to SEI [A.h]"].entries
Q_SEI_cr = sol["Loss of capacity to SEI on cracks [A.h]"].entries
Q_side = sol["Total capacity lost to side reactions [A.h]"].entries
Q_LLI = (
    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_side, label="All side reactions", linestyle=(0, (6, 1)))
plt.xlabel("Throughput capacity [A.h]")
plt.ylabel("Capacity loss [A.h]")
plt.legend()
plt.show()

Also, I have tried to run the example code available here exactly as it is, however it gives me following error message.

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/util.py:59, in FuzzyDict.__getitem__(self, key)
     58 try:
---> 59     return super().__getitem__(key)
     60 except KeyError:

KeyError: 'Loss of capacity to negative SEI [A.h]'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
Cell In[2], line 2
      1 Qt = sol["Throughput capacity [A.h]"].entries
----> 2 Q_SEI = sol["Loss of capacity to negative SEI [A.h]"].entries
      3 Q_SEI_cr = sol["Loss of capacity to negative SEI on cracks [A.h]"].entries
      4 Q_plating = sol["Loss of capacity to negative lithium plating [A.h]"].entries

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/solvers/solution.py:537, in Solution.__getitem__(self, key)
    534     return self._variables[key]
    535 else:
    536     # otherwise create it, save it and then return it
--> 537     self.update(key)
    538     return self._variables[key]

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/solvers/solution.py:476, in Solution.update(self, variables)
    474 cumtrapz_ic = None
    475 pybamm.logger.debug("Post-processing {}".format(key))
--> 476 vars_pybamm = [model.variables_and_events[key] for model in self.all_models]
    478 # Iterate through all models, some may be in the list several times and
    479 # therefore only get set up once
    480 vars_casadi = []

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/solvers/solution.py:476, in <listcomp>(.0)
    474 cumtrapz_ic = None
    475 pybamm.logger.debug("Post-processing {}".format(key))
--> 476 vars_pybamm = [model.variables_and_events[key] for model in self.all_models]
    478 # Iterate through all models, some may be in the list several times and
    479 # therefore only get set up once
    480 vars_casadi = []

File ~/pybamm-eis/env/lib/python3.9/site-packages/pybamm/util.py:62, in FuzzyDict.__getitem__(self, key)
     60 except KeyError:
     61     best_matches = self.get_best_matches(key)
---> 62     raise KeyError(f"'{key}' not found. Best matches are {best_matches}")

KeyError: "'Loss of capacity to negative SEI [A.h]' not found. Best matches are ['Loss of capacity to SEI [A.h]', 'Loss of capacity to SEI on cracks [A.h]', 'Loss of capacity to lithium plating [A.h]']"