Include contact resistance heating in lumped total heating

CHANGELOG.md

@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Bug fixes
+
+- Contact resistance heating is now included in the total heating of every thermal model, instead of only in the `lumped` temperature equation. ([#5890](https://github.com/pybamm-team/PyBaMM/pull/5890))
+
 # [v26.10.0.0](https://github.com/pybamm-team/PyBaMM/tree/pybamm-v26.10.0.0) - 2026-10-09
 
 ## Breaking changes

packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py

@@ -243,6 +243,20 @@ def _get_standard_coupled_variables(self, variables):
         # Total heating
         Q = Q_ohm + Q_rxn + Q_rev + Q_mix + Q_hys
 
+        # Contact resistance heating (I^2*R), added here so every thermal
+        # model includes it in the total heating variables
+        if self.options["contact resistance"] == "true":
+            I = variables["Current [A]"]
+            T_vol_av = variables["Volume-averaged cell temperature [K]"]
+            V_stack = (
+                self.param.n_electrodes_parallel
+                * self.param.L_y
+                * self.param.L_z
+                * self.param.L_x
+            )
+            Q_cr = I**2 * self.param.R_contact(T_vol_av) / V_stack
+            Q = Q + pybamm.FullBroadcast(Q_cr, Q.domain, "current collector")
+
         # Compute the X-average over the entire cell, including current collectors
         # Note: this can still be a function of y and z for higher-dimensional pouch
         # cell models

packages/pybamm/src/pybamm/models/submodels/thermal/lumped.py

@@ -79,14 +79,15 @@ def get_coupled_variables(self, variables):
 
     def set_rhs(self, variables):
         T_vol_av = variables["Volume-averaged cell temperature [K]"]
+        # "Volume-averaged total heating [W.m-3]" already includes contact
+        # resistance heating (added in BaseThermal)
         Q_vol_av = variables["Volume-averaged total heating [W.m-3]"]
         Q_cool_vol_av = variables["Surface total cooling [W.m-3]"]
-        Q_cr_vol_av = variables["Lumped contact resistance heating [W.m-3]"]
         rho_c_p_eff_av = variables[
             "Volume-averaged effective heat capacity [J.K-1.m-3]"
         ]
 
-        self.rhs = {T_vol_av: (Q_vol_av + Q_cr_vol_av + Q_cool_vol_av) / rho_c_p_eff_av}
+        self.rhs = {T_vol_av: (Q_vol_av + Q_cool_vol_av) / rho_c_p_eff_av}
 
     def set_initial_conditions(self, variables):
         T_vol_av = variables["Volume-averaged cell temperature [K]"]

packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py

@@ -212,8 +212,10 @@ def R_contact(T):
         )
         sol_fn = sim_fn.solve([0, 3600])
 
-        dphi_const = sol_const["Contact overpotential [V]"].entries
-        dphi_fn = sol_fn["Contact overpotential [V]"].entries
+        # Compare on a common time grid: the two solves take different steps
+        t = np.linspace(0, 3600, 100)
+        dphi_const = sol_const["Contact overpotential [V]"](t)
+        dphi_fn = sol_fn["Contact overpotential [V]"](t)
 
         if thermal == "isothermal":
             # temperature is fixed at the reference, so the two agree
@@ -222,3 +224,63 @@ 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(
+        "thermal, extra_options",
+        [
+            ("lumped", {}),
+            ("x-full", {}),
+            (
+                "x-lumped",
+                {"current collector": "potential pair", "dimensionality": 1},
+            ),
+            (
+                "x-lumped",
+                {"current collector": "potential pair", "dimensionality": 2},
+            ),
+        ],
+    )
+    def test_contact_resistance_energy_balance(self, thermal, extra_options):
+        # With the temperature held at ambient, the extra "Total heating [W]"
+        # when the contact resistance is switched on must equal the electrical
+        # energy dissipated in the contact resistance
+        R_c = 0.05
+        t_end = 600
+        t = np.linspace(0, t_end, 61)
+
+        parameter_values = pybamm.ParameterValues("NCA_Kim2011")
+        parameter_values.update(
+            {
+                "Current function [A]": 1,
+                # Hold the temperature at ambient for every thermal option
+                "Total heat transfer coefficient [W.m-2.K-1]": 1e6,
+                "Negative current collector surface heat transfer coefficient [W.m-2.K-1]": 1e6,
+                "Positive current collector surface heat transfer coefficient [W.m-2.K-1]": 1e6,
+                "Negative tab heat transfer coefficient [W.m-2.K-1]": 1e6,
+                "Positive tab heat transfer coefficient [W.m-2.K-1]": 1e6,
+                "Edge heat transfer coefficient [W.m-2.K-1]": 1e6,
+            }
+        )
+
+        heating, current = {}, {}
+        for R in [0, R_c]:
+            params = parameter_values.copy()
+            params.update({"Contact resistance [Ohm]": R})
+            sim = pybamm.Simulation(
+                pybamm.lithium_ion.SPMe(
+                    {
+                        "thermal": thermal,
+                        "contact resistance": "true",
+                        **extra_options,
+                    }
+                ),
+                parameter_values=params,
+            )
+            sol = sim.solve([0, t_end])
+            heating[R] = sol["Total heating [W]"](t)
+            current[R] = sol["Current [A]"](t)
+
+        extra_heat = np.trapezoid(heating[R_c] - heating[0], t)
+        expected = np.trapezoid(current[R_c] ** 2 * R_c, t)
+        assert extra_heat > 0
+        np.testing.assert_allclose(extra_heat, expected, rtol=1e-2)