6 - Capturing Hysteresis of LFP in PyBaMM

Agenda: LiionDB database to extract OCP function: https://liiondb.com/ Importing OCP funtion in PyBaMM OCP update for charge and discharge Hysteresis results (open loop) Code: import pybamm dfn = pybamm.lithium_ion.DFN() exp = pybamm.Experiment([ "Discharge at C/50 until 2V", "Charge at C/50 until 3.65V", ]) param = pybamm.ParameterValues("Prada2013") sim = pybamm.Simulation(dfn, experiment=exp, parameter_values=param) sol = sim.solve() t = sol["Time [s]"].entries v = sol["Terminal voltage [V]"].entries Q = sol["Discharge capacity [A.h]"].entries import matplotlib.pyplot as plt plt.plot(Q,v) plt.xlabel("Time [s]") plt.ylabel("Voltage [V]") plt.legend() plt.show() def cOCP(x): import numpy as np FITTING: Equations taken directly from paper table 1 import numpy as np A0 = 3.4510 B0 = -8.8E-3 A1 = 0.6678 B1 = -81.002 C1 = 1.1776 A2 = 2.3738E-9 B2 = 25.222 C2 = 3.4801 A3 = -2.6367E-9 B3 = 25.12 C3 = 3.4879 U = A0 + B0*x + A1*np.exp(B1*x**C1) + A2*np.exp(B2*x**C2) + A3*np.exp(B3*x**C3) return U def dOCP(x): import numpy as np FITTING: Equations taken directly from paper table 1 import numpy as np A0 = 3.4227 B0 = -2.0269E-2 A1 = 0.5087 B1 = -81.163 C1 = 1.0138 A2 = 7.6445E-8 B2 = 25.361 C2 = 3.2983 A3 = -8.4410E-8 B3 = 25.262 C3 = 3.3111 U = A0 + B0*x + A1*np.exp(B1*x**C1) + A2*np.exp(B2*x**C2) + A3*np.exp(B3*x**C3) return U def avgOCP(x): return (cOCP(x) + dOCP(x))/2 model_hysteresis = pybamm.lithium_ion.DFN( options={"open-circuit potential": ("single", "current sigmoid")} ) param_hysteresis = pybamm.ParameterValues("Prada2013") param_hysteresis.update({ "Positive electrode OCP [V]": avgOCP, "Positive electrode lithiation OCP [V]": dOCP, "Positive electrode delithiation OCP [V]": cOCP, }) sim = pybamm.Simulation(model_hysteresis, experiment=exp, parameter_values=param_hysteresis) sol = sim.solve() t = sol["Time [s]"].entries v = sol["Terminal voltage [V]"].entries Q = sol["Discharge capacity [A.h]"].entries import matplotlib.pyplot as plt plt.plot(Q,v) plt.xlabel("Time [s]") plt.ylabel("Voltage [V]") plt.legend() plt.show()