Errors + Long Runtimes When Adding Non-Linear Diffusivity Dependency in Aging Analysis

[MirAbbasAli2A]

Hello PyBamm team and Community -

I am modelling coupled degradation using the OKane2022 parameter set. I have extended the model to include the non-linear diffusivity fitted expressions used by ORegan2022. However, I run into 2 problems. First, each cycle takes extensively long to complete (>4 min). Second, I keep getting the following error messages:

UserWarning: Q_Li=9.0975 Ah is greater than Q_p=8.7299 Ah.
warnings.warn(f"Q_Li={Q_Li:.4f} Ah is greater than Q_p={Q_p:.4f} Ah.")
Traceback (most recent call last):

File ~\anaconda3\lib\site-packages\pybamm\models\full_battery_models\lithium_ion\electrode_soh.py:260 in solve
sol = self._solve_full(inputs, ics)

File ~\anaconda3\lib\site-packages\pybamm\models\full_battery_models\lithium_ion\electrode_soh.py:337 in _solve_full
sol = sim.solve([0], inputs=inputs)

File ~\anaconda3\lib\site-packages\pybamm\simulation.py:643 in solve
self._solution = solver.solve(self.built_model, t_eval, **kwargs)

File ~\anaconda3\lib\site-packages\pybamm\solvers\base_solver.py:839 in solve
new_solution = self._integrate(

File ~\anaconda3\lib\site-packages\pybamm\solvers\algebraic_solver.py:223 in _integrate
raise pybamm.SolverError(

SolverError: Could not find acceptable solution: solver terminated successfully, but maximum solution error (0.5018226999198685) above tolerance (1e-06)

During handling of the above exception, another exception occurred:

Traceback (most recent call last):

File ~\anaconda3\lib\site-packages\pybamm\solvers\solution.py:873 in _get_cycle_summary_variables
esoh_sol = esoh_solver.solve(inputs)

File ~\anaconda3\lib\site-packages\pybamm\models\full_battery_models\lithium_ion\electrode_soh.py:269 in solve
raise split_error

File ~\anaconda3\lib\site-packages\pybamm\models\full_battery_models\lithium_ion\electrode_soh.py:264 in solve
sol = self._solve_split(inputs, ics)

File ~\anaconda3\lib\site-packages\pybamm\models\full_battery_models\lithium_ion\electrode_soh.py:350 in _solve_split
x0_sol = x0_sim.solve([0], inputs=inputs)

File ~\anaconda3\lib\site-packages\pybamm\simulation.py:643 in solve
self._solution = solver.solve(self.built_model, t_eval, **kwargs)

File ~\anaconda3\lib\site-packages\pybamm\solvers\base_solver.py:839 in solve
new_solution = self._integrate(

File ~\anaconda3\lib\site-packages\pybamm\solvers\algebraic_solver.py:223 in _integrate
raise pybamm.SolverError(

SolverError: Could not find acceptable solution: solver terminated successfully, but maximum solution error (0.5315554636851472) above tolerance (1e-06)

During handling of the above exception, another exception occurred:

Traceback (most recent call last):

File ~\anaconda3\lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
exec(code, globals, locals)

File c:\users\abbas\documents\phd docs\pybamm\agingmodel\brosaplanella-tec-reduced-model-978f8c3\agingmodelwithmechanismsanddynamiccycle.py:246
sol = model_solver()

File ~\anaconda3\lib\site-packages\pybamm\simulation.py:826 in solve
cycle_sol = pybamm.make_cycle_solution(

File ~\anaconda3\lib\site-packages\pybamm\solvers\solution.py:809 in make_cycle_solution
cycle_summary_variables = _get_cycle_summary_variables(cycle_solution, esoh_solver)

File ~\anaconda3\lib\site-packages\pybamm\solvers\solution.py:875 in _get_cycle_summary_variables
raise pybamm.SolverError(

SolverError: Could not solve for summary variables, run sim.solve(calc_esoh=False) to skip this step

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 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",
)

param = pybamm.ParameterValues("OKane2022")

Define function to update cracking rate

T = model.variables["Volume-averaged cell temperature [K]"]
sto = model.variables["Negative electrode stoichiometry"]
sto_pos = model.variables["Positive electrode stoichiometry"]

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

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

param["Ambient temperature [K]"] = 298.15
param["Initial temperature [K]"] = 298.15

def solvent_diffusivity_sei(T):
Ea = 37000
arrhenius_v2 = exp(Ea / constants.R * (1 / 298.15 - 1 / T))
return 2.5000000000000002e-22 * arrhenius_v2

def graphite_LGM50_diffusivity_ORegan2022(sto, T):

a0 = 11.17
a1 = -1.553
a2 = -6.136
a3 = -9.725
a4 = 1.85
b1 = 0.2031
b2 = 0.5375
b3 = 0.9144
b4 = 0.5953
c0 = -15.11
c1 = 0.0006091
c2 = 0.06438
c3 = 0.0578
c4 = 0.001356
d = 2092

D_ref = (
    10
    ** (
        a0 * sto
        + c0
        + a1 * exp(-((sto - b1) ** 2) / c1)
        + a2 * exp(-((sto - b2) ** 2) / c2)
        + a3 * exp(-((sto - b3) ** 2) / c3)
        + a4 * exp(-((sto - b4) ** 2) / c4)
    )
    * 3.0321  # correcting factor (see O'Regan et al 2021)
)

E_D_s = d * pybamm.constants.R
arrhenius_v3 = exp(E_D_s / pybamm.constants.R * (1 / 298.15 - 1 / T))

return D_ref * arrhenius_v3

def nmc_LGM50_diffusivity_ORegan2022(sto_pos, T):

a1 = -0.9231
a2 = -0.4066
a3 = -0.993
b1 = 0.3216
b2 = 0.4532
b3 = 0.8098
c0 = -13.96
c1 = 0.002534
c2 = 0.003926
c3 = 0.09924
d = 1449

D_ref_v2 = (
    10
    ** (
        c0
        + a1 * exp(-((sto_pos - b1) ** 2) / c1)
        + a2 * exp(-((sto_pos - b2) ** 2) / c2)
        + a3 * exp(-((sto_pos - b3) ** 2) / c3)
    )
    * 2.7  # correcting factor (see O'Regan et al 2021)
)

E_D_s_v2 = d * pybamm.constants.R
arrhenius_v4 = exp(E_D_s_v2 / pybamm.constants.R * (1 / 298.15 - 1 / T))

return D_ref_v2 * arrhenius_v4

param.update({"Outer SEI solvent diffusivity [m2.s-1]": solvent_diffusivity_sei(T)})
param.update({"Negative electrode diffusivity [m2.s-1]": graphite_LGM50_diffusivity_ORegan2022})
param.update({"Positive electrode diffusivity [m2.s-1]": nmc_LGM50_diffusivity_ORegan2022})

param.update({"Negative electrode cracking rate":anode_cracking_rate_Ai2020})

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
)

var_pts = {
"x_n" : 30, # x-direction, length, negative electrode
"x_s" : 30, # x-direction, length, separator
"x_p" : 30, # 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
}
submesh_types = model.default_submesh_types
submesh_types["negative particle"] = pybamm.MeshGenerator(
pybamm.Exponential1DSubMesh, submesh_params={"side": "right"}
)
submesh_types["positive particle"] = pybamm.MeshGenerator(
pybamm.Exponential1DSubMesh, submesh_params={"side": "right"}
)

#Run model
solver = pybamm.CasadiSolver(mode="safe", atol=1e-6, rtol=1e-3)
sim = pybamm.Simulation(model, experiment=cccv_experiment,parameter_values=param, var_pts=var_pts, solver=solver, submesh_types=submesh_types)

sol = sim.solve(save_at_cycles=10)

Could anyone please tell me what could be causing this error in the model and for that matter, if it is important to consider concentration-dependent diffusivity when modelling aging? Additionally, assuming this is required, any ideas on how I can make the sim run significantly faster than it currently does with concentration-dependent diffusivities?

Thank you for your help.