Can we add charge code in JellyBaMM code ?

[PurvaWarke28]

I am trying to learn JellyBaMM. I saw the paper " Probing Heterogeneity in Li-Ion Batteries with Coupled Multiscale Models of Electrochemistry and Thermal Transport using Tomographic Domains" and I'm trying to replicate the results from this paper with the help of the code you provided as an example @TomTranter

I want to know if I can add PybaMM different submodels in the file simulations.py for capacity fade in model=pybamm.lithium_ion.SPMe

Also, can we add a Charge code to the original code?
I tried to add 10 cycles, but the simulation stopped halfway through and gave an error.
what is the significance of step=0 in liionsolve.py

I am writing the code below
plt.close("all")
wrk = op.Workspace()
wrk.clear()

Geometry of spiral

Nlayers = 19
dtheta = 10
spacing = 209e-6
inner_r = 0.002
pos_tabs = [-1]
neg_tabs = [0]
length_3d = 0.06
tesla_tabs = False
I_app = 3.5
dt = 40
Nsteps = 270
hours = dt * Nsteps / 3600

experiment = pybamm.Experiment(
[
(f"Discharge at {I_app} A for 1 hours", # f"Charge at {I_app} A for {hours} hours"
f"Charge at {I_app} A for 1 hours"),
#"Rest for 10 minutes",
# f"Charge at {I_app} A for {hours} hours or until 4.2 V",
# "Rest for 10 minutes",
],
period=f"{dt} seconds",
)

Create the spiral network and plot it

project, arc_edges = jellybamm.make_spiral_net(
Nlayers, dtheta, spacing, inner_r, pos_tabs, neg_tabs, length_3d, tesla_tabs
)

jellybamm.plot_topology(project.network)

Load parameter values and modify them

param = pybamm.ParameterValues("Chen2020")
param["Electrode width [m]"] = length_3d
param["Electrode height [m]"] = 0.06 # Casing width = 2mm
param["Maximum concentration in negative electrode [mol.m-3]"]=33429
param["Negative current collector conductivity [S.m-1]"]=6e7
param["Positive current collector conductivity [S.m-1]"]=3e7
param["Positive electrode thickness [m]"] = 7.28e-5
param["Negative electrode thickness [m]"] = 8.32e-5
param["Separator thickness [m]"] = 2.08e-5
param["Cation transference number"] = 0.4
param["Negative electrode specific heat capacity [J.kg-1.K-1]"]=1437.4
param["Positive electrode specific heat capacity [J.kg-1.K-1]"]=1269.21
param[ "Separator specific heat capacity [J.kg-1.K-1]"] = 1978.16
param["Negative current collector thermal conductivity [W.m-1.K-1]"]=398
param["Positive current collector thermal conductivity [W.m-1.K-1]"]=238
param["Negative electrode thermal conductivity [W.m-1.K-1]"]=1.04
param["Positive electrode thermal conductivity [W.m-1.K-1]"] =1.58
param["Separator thermal conductivity [W.m-1.K-1]"] = 0.344
param["Negative current collector density [kg.m-3]"] = 8933
param["Positive current collector density [kg.m-3]"] = 2702
param["Negative electrode density [kg.m-3]"] = 1347.33
param["Positive electrode density [kg.m-3]"] = 2328.5
param["Separator density [kg.m-3]"]=1008.98
param["Negative current collector thickness [m]"]=1.04e-5
param["Negative electrode thickness [m]"]=8.32e-5
param["Separator thickness [m]"] = 2.08e-5
param["Positive electrode thickness [m]"]=7.28e-5
param["Positive current collector thickness [m]"] = 1.04e-5
param["Total heat transfer coefficient [W.m-2.K-1]"]= 5 # One of the parameter used for variation
param["Nominal cell capacity [A.h]"]=3.5
param["Current function [A]"]= I_app
param["Open-circuit voltage at 0% SOC [V]"]=2.5
param["Negative electrode OCP entropic change [V.K-1]"]= 0.0004 #dU/dT
param["Positive electrode OCP entropic change [V.K-1]"] = 7e-5

Initial state of charge

initial_soc = None

Run the simulation with jellybamm

project, output = jellybamm.run_simulation_lp(
parameter_values=param,
experiment=experiment,
initial_soc=initial_soc,
project=project,
)

Extract simulation results

time_values = np.array(output['Time [s]'])
cell_current_values = np.array(output['Cell current [A]']).flatten()
terminal_voltage_values = np.array(output['Terminal voltage [V]']).flatten()
temperature = np.array(output['Volume-averaged cell temperature [K]']).flatten()
current_den = np.array(output['Current collector current density [A.m-2]']).flatten()
Vol_avg_tot_heat = np.array(output['Volume-averaged total heating [W.m-3]']).flatten()
capacity = np.cumsum(cell_current_values) * (time_values[1] - time_values[0])/3600

Filtered values

Ter_vol_fil = [x for x in terminal_voltage_values if x!=0]
current_den_fil = [x for x in current_den if x!=0]
xavg_heat_fil = [x for x in Vol_avg_tot_heat if x!=0]
temp_fil = [x for x in temperature if x!=0]

Plot Voltage vs. Capacity

plt.figure(figsize=(8, 6))
plt.plot(capacity[:len(Ter_vol_fil)],Ter_vol_fil , marker='o', linestyle='-', color='blue')

plt.plot(capacity, terminal_voltage_values[:len(capacity)], marker='o', linestyle='-', color='blue')

plt.xlabel('Capacity [Ah]')
plt.ylabel('Terminal voltage [V]')
plt.title('Terminal voltage vs Capacity')
plt.grid(True)