From 32a94b16cb813d20c8353038ce47af9aa1603368 Mon Sep 17 00:00:00 2001 From: Shaurya Shukla Date: Sun, 4 Oct 2026 21:52:06 -0400 Subject: [PATCH 1/9] fix: include SEI film resistive loss in irreversible heating With "SEI film resistance" on, the voltage drop across the film is added to the reaction overpotential eta_r, which lowers the terminal voltage. The thermal model computed irreversible heating as a*j*eta_r, so the energy lost across the film never appeared in "Total heating [W]" or the cell temperature. Use a*j*(eta_r - eta_sei), i.e. the full departure from the open-circuit potential, for both electrodes, every particle phase, and the planar lithium metal electrode. Add an energy-balance test: at fixed temperature, the extra heat from switching on the film resistance equals the extra electrical loss. Fixes #5820 Co-Authored-By: Claude Opus 5.5 --- .../models/submodels/thermal/base_thermal.py | 19 +++-- .../test_lithium_ion/test_thermal_models.py | 72 +++++++++++++++++++ 2 files changed, 86 insertions(+), 5 deletions(-) diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py index fa5e0b590f..163c6da24c 100644 --- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py +++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py @@ -146,8 +146,12 @@ def _get_standard_coupled_variables(self, variables): f"Positive electrode {phase}volumetric interfacial current density [A.m-3]" ] eta_r_p = variables[f"Positive electrode {phase}reaction overpotential [V]"] - # Irreversible electrochemical heating - Q_rxn_p += a_j_p * eta_r_p + eta_sei_p = variables[ + f"Positive electrode {phase}SEI film overpotential [V]" + ] + # Irreversible electrochemical heating, including the resistive loss + # across the SEI film, which eta_r excludes + Q_rxn_p += a_j_p * (eta_r_p - eta_sei_p) # Reversible electrochemical heating dUdT_p = variables[f"Positive electrode {phase}entropic change [V.K-1]"] Q_rev_p += a_j_p * T_p * dUdT_p @@ -178,8 +182,9 @@ def _get_standard_coupled_variables(self, variables): if self.options.electrode_types["negative"] == "planar": i_n = variables["Lithium metal total interfacial current density [A.m-2]"] eta_r_n = variables["Lithium metal interface reaction overpotential [V]"] + eta_sei_n = variables["Negative electrode SEI film overpotential [V]"] Q_rxn_n = pybamm.PrimaryBroadcast( - i_n * eta_r_n / self.param.n.L, + i_n * (eta_r_n - eta_sei_n) / self.param.n.L, ["negative electrode"], "current collector", ) @@ -199,8 +204,12 @@ def _get_standard_coupled_variables(self, variables): eta_r_n = variables[ f"Negative electrode {phase}reaction overpotential [V]" ] - # Irreversible electrochemical heating - Q_rxn_n += a_j_n * eta_r_n + eta_sei_n = variables[ + f"Negative electrode {phase}SEI film overpotential [V]" + ] + # Irreversible electrochemical heating, including the resistive loss + # across the SEI film, which eta_r excludes + Q_rxn_n += a_j_n * (eta_r_n - eta_sei_n) # Reversible electrochemical heating dUdT_n = variables[f"Negative electrode {phase}entropic change [V.K-1]"] Q_rev_n += a_j_n * T_n * dUdT_n diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py index 7f2609c4e8..f3f52d7ef4 100644 --- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py +++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py @@ -222,3 +222,75 @@ def R_contact(T): # cell heats up above the reference, so R (and hence the contact # overpotential magnitude) is larger than the constant case np.testing.assert_array_less(np.abs(dphi_const[1:]), np.abs(dphi_fn[1:])) + + @pytest.mark.parametrize( + "model, options, parameter_set, current", + [ + ( + pybamm.lithium_ion.SPMe, + {"SEI film resistance": "distributed"}, + "OKane2022", + 5, + ), + ( + pybamm.lithium_ion.SPMe, + {"SEI film resistance": "average"}, + "OKane2022", + 5, + ), + ( + pybamm.lithium_ion.DFN, + {"SEI film resistance": "distributed"}, + "OKane2022", + 5, + ), + ( + pybamm.lithium_ion.SPMe, + {"working electrode": "positive", "SEI film resistance": "distributed"}, + "Xu2019", + 2.4e-3, + ), + ], + ) + def test_sei_film_heating(self, model, options, parameter_set, current): + # With the temperature held at ambient, the extra heat generated when the SEI + # film resistance is switched on must match the extra electrical energy lost + parameter_values = pybamm.ParameterValues(parameter_set) + # Xu2019 has no current collector or thermal parameters, so take them from + # Chen2020 + chen2020 = pybamm.ParameterValues("Chen2020") + parameter_values.update( + {k: chen2020[k] for k in chen2020 if k not in parameter_values}, + check_already_exists=False, + ) + parameter_values.update( + { + "Current function [A]": current, + "Initial SEI thickness [m]": 5e-8, + "Total heat transfer coefficient [W.m-2.K-1]": 1e6, + } + ) + t = np.linspace(0, 600, 61) + + voltage, heating = {}, {} + for film in ["none", options["SEI film resistance"]]: + sim = pybamm.Simulation( + model( + { + **options, + "SEI": "constant", + "SEI film resistance": film, + "thermal": "lumped", + } + ), + parameter_values=parameter_values, + ) + sol = sim.solve([0, 600], t_interp=t) + voltage[film] = sol["Voltage [V]"].entries + heating[film] = sol["Total heating [W]"].entries + + film = options["SEI film resistance"] + electrical_loss = np.trapezoid(current * (voltage["none"] - voltage[film]), t) + extra_heat = np.trapezoid(heating[film] - heating["none"], t) + assert electrical_loss > 0 + np.testing.assert_allclose(extra_heat, electrical_loss, rtol=1e-2) From 4acb518733964614936378d8b627f22fb7fbe33b Mon Sep 17 00:00:00 2001 From: Shaurya Shukla Date: Sun, 4 Oct 2026 21:54:43 -0400 Subject: [PATCH 2/9] add changelog entry for #5841 Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 63133f089a..940b2affe5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ - Discretising a model whose differential variables are all scalars or on zero-dimensional domains (such as a lumped current collector) no longer raises SciPy 1.18's `block_diag` `DeprecationWarning`, and its mass matrix stays a sparse matrix rather than becoming a sparse array once SciPy changes `block_diag`'s return type. ([#5782](https://github.com/pybamm-team/PyBaMM/pull/5782)) - The CI test matrix now installs every dependency from `uv.lock`, as local nox sessions do, instead of resolving the latest releases, so only the prebuilt `pybammsolvers` wheel is installed outside the lock. ([#5782](https://github.com/pybamm-team/PyBaMM/pull/5782)) - An editable `pybammsolvers` install, such as the one `uv sync` creates, can be imported on its own: `import pybammsolvers` no longer fails to load `libcasadi` unless `casadi` was imported first. ([#5827](https://github.com/pybamm-team/PyBaMM/pull/5827)) +- Irreversible heating now includes the resistive loss across the SEI film when `"SEI film resistance"` is `"distributed"` or `"average"`. Before, that loss lowered the voltage but never appeared in `"Total heating [W]"` or the cell temperature. ([#5841](https://github.com/pybamm-team/PyBaMM/pull/5841)) ## Optimizations From 2a020dcaa6416be65752674f42476b42d4b77445 Mon Sep 17 00:00:00 2001 From: Shaurya Shukla Date: Sun, 4 Oct 2026 22:35:35 -0400 Subject: [PATCH 3/9] test: read solution values through the interpolating call interface Co-Authored-By: Claude Opus 5.5 --- .../test_lithium_ion/test_thermal_models.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py index f3f52d7ef4..76e85ccb92 100644 --- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py +++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py @@ -285,9 +285,9 @@ def test_sei_film_heating(self, model, options, parameter_set, current): ), parameter_values=parameter_values, ) - sol = sim.solve([0, 600], t_interp=t) - voltage[film] = sol["Voltage [V]"].entries - heating[film] = sol["Total heating [W]"].entries + sol = sim.solve([0, 600]) + voltage[film] = sol["Voltage [V]"](t) + heating[film] = sol["Total heating [W]"](t) film = options["SEI film resistance"] electrical_loss = np.trapezoid(current * (voltage["none"] - voltage[film]), t) From fb2a00bdb095671b609fd9520e633aba899bb7d7 Mon Sep 17 00:00:00 2001 From: rtimms Date: Thu, 8 Oct 2026 16:54:26 +0100 Subject: [PATCH 4/9] fix: include side-reaction heat from SEI growth and lithium plating Irreversible heating covered only the main intercalation reaction, so the electrical energy lost to SEI growth, SEI on cracks and lithium plating never appeared in the heat source. Add each side-reaction current times (delta_phi - U) for that reaction, and an energy-balance test. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 1 + .../models/submodels/thermal/base_thermal.py | 31 ++++++++++ .../test_thermal/test_energy_balance.py | 61 +++++++++++++++++++ 3 files changed, 93 insertions(+) create mode 100644 packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 330c1f4385..1e097fd5ab 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,7 @@ - An editable `pybammsolvers` install, such as the one `uv sync` creates, can be imported on its own: `import pybammsolvers` no longer fails to load `libcasadi` unless `casadi` was imported first. ([#5827](https://github.com/pybamm-team/PyBaMM/pull/5827)) - `FiniteVolume2D` now raises a `DiscretisationError` when a 2D mesh has fewer than 2 nodes in either direction (for example `BasicDFN2D` with 1 point in `z_2d`), instead of an opaque scipy `ValueError` from the harmonic mean. Use a one-dimensional model for problems that do not vary in that direction. ([#5838](https://github.com/pybamm-team/PyBaMM/pull/5838)) - Irreversible heating now includes the resistive loss across the SEI film when `"SEI film resistance"` is `"distributed"` or `"average"`. Before, that loss lowered the voltage but never appeared in `"Total heating [W]"` or the cell temperature. ([#5841](https://github.com/pybamm-team/PyBaMM/pull/5841)) +- Irreversible heating now includes the heat released by the SEI growth, SEI on cracks and lithium plating side reactions, each as its current times the departure of the surface potential difference from that reaction's open-circuit potential. Before, the electrical energy lost to these reactions never appeared in `"Total heating [W]"` or the cell temperature. ([#PRNUM](https://github.com/pybamm-team/PyBaMM/pull/PRNUM)) ## Optimizations diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py index 163c6da24c..16d28d6c86 100644 --- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py +++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py @@ -152,6 +152,7 @@ def _get_standard_coupled_variables(self, variables): # Irreversible electrochemical heating, including the resistive loss # across the SEI film, which eta_r excludes Q_rxn_p += a_j_p * (eta_r_p - eta_sei_p) + Q_rxn_p += self._side_reaction_heating(variables, "positive", phase) # Reversible electrochemical heating dUdT_p = variables[f"Positive electrode {phase}entropic change [V.K-1]"] Q_rev_p += a_j_p * T_p * dUdT_p @@ -210,6 +211,7 @@ def _get_standard_coupled_variables(self, variables): # Irreversible electrochemical heating, including the resistive loss # across the SEI film, which eta_r excludes Q_rxn_n += a_j_n * (eta_r_n - eta_sei_n) + Q_rxn_n += self._side_reaction_heating(variables, "negative", phase) # Reversible electrochemical heating dUdT_n = variables[f"Negative electrode {phase}entropic change [V.K-1]"] Q_rev_n += a_j_n * T_n * dUdT_n @@ -392,6 +394,35 @@ def _current_collector_heating(self, variables): Q_s_cp = self.param.p.sigma_cc * pybamm.grad_squared(phi_s_cp) return Q_s_cn, Q_s_cp + def _side_reaction_heating(self, variables, domain, phase_name): + """ + Irreversible heating from SEI growth and lithium plating in a porous + electrode, as each side-reaction current times ``delta_phi - U`` for that + reaction, which includes any SEI film drop. + """ + Domain = domain.capitalize() + phase = phase_name.strip() or "primary" + domain_options = getattr(self.options, domain) + phase_options = getattr(domain_options, phase) + delta_phi = variables[f"{Domain} electrode surface potential difference [V]"] + prefix = f"{Domain} electrode {phase_name}" + suffix = "volumetric interfacial current density [A.m-3]" + + Q = 0 + if phase_options["SEI"] not in ["none", "constant"]: + U_sei = self.param.domain_params[domain].phase_params[phase].U_sei + reactions = ["SEI"] + if domain_options["SEI on cracks"] == "true": + reactions.append("SEI on cracks") + for reaction in reactions: + a_j_sei = variables[f"{prefix}{reaction} {suffix}"] + Q += a_j_sei * (delta_phi - U_sei) + if phase_options["lithium plating"] != "none": + # Lithium plating has zero open-circuit potential + a_j_pl = variables[f"{prefix}lithium plating {suffix}"] + Q += a_j_pl * delta_phi + return Q + def _heat_of_mixing(self, variables): """Compute heat of mixing source terms.""" diff --git a/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py b/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py new file mode 100644 index 0000000000..2203a56eba --- /dev/null +++ b/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py @@ -0,0 +1,61 @@ +import numpy as np + +import pybamm + + +class TestEnergyBalance: + def test_side_reaction_heat_matches_electrical_loss(self): + # Ohmic plus irreversible heat must equal the power drawn by all reactions + # at their open-circuit potentials minus the terminal power + model = pybamm.lithium_ion.DFN( + { + "calculate heat source for isothermal models": "true", + "SEI": "solvent-diffusion limited", + "SEI on cracks": "true", + "particle mechanics": "swelling and cracking", + "lithium plating": "partially reversible", + } + ) + param = pybamm.ParameterValues("OKane2022") + param.update( + { + "Current function [A]": -10.0, + "Ambient temperature [K]": 268.15, + "Initial temperature [K]": 268.15, + "Initial SEI thickness [m]": 50e-9, + } + ) + sim = pybamm.Simulation(model, parameter_values=param) + solution = sim.solve([0, 1800], initial_soc=0.1) + t = solution.t + + area = param["Electrode width [m]"] * param["Electrode height [m]"] + U_sei = param["SEI open-circuit potential [V]"] + stored_power = 0 + for domain in ["Negative", "Positive"]: + x = sim.mesh[f"{domain.lower()} electrode"].nodes + + def field(name, domain=domain, x=x): + return solution[f"{domain} electrode {name}"](t=t, x=x) + + power_density = field( + "volumetric interfacial current density [A.m-3]" + ) * field("open-circuit potential [V]") + for reaction in ["SEI ", "SEI on cracks "]: + power_density += ( + field(f"{reaction}volumetric interfacial current density [A.m-3]") + * U_sei + ) + # The default mesh is uniform within each electrode + thickness = param[f"{domain} electrode thickness [m]"] + stored_power += power_density.mean(axis=0) * thickness * area + + heat = solution["Ohmic heating [W]"](t) + solution[ + "Irreversible electrochemical heating [W]" + ](t) + electrical_loss = -stored_power - solution["Current [A]"](t) * solution[ + "Voltage [V]" + ](t) + np.testing.assert_allclose( + np.trapezoid(heat, t), np.trapezoid(electrical_loss, t), rtol=1e-3 + ) From 4658454535a3952dcf3785229531b0e90fa56695 Mon Sep 17 00:00:00 2001 From: rtimms Date: Thu, 8 Oct 2026 19:06:53 +0100 Subject: [PATCH 5/9] add changelog link for #5875 Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1e097fd5ab..ee93ba809a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,7 +25,7 @@ - An editable `pybammsolvers` install, such as the one `uv sync` creates, can be imported on its own: `import pybammsolvers` no longer fails to load `libcasadi` unless `casadi` was imported first. ([#5827](https://github.com/pybamm-team/PyBaMM/pull/5827)) - `FiniteVolume2D` now raises a `DiscretisationError` when a 2D mesh has fewer than 2 nodes in either direction (for example `BasicDFN2D` with 1 point in `z_2d`), instead of an opaque scipy `ValueError` from the harmonic mean. Use a one-dimensional model for problems that do not vary in that direction. ([#5838](https://github.com/pybamm-team/PyBaMM/pull/5838)) - Irreversible heating now includes the resistive loss across the SEI film when `"SEI film resistance"` is `"distributed"` or `"average"`. Before, that loss lowered the voltage but never appeared in `"Total heating [W]"` or the cell temperature. ([#5841](https://github.com/pybamm-team/PyBaMM/pull/5841)) -- Irreversible heating now includes the heat released by the SEI growth, SEI on cracks and lithium plating side reactions, each as its current times the departure of the surface potential difference from that reaction's open-circuit potential. Before, the electrical energy lost to these reactions never appeared in `"Total heating [W]"` or the cell temperature. ([#PRNUM](https://github.com/pybamm-team/PyBaMM/pull/PRNUM)) +- Irreversible heating now includes the heat released by the SEI growth, SEI on cracks and lithium plating side reactions, each as its current times the departure of the surface potential difference from that reaction's open-circuit potential. Before, the electrical energy lost to these reactions never appeared in `"Total heating [W]"` or the cell temperature. ([#5875](https://github.com/pybamm-team/PyBaMM/pull/5875)) ## Optimizations From 06737acf75d5e163403c666146efa4bd346aba19 Mon Sep 17 00:00:00 2001 From: rtimms Date: Thu, 8 Oct 2026 19:15:40 +0100 Subject: [PATCH 6/9] fix: split lithium metal heat between plating and SEI growth In a half cell the total interfacial current is plating plus SEI growth, but irreversible heating applied the plating overpotential to all of it. Use the plating current for the main term and add SEI growth current times (delta_phi - U_sei). Extend the energy-balance test to composite and half-cell models. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 2 +- .../models/submodels/thermal/base_thermal.py | 14 ++- .../test_thermal/test_energy_balance.py | 109 +++++++++++++----- 3 files changed, 95 insertions(+), 30 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ee93ba809a..5a0447bd08 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,7 +25,7 @@ - An editable `pybammsolvers` install, such as the one `uv sync` creates, can be imported on its own: `import pybammsolvers` no longer fails to load `libcasadi` unless `casadi` was imported first. ([#5827](https://github.com/pybamm-team/PyBaMM/pull/5827)) - `FiniteVolume2D` now raises a `DiscretisationError` when a 2D mesh has fewer than 2 nodes in either direction (for example `BasicDFN2D` with 1 point in `z_2d`), instead of an opaque scipy `ValueError` from the harmonic mean. Use a one-dimensional model for problems that do not vary in that direction. ([#5838](https://github.com/pybamm-team/PyBaMM/pull/5838)) - Irreversible heating now includes the resistive loss across the SEI film when `"SEI film resistance"` is `"distributed"` or `"average"`. Before, that loss lowered the voltage but never appeared in `"Total heating [W]"` or the cell temperature. ([#5841](https://github.com/pybamm-team/PyBaMM/pull/5841)) -- Irreversible heating now includes the heat released by the SEI growth, SEI on cracks and lithium plating side reactions, each as its current times the departure of the surface potential difference from that reaction's open-circuit potential. Before, the electrical energy lost to these reactions never appeared in `"Total heating [W]"` or the cell temperature. ([#5875](https://github.com/pybamm-team/PyBaMM/pull/5875)) +- Irreversible heating now includes the heat released by the SEI growth, SEI on cracks and lithium plating side reactions, including SEI growth on a lithium metal electrode, each as its current times the departure of the surface potential difference from that reaction's open-circuit potential. Before, the electrical energy lost to these reactions never appeared in `"Total heating [W]"` or the cell temperature. ([#5875](https://github.com/pybamm-team/PyBaMM/pull/5875)) ## Optimizations diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py index 16d28d6c86..0da9dd4a10 100644 --- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py +++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py @@ -181,11 +181,21 @@ def _get_standard_coupled_variables(self, variables): ocp_options = (ocp_options,) * len(phase_names) if self.options.electrode_types["negative"] == "planar": - i_n = variables["Lithium metal total interfacial current density [A.m-2]"] + # The total current splits between lithium plating and SEI growth + j_pl = variables["Lithium metal plating current density [A.m-2]"] eta_r_n = variables["Lithium metal interface reaction overpotential [V]"] eta_sei_n = variables["Negative electrode SEI film overpotential [V]"] + Q_rxn_n = j_pl * (eta_r_n - eta_sei_n) + if self.options.negative["SEI"] not in ["none", "constant"]: + j_sei = variables[ + "Negative electrode SEI interfacial current density [A.m-2]" + ] + delta_phi_n = variables[ + "Lithium metal interface surface potential difference [V]" + ] + Q_rxn_n += j_sei * (delta_phi_n - self.param.n.prim.U_sei) Q_rxn_n = pybamm.PrimaryBroadcast( - i_n * (eta_r_n - eta_sei_n) / self.param.n.L, + Q_rxn_n / self.param.n.L, ["negative electrode"], "current collector", ) diff --git a/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py b/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py index 2203a56eba..07bdd11743 100644 --- a/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py +++ b/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py @@ -1,54 +1,109 @@ import numpy as np +import pytest import pybamm class TestEnergyBalance: - def test_side_reaction_heat_matches_electrical_loss(self): + @pytest.mark.parametrize( + "options, parameter_set, current, initial_soc", + [ + ( + { + "SEI": "solvent-diffusion limited", + "SEI on cracks": "true", + "particle mechanics": "swelling and cracking", + "lithium plating": "partially reversible", + }, + "OKane2022", + -10.0, + 0.1, + ), + ( + { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + "SEI": "solvent-diffusion limited", + "lithium plating": "partially reversible", + }, + "Bonkile2024", + -10.0, + 0.1, + ), + ( + {"working electrode": "positive", "SEI": "reaction limited"}, + "Xu2019", + 2.4e-3, + None, + ), + ], + ) + def test_side_reaction_heat_matches_electrical_loss( + self, options, parameter_set, current, initial_soc + ): # Ohmic plus irreversible heat must equal the power drawn by all reactions # at their open-circuit potentials minus the terminal power model = pybamm.lithium_ion.DFN( - { - "calculate heat source for isothermal models": "true", - "SEI": "solvent-diffusion limited", - "SEI on cracks": "true", - "particle mechanics": "swelling and cracking", - "lithium plating": "partially reversible", - } + {"calculate heat source for isothermal models": "true", **options} + ) + param = pybamm.ParameterValues(parameter_set) + # Xu2019 has no SEI growth, current collector or thermal parameters + okane2022 = pybamm.ParameterValues("OKane2022") + param.update( + {k: okane2022[k] for k in okane2022 if k not in param}, + check_already_exists=False, ) - param = pybamm.ParameterValues("OKane2022") param.update( { - "Current function [A]": -10.0, + "Current function [A]": current, "Ambient temperature [K]": 268.15, "Initial temperature [K]": 268.15, "Initial SEI thickness [m]": 50e-9, } ) + if parameter_set == "Xu2019": + # Grow SEI fast enough on the lithium metal for its heat to matter + param["SEI reaction exchange current density [A.m-2]"] = 1e-3 sim = pybamm.Simulation(model, parameter_values=param) - solution = sim.solve([0, 1800], initial_soc=0.1) + solution = sim.solve([0, 1800], initial_soc=initial_soc) t = solution.t area = param["Electrode width [m]"] * param["Electrode height [m]"] - U_sei = param["SEI open-circuit potential [V]"] stored_power = 0 - for domain in ["Negative", "Positive"]: - x = sim.mesh[f"{domain.lower()} electrode"].nodes + for domain in ["negative", "positive"]: + Domain = domain.capitalize() + if model.options.electrode_types[domain] == "planar": + # Lithium plating has zero open-circuit potential + j_sei = solution[ + "Negative electrode SEI interfacial current density [A.m-2]" + ] + stored_power += ( + j_sei(t) * param["SEI open-circuit potential [V]"] * area + ) + continue + x = sim.mesh[f"{domain} electrode"].nodes + phases = model.options.phases[domain] + for phase in phases: + phase_name = f"{phase} " if len(phases) > 1 else "" + prefix = f"{phase.capitalize()}: " if len(phases) > 1 else "" + U_sei = param[f"{prefix}SEI open-circuit potential [V]"] - def field(name, domain=domain, x=x): - return solution[f"{domain} electrode {name}"](t=t, x=x) + def field(name, prefix=f"{Domain} electrode {phase_name}", x=x): + return solution[f"{prefix}{name}"](t=t, x=x) - power_density = field( - "volumetric interfacial current density [A.m-3]" - ) * field("open-circuit potential [V]") - for reaction in ["SEI ", "SEI on cracks "]: - power_density += ( - field(f"{reaction}volumetric interfacial current density [A.m-3]") - * U_sei - ) - # The default mesh is uniform within each electrode - thickness = param[f"{domain} electrode thickness [m]"] - stored_power += power_density.mean(axis=0) * thickness * area + power_density = field( + "volumetric interfacial current density [A.m-3]" + ) * field("open-circuit potential [V]") + for reaction in ["SEI ", "SEI on cracks "]: + power_density += ( + field( + f"{reaction}volumetric interfacial current density [A.m-3]" + ) + * U_sei + ) + # The default mesh is uniform within each electrode + thickness = param[f"{Domain} electrode thickness [m]"] + stored_power += power_density.mean(axis=0) * thickness * area heat = solution["Ohmic heating [W]"](t) + solution[ "Irreversible electrochemical heating [W]" From d15d8d0dc9fe7a18de2a5d7b66859e3778c88785 Mon Sep 17 00:00:00 2001 From: Shaurya Shukla Date: Thu, 8 Oct 2026 15:39:53 -0400 Subject: [PATCH 7/9] size-average the SEI film overpotential and add an MPM heating test Co-Authored-By: Claude Opus 5.5 --- .../pybamm/models/submodels/interface/base_interface.py | 4 ++++ .../src/pybamm/models/submodels/thermal/base_thermal.py | 8 ++++---- .../test_lithium_ion/test_thermal_models.py | 8 ++++++++ 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py b/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py index 6624b646f8..cbe7d2851c 100644 --- a/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py +++ b/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py @@ -391,6 +391,10 @@ def _get_standard_sei_film_overpotential_variables(self, eta_sei): } return variables + # Size average, as for eta_r in "_get_standard_overpotential_variables" + if eta_sei.domain and eta_sei.domain[0].endswith("particle size"): + eta_sei = pybamm.size_average(eta_sei) + # Average, and broadcast if necessary eta_sei_av = pybamm.x_average(eta_sei) if eta_sei.domain == []: diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py index 163c6da24c..b7956e004c 100644 --- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py +++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py @@ -149,8 +149,8 @@ def _get_standard_coupled_variables(self, variables): eta_sei_p = variables[ f"Positive electrode {phase}SEI film overpotential [V]" ] - # Irreversible electrochemical heating, including the resistive loss - # across the SEI film, which eta_r excludes + # Irreversible heating. eta_r is the kinetic overpotential net of the + # film drop, so eta_r - eta_sei adds the resistive loss in the SEI film Q_rxn_p += a_j_p * (eta_r_p - eta_sei_p) # Reversible electrochemical heating dUdT_p = variables[f"Positive electrode {phase}entropic change [V.K-1]"] @@ -207,8 +207,8 @@ def _get_standard_coupled_variables(self, variables): eta_sei_n = variables[ f"Negative electrode {phase}SEI film overpotential [V]" ] - # Irreversible electrochemical heating, including the resistive loss - # across the SEI film, which eta_r excludes + # Irreversible heating. eta_r is the kinetic overpotential net of the + # film drop, so eta_r - eta_sei adds the resistive loss in the SEI film Q_rxn_n += a_j_n * (eta_r_n - eta_sei_n) # Reversible electrochemical heating dUdT_n = variables[f"Negative electrode {phase}entropic change [V.K-1]"] diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py index 76e85ccb92..01b3aa4454 100644 --- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py +++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py @@ -244,6 +244,12 @@ def R_contact(T): "OKane2022", 5, ), + ( + pybamm.lithium_ion.MPM, + {"SEI film resistance": "distributed"}, + "OKane2022", + 5, + ), ( pybamm.lithium_ion.SPMe, {"working electrode": "positive", "SEI film resistance": "distributed"}, @@ -263,6 +269,8 @@ def test_sei_film_heating(self, model, options, parameter_set, current): {k: chen2020[k] for k in chen2020 if k not in parameter_values}, check_already_exists=False, ) + if model is pybamm.lithium_ion.MPM: + parameter_values = pybamm.get_size_distribution_parameters(parameter_values) parameter_values.update( { "Current function [A]": current, From c412b242971e389171293ad9451312abe2535431 Mon Sep 17 00:00:00 2001 From: rtimms Date: Mon, 12 Oct 2026 07:43:55 +0100 Subject: [PATCH 8/9] test: cover side-reaction heating in unit tests and pass defaulted options explicitly Add well-posedness tests for lumped thermal with SEI growth, SEI on cracks and plating, and for a half cell with SEI growth, so the new heating branches are covered by the unit suite. Pass the options that #5818 deprecates relying on. Co-Authored-By: Claude Opus 5.5 --- .../test_thermal/test_energy_balance.py | 9 ++++++++- .../base_lithium_ion_half_cell_tests.py | 9 +++++++++ .../test_lithium_ion/base_lithium_ion_tests.py | 13 +++++++++++++ 3 files changed, 30 insertions(+), 1 deletion(-) diff --git a/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py b/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py index 07bdd11743..6ef7fdcf3b 100644 --- a/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py +++ b/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py @@ -13,6 +13,7 @@ class TestEnergyBalance: "SEI": "solvent-diffusion limited", "SEI on cracks": "true", "particle mechanics": "swelling and cracking", + "stress-induced diffusion": "true", "lithium plating": "partially reversible", }, "OKane2022", @@ -25,6 +26,7 @@ class TestEnergyBalance: "open-circuit potential": (("single", "current sigmoid"), "single"), "SEI": "solvent-diffusion limited", "lithium plating": "partially reversible", + "surface form": "algebraic", }, "Bonkile2024", -10.0, @@ -44,7 +46,12 @@ def test_side_reaction_heat_matches_electrical_loss( # Ohmic plus irreversible heat must equal the power drawn by all reactions # at their open-circuit potentials minus the terminal power model = pybamm.lithium_ion.DFN( - {"calculate heat source for isothermal models": "true", **options} + { + "calculate heat source for isothermal models": "true", + "SEI film resistance": "distributed", + "total interfacial current density as a state": "true", + **options, + } ) param = pybamm.ParameterValues(parameter_set) # Xu2019 has no SEI growth, current collector or thermal parameters diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py index e8fb3d9eb6..66f7885a81 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py @@ -58,6 +58,15 @@ def test_well_posed_reaction_limited_sei(self): options = {"SEI": "reaction limited"} self.check_well_posedness(options) + def test_well_posed_reaction_limited_sei_lumped_thermal(self): + options = { + "SEI": "reaction limited", + "thermal": "lumped", + "SEI film resistance": "distributed", + "total interfacial current density as a state": "true", + } + self.check_well_posedness(options) + def test_well_posed_asymmetric_reaction_limited_sei(self): options = {"SEI": "reaction limited (asymmetric)"} self.check_well_posedness(options) diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py index 965f5b8949..1a0cc7b534 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py @@ -372,6 +372,19 @@ def test_well_posed_sei_reaction_limited_on_cracks(self): } self.check_well_posedness(options) + def test_well_posed_side_reactions_lumped_thermal(self): + options = { + "SEI": "reaction limited", + "SEI on cracks": "true", + "particle mechanics": "swelling and cracking", + "lithium plating": "reversible", + "thermal": "lumped", + "SEI film resistance": "distributed", + "stress-induced diffusion": "true", + "total interfacial current density as a state": "true", + } + self.check_well_posedness(options) + def test_well_posed_sei_solvent_diffusion_limited_on_cracks(self): options = { "SEI": "solvent-diffusion limited", From b71f139776c1837cefa0aa779e6f78c03475398f Mon Sep 17 00:00:00 2001 From: Robert Timms <43040151+rtimms@users.noreply.github.com> Date: Mon, 2 Nov 2026 12:48:52 +0000 Subject: [PATCH 9/9] docs: shorten changelog entries to the 40-word limit --- CHANGELOG.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ceeb89e648..216453b669 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -36,8 +36,8 @@ - `RegPower` now evaluates the same operations in `evaluate` and CasADi, taking a literal exponent of 0.5 or 2 as square roots, so they agree bitwise; CasADi used a `hypot` form. With `delta=0`, `evaluate` and the Jacobian use `base ** exponent`, as CasADi did. ([#5955](https://github.com/pybamm-team/PyBaMM/pull/5955)) - CasADi linear and `pchip` interpolants on a uniform grid now find a point's segment by binary search. The `"exact"` lookup could pick the neighbouring segment at or one ulp beside a knot, giving the Jacobian that segment's slope. ([#5955](https://github.com/pybamm-team/PyBaMM/pull/5955)) - `QuickPlot`, and so `Simulation.plot`, now uses matplotlib's color cycle (`axes.prop_cycle`) for default colors when a style sheet or rcParams changes it; otherwise it keeps its red, blue, black, green, magenta and cyan colors. ([#5953](https://github.com/pybamm-team/PyBaMM/pull/5953)) -- Irreversible heating now includes the resistive loss across the SEI film when `"SEI film resistance"` is `"distributed"` or `"average"`. Before, that loss lowered the voltage but never appeared in `"Total heating [W]"` or the cell temperature. ([#5841](https://github.com/pybamm-team/PyBaMM/pull/5841)) -- Irreversible heating now includes the heat released by the SEI growth, SEI on cracks and lithium plating side reactions, including SEI growth on a lithium metal electrode, each as its current times the departure of the surface potential difference from that reaction's open-circuit potential. Before, the electrical energy lost to these reactions never appeared in `"Total heating [W]"` or the cell temperature. ([#5875](https://github.com/pybamm-team/PyBaMM/pull/5875)) +- Irreversible heating now includes the resistive loss across the SEI film when `"SEI film resistance"` is `"distributed"` or `"average"`, so it reaches `"Total heating [W]"` and the cell temperature. ([#5841](https://github.com/pybamm-team/PyBaMM/pull/5841)) +- Irreversible heating now includes the heat from SEI growth, SEI on cracks and lithium plating side reactions, so their electrical losses reach `"Total heating [W]"` and the cell temperature. ([#5875](https://github.com/pybamm-team/PyBaMM/pull/5875)) ## Optimizations