Issue with EIS Implementation in PyBaMM for LGM50T 21700 Cell (Chen2020)

[mlnsai]

Hi Team,

I'm currently working on implementing EIS (Electrochemical Impedance Spectroscopy) in PyBaMM for a Chenn2020 Parameters-based LGM50T 21700 cell.

Below is the code snippet I've used:

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

# Set up
model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}, name="SPM")
parameter_values = pybamm.ParameterValues("Chen2020")
frequencies = np.logspace(-4, 3, 30)

# EIS parameters
I = 1 * 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)


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

sim = pybamm.Simulation(
    model, parameter_values=parameter_values, 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)

# Nyquist Plot
_, ax = plt.subplots()
ax = pbeis.nyquist_plot(impedances_time, ax=ax, linestyle='-', alpha=0.7)

ax.legend()
plt.suptitle(f"{model.name}")
plt.show()


# Convert impedance data to arrays for real and imaginary parts
impedance_array = np.array(impedances_time)

# Calculate the magnitude (resultant) of the impedance
resultant_impedance = np.abs(impedance_array)


# Convert impedance to milliohms
impedance_milliohms = resultant_impedance * 1000  # Ohms to milliohms
# Print impedance values in milliohms at specific frequencies
target_frequencies = [0.01, 1000]
for freq in target_frequencies:
    index = np.argmin(np.abs(frequencies - freq))
    impedance_at_freq_milliohms = impedance_milliohms[index]
    print(f'Impedance at {freq} Hz: {impedance_at_freq_milliohms:.2f} mΩ')


# Create the line plot
plt.figure()
plt.plot(frequencies, impedance_milliohms, linestyle='-')  # Line plot
plt.xlabel('Frequency (Hz)')
plt.ylabel('Impedance (mΩ)')

# Set y-axis to have at least 5-6 major points
ax = plt.gca()  # Get current axis
ax.yaxis.set_major_locator(ticker.MaxNLocator(6))

plt.title('Impedance vs Frequencies')
plt.xscale('log')  # Log scale for frequencies
plt.grid(True)
plt.show()

I've attempted to extract impedances at various frequencies and plotted the results. However, at 1kHz frequency, the impedance I'm obtaining is not aligning with the expected value of around 25mOhm, as per the datasheet

Here are the corresponding Nyquist and Impedance Plots.

Any guidance on additional steps I might need to take or checks I should perform would be greatly appreciated.
@brosaplanella @rtimms @tinosulzer

Thank you!

[mlnsai]

Hey team!! A gentle reminder. Any leads would be greatly appreciated :)

I still couldn't figure out the issue.

@brosaplanella @rtimms @tinosulzer

[mlnsai]

Hi @Akila1993 , thanks for checking in! I did make some progress. Using DFN and SPMe models didn't show significant variation, but tweaking the Contact resistance did. However, determining the right value remains a challenge. Currently, I'm using the initial contact resistance to match the expected 25mOhm at 1 kHz. Open to any suggestions!

Before:
"Contact resistance [Ohm]": 0,
Impedance at 0.01 Hz: 54.84 mΩ
Impedance at 1000 Hz: 2.48 mΩ

After:
"Contact resistance [Ohm]": 0.023,
Impedance at 0.01 Hz: 77.31 mΩ
Impedance at 1000 Hz: 25.47 mΩ

[Akila1993]

Hi, @mlnsai

Thank you for your reply,

However, I tried to change the "Contact resistance [Ohm]" in the Chen2020 parameter set but in the list of parameters, it is not listed. @valentinsulzer is it due to the version of my PyBaMM or anything else? When I am running it with Google colab it shows me the Contact resistance [Ohm] parameter.

Please advice.

My snipped code is below.

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

# Set up
model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}, name="SPM")
parameter_values = pybamm.ParameterValues("Chen2020")

parameter_values 

