Zitao Ni

Ph.D. in Materials Science and Engineering

Study Notes · Part 4: Full SPM Implementation & Voltage Coupling

P2D Model Study Notes · Part 4: Full SPM Implementation & Voltage Coupling

Date: 2026-05-28
Task: Connect OCV, Diffusion, and BV into a complete SPM engine; understand voltage coupling and run simulations
Progress: Full SPM implementation · step() four-step dataflow + voltage coupling + discharge curves + C-rate analysis
Learning Mode: Guided Q&A (question first → think → hand-derive → confirm answer)

Database Scan Results (Pre-Teaching)

1
2
3
4
5
📚 Database Scan Results (Pre-Ch4):
  [Brett]        Ch2 §2.2-2.4  — cell potential calculation
  [Brett]        Ch1 §1.3      — thermodynamics & kinetics, overpotential introduction
  [Rechargeable]  Ch1 §1.2.2   — electrode kinetics & polarization losses
                               — actual voltage = theoretical − polarization − ohmic drop

No new books/chapters added; existing index already covers this chapter.

Learning Method Overview

This chapter connects the knowledge from the first three Sessions (OCV, Diffusion, BV) into one complete working engine.

Session Topic Core Activity Output
4.1 step() four-step dataflow Hand-derive current→flux→overpotential→diffusion→voltage step_dataflow.png
4.2 Voltage coupling + code practice Expand voltage formula → write spm_my.py from skeleton → validate spm_my.py (0 nV diff)
4.3 Discharge curves + C-rate Voltage breakdown chart + C-rate comparison chart spm_discharge.png + crate_comparison.png
4.4 Fully-implicit coupling Derive why diffusion is independent of η → intra-step iteration unnecessary Key insight

Where the Four SPM Equations Sit in the Module

1
2
3
4
5
6
7
8
9
Current i_app
    ↓
Step 1: j = ±i_app/(a_s·L·F)    ← current→flux (Ch1 OCV basics)
    ↓
Step 2: η = BV(j, c_s_surf)      ← flux→overpotential (Ch3 BV)
    ↓        ↑
c_s_surf ← Step 3: diffusion updates c(r) ← concentration←diffusion (Ch2 diffusion)
    ↓
Step 4: V = U_p + η_p − U_n − η_n ← voltage coupling (this chapter's core)

Coupling mode: cross-step coupling — the new c_s_surf produced by Step 3 is used in Step 2 of the next time step.

SPM Dataflow

Session 4.1: step() Four-Step Dataflow

4.1.1 Warm-Up: Review of the Four SPM Equations

# Equation Core Formula Role
① OCV curve $U(\theta) = f(\mathrm{SOC})$ Equilibrium electrode potential
② Diffusion equation $\partial c / \partial t = D \cdot \nabla^2 c$ Spatiotemporal lithium distribution inside spherical particles
③ Butler-Volmer $j = j_0 \cdot [e^{\alpha F\eta/RT} - e^{-(1-\alpha)F\eta/RT}]$ Overpotential needed for interface reaction
④ Voltage coupling $V = U_p + \eta_p - U_n - \eta_n$ First three assembled into terminal voltage

4.1.2 Causal Order of step()

Guided Question: During charging, four things happen — display SOC (A), current enters (B), interface reaction (C), lithium diffusion (D) — which happens first?

Intuitive ordering: A→B→C→D (voltage→current→reaction→diffusion)

Correction: Physically the causality is reversed. The power source forces current → current becomes flux → interface reaction → diffusion → voltage change. Voltage is the result, not the cause. The correct order B→C→D→A corresponds exactly to step()’s four steps:

step() Step Physical Meaning
Step 1: $i_{app} \rightarrow j$ Current converted to reaction flux
Step 2: $j \rightarrow \eta$ BV inversion of overpotential
Step 3: $j$ drives diffusion Update particle interior concentration
Step 4: $c \rightarrow U \rightarrow V$ Compute voltage from new concentrations

4.1.3 Step 1 Derivation: Current → Reaction Flux j

Unit Derivation

Guided Question: What must $i_{app}$ (C/s/m²) be divided by to obtain $j$ (mol/s/m²)?

Answer: Divide by $F$ (C/mol)

Verification:

$$\frac{i_{app}}{F} = \frac{\mathrm{C/(s \cdot m^2)}}{\mathrm{C/mol}} = \frac{\mathrm{mol}}{\mathrm{s \cdot m^2}}$$

Why is $a_s \times L$ in the Denominator?

Guided Question: What is the physical intuition for placing $a_s \times L$ in the denominator?

Answer: $i_{app}$ is constant; the more active material (larger surface area), the smaller the reaction flux per unit area.

Sign Convention

Guided Question: During discharge, the anode oxidizes (j < 0) and the cathode reduces (j > 0). What are the signs of $j_n$ and $j_p$?

Answer: $j_n$ uses $-i_{app}$, $j_p$ uses $+i_{app}$

Final Formula (discharge, $i_{app} > 0$):

$$j_n = -\frac{i_{app}}{a_{s,n} \cdot L_n \cdot F}, \quad j_p = +\frac{i_{app}}{a_{s,p} \cdot L_p \cdot F}$$

Parameter Meaning Unit
$i_{app}$ Applied current density A/m²
$a_s$ Active specific surface area $= 3\varepsilon_s/R$ m⁻¹
$L$ Electrode thickness m
$F$ Faraday constant $= 96485$ C/mol
$j$ Reaction flux mol/(m²·s)

4.1.4 Step 3 Inputs and Initial Concentration

Guided Question: Each time we run one diffusion step, besides fixed parameters, what two varying quantities are needed?

Answer: The current concentration array $c_{old}$ (all grid points) and the surface flux $j$ (from Step 1, as boundary condition). The diffusion solver does not need $\eta$.

Guided Question: At t=0, what is the concentration distribution inside the particle? A) Surface high, center low / B) Center high, surface low / C) Uniform throughout?

