How should crack model parameters, mechanical properties parameters and LAM parameters be implemented from Ai2020 into another parameter set?

[Mykle99]

Hi! I am trying to use li-plating and LAM submodels with Mohtat2020 parameters, but I am getting errors. Ai2020 params have been inserted into Mohtat2020 to account for the "missing" parameters. Also note that the code worked as intended before introducing LAM submodels and LAM parameters (now they are tested alone), but I have solved the model with SEI and Li-plating, which I want to do together with LAM, but I cannot even get LAM to work by itself. Do I need to use a particular discretisation, and thereafter a particular solver?

############################################################################
# Required functions

# https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/
# lithium_ion/negative_electrodes/graphite_Ai2020/graphite_volume_change_Ai2020.py
# ingen imports
def graphite_volume_change_Ai2020(sto):
    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 ** 9
        + p2 * sto ** 8
        + p3 * sto ** 7
        + p4 * sto ** 6
        + p5 * sto ** 5
        + p6 * sto ** 4
        + p7 * sto ** 3
        + p8 * sto ** 2
        + p9 * sto
        + p10
    )
    return t_change


# https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/
# lithium_ion/negative_electrodes/graphite_Ai2020/graphite_cracking_rate_Ai2020.py
from pybamm import Parameter, constants, exp

def graphite_cracking_rate_Ai2020(T_dim):
    k_cr = 3.9e-20
    T_ref = Parameter("Reference temperature [K]")
    Eac_cr = Parameter(
        "Negative electrode activation energy for cracking rate [J.mol-1]"
    )
    arrhenius = exp(Eac_cr / constants.R * (1 / T_dim - 1 / T_ref))
    return k_cr * arrhenius

# https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/
# lithium_ion/positive_electrodes/lico2_Ai2020/lico2_volume_change_Ai2020.py
from pybamm import Parameter

def lico2_volume_change_Ai2020(sto):
    omega = Parameter("Positive electrode partial molar volume [m3.mol-1]")
    c_p_max = Parameter("Maximum concentration in positive electrode [mol.m-3]")
    t_change = omega * c_p_max * sto
    return t_change

# https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/
# lithium_ion/positive_electrodes/lico2_Ai2020/lico2_cracking_rate_Ai2020.py
from pybamm import Parameter, constants, exp

def lico2_cracking_rate_Ai2020(T_dim):
    k_cr = 3.9e-20
    T_ref = Parameter("Reference temperature [K]")
    Eac_cr = Parameter(
        "Positive electrode activation energy for cracking rate [J.mol-1]"
    )
    arrhenius = exp(Eac_cr / constants.R * (1 / T_dim - 1 / T_ref))
    return k_cr * arrhenius

############################################################################

# My script:
# setting the logging level to notice for visualization
pybamm.set_logging_level("NOTICE")

# defining important variables for the experiments which are conducted
# how many cycles for each experiment + which SoC ranges to test in the
# experiments, C-rate for charge and discharge is also set. (1C = 5A)
N_cycles = 2
SoC_ranges = np.array([[0.1, 0.4], [0.2, 0.5], [0.3, 0.6], 
                       [0.4, 0.7], [0.5, 0.8], [0.6, 0.9]])
charging_current, discharging_current = 5, 5

LAM_options = {
        "particle mechanics": "swelling and cracking",
        "loss of active material" : "stress-driven"
}

mohtat2020_params = pybamm.parameter_sets.Mohtat2020

# Lithium plating parameters inserted from Okane2020 parameter set:
# https://github.com/pybamm-team/PyBaMM/tree/develop/pybamm/input/parameters/
# lithium_ion/lithium_platings/okane2020_Li_plating
mohtat2020_params["lithium plating"] = "okane2020_Li_plating"

dfn = pybamm.lithium_ion.DFN(options = LAM_options)

params = pybamm.ParameterValues(chemistry = mohtat2020_params)

###################################################################################################
# Negative electrode parameters inserted from Ai2020 parameter set:
# https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/
# lithium_ion/negative_electrodes/graphite_Ai2020/parameters.csv