{'1 + dlnf/dlnc': 1.0,
 'Ambient temperature [K]': 298.15,
 'Bulk solvent concentration [mol.m-3]': 2636.0,
 'Cation transference number': 0.2594,
 'Cell cooling surface area [m2]': 0.00531,
 'Cell thermal expansion coefficient [m.K-1]': 1.1e-06,
 'Cell volume [m3]': 2.42e-05,
 'Current function [A]': 5.0,
 'EC diffusivity [m2.s-1]': 2e-18,
 'EC initial concentration in electrolyte [mol.m-3]': 4541.0,
 'Electrode height [m]': 0.065,
 'Electrode width [m]': 1.58,
 'Electrolyte conductivity [S.m-1]': <function electrolyte_conductivity_Nyman2008 at 0x14bf8b310>,
 'Electrolyte diffusivity [m2.s-1]': <function electrolyte_diffusivity_Nyman2008 at 0x14bf8b280>,
 'Initial concentration in electrolyte [mol.m-3]': 1000.0,
 'Initial concentration in negative electrode [mol.m-3]': 29866.0,
 'Initial concentration in positive electrode [mol.m-3]': 17038.0,
 'Initial inner SEI thickness [m]': 2.5e-09,
 'Initial outer SEI thickness [m]': 2.5e-09,
 'Initial temperature [K]': 298.15,
 'Inner SEI electron conductivity [S.m-1]': 8.95e-14,
 'Inner SEI lithium interstitial diffusivity [m2.s-1]': 1e-20,
 'Inner SEI open-circuit potential [V]': 0.1,
 'Inner SEI partial molar volume [m3.mol-1]': 9.585e-05,
 'Inner SEI reaction proportion': 0.5,
 'Lithium interstitial reference concentration [mol.m-3]': 15.0,
 'Lower voltage cut-off [V]': 2.5,
 'Maximum concentration in negative electrode [mol.m-3]': 33133.0,
 'Maximum concentration in positive electrode [mol.m-3]': 63104.0,
 'Negative current collector conductivity [S.m-1]': 58411000.0,
 'Negative current collector density [kg.m-3]': 8960.0,
 'Negative current collector specific heat capacity [J.kg-1.K-1]': 385.0,
 'Negative current collector thermal conductivity [W.m-1.K-1]': 401.0,
 'Negative current collector thickness [m]': 1.2e-05,
 'Negative electrode Bruggeman coefficient (electrode)': 1.5,
 'Negative electrode Bruggeman coefficient (electrolyte)': 1.5,
 'Negative electrode OCP [V]': <function graphite_LGM50_ocp_Chen2020 at 0x14b7b5a60>,
 'Negative electrode OCP entropic change [V.K-1]': 0.0,
 'Negative electrode active material volume fraction': 0.75,
 'Negative electrode cation signed stoichiometry': -1.0,
 'Negative electrode charge transfer coefficient': 0.5,
 'Negative electrode conductivity [S.m-1]': 215.0,
 'Negative electrode density [kg.m-3]': 1657.0,
 'Negative electrode diffusivity [m2.s-1]': 3.3e-14,
 'Negative electrode double-layer capacity [F.m-2]': 0.2,
 'Negative electrode electrons in reaction': 1.0,
 'Negative electrode exchange-current density [A.m-2]': <function graphite_LGM50_electrolyte_exchange_current_density_Chen2020 at 0x14bf8b040>,
 'Negative electrode porosity': 0.25,
 'Negative electrode reaction-driven LAM factor [m3.mol-1]': 0.0,
 'Negative electrode specific heat capacity [J.kg-1.K-1]': 700.0,
 'Negative electrode thermal conductivity [W.m-1.K-1]': 1.7,
 'Negative electrode thickness [m]': 8.52e-05,
 'Negative particle radius [m]': 5.86e-06,
 'Nominal cell capacity [A.h]': 5.0,
 'Number of cells connected in series to make a battery': 1.0,
 'Number of electrodes connected in parallel to make a cell': 1.0,
 'Outer SEI open-circuit potential [V]': 0.8,
 'Outer SEI partial molar volume [m3.mol-1]': 9.585e-05,
 'Outer SEI solvent diffusivity [m2.s-1]': 2.5000000000000002e-22,
 'Positive current collector conductivity [S.m-1]': 36914000.0,
 'Positive current collector density [kg.m-3]': 2700.0,
 'Positive current collector specific heat capacity [J.kg-1.K-1]': 897.0,
 'Positive current collector thermal conductivity [W.m-1.K-1]': 237.0,
 'Positive current collector thickness [m]': 1.6e-05,
 'Positive electrode Bruggeman coefficient (electrode)': 1.5,
 'Positive electrode Bruggeman coefficient (electrolyte)': 1.5,
 'Positive electrode OCP [V]': <function nmc_LGM50_ocp_Chen2020 at 0x14bf8b160>,
 'Positive electrode OCP entropic change [V.K-1]': 0.0,
 'Positive electrode active material volume fraction': 0.665,
 'Positive electrode cation signed stoichiometry': -1.0,
 'Positive electrode charge transfer coefficient': 0.5,
 'Positive electrode conductivity [S.m-1]': 0.18,
 'Positive electrode density [kg.m-3]': 3262.0,
 'Positive electrode diffusivity [m2.s-1]': 4e-15,
 'Positive electrode double-layer capacity [F.m-2]': 0.2,
 'Positive electrode electrons in reaction': 1.0,
 'Positive electrode exchange-current density [A.m-2]': <function nmc_LGM50_electrolyte_exchange_current_density_Chen2020 at 0x14bf8b1f0>,
 'Positive electrode porosity': 0.335,
 'Positive electrode reaction-driven LAM factor [m3.mol-1]': 0.0,
 'Positive electrode specific heat capacity [J.kg-1.K-1]': 700.0,
 'Positive electrode thermal conductivity [W.m-1.K-1]': 2.1,
 'Positive electrode thickness [m]': 7.56e-05,
 'Positive particle radius [m]': 5.22e-06,
 'Ratio of lithium moles to SEI moles': 2.0,
 'Reference temperature [K]': 298.15,
 'SEI growth activation energy [J.mol-1]': 0.0,
 'SEI kinetic rate constant [m.s-1]': 1e-12,
 'SEI open-circuit potential [V]': 0.4,
 'SEI reaction exchange current density [A.m-2]': 1.5e-07,
 'SEI resistivity [Ohm.m]': 200000.0,
 'Separator Bruggeman coefficient (electrolyte)': 1.5,
 'Separator density [kg.m-3]': 397.0,
 'Separator porosity': 0.47,
 'Separator specific heat capacity [J.kg-1.K-1]': 700.0,
 'Separator thermal conductivity [W.m-1.K-1]': 0.16,
 'Separator thickness [m]': 1.2e-05,
 'Total heat transfer coefficient [W.m-2.K-1]': 10.0,
 'Typical current [A]': 5.0,
 'Typical electrolyte concentration [mol.m-3]': 1000.0,
 'Upper voltage cut-off [V]': 4.2,
 'citations': ['Chen2020']}

