From 83bf4e4f2ae062b7e060cec650c468a67901eb50 Mon Sep 17 00:00:00 2001 From: bradyplanden Date: Tue, 22 Sep 2026 10:20:41 +0100 Subject: [PATCH] fix: sum the electrochemical heat sources over every particle phase MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `BaseThermal._get_standard_coupled_variables` builds `Q_rxn`, `Q_rev` and `Q_hys` by zipping the phase names against the `"open-circuit potential"` option. That option stays a plain string when it is left at its default, so `[ocp_options]` had length one and zip(["primary ", "secondary "], ["single"], strict=False) silently truncated the loop after the primary phase. On a composite electrode the secondary phase contributed no irreversible, reversible or hysteresis heat at all. A Chen2020 graphite electrode split into two identical half-volume-fraction phases reported 13657 W.m-3 of volume-averaged irreversible heating against the 24371 W.m-3 of the equivalent single-phase electrode, and its lumped cell temperature ran 2 K cold over 1200 s. Passing the option once per phase, as `(("single", "single"), "single")`, was a workaround and gave the right answer. Padding a scalar option out to the phase count is the whole fix; the zips stay non-strict, so the one remaining length mismatch — two open-circuit potential options declared against a single particle phase, which the options validator accepts today — keeps being ignored rather than raising from inside the thermal submodel. The bug dates to #4893. Found while reviewing #5769, which fixes the same class of omission in the heat of mixing. --- CHANGELOG.md | 1 + .../models/submodels/thermal/base_thermal.py | 6 +- .../test_compare_outputs_two_phase.py | 61 +++++++++++++++++++ 3 files changed, 66 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index aafbacc4d0..2f837dd66b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ ## Bug fixes +- The irreversible, reversible and hysteresis heat sources now sum over every particle phase. `"open-circuit potential"` stays a plain string when it is left at its default, so on a composite electrode it was zipped against the two phase names and silently truncated the loop after the primary phase: a graphite electrode split into two identical halves reported 56% of the irreversible heating of the equivalent single-phase electrode. Passing the option once per phase was a workaround. ([#5770](https://github.com/pybamm-team/PyBaMM/pull/5770)) - `pybamm.Function(scipy.special.erf, ...)` converts to CasADi again; CasADi 3.8 stopped dispatching scipy's `erf` ufunc onto symbolic values, so the generic `Function` path now maps it to `casadi.erf` explicitly. `pybamm.erf` was unaffected. ([#5761](https://github.com/pybamm-team/PyBaMM/pull/5761)) - The dimensional minimum particle surface concentration is registered as `"Minimum particle surface concentration [mol.m-3]"`. An implicit string concatenation repeated its first line, so the variable was only reachable under the doubled name `"Minimum negative particle Minimum negative particle surface concentration [mol.m-3]"`, and the intended name resolved to nothing. The `Maximum` sibling was unaffected. ([#5759](https://github.com/pybamm-team/PyBaMM/pull/5759)) - `FiniteVolume` node-to-edge shifts now use the true node spacing instead of weights that assume a uniform mesh, so the exterior edge values agree with `boundary_value`. On a non-uniform mesh they did not, which broke lithium conservation wherever a flux is imposed as a gradient divided by a surface value, as the Fickian particle does: with a concentration-dependent particle diffusivity on an `Exponential1DSubMesh`, `ORegan2022` gained 26% of its particle lithium over three 1C cycles. Uniform-mesh results are unchanged. ([#5755](https://github.com/pybamm-team/PyBaMM/pull/5755)) 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 6a00e1ef70..fa5e0b590f 100644 --- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py +++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py @@ -136,7 +136,8 @@ def _get_standard_coupled_variables(self, variables): phase_names = ["primary ", "secondary "] ocp_options = self.options.positive["open-circuit potential"] if isinstance(ocp_options, str): - ocp_options = [ocp_options] + # a single option applies to every phase + ocp_options = (ocp_options,) * len(phase_names) Q_rxn_p, Q_rev_p, Q_hys_p = 0, 0, 0 T_p = variables["Positive electrode temperature [K]"] @@ -171,7 +172,8 @@ def _get_standard_coupled_variables(self, variables): phase_names = ["primary ", "secondary "] ocp_options = self.options.negative["open-circuit potential"] if isinstance(ocp_options, str): - ocp_options = [ocp_options] + # a single option applies to every phase + 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]"] diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py index fdd2823f0f..dea71aca7a 100644 --- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py +++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_compare_outputs_two_phase.py @@ -176,3 +176,64 @@ def test_compare_SPMe_silicon_graphite(self): def test_compare_DFN_silicon_graphite(self): model_class = pybamm.lithium_ion.DFN self.compare_outputs_two_phase_silicon_graphite(model_class) + + def compare_heat_sources_two_phase_graphite_graphite(self, model_class): + """ + Every heat source is a sum over the particle phases, so splitting one graphite + phase into two identical halves must reproduce the one-phase heat sources + """ + options = { + "thermal": "isothermal", + "calculate heat source for isothermal models": "true", + } + t_eval = [0, 3600] + t_interp = np.linspace(0, 3600) + + parameter_values = pybamm.ParameterValues("Chen2020") + sol = pybamm.Simulation( + model_class(options), parameter_values=parameter_values + ).solve(t_eval=t_eval, t_interp=t_interp) + + parameter_values_two_phase = pybamm.ParameterValues("Chen2020") + for parameter in [ + "Negative electrode OCP [V]", + "Negative electrode OCP entropic change [V.K-1]", + "Maximum concentration in negative electrode [mol.m-3]", + "Initial concentration in negative electrode [mol.m-3]", + "Negative particle radius [m]", + "Negative particle diffusivity [m2.s-1]", + "Negative electrode exchange-current density [A.m-2]", + "Negative electrode active material volume fraction", + ]: + value = parameter_values_two_phase[parameter] + if parameter.endswith("active material volume fraction"): + value = value / 2 + parameter_values_two_phase.update( + {f"Primary: {parameter}": value, f"Secondary: {parameter}": value} + ) + del parameter_values_two_phase[parameter] + + sol_two_phase = pybamm.Simulation( + model_class({"particle phases": ("2", "1"), **options}), + parameter_values=parameter_values_two_phase, + ).solve(t_eval=t_eval, t_interp=t_interp) + + for variable in [ + "Volume-averaged irreversible electrochemical heating [W.m-3]", + "Volume-averaged reversible heating [W.m-3]", + "Volume-averaged hysteresis electrochemical heating [W.m-3]", + "Volume-averaged total heating [W.m-3]", + "Voltage [V]", + ]: + np.testing.assert_allclose( + sol[variable](t_interp), + sol_two_phase[variable](t_interp), + rtol=1e-2, + atol=1e-8, + ) + + def test_compare_heat_sources_SPM_graphite_graphite(self): + self.compare_heat_sources_two_phase_graphite_graphite(pybamm.lithium_ion.SPM) + + def test_compare_heat_sources_DFN_graphite_graphite(self): + self.compare_heat_sources_two_phase_graphite_graphite(pybamm.lithium_ion.DFN)