Skip to content

Request for help with the implementation of the SPM model using an implicit or explicit scheme. #5059

New issue
New issue
Open
Open
Request for help with the implementation of the SPM model using an implicit or explicit scheme.#5059
@simdieudoboinzemwende

Description

@simdieudoboinzemwende
simdieudoboinzemwende
opened on Jun 13, 2025

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.
I need the code structure for developing a PINN framework. Currently, I have a pipeline, but I notice that the resolution is not good. I compared the results with those given by the SPM model in PyBaMM, and the discrepancies are significant. My code is below, but I don’t know what is wrong.

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

1) Instanciez d’abord vos paramètres SPM

params = SPMParameters(
I_app= -c, # courant appliqué (A), négatif pour décharge
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
)

les solveurs pour chaque électrode

solver_pos = SPMSolver(params, electrode='positive', method='implicit')
solver_neg = SPMSolver(params, electrode='negative', method='implicit')

visualisation

viz = SPMVisualization(solver_pos, solver_neg, params)

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() # surpotentiels η
viz.plot_cell_voltage() # tension terminale

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions