Unusual reverse capacity trend during aging simulation despite changing initial concentrations

[MirAbbasAli2A]

Hello Pybamm Team -

I was running an aging experiment to reproduce Figure 5 in this paper written by @tinosulzer. The only differences are that I was using a DFN model and using a coupled degradation analysis. I get a very unusual trend where my discharge capacity increases with cycle number.

When analyzing the problem further, I find that the porosity is just changing at the separator and not anywhere else in the negative electrode.

Additionally, I kept getting the warning of Q_Li>Q_p. I have read over previous discussions and have found that changing the initial concentrations of the electrodes was a potential solution, but this didn't help (at least didn't help beyond a certain point).

Here is my code:

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"]=2.0
param["Initial inner SEI thickness [m]"]=0.0
param["Initial outer SEI thickness [m]"] = 2.5e-9
param["SEI growth activation energy [J.mol-1]"]=5000
param["Inner SEI partial molar volume [m3.mol-1]"] = 5.22e-5
param["Outer SEI partial molar volume [m3.mol-1]"] = 5.22e-5
# param["Initial concentration in negative electrode [mol.m-3]"]= 48.8682*5
# param["Initial concentration in positive electrode [mol.m-3]"]=31513.0/70
# 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-10}, 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-07}, 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.84e-09)}, 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.98e-18)}, 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 = (5.29e-25)
    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 = 5000
    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/20",
    "Rest for 1 hour",
    "Discharge at C/20 until 2.8V",
    "Hold at 2.8V until C/20"
    "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 = 5
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)

Can someone please help explain to me where I am going wrong?

Thank you for your assistance.

[MirAbbasAli2A]

Pybamm team, does anyone have any thoughts on how I can fix this behaviour?

@tinosulzer , @brosaplanella , @TomTranter , @DrSOKane , @rtimms

[valentinsulzer]

It's not immediately obvious what is going wrong, you will need to keep debugging by looking at the outputs like you are doing.
For a start you could try setting "SEI porosity change" and/or "lithium plating porosity change" to zero to see if one of those is driving the behavior. Also, limit to 50 cycles to debug so you don't have to wait for 250 cycles

[MirAbbasAli2A]

@tinosulzer , thank you for your reply. yes I did try turning the porosity changes off but it was yielding similar behaviour. I will take your advice and proceed with debugging.