fix: include side-reaction heat from SEI growth and lithium plating

CHANGELOG.md

@@ -36,6 +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"`, 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
 

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 == []:

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

@@ -146,8 +146,13 @@ 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 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)
+            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
@@ -176,10 +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 / self.param.n.L,
+                Q_rxn_n / self.param.n.L,
                 ["negative electrode"],
                 "current collector",
             )
@@ -199,8 +215,13 @@ 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 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)
+                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
@@ -383,6 +404,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."""
 

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

@@ -222,3 +222,83 @@ 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.MPM,
+                {"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,
+        )
+        if model is pybamm.lithium_ion.MPM:
+            parameter_values = pybamm.get_size_distribution_parameters(parameter_values)
+        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])
+            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)
+        extra_heat = np.trapezoid(heating[film] - heating["none"], t)
+        assert electrical_loss > 0
+        np.testing.assert_allclose(extra_heat, electrical_loss, rtol=1e-2)

packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py

@@ -0,0 +1,123 @@
+import numpy as np
+import pytest
+
+import pybamm
+
+
+class TestEnergyBalance:
+    @pytest.mark.parametrize(
+        "options, parameter_set, current, initial_soc",
+        [
+            (
+                {
+                    "SEI": "solvent-diffusion limited",
+                    "SEI on cracks": "true",
+                    "particle mechanics": "swelling and cracking",
+                    "stress-induced diffusion": "true",
+                    "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",
+                    "surface form": "algebraic",
+                },
+                "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 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
+        okane2022 = pybamm.ParameterValues("OKane2022")
+        param.update(
+            {k: okane2022[k] for k in okane2022 if k not in param},
+            check_already_exists=False,
+        )
+        param.update(
+            {
+                "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=initial_soc)
+        t = solution.t
+
+        area = param["Electrode width [m]"] * param["Electrode height [m]"]
+        stored_power = 0
+        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, 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
+
+        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
+        )

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)

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",