Cyclic experiment with degradation mechanism and EIS

[Akila1993]

Hello Everyone,

As per the code in this discussion (ignore the errors in the discussion as it is solved), I am getting the EIS for different C-rates and different cycle numbers. Please see the below EIS figure for SEI growth enabled model.

Model options:
SEI – Solvent Diffusion Limited
Porosity change – True

I understand the behaviour for C-rate is reasonable because the higher C-rates have the impact on charge-transfer resistance and diffusion time constant.

But still I have the doubt about why increasing number of cycles does not reflect any significant changes in EIS.

Please provide your insights in this regard.

Thank you.

[Akila1993]

@brosaplanella @tinosulzer Would appreciate it if you could give your feedback in this regard.

[Akila1993]

Hi Everyone @brosaplanella @tinosulzer @rtimms ,

Does anyone have a comment on this? Would appreciate it if you could give your thoughts.

I still couldn't figure out the issue.

Thank you.

[brosaplanella]

Do you have a minimum working example? Changing the number of cycles should definitely affect EIS, so maybe you are taking the initial states for all cases rather than the final states.

[brosaplanella]

Hi! I meant adjusting the SEI reaction rate parameter, not the diffusion one, as this will be more likely to impact the degradation. I have tried your code but the Nyquist plot comes out as a flat line. Can you check that works at your end? It would also help to simplify the code as much as possible, while retaining the issues, so it is easier to debug.

[Akila1993]

Hi @brosaplanella ,

Yes. the code is working for me I think the current amplitude should be changed to 1mA in the previous code., Now I tried with the SEI reaction rate changing and ran for 5 cycles and 25 cycles but the EIS is the same, please look at the attached for the exact code I used.

I also have attached the impedance that I have generated from this code.

EIS with 10X SEI rate 25 cycle.csv
EIS with 10X SEI rate 5 cycle.csv

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

pp = pybamm.ParameterValues("OKane2022")

model = pybamm.lithium_ion.DFN(
    {   
        "thermal": "lumped", #lumped thermal submodel 
        "surface form": "differential", #Adding capacitance to the model 
        "SEI":"solvent-diffusion limited", #"reaction 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",
        
    }
)

var_pts = {
    "x_n": 5,  # negative electrode
    "x_s": 5,  # separator
    "x_p": 5,  # positive electrode
    "r_n": 30,  # negative particle
    "r_p": 30,  # positive particle
}

def get_chen2020_parameters():
    parameter_values = pybamm.ParameterValues("Chen2020")
    parameters = parameter_values
    return parameters

def get_okane2022_degradation_parameters():

 degradation_parameters = {
 'Dead lithium decay constant [s-1]': 1e-06,
 'Dead lithium decay rate [s-1]':pp["Dead lithium decay rate [s-1]"],
 'Initial plated lithium concentration [mol.m-3]': 0.0,
 'Lithium metal partial molar volume [m3.mol-1]': 1.3e-05,
 'Lithium plating kinetic rate constant [m.s-1]': 1e-09,
 'Lithium plating transfer coefficient': 0.65,
 'Negative electrode LAM constant exponential term': 2.0,
 'Negative electrode LAM constant proportional term [s-1]': 2.7778e-07,
 "Negative electrode Paris' law constant b": 1.12,
 "Negative electrode Paris' law constant m": 2.2,
 "Negative electrode Poisson's ratio": 0.3,
 "Negative electrode Young's modulus [Pa]": 15000000000.0,
 'Negative electrode cracking rate':pp["Negative electrode cracking rate"],
 'Negative electrode critical stress [Pa]': 60000000.0,
 'Negative electrode initial crack length [m]': 2e-08,
 'Negative electrode initial crack width [m]': 1.5e-08,
 'Negative electrode number of cracks per unit area [m-2]': 3180000000000000.0,
 'Negative electrode partial molar volume [m3.mol-1]': 3.1e-06,
 'Negative electrode reference concentration for free of deformation [mol.m-3]': 0.0,
 'Negative electrode volume change':pp["Negative electrode volume change"],
 'Positive electrode LAM constant exponential term': 2.0,
 'Positive electrode LAM constant proportional term [s-1]': 2.7778e-07,
 'Positive electrode OCP entropic change [V.K-1]': 0.0,
 "Positive electrode Paris' law constant b": 1.12,
 "Positive electrode Paris' law constant m": 2.2,
 "Positive electrode Poisson's ratio": 0.2,
 "Positive electrode Young's modulus [Pa]": 375000000000.0,
 'Positive electrode active material volume fraction': 0.665,
 'Positive electrode cracking rate':pp["Positive electrode cracking rate"],
 'Positive electrode critical stress [Pa]': 375000000.0,
 'Positive electrode initial crack length [m]': 2e-08,
 'Positive electrode initial crack width [m]': 1.5e-08,
 'Positive electrode number of cracks per unit area [m-2]': 3180000000000000.0,
 'Positive electrode partial molar volume [m3.mol-1]': 1.25e-05,
 'Positive electrode reference concentration for free of deformation [mol.m-3]': 0.0,
 'Positive electrode volume change':pp["Positive electrode volume change"],
 'Typical plated lithium concentration [mol.m-3]': 1000.0,
 'Exchange-current density for stripping [A.m-2]':pp["Exchange-current density for stripping [A.m-2]"],
 'Exchange-current density for plating [A.m-2]':pp["Exchange-current density for plating [A.m-2]"],
      
    }
 return degradation_parameters

