Hello dear members, I am Dieudonné SIMPORE, a PhD student at Joseph Ki-ZERBO University. I use PyBaMM extensively for my work and find the content truly fascinating. My issue is that I am trying to implement the SPM model using a finite difference method with either an implicit or explicit scheme. I am using the Chen2020 dataset, but I am getting results that are quite different from what one would expect. For example, starting from 100% SOC and after 1 hour of discharge, the SOC does not fully decrease, and the voltage evolves in a rather incorrect manner.
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import csc_matrix
from scipy.sparse.linalg import spsolve
Cell 1: définition de la classe SPMParameters
class SPMParameters:
def init(self,
I_app=1.0, # courant appliqué [A]
total_time=100, # durée totale [s]
dt=0.5, # pas de temps [s]
Nx=50, # points spatiaux
soc_initiale = 1
):
# constantes universelles
self.F = 96485.0 # [C/mol]
self.R_gas = 8.314 # [J/(mol·K)]
self.T = 298.15 # [K]
# électrode positive (NMC)
self.A_pos = 10.27e-2
self.L_pos = 75.6e-6
self.Rp_pos = 5.22e-6
self.Ds_pos = 1.48e-15
self.Cs_max_pos = 51765.0
self.epsilon_act_pos = 0.665
self.k_pos = 3.42e-6
self.Rct_pos = 1.80
self.c_pos_min = 0.2661
self.c_pos_max = 0.9084
# électrode négative (Graphite-SiOx)
self.A_neg = 10.27e-2
self.L_neg = 85.2e-6
self.Rp_neg = 5.86e-6
self.Ds_neg = 1.74e-15
self.Cs_max_neg = 29583.0
self.epsilon_act_neg = 0.75
self.k_neg = 6.48e-7
self.Rct_neg = 14.72
self.c_neg_max = 0.9014
self.c_neg_min = 0.0279
# électrolyte commun
self.ce = 1000.0 # [mol/m³]
# paramètres numériques
self.dt = dt
self.Nx = Nx
self.total_time = total_time
self.I_app = I_app
self.soc_initiale = soc_initiale
Cell 2: définition de la classe SPMUtilities
class SPMUtilities:
@staticmethod
def compute_flux(I_app, F, A, L, epsilon_act, Rp):
a_s = 3.0 * epsilon_act / Rp
return I_app / (F * A * L * a_s) # [mol/(m²·s)]
@staticmethod
def compute_stoichiometry(Cs, Cs_max):
return Cs / Cs_max
@staticmethod
def compute_soc(c, c_min, c_max):
return (c - c_min) / (c_max - c_min)
@staticmethod
def compute_ocv(Cs, Cs_max, electrode='positive'):
c = Cs / Cs_max
if electrode.lower() == 'positive':
return (
-0.8090 * c + 4.4875
- 0.0428 * np.tanh(18.5138 * (c - 0.5542))
- 17.7326 * np.tanh(15.7890 * (c - 0.3117))
+ 17.5842 * np.tanh(15.9308 * (c - 0.3120))
)
elif electrode.lower() == 'negative':
return (
1.9739 * np.exp(-39.3631 * c) + 0.2482
- 0.0909 * np.tanh(29.8538 * (c - 0.1234))
- 0.04478 * np.tanh(14.9159 * (c - 0.2769))
- 0.0205 * np.tanh(30.4444 * (c - 0.6103))
)
else:
raise ValueError("Electrode must be 'positive' or 'negative'.")
@staticmethod
def construct_implicit_matrix(N, dt, dr, D, r, j_surface):
A = np.zeros((N, N))
for i in range(1, N - 1):
alpha = D * dt / dr**2
beta = D * dt / (r[i] * dr)
A[i, i-1] = - (alpha + beta)
A[i, i] = 1 + 2 * alpha
A[i, i+1] = - (alpha - beta)
A[0, 0], A[0, 1] = 1, -1
A[-1, :] = 0
A[-1, -2], A[-1, -1] = -1, 1
return csc_matrix(A)
@staticmethod
def compute_j0(k, ce, c_surf, Cs_max):
return k * np.sqrt(ce * c_surf * (Cs_max - c_surf ))
@staticmethod
def compute_eta(j, j0, T, F, R=8.314):
return (2 * R * T / F) * np.arcsinh(j * F/ (2 * j0))
@staticmethod
def compute_Rct_term(Rct_pos, a_s_pos, L_pos, Rct_neg, a_s_neg, L_neg):
return (Rct_pos + Rct_neg )
@staticmethod
def compute_cell_voltage(ocv_pos, ocv_neg, eta_pos, eta_neg, Rct_term, I_app):
return (ocv_pos - ocv_neg) + (eta_pos - eta_neg) - 0*(Rct_term * I_app)
##########################
CLASSE 3: SPMSolver
##########################
class SPMSolver:
def init(self, params: SPMParameters, electrode=‘positive’, method=‘implicit’):
“”"
Initialise le solveur du modèle SPM pour une électrode.
Arguments:
params : instance de SPMParameters
electrode: 'positive' ou 'negative'
method : 'implicit' ou 'explicit'
"""
self.params = params
self.electrode = electrode.lower()
self.method = method.lower()
# Sélection des paramètres selon l'électrode choisie
if self.electrode == 'positive':
self.Ds = params.Ds_pos
self.Rp = params.Rp_pos
self.Cs_max = params.Cs_max_pos
self.epsilon_act = params.epsilon_act_pos
self.A = params.A_pos
self.L = params.L_pos
self.c0 = ((1-params.soc_initiale) * (params.c_pos_max - params.c_pos_min) + params.c_pos_min )* params.Cs_max_pos
elif self.electrode == 'negative':
self.Ds = params.Ds_neg
self.Rp = params.Rp_neg
self.Cs_max = params.Cs_max_neg
self.epsilon_act = params.epsilon_act_neg
self.A = params.A_neg
self.L = params.L_neg
self.c0 = (params.soc_initiale * (params.c_neg_max - params.c_neg_min) + params.c_neg_min )* params.Cs_max_neg
else:
raise ValueError("L'électrode doit être 'positive' ou 'negative'.")
# Discrétisation spatiale de la particule
self.Nx = params.Nx
self.dr = self.Rp / (self.Nx - 1)
self.r = np.linspace(0, self.Rp, self.Nx)
self.dt = params.dt
self.num_steps = int(params.total_time / self.dt)
# Condition initiale : concentration uniforme
self.c_init = np.ones(self.Nx) * self.c0
# Calcul du flux interfacial (via SPMUtilities)
flux = SPMUtilities.compute_flux(
I_app=params.I_app,
F=params.F,
A=self.A,
L=self.L,
epsilon_act=self.epsilon_act,
Rp=self.Rp
)
# Pour l'électrode négative, le flux doit être multiplié par -1
if self.electrode == 'negative':
flux = -flux
self.j_surface = flux
# Construction de la matrice implicite
self.A_matrix = SPMUtilities.construct_implicit_matrix(
N=self.Nx,
dt=self.dt,
dr=self.dr,
D=self.Ds,
r=self.r,
j_surface=self.j_surface
)
def simulate_implicit(self):
"""
Résout le modèle par schéma implicite.
Renvoie (times, r, sol) où sol est la matrice solution (profil de concentration à chaque instant).
"""
c = self.c_init.copy()
sol = np.zeros((self.num_steps + 1, self.Nx))
sol[0, :] = c.copy()
times = np.linspace(0, self.params.total_time, self.num_steps + 1)
for n in range(1, self.num_steps + 1):
b = c.copy()
b[0] = 0.0 # symétrie au centre
b[-1] = - (self.j_surface * self.dr) / self.Ds # condition flux à la surface
c = spsolve(self.A_matrix, b)
sol[n, :] = c.copy()
return times, self.r, sol
def simulate_explicit(self):
"""
Résout le modèle par schéma explicite.
Renvoie (times, r, sol) de la même manière.
"""
c = self.c_init.copy()
sol = np.zeros((self.num_steps + 1, self.Nx))
sol[0, :] = c.copy()
times = np.linspace(0, self.params.total_time, self.num_steps + 1)
for n in range(1, self.num_steps + 1):
c_new = c.copy()
for i in range(1, self.Nx - 1):
d2c = (c[i+1] - 2 * c[i] + c[i-1]) / self.dr**2
dc = (c[i+1] - c[i-1]) / (2 * self.dr)
spherical_term = (2 * dc / self.r[i]) if self.r[i] > 0 else 0.0
c_new[i] = c[i] + self.dt * self.Ds * (d2c + spherical_term)
# Centre (r=0)
c_new[0] = c[0] + self.dt * self.Ds * (2 * (c[1] - c[0]) / self.dr**2)
# Surface (r=Rp)
c_new[-1] = c[-2] - (self.j_surface * self.dr) / self.Ds
c = c_new.copy()
sol[n, :] = c.copy()
return times, self.r, sol
def simulate(self):
"""
Sélectionne la méthode de résolution selon self.method.
"""
if self.method == 'implicit':
return self.simulate_implicit()
elif self.method == 'explicit':
return self.simulate_explicit()
else:
raise ValueError("Méthode inconnue. Choisissez 'implicit' ou 'explicit'.")
##########################
CLASSE 4: SPMVisualization
##########################
class SPMVisualization:
def init(self, solver_pos: SPMSolver, solver_neg: SPMSolver, params: SPMParameters):
“”"
Initialise la visualisation pour les deux électrodes.
Arguments:
solver_pos: instance de SPMSolver pour l'électrode positive.
solver_neg: instance de SPMSolver pour l'électrode négative.
params : instance de SPMParameters.
"""
self.solver_pos = solver_pos
self.solver_neg = solver_neg
self.params = params
def plot_concentration_evolution(self):
"""
Affiche l'évolution des concentrations de surface pour chaque électrode
sur deux axes y distincts (similaire à Matlab yyaxis).
"""
times_pos, _, sol_pos = self.solver_pos.simulate()
times_neg, _, sol_neg = self.solver_neg.simulate()
Cs_pos = sol_pos[:, -1] # concentration de surface positive
Cs_neg = sol_neg[:, -1] # concentration de surface négative
fig, ax1 = plt.subplots(figsize=(8, 5))
ax1.plot(times_pos, Cs_pos, 'r-', label="Positive")
ax1.set_xlabel("Temps (s)")
ax1.set_ylabel("Concentration positive (mol/m³)", color='r')
ax1.tick_params(axis='y', labelcolor='r')
ax2 = ax1.twinx()
ax2.plot(times_neg, Cs_neg, 'b-', label="Negative")
ax2.set_ylabel("Concentration négative (mol/m³)", color='b')
ax2.tick_params(axis='y', labelcolor='b')
plt.title("Évolution des concentrations de surface")
fig.tight_layout()
def plot_stoichiometry(self):
"""
Affiche l'évolution de la stœchiométrie pour chaque électrode sur deux axes y distincts.
"""
times_pos, _, sol_pos = self.solver_pos.simulate()
times_neg, _, sol_neg = self.solver_neg.simulate()
stoich_pos = SPMUtilities.compute_stoichiometry(sol_pos[:, -1], self.params.Cs_max_pos)
stoich_neg = SPMUtilities.compute_stoichiometry(sol_neg[:, -1], self.params.Cs_max_neg)
fig, ax1 = plt.subplots(figsize=(8, 5))
ax1.plot(times_pos, stoich_pos, 'r-', label="Positive")
ax1.set_xlabel("Temps (s)")
ax1.set_ylabel("Stœchiométrie positive", color='r')
ax1.tick_params(axis='y', labelcolor='r')
ax2 = ax1.twinx()
ax2.plot(times_neg, stoich_neg, 'b-', label="Negative")
ax2.set_ylabel("Stœchiométrie négative", color='b')
ax2.tick_params(axis='y', labelcolor='b')
plt.title("Évolution de la stœchiométrie")
fig.tight_layout()
plt.show()
def plot_soc(self):
"""
Affiche l'évolution de la stœchiométrie pour chaque électrode sur deux axes y distincts.
"""
times_neg, _, sol_neg = self.solver_neg.simulate()
stoich_neg = SPMUtilities.compute_stoichiometry(sol_neg[:, -1], self.params.Cs_max_neg)
soc = SPMUtilities.compute_soc(stoich_neg, self.params.c_neg_min, self.params.c_neg_max)
plt.figure(figsize=(8,5))
plt.plot(times_neg, soc * 100, color='green')
plt.ylabel('SOC (%)')
plt.title('State of charge')
plt.tight_layout()
plt.grid()
plt.show()
def plot_ocP(self):
"""
Affiche l'évolution de l'OCP pour chaque électrode sur deux axes y distincts.
"""
times_pos, _, sol_pos = self.solver_pos.simulate()
times_neg, _, sol_neg = self.solver_neg.simulate()
Cs_pos = sol_pos[:, -1]
Cs_neg = sol_neg[:, -1]
ocP_pos = SPMUtilities.compute_ocv(Cs_pos, self.params.Cs_max_pos, electrode="positive")
ocP_neg = SPMUtilities.compute_ocv(Cs_neg, self.params.Cs_max_neg, electrode="negative")
fig, ax1 = plt.subplots(figsize=(8, 5))
ax1.plot(times_pos, ocP_pos, 'r-', label="OCP Positive")
ax1.set_xlabel("Temps (s)")
ax1.set_ylabel("OCP Positive (V)", color='r')
ax1.tick_params(axis='y', labelcolor='r')
ax2 = ax1.twinx()
ax2.plot(times_neg, ocP_neg, 'b-', label="OCP Négative")
ax2.set_ylabel("OCP Négative (V)", color='b')
ax2.tick_params(axis='y', labelcolor='b')
plt.title("Évolution de l'OCP")
fig.tight_layout()
plt.show()
def plot_overpotential(self):
# calcul du j et j0, et des η comme vu précédemment...
times_pos, _, sol_pos = self.solver_pos.simulate()
times_neg, _, sol_neg = self.solver_neg.simulate()
J_surface_pos = SPMUtilities.compute_flux(self.params.I_app, self.params.F, self.params.A_pos, self.params.L_pos, self.params.epsilon_act_pos, self.params.Rp_pos)
J_surface_neg = -SPMUtilities.compute_flux(self.params.I_app, self.params.F, self.params.A_neg, self.params.L_neg, self.params.epsilon_act_neg, self.params.Rp_neg)
jpos = self.params.F * J_surface_pos
jneg = self.params.F * J_surface_neg
cspos = sol_pos[:,-1]
csneg = sol_neg[:,-1]
j0pos = SPMUtilities.compute_j0(self.params.k_pos, self.params.ce, cspos, self.params.Cs_max_pos)
j0neg = SPMUtilities.compute_j0(self.params.k_neg, self.params.ce, csneg, self.params.Cs_max_neg)
etapos = SPMUtilities.compute_eta(jpos, j0pos, self.params.T, self.params.F)
etaneg = SPMUtilities.compute_eta(jneg, j0neg, self.params.T, self.params.F)
fig, ax1 = plt.subplots(figsize=(8, 5))
ax1.plot(times_pos, etapos, 'r-', label="η positive")
ax1.set_xlabel("Temps (s)")
ax1.set_ylabel("Over potentielle Positive (V)", color='r')
ax1.tick_params(axis='y', labelcolor='r')
ax2 = ax1.twinx()
ax2.plot(times_neg, etaneg, 'b-', label="η négative")
ax2.set_ylabel("Over potentielle Négative (V)", color='b')
ax2.tick_params(axis='y', labelcolor='b')
plt.title("Évolution de l'Overtielle")
fig.tight_layout()
plt.show()
def plot_ocv(self):
"""
Affiche l'évolution de l'OCP pour chaque électrode sur deux axes y distincts.
"""
times_pos, _, sol_pos = self.solver_pos.simulate()
times_neg, _, sol_neg = self.solver_neg.simulate()
Cs_pos = sol_pos[:, -1]
Cs_neg = sol_neg[:, -1]
ocP_pos = SPMUtilities.compute_ocv(Cs_pos, self.params.Cs_max_pos, electrode="positive")
ocP_neg = SPMUtilities.compute_ocv(Cs_neg, self.params.Cs_max_neg, electrode="negative")
ocv=ocP_pos-ocP_neg
plt.figure(figsize=(8,5))
plt.plot(times_pos, ocv, color='green')
plt.ylabel('OCV (V)')
plt.title('Tension à circuit ouvert')
plt.tight_layout()
plt.grid()
plt.show()
def plot_cell_voltage(self):
# OCV
times_pos, _, sol_pos = self.solver_pos.simulate()
times_neg, _, sol_neg = self.solver_neg.simulate()
ocv_p = SPMUtilities.compute_ocv(sol_pos[:,-1], self.params.Cs_max_pos, 'positive')
ocv_n = SPMUtilities.compute_ocv(sol_neg[:,-1], self.params.Cs_max_neg, 'negative')
# η
J_surface_pos = SPMUtilities.compute_flux(self.params.I_app, self.params.F, self.params.A_pos, self.params.L_pos, self.params.epsilon_act_pos, self.params.Rp_pos)
J_surface_neg = -SPMUtilities.compute_flux(self.params.I_app, self.params.F, self.params.A_neg, self.params.L_neg, self.params.epsilon_act_neg, self.params.Rp_neg)
jpos = self.params.F * J_surface_pos
jneg = self.params.F * J_surface_neg
cspos = sol_pos[:,-1]
csneg = sol_neg[:,-1]
j0pos = SPMUtilities.compute_j0(self.params.k_pos, self.params.ce, cspos, self.params.Cs_max_pos)
j0neg = SPMUtilities.compute_j0(self.params.k_neg, self.params.ce, csneg, self.params.Cs_max_neg)
etapos = SPMUtilities.compute_eta(jpos, j0pos, self.params.T, self.params.F)
etaneg = SPMUtilities.compute_eta(jneg, j0neg, self.params.T, self.params.F)
# Rct_term
a_s_pos = 3*self.params.epsilon_act_pos/self.params.Rp_pos
a_s_neg = 3*self.params.epsilon_act_neg/self.params.Rp_neg
Rct_term = SPMUtilities.compute_Rct_term(self.params.Rct_pos, a_s_pos, self.params.L_pos,
self.params.Rct_neg, a_s_neg, self.params.L_neg)
Vcell = SPMUtilities.compute_cell_voltage(ocv_p, ocv_n, etapos, etaneg, Rct_term, self.params.I_app)
print(Vcell)
plt.figure(figsize=(8,5))
plt.plot(times_pos, Vcell, color='purple')
plt.ylabel('Tension terminale (V)')
plt.title('Tension terminale')
plt.tight_layout()
plt.grid()
plt.show()
Cell 5: Exemple d’utilisation
c = 5.20689
params = SPMParameters(
I_app= -c, # courant
total_time=3600, # durée totale (s)
dt=0.5, # pas de temps (s)
Nx=50, # nombre de points spatiaux
soc_initiale = 1.0 # l’etat de charge initiale
)
solver_pos = SPMSolver(params, electrode=‘positive’, method=‘implicit’)
solver_neg = SPMSolver(params, electrode=‘negative’, method=‘implicit’)
viz = SPMVisualization(solver_pos, solver_neg, params)
resultats
viz.plot_concentration_evolution() # concentrations de surface
#viz.plot_stoichiometry() # stœchiométries
viz.plot_soc() # soc
#viz.plot_ocP() # potentiels OCP
viz.plot_ocv() # potentiels OCV
#viz.plot_overpotential() # surperpotentiels η
viz.plot_cell_voltage() # tension terminale V_cell