# Negative electrode, mechanical properties:
params.update({"Negative electrode Poisson's ratio" : 0.3}, check_already_exists = False)
params.update({"Negative electrode Young's modulus [Pa]" : 15e9}, check_already_exists = False)
params.update({"Negative electrode reference concentration for free of deformation [mol.m-3]" : 0}, check_already_exists = False)
params.update({"Negative electrode partial molar volume [m3.mol-1]" :  3.1e-6}, check_already_exists = False)
params.update({"Negative electrode volume change" : graphite_volume_change_Ai2020}, check_already_exists = False)

# Negative electrode, crack model:
params.update({"Negative electrode initial crack length [m]" : 20e-9}, check_already_exists = False)
params.update({"Negative electrode initial crack width [m]" : 15e-9}, check_already_exists = False)
params.update({"Negative electrode number of cracks per unit area [m-2]" : 3.18e15}, check_already_exists = False)
params.update({"Negative electrode Paris' law constant b" : 1.12}, check_already_exists = False)
params.update({"Negative electrode Paris' law constant m" : 2.2}, check_already_exists = False)
params.update({"Negative electrode cracking rate" : graphite_cracking_rate_Ai2020}, check_already_exists = False)
params.update({"Negative electrode activation energy for cracking rate [J.mol-1]" : 0}, check_already_exists = False)

# Negative electrode, loss of active materials (LAM) model:
params.update({"Negative electrode LAM constant proportional term [s-1]" : 0}, check_already_exists = False)
params.update({"Negative electrode LAM constant exponential term" : 2}, check_already_exists = False)
params.update({"Negative electrode critical stress [Pa]" : 60e6}, check_already_exists = False)

# Positve electrode parameters inserted from Ai2020 parameter set:
# https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/
# lithium_ion/positive_electrodes/lico2_Ai2020/parameters.csv

# Maybe this one is needed? ("Missing" parameter in Mohtat2020, found in Ai2020)
# Microstructure:
# mohtat2020_params["Positive electrode surface area to volume ratio [m-1]"] = 620000

# Positve electrode, mechanical properties:
params.update({"Positive electrode Poisson's ratio" : 0.2}, check_already_exists = False)
params.update({"Positive electrode Young's modulus [Pa]" : 375e9}, check_already_exists = False)
params.update({"Positive electrode reference concentration for free of deformation [mol.m-3]" : 0}, check_already_exists = False)
params.update({"Positive electrode partial molar volume [m3.mol-1]" : -7.28e-7}, check_already_exists = False)
params.update({"Positive electrode volume change" : lico2_volume_change_Ai2020}, check_already_exists = False)

# Positve electrode, crack model:
params.update({"Positive electrode initial crack length [m]" : 20e-9}, check_already_exists = False)
params.update({"Positive electrode initial crack width [m]" : 15e-9}, check_already_exists = False)
params.update({"Positive electrode number of cracks per unit area [m-2]" : 3.18e15}, check_already_exists = False)
params.update({"Positive electrode Paris' law constant b" : 1.12}, check_already_exists = False)
params.update({"Positive electrode Paris' law constant m" : 2.2}, check_already_exists = False)
params.update({"Positive electrode cracking rate" : lico2_cracking_rate_Ai2020}, check_already_exists = False)
params.update({"Positive electrode activation energy for cracking rate [J.mol-1]" : 0}, check_already_exists = False)

# Positve electrode, loss of active materials (LAM) model:
params.update({"Positive electrode LAM constant proportional term [s-1]" : 2.78e-13}, check_already_exists = False)
params.update({"Positive electrode LAM constant exponential term" : 2}, check_already_exists = False)
params.update({"Positive electrode critical stress [Pa]" : 375e6}, check_already_exists = False)

# Cell parameter inserted from Ai2020:
# https://github.com/pybamm-team/PyBaMM/blob/develop/pybamm/input/parameters/
# lithium_ion/cells/Enertech_Ai2020/parameters.csv
params.update({"Cell thermal expansion coefficient [m.K-1]" : 1.1E-6}, check_already_exists = False)
# Maybe Cell emissivity 0.95 has to be set? It is set in Ai2020, but does not exist in Mohtat2020

# 1 + dllnf/dllnc is also different in Ai2020 and Mohtat2020
###################################################################################################

