Inquiry About Parameters Influencing Negative Electrode Thickness Change in PyBaMM

[Yull-1]

PyBaMM Version

24.11.1

Python Version

3.11.9

Describe the bug

Hi PyBaMM Team,

I am currently using PyBaMM to simulate a battery with silicon-carbon as the negative electrode material. In my experiment, when the battery charges from 0% SOC to 100% SOC, the thickness of the negative electrode expands by approximately 60 μm. However, in the simulation, the negative electrode thickness change is only around 1.4 μm.

Could you advise which parameters primarily influence the "Negative electrode thickness change [m]" in the model? Additionally, how can I adjust the model to better reflect the experimental observation of significant electrode expansion?

Thank you for your support and guidance

Steps to Reproduce

`

import pybamm
import matplotlib.pyplot as plt
import os
import importlib.util
import sys
import numpy as np

pybamm.set_logging_level('NOTICE')
model = pybamm.lithium_ion.DFN(
{
"SEI": "solvent-diffusion limited",
"SEI porosity change": "true",
"lithium plating": "partially reversible",
"lithium plating porosity change": "false", # 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": "false",
}
)

parameter_values = pybamm.ParameterValues('OKane2022')

path, _ = os.path.split(os.path.abspath(file))
graphite_LGM50_ocp_Chen2020_data = pybamm.parameters.process_1D_data(
"my40_soc_ocp.csv", path=path
)

def graphite_LGM50_ocp_Chen2020(sto):
name, (x, y) = graphite_LGM50_ocp_Chen2020_data
return pybamm.Interpolant(x, y, sto, name=name, interpolator="cubic")

path, _ = os.path.split(os.path.abspath(file))
my_ocp_data = pybamm.parameters.process_1D_data("my_ocp_function_P.csv", path=path)

def my_ocp_function_P(sto):
name, (x, y) = my_ocp_data
#print(name,(x,y))
return pybamm.Interpolant(x, y, sto, name=name, interpolator="cubic")
parameter_values["Negative electrode OCP [V]"] = graphite_LGM50_ocp_Chen2020
parameter_values['Positive electrode OCP [V]'] = my_ocp_function_P
print(parameter_values['Positive electrode OCP [V]'], parameter_values["Negative electrode OCP [V]"])

def negative_exchange_current_density(c_e, c_s_surf, c_s_max, T):
m_ref = 2.7e-7
E_r = 35000
arrhenius = np.exp(E_r / pybamm.constants.R * (1 / 298.15 - 1 / T))
return m_ref * arrhenius * c_e0.5 * c_s_surf0.5 * (c_s_max - c_s_surf) ** 0.5

def negative_diffusivity(sto, T):
D_ref = 2.6597e-15
E_D_s = 3.03e4
arrhenius = np.exp(E_D_s / pybamm.constants.R * (1 / 298.15 - 1 / T))
return D_ref * arrhenius

parameter_values.update(
{

"Negative electrode initial crack width [m]": 1.5e-08,
'Negative electrode number of cracks per unit area [m-2]':2e11,
'Negative electrode initial crack length [m]':2e-6,
"Negative electrode Poisson's ratio": 0.28,
"Negative electrode critical stress [Pa]": 3.24e8,#3.24e8,
"Negative electrode LAM constant exponential term": 4,
"Negative electrode LAM constant proportional term [s-1]": 2.7778e-07,
"Negative electrode Paris' law constant b": 0.9,
"Negative electrode Paris' law constant m": 2.15,#2.15,
"Negative electrode Young's modulus [Pa]": 3e10,
"Negative electrode partial molar volume [m3.mol-1]": 6.6e-05,

'SEI resistivity [Ohm.m]': 200000.0,
'SEI growth activation energy [J.mol-1]': 34000.0,
"Initial outer SEI thickness [m]": 3e-7,
"Initial inner SEI thickness [m]": 0,
"Outer SEI partial molar volume [m3.mol-1]": 9.585e-06,
"Inner SEI partial molar volume [m3.mol-1]": 9.585e-06,
"SEI reaction exchange current density [A.m-2]": 1.5e-07,
"Outer SEI solvent diffusivity [m2.s-1]": 2.1000000000000002e-20,
"Inner SEI lithium interstitial diffusivity [m2.s-1]": 1e-19,
"SEI kinetic rate constant [m.s-1]": 1e-12,


"Lithium plating transfer coefficient": 0.65,
"Lithium plating kinetic rate constant [m.s-1]": 1e-09,
"Dead lithium decay constant [s-1]": 1e-06,


'Electrode height [m]': 0.0775,
'Electrode width [m]': 2.97,
'Positive electrode thickness [m]': 6.085e-05,
'Maximum concentration in negative electrode [mol.m-3]': 80644.074,
'Maximum concentration in positive electrode [mol.m-3]': 49765,#49765,
'Separator thickness [m]': 1.5e-05,
'Nominal cell capacity [A.h]': 9.73,
'Current function [A]': 8.5114,
'Initial concentration in positive electrode [mol.m-3]': 14750,#14750,
'Initial concentration in negative electrode [mol.m-3]': 64515.26,

'Positive electrode conductivity [S.m-1]': 3.163,
'Positive electrode porosity': 0.286,
'Positive electrode active material volume fraction': 0.7137,
"Positive particle radius [m]": 5.22e-06,
"Positive electrode Young's modulus [Pa]": 375000000000.0,

'Negative electrode porosity': 0.4543,
'Negative electrode thickness [m]': 3.72389e-05,
'Negative electrode conductivity [S.m-1]': 71.9489,
'Negative electrode active material volume fraction': 0.5457,
"Negative electrode exchange-current density [A.m-2]": negative_exchange_current_density,#
'Separator porosity': 0.38,
"Negative particle diffusivity [m2.s-1]":negative_diffusivity,
"Negative electrode Bruggeman coefficient (electrolyte)": 1.5,
"Negative electrode Bruggeman coefficient (electrode)": 1.5,
 }
)

c_n_init = parameter_values["Initial concentration in negative electrode [mol.m-3]"]
c_n_max = parameter_values["Maximum concentration in negative electrode [mol.m-3]"]

soc = c_n_init / c_n_max
print(f"SOC: {soc}")

var_pts = {
"x_n":10,
'x_s':5,
"x_p":10,
"r_n":30,
"r_p":30
}

exp = pybamm.Experiment(
[(
'Rest for 10 minutes',
'Discharge at 8.5613A until 2.5V',
'Rest for 10 minutes',
'Charge at 8.56131A until 4.2V',
'Hold at 4.2V until 0.05C'
)]*10)

safe_solver = pybamm.CasadiSolver(atol=1e-3, rtol=1e-3, mode="safe")

sim = pybamm.Simulation(model, experiment = exp, parameter_values= parameter_values, var_pts = var_pts, solver = safe_solver)
sim.solve()
#sim.plot()
sol = sim.solution

sol.save('40%1CCCV.pkl')

`