Answer: (Choose B — center concentrated, surface depleted — based on the intuition that lithium is stuffed into the particle interior, so the center has more.)

Correction: In the absence of sustained external driving force, diffusion’s steady state is uniform everywhere. At t=0, the cell has been idle for a long time; no current means no interface reaction, and diffusion has already eliminated all concentration differences. Just like an ink drop in a glass of water will eventually become uniform if you wait long enough.

Key insight: When $\partial c / \partial t = 0$ in the diffusion equation, $\nabla^2 c = 0$, which under spherical symmetry implies $c(r) = \mathrm{const}$.

4.1.5 Step 4 Derivation: Voltage Coupling

Guided Question: $V_{cell} = (U_p + \eta_p) - (U_n + \eta_n)$. What does the expansion look like?

Step 1 (substitute):
$$V_{cell} = (U_p + \eta_p) - (U_n + \eta_n)$$

Step 2 (remove parentheses):
$$V_{cell} = U_p + \eta_p - U_n - \eta_n$$

Term Meaning During Discharge
$U_p$ Cathode equilibrium potential 3.54.0 V
$-U_n$ Subtract anode equilibrium potential −0.1 V
$+\eta_p$ Cathode overpotential $\eta_p < 0$ → penalty
$-\eta_n$ Subtract anode overpotential $\eta_n > 0$ → also penalty

During discharge, both overpotential contributions are negative:

$$V_{cell} = \underbrace{(U_p - U_n)}{\mathrm{ideal\ OCV}} - \underbrace{(|\eta_p| + |\eta_n|)}{\mathrm{polarization\ loss}}$$

Session 4.2: Code Practice — Hand-Writing spm_my.py

4.2.1 Skeleton Construction

Import component functions build_diffusion_matrix() + solve_one_step() from Chapter 2 and bv_overpotential() from Chapter 3, then assemble the SPMModelMy class skeleton containing 8 TODOs.

4.2.2 Code Error Summary

Issue Initial Write Correct Write Reason
TODO 3 parameter order (N, D, dr, r, dt) (N, D, dt, dr, r) Function signature’s third parameter is dt
TODO 4 variable confusion -self._A_n / (...) -i_app / (...) self._A_n is the diffusion matrix, not current
TODO 5 extra def def bv_overpotential(...) bv_overpotential(...) def is only for defining, not calling
TODO 5 alpha path self.alpha self.p.alpha alpha lives in the parameter object
TODO 6 matrix confusion self._A_n (used for cathode) self._A_p Cathode diffusion must use its own matrix
TODO 6 wrong punctuation self.p,c_max_p self.p.c_max_p Comma changed to dot

4.2.3 Core Code Segments (Algorithm Pseudocode)

Step 1 (current→flux):

1
2
3
4
5
6
Function: compute_flux(i_app, a_s_n, L_n, F, a_s_p, L_p)
  Purpose: Convert applied current density to reaction flux for both electrodes
  Algorithm:
    1. j_n = -i_app / (a_s_n * L_n * F)
    2. j_p =  i_app / (a_s_p * L_p * F)
  Returns: (j_n, j_p)

Step 2 (BV→overpotential):

