diff --git a/CHANGELOG.md b/CHANGELOG.md index e9f8ccc197..3027411c7d 100644 --- a/CHANGELOG.md +++ b/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 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..dab4a730f2 100644 --- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py +++ b/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 diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/lumped.py b/packages/pybamm/src/pybamm/models/submodels/thermal/lumped.py index a0305e23cf..3c0e36381a 100644 --- a/packages/pybamm/src/pybamm/models/submodels/thermal/lumped.py +++ b/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]"] 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..d686d3ddd8 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 @@ -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)