I cannot simulate the voltage plateau during relaxation caused by Li plating

[Travis56498]

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!

[github-actions[bot]]

Hi! We have now moved our discussions to Discourse. Please post your question there.