# Have tried to do a bunch of stuff to the solver in order to get the sim to run. Have also tried different solvers implemented in 
# PyBaMM
solver = pybamm.CasadiSolver(atol = 1e-2, rtol = 1e-2, root_tol = 1e-2, extrap_tol = 1e-2, max_step_decrease_count = 10, mode = "fast with events", extra_options_setup = {"max_num_steps" : 100000}) # , "print_stats" : True

# running an experiment to get the initial battery cell capacity 
mapping_experiment = pybamm.Experiment([
        ("Charge at 1 A for 1 minute")]   
)

# extracting "Capacity [A.h]" by solving the model
mapping_sim = pybamm.Simulation(dfn, parameter_values = params, solver = solver, experiment = mapping_experiment)
mapping_sol = mapping_sim.solve()
mapping_solution = mapping_sim.solution
initial_cap = mapping_solution.summary_variables["Capacity [A.h]"][0]

# making empty lists such that information can be appended at the
# end of each experiment
sims, sols, solutions,  = [], [], []
# iterating through defined SoC ranges
for SoC_range in SoC_ranges:
    # label for each experiment
    # calculating values needed for the actual experiment 
    initial_lower_SoC = SoC_range[0] * initial_cap
    initial_upper_SoC =  SoC_range[1] * initial_cap
    print(initial_lower_SoC, initial_upper_SoC, initial_cap)
    charge_time = (initial_upper_SoC - initial_lower_SoC)/charging_current          # t [h] = Ah/A
    discharge_time = (initial_upper_SoC - initial_lower_SoC)/discharging_current    # t [h] = Ah/A

    # setting up the first cycle of the experiment, this is done to ensure 
    # that the initial SoC is set to the correct value
    first_cycle_experiment = pybamm.Experiment([
        (f"Charge at {charging_current} A for {charge_time} hours", 
         f"Discharge at {discharging_current} A for {discharge_time} hours")]      
    )

    sim = pybamm.Simulation(dfn, parameter_values = params, solver = solver, experiment = first_cycle_experiment)
    # initial_soc is set to the correct starting value
    sol = sim.solve(initial_soc = SoC_range[0])
    solution = sim.solution

    for cycle_N in range(1, N_cycles): # 1 to N since the first cycle is conducted outside this for-loop
        # finding the remaining capacity
        remaining_cap = solution.summary_variables["Capacity [A.h]"][cycle_N-1]
        # calculating new capacity range, discharge and charge times, for constant SoC range throughout the experiment
        newCap_lower_SoC, newCap_upper_SoC = SoC_range[0] * remaining_cap, SoC_range[1] * remaining_cap
        charge_time = (newCap_upper_SoC - newCap_lower_SoC)/charging_current
        discharge_time = (newCap_upper_SoC - newCap_lower_SoC)/discharging_current
        remain_cycles_experiment = pybamm.Experiment([
            (f"Charge at {charging_current} A for {charge_time} hours", 
             f"Discharge at {discharging_current} A for {discharge_time} hours")]      
        )
        sim = pybamm.Simulation(dfn, parameter_values = params, solver = solver, experiment = remain_cycles_experiment)
        sol = sim.solve(starting_solution = sol)
        solution = sim.solution

    sims.append(sim)
    sols.append(sol)
    solutions.append(solution)
    
pybamm.dynamic_plot(sols)

Error:

2022-04-28 17:54:32,319 - [NOTICE] simulation.solve(789): Cycle 1/1, step 1/2: Charge at 5 A for 0.2981359607567931 hours
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:849: Calculating Jacobian failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:849: Calculating Jacobian failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:849: Calculating Jacobian failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:849: Calculating Jacobian failed
psetup failed: .../casadi/interfaces/sundials/idas_interface.cpp:849: Calculating Jacobian failed
The residual routine or the linear setup or solve routine had a recoverable error, but IDACalcIC was unable to recover.

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\pybamm\solvers\casadi_solver.py:693, in CasadiSolver._run_integrator(self, model, y0, inputs_dict, inputs, t_eval, use_grid, extract_sensitivities_in_solution)
    692 timer = pybamm.Timer()
--> 693 casadi_sol = integrator(
    694     x0=y0_diff, z0=y0_alg, p=inputs_with_tmin, **self.extra_options_call
    695 )
    696 integration_time = timer.time()

