EIS for coupled-degradation

[Akila1993]

Hi All,

I have tried to obtain EIS using PyBaMM-EIS for the coupled-degradation model proposed here by @DrSOKane . The code for the said scenario is mentioned below.

And it gives the following error: SolverError: Cannot use ODE solver 'Scipy solver (BDF)' to solve DAE model

Can someone guide me in this regard please?

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

Set up

options = {
"thermal": "lumped",
"surface form": "differential",
"SEI": "solvent-diffusion limited",
"SEI porosity change": "true",
"lithium plating": "partially reversible",
"lithium plating porosity change": "true", # alias for "SEI porosity change"
"particle mechanics": ("swelling and cracking", "swelling only"),
"SEI on cracks": "true",
"loss of active material": "stress-driven",
"calculate discharge energy": "true", # for compatibility with older PyBaMM versions
}

model = pybamm.lithium_ion.DFN(options,name="DFN",)
param = pybamm.ParameterValues("OKane2022")

frequencies = np.logspace(-4, 2, 30)

Time domain

I = 50 * 1e-3
number_of_periods = 20
samples_per_period = 16

def current_function(t):
return I * pybamm.sin(2 * np.pi * pybamm.InputParameter("Frequency [Hz]") * t)

param["Current function [A]"] = current_function

start_time = timer.time()

sim = pybamm.Simulation(model, parameter_values=param,solver=pybamm.ScipySolver())

impedances_time = []
for frequency in frequencies:
# Solve
period = 1 / frequency
dt = period / samples_per_period
t_eval = np.array(range(0, 1 + samples_per_period * number_of_periods)) * dt
sol = sim.solve(t_eval, inputs={"Frequency [Hz]": frequency})
# Extract final two periods of the solution
time = sol["Time [s]"].entries[-3 * samples_per_period - 1 :]
current = sol["Current [A]"].entries[-3 * samples_per_period - 1 :]
voltage = sol["Terminal voltage [V]"].entries[-3 * samples_per_period - 1 :]
# FFT
current_fft = fft(current)
voltage_fft = fft(voltage)
# Get index of first harmonic
idx = np.argmax(np.abs(current_fft))
impedance = -voltage_fft[idx] / current_fft[idx]
impedances_time.append(impedance)

end_time = timer.time()
time_elapsed = end_time - start_time
print("Time domain method: ", time_elapsed, "s")

Frequency domain

methods = ["direct"]
impedances_freqs = []
for method in methods:
start_time = timer.time()
eis_sim = pbeis.EISSimulation(model, parameter_values=parameter_values)
impedances_freq = eis_sim.solve(frequencies, method)
end_time = timer.time()
time_elapsed = end_time - start_time
print(f"Frequency domain ({method}): ", time_elapsed, "s")
impedances_freqs.append(impedances_freq)

Compare

_, ax = plt.subplots()
ax = pbeis.nyquist_plot(impedances_time, ax=ax, label="Time", alpha=0.7)
for i, method in enumerate(methods):
ax = pbeis.nyquist_plot(
impedances_freqs[i], ax=ax, label=f"Frequency ({method})", alpha=0.7
)
ax.legend()
plt.suptitle(f"{model.name}")
plt.show()

[TomTranter]

Try removing the solver=pybamm.ScipySolver() from the Simulation. It should default to Casadi and that can solve DFN model

[Akila1993]

Hi @DrSOKane ,

I have included all the changes and it gives me the promising results (Did not include SEI on cracks for now). And also as you mentioned the EIS behaviour is quite interesting with coupled degradation models.

Thanks for your time and invaluable comments and advises.

[Akila1993]

Hello @DrSOKane ,

Is it possible use other parameter sets such as Chen2020 with for the following options.

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",
        "lithium plating": "partially reversible",
        "lithium plating porosity change": "true",  # alias for "SEI porosity change"
        "particle mechanics":"swelling only",
        "SEI on cracks": "true",
        "loss of active material": "stress-driven",
        
    }
)

I am getting errors with lithium plating options and subsequently we turn of plating option we cannot get "particle mechanics", "SEI on cracks": "loss of active material": too.

Please advise.

Thanks in advance.

[brosaplanella]

Some model combinations might be incompatible, can you please paste the full error message to check if that's the case?

[DrSOKane]

OKane2022 is the only ready-made parameter set that works with all degradation options. It was designed for that purpose.

Fortunately, you can add degradation parameters to any parameter set by making a copy of the parameter file and copy-pasting the degradation parameters you need from OKane2022. You can then load the parameters using

from your_parameter_set import get_parameter_values
parameter_values = get_parameter_values()

[Akila1993]

@DrSOKane Thank you for your input. It works!