fix: heat of mixing with polynomial particle profiles

CHANGELOG.md

@@ -30,6 +30,9 @@
 - Fixed copying a function without simplification, as `ParameterValues` does for a parameter set loaded with `Serialise.load_parameters`, returning a generic `pybamm.Function`: `pybamm.exp` and the other specific functions lost their class, and an `Arcsinh2` copy could not be evaluated. ([#5954](https://github.com/pybamm-team/PyBaMM/pull/5954))
 - `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))
+- The heat of mixing now sums over every particle phase, so `"heat of mixing": "true"` works on a composite electrode. It previously read the single-phase variable names (`"Negative particle radius [m]"`) and `param.n.prim` only, so any model with `"particle phases": "2"` raised `ModelError: Missing variable for submodel 'thermal': 'Negative particle radius [m]'` at build time. Single-phase results are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
+- The heat of mixing now uses the particle's effective diffusivity, the one its flux uses, instead of re-evaluating the base diffusivity. With `"diffusivity": "current sigmoid"` or `"stress-induced diffusion": "true"` the two differ, and the heat sources stopped conserving energy: a phase with lithiation and delithiation diffusivities of a third and three times the base value was off by 33% of the heat. A uniform particle profile now has zero heat of mixing, so a DFN with `"particle": "uniform profile"` and `"heat of mixing": "true"` builds again. Results with the default options are unchanged. ([#5769](https://github.com/pybamm-team/PyBaMM/pull/5769))
+- `"heat of mixing": "true"` now works with `"particle": "quadratic profile"` or `"quartic profile"`, which failed to discretise with a `ShapeError`. These profiles reconstruct the concentration as a polynomial in `r`, so the heat of mixing now uses their closed-form flux instead of a finite-volume gradient. ([#5935](https://github.com/pybamm-team/PyBaMM/pull/5935))
 
 ## Optimizations
 

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

@@ -383,65 +383,76 @@ 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 _phase_heat_of_mixing(self, variables, domain, phase_name, phase_param):
+        """
+        Compute the heat of mixing source term for a single particle phase.
+
+        The OCP is the phase's equilibrium OCP, the same one the hysteresis heat is
+        referenced to, not the hysteresis-weighted OCP that sets the voltage. Energy
+        only balances if both heat sources share that reference.
+        """
+        phase_options = getattr(getattr(self.options, domain), phase_param.phase)
+        if phase_options["particle"] == "uniform profile":
+            # no concentration gradient, so nothing to mix
+            return pybamm.FullBroadcast(0, [f"{domain} electrode"], "current collector")
+        Domain = domain.capitalize()
+        F = pybamm.constants.F.value
+        a = variables[
+            f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]"
+        ]
+        R = variables[f"{Domain} {phase_name}particle radius [m]"]
+        N = a / (4 * np.pi * R**2)
+        # the diffusivity the particle flux actually uses, including any current
+        # sigmoid or stress-induced diffusion, so the dissipation matches the flux
+        prefix = f"X-averaged {domain}" if self.x_average else Domain
+        c = variables[f"{prefix} {phase_name}particle concentration [mol.m-3]"]
+        T = variables[f"{prefix} electrode temperature [K]"]
+        D = variables[f"{prefix} {phase_name}particle effective diffusivity [m2.s-1]"]
+        T_part = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle"])
+        U_eq = phase_param.U(c / phase_param.c_max, T_part)
+        if phase_options["particle"] in ["quadratic profile", "quartic profile"]:
+            # c is a polynomial in r on the nodes, so D (dc/dr)^2 dU/dc = -N dU/dr
+            # with the closed-form flux N; a finite-volume gradient doesn't apply
+            N_s = variables[f"{prefix} {phase_name}particle flux [mol.m-2.s-1]"]
+            r = pybamm.SpatialVariable(
+                f"r_{domain[0]}", domains=c.domains, coord_sys="spherical polar"
+            )
+            integrand_r = -N_s * U_eq.diff(r)
+        else:
+            dc_dr2 = pybamm.inner(pybamm.grad(c), pybamm.grad(c))
+            integrand_r = D * dc_dr2 * U_eq.diff(c)
+        integration_variable_r = [
+            pybamm.SpatialVariable("r", domain=integrand_r.domain)
+        ]
+        integral_r = pybamm.Integral(integrand_r, integration_variable_r)
+        return -F * N * integral_r
+
+    def _electrode_heat_of_mixing(self, variables, domain):
+        """Sum the heat of mixing over the particle phases of one electrode."""
+        electrode_param = self.param.domain_params[domain]
+        phases = self.options.phases[domain]
+        names = [""] if len(phases) == 1 else [f"{phase} " for phase in phases]
+        return sum(
+            self._phase_heat_of_mixing(
+                variables, domain, name, electrode_param.phase_params[phase]
+            )
+            for name, phase in zip(names, phases, strict=True)
+        )
+
     def _heat_of_mixing(self, variables):
         """Compute heat of mixing source terms."""
 
         if self.options["heat of mixing"] == "true":
-            F = pybamm.constants.F.value
-            pi = np.pi
-
             # Compute heat of mixing in negative electrode
             if self.options.electrode_types["negative"] == "planar":
                 Q_mix_s_n = pybamm.FullBroadcast(
                     0, ["negative electrode"], "current collector"
                 )
             else:
-                a_n = variables["Negative electrode surface area to volume ratio [m-1]"]
-                R_n = variables["Negative particle radius [m]"]
-                N_n = a_n / (4 * pi * R_n**2)
-                if self.x_average:
-                    c_n = variables[
-                        "X-averaged negative particle concentration [mol.m-3]"
-                    ]
-                    T_n = variables["X-averaged negative electrode temperature [K]"]
-                else:
-                    c_n = variables["Negative particle concentration [mol.m-3]"]
-                    T_n = variables["Negative electrode temperature [K]"]
-                T_n_part = pybamm.PrimaryBroadcast(T_n, ["negative particle"])
-                dc_n_dr2 = pybamm.inner(pybamm.grad(c_n), pybamm.grad(c_n))
-                D_n = self.param.n.prim.D(c_n, T_n_part)
-                dUeq_n = self.param.n.prim.U(
-                    c_n / self.param.n.prim.c_max, T_n_part
-                ).diff(c_n)
-                integrand_r_n = D_n * dc_n_dr2 * dUeq_n
-                integration_variable_r_n = [
-                    pybamm.SpatialVariable("r", domain=integrand_r_n.domain)
-                ]
-                integral_r_n = pybamm.Integral(integrand_r_n, integration_variable_r_n)
-                Q_mix_s_n = -F * N_n * integral_r_n
+                Q_mix_s_n = self._electrode_heat_of_mixing(variables, "negative")
 
             # Compute heat of mixing in positive electrode
-            a_p = variables["Positive electrode surface area to volume ratio [m-1]"]
-            R_p = variables["Positive particle radius [m]"]
-            N_p = a_p / (4 * pi * R_p**2)
-            if self.x_average:
-                c_p = variables["X-averaged positive particle concentration [mol.m-3]"]
-                T_p = variables["X-averaged positive electrode temperature [K]"]
-            else:
-                c_p = variables["Positive particle concentration [mol.m-3]"]
-                T_p = variables["Positive electrode temperature [K]"]
-            T_p_part = pybamm.PrimaryBroadcast(T_p, ["positive particle"])
-            dc_p_dr2 = pybamm.inner(pybamm.grad(c_p), pybamm.grad(c_p))
-            D_p = self.param.p.prim.D(c_p, T_p_part)
-            dUeq_p = self.param.p.prim.U(c_p / self.param.p.prim.c_max, T_p_part).diff(
-                c_p
-            )
-            integrand_r_p = D_p * dc_p_dr2 * dUeq_p
-            integration_variable_r_p = [
-                pybamm.SpatialVariable("r", domain=integrand_r_p.domain)
-            ]
-            integral_r_p = pybamm.Integral(integrand_r_p, integration_variable_r_p)
-            Q_mix_s_p = -F * N_p * integral_r_p
+            Q_mix_s_p = self._electrode_heat_of_mixing(variables, "positive")
             Q_mix_s_s = pybamm.FullBroadcast(0, ["separator"], "current collector")
         else:
             Q_mix_s_n = pybamm.FullBroadcast(

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

@@ -178,24 +178,40 @@ 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):
+    def compare_heat_sources_two_phase_graphite_graphite(
+        self, model_class, thermal_options, parameter_set="Chen2020"
+    ):
         """
         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
         """
+        # composite models default to the algebraic surface form, so pin it for the
+        # one-phase model too and compare like with like
         options = {
-            "thermal": "isothermal",
-            "calculate heat source for isothermal models": "true",
+            **thermal_options,
+            "heat of mixing": "true",
+            "surface form": "algebraic",
         }
         t_eval = [0, 3600]
         t_interp = np.linspace(0, 3600)
 
-        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values = pybamm.ParameterValues(parameter_set)
+        if parameter_set == "Chen2020":
+            # Chen2020 has zero entropic change, which would leave the reversible
+            # heat and the temperature dependence of the heat of mixing untested
+            entropic_change = pybamm.ParameterValues("Ai2020")
+            for domain in ["Negative", "Positive"]:
+                name = f"{domain} electrode OCP entropic change [V.K-1]"
+                parameter_values[name] = entropic_change[name]
+
+        # the default rtol lets a temperature state near 300 K drift by ~0.03 K,
+        # which is comparable to the one- and two-phase differences being tested
+        solver = pybamm.IDAKLUSolver(rtol=1e-6, atol=1e-8)
         sol = pybamm.Simulation(
-            model_class(options), parameter_values=parameter_values
+            model_class(options), parameter_values=parameter_values, solver=solver
         ).solve(t_eval=t_eval, t_interp=t_interp)
 
-        parameter_values_two_phase = pybamm.ParameterValues("Chen2020")
+        parameter_values_two_phase = parameter_values.copy()
         for parameter in [
             "Negative electrode OCP [V]",
             "Negative electrode OCP entropic change [V.K-1]",
@@ -217,27 +233,129 @@ def compare_heat_sources_two_phase_graphite_graphite(self, model_class):
         sol_two_phase = pybamm.Simulation(
             model_class({"particle phases": ("2", "1"), **options}),
             parameter_values=parameter_values_two_phase,
+            solver=solver,
         ).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 heat of mixing [W.m-3]",
             "Volume-averaged total heating [W.m-3]",
             "Voltage [V]",
         ]:
+            one_phase = sol[variable](t_interp)
+            # scale atol to the variable, since the reversible heat changes sign
             np.testing.assert_allclose(
-                sol[variable](t_interp),
                 sol_two_phase[variable](t_interp),
+                one_phase,
+                rtol=1e-2,
+                atol=1e-3 * np.nanmax(np.abs(one_phase)) + 1e-8,
+            )
+
+        if options["thermal"] != "isothermal":
+            temperature = "Volume-averaged cell temperature [K]"
+            temperature_rise = sol[temperature](t_interp) - sol[temperature](0)
+            temperature_rise_two_phase = sol_two_phase[temperature](
+                t_interp
+            ) - sol_two_phase[temperature](0)
+            np.testing.assert_allclose(
+                temperature_rise_two_phase,
+                temperature_rise,
                 rtol=1e-2,
-                atol=1e-8,
+                atol=1e-3 * np.nanmax(temperature_rise),
             )
 
-    def test_compare_heat_sources_SPM_graphite_graphite(self):
-        self.compare_heat_sources_two_phase_graphite_graphite(pybamm.lithium_ion.SPM)
+    @pytest.mark.parametrize(
+        ("model_class", "thermal_options", "parameter_set"),
+        [
+            pytest.param(model_class, options, parameter_set, id=f"{name}-{label}")
+            for name, model_class in [
+                ("SPM", pybamm.lithium_ion.SPM),
+                ("SPMe", pybamm.lithium_ion.SPMe),
+                ("DFN", pybamm.lithium_ion.DFN),
+            ]
+            for label, options, parameter_set in [
+                (
+                    "isothermal",
+                    {
+                        "thermal": "isothermal",
+                        "calculate heat source for isothermal models": "true",
+                    },
+                    "Chen2020",
+                ),
+                ("lumped", {"thermal": "lumped"}, "Chen2020"),
+                # x-full resolves the temperature in x, so the heat of mixing uses
+                # the local rather than the x-averaged electrode temperature
+                (
+                    "x-full",
+                    {"thermal": "x-full", "cell geometry": "pouch"},
+                    "Marquis2019",
+                ),
+                (
+                    "x-lumped-1plus1D",
+                    {
+                        "thermal": "x-lumped",
+                        "cell geometry": "pouch",
+                        "current collector": "potential pair",
+                        "dimensionality": 1,
+                    },
+                    "Marquis2019",
+                ),
+            ]
+        ],
+    )
+    def test_compare_heat_sources_graphite_graphite(
+        self, model_class, thermal_options, parameter_set
+    ):
+        self.compare_heat_sources_two_phase_graphite_graphite(
+            model_class, thermal_options, parameter_set
+        )
+
+    def heat_of_mixing_silicon_graphite(self, model_class):
+        # Check the heat of mixing of two genuinely different phases heats the cell
+        options = {
+            "particle phases": ("2", "1"),
+            "open-circuit potential": (("single", "current sigmoid"), "single"),
+            "thermal": "lumped",
+        }
+        parameter_values = pybamm.ParameterValues("Chen2020_composite")
+        # thermal properties are per layer, but the set only has per-phase densities
+        parameter_values.update(
+            {"Negative electrode density [kg.m-3]": 1657.0},
+            check_already_exists=False,
+        )
+        # the shipped graphite OCP is a cubic Interpolant, whose derivative cannot
+        # yet be converted for the solver, so use the analytic Chen2020 fit
+        parameter_values["Primary: Negative electrode OCP [V]"] = (
+            pybamm.ParameterValues("Chen2020")["Negative electrode OCP [V]"]
+        )
+        t_eval = [0, 3600]
+        t_interp = np.linspace(0, 3600)
+
+        solutions = {
+            heat_of_mixing: pybamm.Simulation(
+                model_class({**options, "heat of mixing": heat_of_mixing}),
+                parameter_values=parameter_values,
+            ).solve(t_eval=t_eval, t_interp=t_interp)
+            for heat_of_mixing in ["false", "true"]
+        }
+
+        heat_of_mixing = solutions["true"]["Volume-averaged heat of mixing [W.m-3]"](
+            t_interp
+        )
+        heat_of_mixing = heat_of_mixing[~np.isnan(heat_of_mixing)]
+        assert np.all(heat_of_mixing > 0)
+        temperature = "Volume-averaged cell temperature [K]"
+        assert solutions["true"][temperature](t_interp[-2]) > solutions["false"][
+            temperature
+        ](t_interp[-2])
+
+    def test_heat_of_mixing_SPM_silicon_graphite(self):
+        self.heat_of_mixing_silicon_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)
+    def test_heat_of_mixing_DFN_silicon_graphite(self):
+        self.heat_of_mixing_silicon_graphite(pybamm.lithium_ion.DFN)
 
     @pytest.mark.parametrize(
         "model_class",