def update_parameters():
    # Load Chen2020 parameters
    chen2020_params = get_chen2020_parameters()

    # Load OKane2022 degradation parameters
    okane2022_degradation_params = get_okane2022_degradation_parameters()

    # Update Chen2020 parameters with OKane2022 degradation parameters
    chen2020_params.update(okane2022_degradation_params, check_already_exists=False)

    return chen2020_params

updated_parameters = update_parameters()

updated_parameters["SEI kinetic rate constant [m.s-1]"] = 1e-12*10


cycle_number = 5
exp = pybamm.Experiment(
    [
        (
            "Discharge at 1C until 2.5 V",  # ageing cycles
            "Charge at 0.6C until 4.2 V (5 minute period)",
            "Hold at 4.2 V until C/100 (5 minute period)",
        )
    ]
    * cycle_number
)

sim = pybamm.Simulation(model, parameter_values=updated_parameters,experiment=exp,solver=pybamm.CasadiSolver(mode='safe without grid'),var_pts=var_pts)

solution = sim.solve()

#EIS Running 

start_time = timer.time()

# Set initial conditions of aged_model from the solution
aged_model = model.set_initial_conditions_from(solution)

#frequency range
frequencies = np.logspace(-4, 2, 60)

#amplitude of the current for EIS 
I = 1 * 1e-3

number_of_periods = 20
samples_per_period = 20

#Current function 

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

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

sim_aged = pybamm.Simulation(aged_model, parameter_values=updated_parameters,solver=pybamm.CasadiSolver(mode='safe without grid'),var_pts=var_pts)


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_aged.solve(t_eval, inputs={"Frequency [Hz]": frequency},initial_soc=0.9)
    # 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")

# Extract real and imaginary parts from the impedance data
real_part = np.real(impedances_time)
imaginary_part = np.imag(impedances_time)

# Plot the Nyquist plot
plt.figure(figsize=(8, 8))
plt.plot(real_part, -imaginary_part, marker='o', linestyle='-')
plt.title('Electrochemical Impedance Spectroscopy (EIS) - Nyquist Plot')
plt.xlabel('Real(Z) ohms')
plt.ylabel('-Imaginary(Z) ohms')


plt.grid(True)
plt.show()

[Akila1993]

@brosaplanella Please also note that with 10 times higher SEI creation rate than the default value the EIS has shifted to right side showing a real impedance increase, However with the number of cycle it does not change as indicated in previous comment.

[Akila1993]

Hi Everyone,

As per this discussion the cycle number affects the particle cracking but not SEI Layer growth and Li Plating model.

Accordingly, can we expect no changes in EIS for different cycle numbers when only SEI layer growth or the Li plating model is activated?

Please help with your comments or opinions.

Thanks.

[rtimms]

Degradation due to SEI growth is dominated by calendar ageing, not cycling, which is why you aren't seeing much difference with cycles but you are with time.