In many studies, voltage-time plots obtained from P2D model simulations after CCCV charging show the occurrence of a voltage plateau during the relaxation phase, and a characteristic peak can be observed after taking the second derivative. This feature is used to determine whether lithium plating (Li plating) has occurred. However, I am unable to simulate this phenomenon in PyBaMM. Below is an example.
Figure 2 from Brodsky Ringler et al. (2023) illustrates the voltage plateau observed during the relaxation phase and the characteristic peak after taking the second derivative. This image shows the lithium plating and stripping dynamics during fast charging, and is a key reference for understanding the behavior of Li plating under such conditions (Brodsky Ringler, P., Wise, M., Ramesh, P., Kim, J. H., Canova, M., Bae, C., β¦ & Park, H. (2023). Modeling of lithium plating and stripping dynamics during fast charging. Batteries, 9(7), 337).
This is my result, including the voltage-time and voltage-time derivative plots:
This is my code:
import pybamm
import matplotlib.pyplot as plt
import pandas as pdvar_pts = {
"x_n": 30, # negative electrode
"x_s": 30, # separator
"x_p": 30, # positive electrode
"r_n": 50, # negative particle
"r_p": 50, # positive particle
}
model = pybamm.lithium_ion.DFN(
{
"SEI": "solvent-diffusion limited",
"SEI porosity change": "true",
"calculate discharge energy": "true",
"lithium plating": "reversible",
"lithium plating porosity change": "true",
}
)
solver = pybamm.CasadiSolver(dt_max=5)
cycle_number = 1
param = pybamm.ParameterValues("OKane2022")
param["Ambient temperature [K]"] = 263.15
param["Lithium plating kinetic rate constant [m.s-1]"] = 5.5E-9
param["Positive electrode charge transfer coefficient"] = 0.40
param["Negative electrode charge transfer coefficient"] = 0.60
exp2 = pybamm.Experiment(
[
"Discharge at C/20 until 2.5V (10 minutes period)",
"Charge at 1C until 4.2 V (50 second period)",
"Hold at 4.2 V until C/20 (50 second period)"
] * cycle_number
+
["Rest for 720 minutes (50 second period)"])
exp3 = pybamm.Experiment(
[
"Discharge at C/20 until 2.5V (10 minutes period)",
"Charge at C/2 until 4.2 V (50 second period)",
"Hold at 4.2 V until C/20 (50 second period)"
] * cycle_number
+
["Rest for 720 minutes (50 second period)"])
exp4 = pybamm.Experiment(
[
"Discharge at C/20 until 2.5V (10 minutes period)",
"Charge at C/4 until 4.2 V (50 second period)",
"Hold at 4.2 V until C/20 (50 second period)"
] * cycle_number
+
["Rest for 720 minutes (50 second period)"])
exp5 = pybamm.Experiment(
[
"Discharge at C/20 until 2.5V (10 minutes period)",
"Charge at C/8 until 4.2 V (50 second period)",
"Hold at 4.2 V until C/20 (50 second period)"
] * cycle_number
+
["Rest for 720 minutes (50 second period)"])sim1 = pybamm.Simulation(model, parameter_values=param, experiment=exp1, var_pts=var_pts, solver=solver)
sim2 = pybamm.Simulation(model, parameter_values=param, experiment=exp2, var_pts=var_pts, solver=solver)
sim3 = pybamm.Simulation(model, parameter_values=param, experiment=exp3, var_pts=var_pts, solver=solver)
sim4 = pybamm.Simulation(model, parameter_values=param, experiment=exp4, var_pts=var_pts, solver=solver)
sim5 = pybamm.Simulation(model, parameter_values=param, experiment=exp5, var_pts=var_pts, solver=solver)
sol2 = sim2.solve()
sol3 = sim3.solve()
sol4 = sim4.solve()
sol5 = sim5.solve()V1 = sol1["Voltage [V]"].entries
V2 = sol2["Voltage [V]"].entries
V3 = sol3["Voltage [V]"].entries
V4 = sol4["Voltage [V]"].entries
V5 = sol5["Voltage [V]"].entries
T2 = sol2["Time [min]"].entries
T3 = sol3["Time [min]"].entries
T4 = sol4["Time [min]"].entries
T5 = sol5["Time [min]"].entries
T2_zero = T2[360:1225] - T2[360]
T3_zero = T3[382:1248] - T3[382]
T4_zero = T4[469:1335] - T4[469]
T5_zero = T5[696:1562] - T5[696]
plt.figure(figsize=(10, 6))
plt.plot(T2_zero, V2[360:1225], "-r", markersize=5, label="1C")
plt.plot(T3_zero, V3[382:1248], "-b", markersize=5, label="C/2")
plt.plot(T4_zero, V4[469:1335], "-y", markersize=5, label="C/4")
plt.plot(T5_zero, V5[696:1562], "-k", markersize=5, label="C/8")
plt.xlabel("Time [min]", fontsize=16)
plt.ylabel("Voltage [V]", fontsize=16)
plt.title("Voltage [V] vs Time [min] reversible", fontsize=18)
plt.tick_params(axis='y', labelsize=14)
plt.tick_params(axis='x', labelsize=14)
plt.legend(loc="upper right", fontsize=12)
plt.ylim(4.11,4.18)
plt.xlim(-0,35)
plt.grid(True)
plt.show()import numpy as np
import matplotlib.pyplot as plt
def calculate_dvdt(time, voltage):
dv = np.diff(voltage)
dt = np.diff(time)
dvdt = dv / dt
# θΏεδΈιι»ηζιοΌιΏε
ι·εΊ¦δΈδΈθ΄
time_mid = (time[:-1] + time[1:]) / 2
return time_mid, dvdt
T2_mid, dvdt_T2 = calculate_dvdt(T2_zero, V2[360:1225])
T3_mid, dvdt_T3 = calculate_dvdt(T3_zero, V3[382:1248])
T4_mid, dvdt_T4 = calculate_dvdt(T4_zero, V4[469:1335])
T5_mid, dvdt_T5 = calculate_dvdt(T5_zero, V5[696:1562])
plt.figure(figsize=(10, 6))
plt.plot(T2_mid, dvdt_T2, "-r", label="1C")
plt.plot(T3_mid, dvdt_T3, "-b", label="C/2")
plt.plot(T4_mid, dvdt_T4, "-y", label="C/4")
plt.plot(T5_mid, dvdt_T5, "-k", label="C/8")
plt.xlabel("Time [min]", fontsize=16)
plt.ylabel("dV/dt [V/min]", fontsize=16)
plt.title("dV/dt vs Time [min]", fontsize=18)
plt.tick_params(axis='y', labelsize=14)
plt.tick_params(axis='x', labelsize=14)
plt.legend(loc="upper right", fontsize=12)
plt.ylim(-0.005, 0.0)
plt.xlim(1,35)
plt.grid(True)
plt.show()```
Iβm not sure if there are any programming details that I might have missed. If possible, please help me simulate the correct feature. Thank you!