File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\casadi\casadi.py:13453, in Function.__call__(self, *args, **kwargs)
  13451   else:
  13452 # Named inputs -> return dictionary
> 13453     return self.call(kwargs)

File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\casadi\casadi.py:12324, in Function.call(self, *args)
  12257 """
  12258   Generate a Jacobian function of output oind with respect to input iind.
  12259 
   (...)
  12322 
  12323 """
> 12324 return _casadi.Function_call(self, *args)

RuntimeError: .../casadi/interfaces/sundials/idas_interface.cpp:591: IDACalcIC returned "IDA_NO_RECOVERY". Consult IDAS documentation.

During handling of the above exception, another exception occurred:

SolverError                               Traceback (most recent call last)
Input In [29], in <cell line: 191>()
    207 sim = pybamm.Simulation(dfn, parameter_values = params, solver = solver, experiment = first_cycle_experiment)
    208 # initial_soc is set to the correct starting value
--> 209 sol = sim.solve(initial_soc = SoC_range[0])
    210 solution = sim.solution
    212 for cycle_N in range(1, N_cycles): # 1 to N since the first cycle is conducted outside this for-loop
    213     # finding the remaining capacity

File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\pybamm\simulation.py:802, in Simulation.solve(self, t_eval, solver, check_model, save_at_cycles, calc_esoh, starting_solution, initial_soc, **kwargs)
    800 # Make sure we take at least 2 timesteps
    801 npts = max(int(round(dt / exp_inputs["period"])) + 1, 2)
--> 802 step_solution = solver.step(
    803     current_solution,
    804     model,
    805     dt,
    806     npts=npts,
    807     save=False,
    808     **kwargs,
    809 )
    810 steps.append(step_solution)
    811 current_solution = step_solution

File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\pybamm\solvers\base_solver.py:1261, in BaseSolver.step(self, old_solution, model, dt, npts, external_variables, inputs, save)
   1254 pybamm.logger.verbose(
   1255     "Stepping for {:.0f} < t < {:.0f}".format(
   1256         t * model.timescale_eval,
   1257         (t + dt_dimensionless) * model.timescale_eval,
   1258     )
   1259 )
   1260 timer.reset()
-> 1261 solution = self._integrate(model, t_eval, ext_and_inputs)
   1262 solution.solve_time = timer.time()
   1264 # Check if extrapolation occurred

File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\pybamm\solvers\casadi_solver.py:177, in CasadiSolver._integrate(self, model, t_eval, inputs_dict)
    173 # Create an integrator with the grid (we just need to do this once)
    174 self.create_integrator(
    175     model, inputs, t_eval, use_event_switch=use_event_switch
    176 )
--> 177 solution = self._run_integrator(
    178     model, model.y0, inputs_dict, inputs, t_eval
    179 )
    180 # Check if the sign of an event changes, if so find an accurate
    181 # termination point and exit
    182 solution = self._solve_for_event(solution, init_event_signs)

File ~\AppData\Local\Programs\Python\Python39\lib\site-packages\pybamm\solvers\casadi_solver.py:743, in CasadiSolver._run_integrator(self, model, y0, inputs_dict, inputs, t_eval, use_grid, extract_sensitivities_in_solution)
    740     return sol
    741 except RuntimeError as e:
    742     # If it doesn't work raise error
--> 743     raise pybamm.SolverError(e.args[0])

SolverError: .../casadi/interfaces/sundials/idas_interface.cpp:591: IDACalcIC returned "IDA_NO_RECOVERY". Consult IDAS documentation.

[Mykle99]

Working solution found:

dfn = pybamm.lithium_ion.DFN()
print(dfn.default_var_pts)

-> {'x_n': 20, 'x_s': 20, 'x_p': 20, 'r_n': 20, 'r_p': 20, 'y': 10, 'z': 10, 'R_n': 30, 'R_p': 30}

The default mesh is not fine enough. Manually set the mesh to be finer with var_pts by increasing the number of points.

# most values increased to 40, up from 20. y, z, R_n, and R_p kept at original values
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
}

# insert new mesh resolution into sim
sim = pybamm.Simulation(dfn, var_pts = var_pts)
# model can now be solved
sol = sim.solve()

[RobinBluehub]

How many cycles can you cycle through? Can it be cycled for thousands of cycles?