• plotting code
  `

import pybamm
import matplotlib.pyplot as plt

sol = pybamm.load('40%1CCCV.pkl')

#print(sol.summary_variables['Negative total SEI thickness [m]'].shape)
#sol.plot(['Negative total SEI thickness [m]'])
sei_thick = sol['Negative total SEI thickness [m]'].entries[5]
sei_thick_1 = sol['Negative total SEI thickness [m]'].entries[0]
sei_concentration = sol['Negative SEI concentration [mol.m-3]'].entries[5]
sei_loss = sol['Loss of lithium to negative SEI [mol]'].entries
neg_thick_change_1 = sol.cycles[1]['Negative electrode thickness change [m]'].entries
#print(neg_thick_change, neg_thick_change.shape)

print(sei_concentration, sei_concentration.shape)
print(sei_thick*1e9, sei_thick.shape)

t = sol['Time [h]'].entries
print(t, t.shape)
V = sol.cycles[1]['Terminal voltage [V]'].entries
t_1 = sol.cycles[1]['Time [s]'].entries

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(18, 4))
ax1.plot(t, sei_thick*1e9, label="sei")
ax1.set_xlabel("Time [h]")
ax1.set_ylabel("sei")
ax1.legend()
ax2.plot(t, sei_concentration, label="sei")
ax2.set_xlabel("Time [h]")
ax2.set_ylabel("sei")
ax2.legend()

ax3.plot(t_1, neg_thick_change_1 * -1e6, label="Thickness Change", color="tab:blue")
ax3.set_xlabel("Time [s]")
ax3.set_ylabel("Negative Electrode Thickness Change [μm]", color="tab:blue")
ax3.tick_params(axis='y', labelcolor="tab:blue")
ax3.legend(loc="upper left")

ax4 = ax3.twinx()
ax4.plot(t_1, V, label="Terminal Voltage", color="tab:red")
ax4.set_ylabel("Terminal Voltage [V]", color="tab:red")
ax4.tick_params(axis='y', labelcolor="tab:red")
ax4.legend(loc="upper right")

plt.tight_layout()
plt.show()

pybamm.plot_summary_variables(sol)
`

-result

Relevant log output