Errors in aging model + discrepancy in OKane2022 beta parameter

[MirAbbasAli2A]

Hello PyBamm Team -

I am running an aging model considering all aging mechanisms. I have extended the solvent diffusivity within the SEI to be temperature dependent with the arrhenius relationship and multiplied the LAM proportionality constant (beta) by 20 at the anode and by 1 at the cathode. When I perform this setup, I get an IDAconv fail error. Below is a working example of my code:

#Import necessary libraries
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

model = pybamm.lithium_ion.DFN(options = {
    "thermal":"lumped",
    "dimensionality":0,
    "cell geometry": "arbitrary", 
    "SEI": "solvent-diffusion limited",
    "SEI porosity change": "true",
    "particle mechanics": ("swelling and cracking", "swelling and cracking"),
    "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",
    )

param = pybamm.ParameterValues("OKane2022")

param = set_thermal_parameters(param, 16, 2.32e6, temperature)

param["Ambient temperature [K]"] = 298.15
param["Initial temperature [K]"] = 298.15
param["Negative electrode LAM constant proportional term [s-1]"] = 20*2.7778e-07
param["Positive electrode LAM constant proportional term [s-1]"] = 2.7778e-07

T = model.variables["Volume-averaged cell temperature [K]"]
 
def solvent_diffusivity_sei(T):
    Ea = 37000
    arrhenius = exp(Ea / constants.R * (1 / 298.15 - 1 / T))
    return 2.5000000000000002e-22 * arrhenius
 
param.update({"Outer SEI solvent diffusivity [m2.s-1]": solvent_diffusivity_sei(T)})

pybamm.set_logging_level("NOTICE")

#Simple cycling experiment
N=550
 
cccv_experiment = pybamm.Experiment([
    ("Charge at 1C until 4.2V",
    "Hold at 4.2V until C/20",
    "Discharge at 1C until 2.5V",
    "Hold at 2.5V until C/20",
    )
    ]*N
    )

#Refine mesh
var_pts = {
    "x_n" : 20,  # x-direction, length, negative electrode
    "x_s" : 20,  # x-direction, length, separator
    "x_p" : 20,  # x-direction, length, positive electrode
    "r_n" : 30,  # number of volumes in the radial direction, negative particle
    "r_p" : 30   # number of volumes in the radial direction, positive particle
}

#Run model
sim = pybamm.Simulation(model, experiment=cccv_experiment,parameter_values=param, var_pts=var_pts)
sol = sim.solve(save_at_cycles= 10)

It is worth noting that without multiplying beta by 20, the simulation runs just fine. I also wanted to know what the actual value of beta is? In the OKane2022 paper it is 0.001 but in the paramater set it is 2.7778e-07. When I run with 0.001, the code fails after 20 cycles with a Linear solver did not converge error.

FYI, here is the error code I am getting:

As you can see, it fails at 307 cycles.

@DrSOKane , @tinosulzer, @brosaplanella I would be truly grateful if you can suggest anything to help fix this problem I am having with my model.

Thanks in advance.

[MirAbbasAli2A]

Update: I refined the mesh to have 30 pts for x_n, x_s and x_p as well and it ran till 449 cycles. However, it then failed with the same exact error as before. Any advice/suggestions on how I can fix this?

[MirAbbasAli2A]

Another update: it runs successfully if I remove the cracking option for the cathode per the coupled degradation notebook. Are there any issues with trying to model particle cracking in the cathode? Wouldn't it just be calculating the loss in solid phase volume fraction just like at the anode?

[DrSOKane]

Hi @MirAbbasAli2A, I'm pleased to see you managed to resolve the problem on your own. Now I'm here, I can give you a detailed explanation of why that solution worked.

There is no point using the "swelling and cracking" option in the positive electrode if you're modelling degradation, because there is no SEI to grow on the crack. The option is available if you're interested in studying how the crack propagates within the positive electrode particles, but it does not affect the degradation behaviour. (It also allows SEI on cracks to be enabled in graphite half-cells, where the graphite is the positive electrode, but you don't specify it as a tuple for that case anyway, because the particles mechanics model isn't applicable to lithium metal electrodes, which has no particles.)

The crack model solves Paris' Law of fatigue crack propagation. The solution to Paris' Law is the crack length grows exponentially with time. Exponential growth poses numerical difficulties, so Paris' Law should only be solved where necessary.

The "swelling only" model is much easier to solve, and still gives you the loss of solid phase volume fraction, which is what you're interested in for the positive electrode.

[DrSOKane]

According to the poster, cycling at the same rate but with only 60% DoD means the knee point is avoided and the cell lasts a long time. Which SoC range is avoided? The high one? Or the low one?

[MirAbbasAli2A]

@DrSOKane , they actually reduce the discharge window from both ends. So instead of discharging from 4.2V to 2.5V, they discharge from 3.98V to 3.42V.

See this table from the paper I attached in the previous reply:

[MirAbbasAli2A]

@DrSOKane, I thought I would let you know that I got good correlation to the data set at a 1C discharge rate at 80pDOD and 100pDOD:

This was done by multiplying the cracking rate by 50, applying an arrhenius correlation to the SEI solvent diffusivity of 2.5E-21, a decay constant of 2.5E-6 and MOST importantly, changing the LAM beta constant of the anode to 80 for the 100pDOD case and 70 for the 80PDOD case. I found that the LAM beta factor makes the most significant change.

However when I try a beta factor of 50 for the anode (because all other values seem to fail) for a 60pDOD case (keeping everything else constant from the previous runs), then the correlation is not so great:

Are the model parameters meant to be adjusted for every test condition or are they meant to stay constant? If they are meant to be adjusted, in the absence of experimental data, is this just meant to be done intuitively (i.e. low c-rates and depth of discharge means limited plating and LAM? Or is there a general guideline to tune these parameters according to the input conditions?

[DrSOKane]

You should not adjust the model parameters for each test condition. The test of a truly accurate parameter set is its ability to predict correct results for the same cell under different conditions.

When you stop to consider that lithium-ion batteries are known to have more degradation mechanisms than the ones currently included in PyBaMM, the quality of the results you have shown is impressive.

Another limitation of the OKane2022 parameter set is that the solid-phase diffusivities are functions of temperature only, not concentration. This will have a large impact on the particle mechanics models. Adding concentration dependence will increase the accuracy of the model, but also add significant computational expense.

[MirAbbasAli2A]

@DrSOKane I do see where you are coming from. However, the part that troubles me is how can we assume the parameter set is the same when the same degree of each degradation mechanism is not observed for each input condition. For example, how can we assume that the beta value is the same for a 1C 100% DOD vs. a 0.5C 60% DOD. The former situation has way more stress-driven LAM vs. the latter? Is this difference just supposed to be accounted for by the hydrostatic stress? This also applies to any of the other fitted parameters (e.g how can the plating rate be the same for a charging rate that is 0.5C vs 2C?)

I tried doing the exercise of extending the model to have concentration dependent solid phase diffusivities over here with the ORegan2022 definition. However, like you said, it takes about 20 min to complete 1 cycle despite me meshing it just like Oregan2022 did and even then, it runs for 20 cycles before failing.