1
2
3
4
5
6
Function: compute_overpotential(j_n, j_p, params)
  Purpose: Invert Butler-Volmer to obtain overpotentials
  Algorithm:
    1. eta_n = bv_overpotential(j_n, k_n, c_e, c_s_surf_n(), c_s_max_n, alpha, T)
    2. eta_p = bv_overpotential(j_p, k_p, c_e, c_s_surf_p(), c_s_max_p, alpha, T)
  Returns: (eta_n, eta_p)

Step 3 (diffusion updates concentration):

1
2
3
4
5
6
Function: update_concentration(c_s_n, c_s_p, A_n, A_p, j_n, j_p, D_s_n, D_s_p, dr_n, dr_p)
  Purpose: Advance diffusion by one time step for both electrodes
  Algorithm:
    1. c_s_n_new = solve_one_step(c_s_n, A_n, j_n, D_s_n, dr_n)
    2. c_s_p_new = solve_one_step(c_s_p, A_p, j_p, D_s_p, dr_p)
  Returns: (c_s_n_new, c_s_p_new)

Step 4 (voltage coupling):

1
2
3
4
5
Function: compute_voltage(U_n, U_p, eta_n, eta_p)
  Purpose: Assemble terminal voltage from OCV and overpotentials
  Algorithm:
    1. V_cell = U_p + eta_p - U_n - eta_n
  Returns: V_cell

simulate() (time loop):

1
2
3
4
5
6
7
8
9
10
11
Function: simulate(t_end, i_app, dt)
  Purpose: Run full discharge/charge simulation
  Algorithm:
    1. n_steps = int(t_end / dt)
    2. Initialize time_arr[n_steps+1], V_arr[n_steps+1]
    3. time_arr[0] = 0, V_arr[0] = OCV_p - OCV_n  (t=0: no overpotential, V = pure OCV difference)
    4. For k = 1 to n_steps:
       a. V = step(i_app)
       b. time_arr[k] = current_time
       c. V_arr[k] = V
  Returns: {'time': time_arr, 'voltage': V_arr}

4.2.4 Numerical Verification

1
2
3
4
model.py: V=3.580372V, eta_n=-10.7475mV, eta_p=7.8621mV
spm_my  : V=3.580372V, eta_n=-10.7475mV, eta_p=7.8621mV
V diff = 0.00 nV
c_s_n diff = 0.00e+00

Conclusion: Perfect match — not even a single nanovolt of difference.

Session 4.3: Discharge Curves and Voltage Breakdown

4.3.1 1C Discharge Voltage Breakdown

Metric 1C Discharge (3600s)
$V_{start}$ (ideal OCV) 3.5621 V
$V_{end}$ 3.2470 V
Total drop 315.2 mV
$\Delta(U_p-U_n)$ (OCV change) 348.2 mV (dominant, ~90%)
$\eta_p-\eta_n$ (overpotential contribution) 33.0 mV (secondary, ~10%)

Key insight: The primary driver of voltage drop is the shrinking OCV difference — during discharge, lithium moves from anode to cathode, both electrode SOCs change, and the OCV difference steadily decreases. The overpotential contribution at 1C is relatively small.

Discharge curve four subplots

4.3.2 C-rate Effect

Guided Question: Going from 1C to 2C, does the overpotential contribution increase or decrease?

Answer: It increases

Sub-linearity of arcsinh

$\eta = (2RT/F) \cdot \mathrm{arcsinh}(j/2j_0)$. Properties of arcsinh:

  • Small $x$: $\mathrm{arcsinh}(x) \approx x$ (nearly linear)
  • Large $x$: $\mathrm{arcsinh}(x) \approx \ln(2x)$ (logarithmic growth, slower than linear)

First thought: arcsinh is super-linear → Corrected: arcsinh is sub-linear — doubling the input yields less than double the output. In the operating range (10~20 A/m²), it’s still close to linear: 63.3 ÷ 33.0 ≈ 1.92×.

Why Does the Fraction Also Increase?

Guided Question: For the same SOC change, is $\Delta(U_p-U_n)$ the same at C/5 and at 2C?

Answer: Yes — OCV depends only on SOC, not on current.

Therefore: OCV change is fixed, η grows with current, so its share naturally rises.

Final Data

C-rate Current (A/m²) $\eta_p-\eta_n$ (mV) η Share
C/5 2 6.6 ~2%
C/2 5 16.6 ~5%
1C 10 33.0 ~10%
2C 20 63.3 ~19%

Core conclusion: The higher the C-rate, the larger the overpotential’s share of total voltage drop.

Note: When comparing voltages at different C-rates, comparison must be at the same SOC, not the same time.

C-rate comparison four subplots