[mlnsai]

Hi @Akila1993 , yes, it has do with the version of pybamm. This parameter is added as input variable in recent versions. I had to upgrade my version from 22.3 to 24.1 (recent) to tune this parameter. Hope this helps.

[Akila1993]

Hi @mlnsai,

I tried to change the contact resistance and the EIS remains the same even in PyBaMM EIS. It is possible to share your code which has obtained the following results.

Before:
"Contact resistance [Ohm]": 0,
Impedance at 0.01 Hz: 54.84 mΩ
Impedance at 1000 Hz: 2.48 mΩ

After:
"Contact resistance [Ohm]": 0.023,
Impedance at 0.01 Hz: 77.31 mΩ
Impedance at 1000 Hz: 25.47 mΩ

Thank you for your consideration.

[rtimms]

You need to set the option “contact resistance” to “true” to see this effect.

[alirahdarian]

Hello,
I am using the same code as you but I do not reach the same result. The problem is that I do not have the semicircle in my plot. Can you please let me know what is the problem of my code?
import pybamm
import numpy as np
import matplotlib.pyplot as plt
import time as timer
from scipy.fft import fft
import matplotlib.ticker as ticker

#Set up
model = pybamm.lithium_ion.SPM(options={"surface form": "differential"}, name="SPM")
parameter_values = pybamm.ParameterValues("Chen2020")
frequencies = np.logspace(-4, 3, 100)

#EIS parameters
I = 1 * 1e-3
number_of_periods = 20
samples_per_period = 100

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

parameter_values["Current function [A]"] = current_function
#parameter_values["Contact resistance [Ohm]"] = 0.023
sim = pybamm.Simulation(
model, parameter_values=parameter_values, 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[-5 * samples_per_period - 1 :]
#print('time=',sol["Time [s]"])
#print('stable time=',time)
current = sol["Current [A]"].entries[-5 * samples_per_period - 1 :]
voltage = sol["Terminal voltage [V]"].entries[-5 * 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)

#Nyquist Plot
print('impedances:', impedances_time)
data = np.array(impedances_time)
print('data after array=', data)
plt.plot(data.real, -data.imag,linestyle='None' , marker='o')
plt.xlabel(r"$Z_\mathrm{Re}$ [Ohm]")
plt.ylabel(r"$-Z_\mathrm{Im}$ [Ohm]")
plt.show()

sim.plot(output_variables=["Current [A]",
"Terminal voltage [V]", 'Battery open-circuit voltage [V]'])

[rtimms]

This is a bug. See #4139

[alirahdarian]

How can we specify SOC (state of charge) in the code of EIS?
can we do it as follows:
sol = sim.solve(t_eval, inputs={"Frequency [Hz]": frequency}, initial_soc=0.1)

[Akila1993]

Hi @alirahdarian ,

Please see this discussion: #3572 (comment).