Session 4.4: Fully-Implicit Coupling — Step-by-Step Derivation

4.4.1 Does Diffusion Need η?

Guided Question: Does solve_one_step(c_old, ab, j, D, dr) take η as an input?

Answer: No

4.4.2 Does BV Need the Entire c Array?

Guided Question: Does bv_overpotential(j, k, c_e, c_s_surf, ...) need the concentration at every point inside the particle?

Answer: No — only $c_{s,surf}$ is needed.

4.4.3 Piecing the Two Tables Together

Diffusion BV
Input $c_{old}$, $j$ $c_{s,surf}$, $j$
Needs η? No —
Needs entire c? — No

4.4.4 Can a New η Change the Diffusion Result?

Guided Question: Diffusion produces a new $c_{s,surf}$ → BV uses it to compute a new $\eta$ → take that new $\eta$ and re-run diffusion — does the result change?

Answer: No — diffusion does not depend on $\eta$. $c_{old}$ and $j$ are unchanged, so the output is always the same.

4.4.5 Where Is the Coupling?

Guided Question: Where does the coupling between diffusion and BV occur?

Answer: Across time steps — this step’s $c_{s,surf}$ → next step’s $\eta$.

4.4.6 Summary

1
2
3
4
5
6
7
Diffusion depends only on c_old and j (fixed within a step)
    ↓
A new η cannot change the diffusion result
    ↓
Iterating BV and diffusion within the same step is pointless
    ↓
SPM coupling is cross-step: this step's c_surf → next step's η

Key insight: “Fully-implicit coupling” iteration is only needed when diffusion and BV mutually depend on each other as inputs within the same time step. In SPM, they do not directly depend on each other. Staggering by one step is not just an engineering simplification — it is mathematically correct. What truly requires iteration is the P2D model.

Summary — Knowledge Graph

All Formulas Learned

Formula Meaning Source
$j_n = -i_{app}/(a_{s,n} L_n F)$ Anode flux (oxidation during discharge, j<0) §4.1
$j_p = +i_{app}/(a_{s,p} L_p F)$ Cathode flux (reduction during discharge, j>0) §4.1
$\eta = (2RT/F) \cdot \mathrm{arcsinh}(j/2j_0)$ BV inversion (Chapter 3) §3.2
$V_{cell} = U_p + \eta_p - U_n - \eta_n$ Voltage coupling §4.1
$V_{cell} = (U_p-U_n) - ( \eta_p +
Coupling: cross-step conduction, no intra-step iteration SPM coupling nature §4.4

All Code Components

Function/Class Purpose File Status
SPMModelMy.__init__() Initialize grid, concentrations, diffusion matrices spm_my.py ✅ Self-written
step() Advance dt seconds through four steps spm_my.py ✅ Self-written
simulate() Time loop + result recording spm_my.py ✅ Self-written
build_diffusion_matrix() Diffusion matrix construction (Ch2 component) diffusion_solver_my.py ✅
solve_one_step() One-step diffusion (Ch2 component) diffusion_solver_my.py ✅
bv_overpotential() BV inversion (Ch3 component) bv_my.py ✅

Pitfalls Encountered

  1. build_diffusion_matrix parameter order is (N, D, dt, dr, r), not (N, D, dr, r, dt)
  2. Step 1 variable is i_app (function parameter), not self._A_n (diffusion matrix)
  3. Calling a function does not require def
  4. Parameters live in self.p: self.p.alpha, not self.alpha
  5. Cathode diffusion matrix is self._A_p — must not mix with anode’s

Position in the Module

1
2
3
4
5
6
7
8
9
10
┌─────────────────────────────────────────────────────────────┐
│                    The Four SPM Equations                     │
│                                                             │
│  Ch1 OCV    →  U(θ) = f(SOC)         Thermodynamic equilibrium potential │
│  Ch2 Diffusion → ∂c/∂t = D·∇²c        Lithium distribution in particles    │
│  Ch3 BV     →  j ↔ η                 Interface reaction rate               │
│  Ch4 Coupling → V = U + η_p − η_n    Voltage = ideal − loss                │
│                                                             │
│  Missing (SPMe/P2D additions): electrolyte concentration gradient, ohmic drop │
└─────────────────────────────────────────────────────────────┘

Task Completion Status

Task Status
Hand-derive Step 1 (current→flux) sign convention and formula ✅
Hand-derive Step 4 (voltage coupling) expansion and sign verification ✅
Write spm_my.py from skeleton (8 TODOs) ✅
Numerical verification vs model.py (0 nV difference) ✅
Discharge curve four subplots ✅
C-rate comparison four subplots ✅
0%