Repository: pybamm-team/PyBaMM (python)

Issue:
Add option to include heat of mixing
Add an option to include heat of mixing (irreversible) and use the expression for the reversible heating from [this paper](https://arxiv.org/abs/2104.01590)

Pull request description:
Issue 1778 - heat of mixing
# Description
implemented heat of mixing

Please include a summary of the change and which issue is fixed. Please also include relevant motivation and context. List any dependencies that are required for this change.

Fixes #1778 

## Type of change

Please add a line in the relevant section of [CHANGELOG.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CHANGELOG.md) to document the change (include PR #) - note reverse order of PR #s. If necessary, also add to the list of breaking changes.

- [x] New feature (non-breaking change which adds functionality)
- [ ] Optimization (back-end change that speeds up the code)
- [ ] Bug fix (non-breaking change which fixes an issue)

# Key checklist:

- [ ] No style issues: `$ pre-commit run` (see [CONTRIBUTING.md](https://github.com/pybamm-team/PyBaMM/blob/develop/CONTRIBUTING.md#installing-and-using-pre-commit) for how to set this up to run automatically when committing locally, in just two lines of code)
- [ ] All tests pass: `$ python run-tests.py --all`
- [ ] The documentation builds: `$ python run-tests.py --doctest`

You can run unit and doctests together at once, using `$ python run-tests.py --quick`.

## Further checks:

- [ ] Code is commented, particularly in hard-to-understand areas
- [ ] Tests added that prove fix is effective or that feature works


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8dfb4a7a1..4e4638dff 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,7 +4,7 @@
 
 - Added new parameters `"f{pref]Initial inner SEI on cracks thickness [m]"` and `"f{pref]Initial outer SEI on cracks thickness [m]"`, instead of hardcoding these to `L_inner_0 / 10000` and `L_outer_0 / 10000`. ([#4168](https://github.com/pybamm-team/PyBaMM/pull/4168))
 - Added `pybamm.DataLoader` class to fetch data files from [pybamm-data](https://github.com/pybamm-team/pybamm-data/releases/tag/v1.0.0) and store it under local cache. ([#4098](https://github.com/pybamm-team/PyBaMM/pull/4098))
-- Transport efficiency submodel has new options from the literature relating to different tortuosity factor models and also a new option called "tortuosity factor" for specifying the value or function directly as parameters ([#3437](https://github.com/pybamm-team/PyBaMM/pull/3437))
+- Added `time` as an option for `Experiment.termination`. Now allows solving up to a user-specified time while also allowing different cycles and steps in an experiment to be handled normally. ([#4073](https://github.com/pybamm-team/PyBaMM/pull/4073))
 - Added `plot_thermal_components` to plot the contributions to the total heat generation in a battery ([#4021](https://github.com/pybamm-team/PyBaMM/pull/4021))
 - Added functions for normal probability density function (`pybamm.normal_pdf`) and cumulative distribution function (`pybamm.normal_cdf`) ([#3999](https://github.com/pybamm-team/PyBaMM/pull/3999))
 - "Basic" models are now compatible with experiments ([#3995](https://github.com/pybamm-team/PyBaMM/pull/3995))
@@ -20,7 +20,8 @@
 - Add support for BPX version 0.4.0 which allows for blended electrodes and user-defined parameters in BPX([#3414](https://github.com/pybamm-team/PyBaMM/pull/3414))
 - Added `by_submodel` feature in `print_parameter_info` method to allow users to print parameters and types of submodels in a tabular and readable format ([#3628](https://github.com/pybamm-team/PyBaMM/pull/3628))
 - Added `WyciskOpenCircuitPotential` for differential capacity hysteresis state open-circuit potential submodel ([#3593](https://github.com/pybamm-team/PyBaMM/pull/3593))
-- Added `time` as an option for `Experiment.termination`. Now allows solving up to a user-specified time while also allowing different cycles and steps in an experiment to be handled normally. ([#4073](https://github.com/pybamm-team/PyBaMM/pull/4073))
+- Transport efficiency submodel has new options from the literature relating to different tortuosity factor models and also a new option called "tortuosity factor" for specifying the value or function directly as parameters ([#3437](https://github.com/pybamm-team/PyBaMM/pull/3437))
+- Heat of mixing source term can now be included into thermal models ([#2837](https://github.com/pybamm-team/PyBaMM/pull/2837))
 
 ## Bug Fixes
 
diff --git a/examples/scripts/compare_lithium_ion_heat_of_mixing.py b/examples/scripts/compare_lithium_ion_heat_of_mixing.py
new file mode 100644
index 000000000..455f5b1e8
--- /dev/null
+++ b/examples/scripts/compare_lithium_ion_heat_of_mixing.py
@@ -0,0 +1,112 @@
+#
+# Compare SPMe model with and without heat of mixing
+#
+import pybamm
+import numpy as np
+import matplotlib.pyplot as plt
+
+pybamm.set_logging_level("INFO")
+
+# load models
+models = [
+    pybamm.lithium_ion.SPMe(
+        {"heat of mixing": "true", "thermal": "lumped"}, name="SPMe with heat of mixing"
+    ),
+    pybamm.lithium_ion.SPMe({"thermal": "lumped"}, name="SPMe without heat of mixing"),
+]
+
+# set parametrisation (Ecker et al., 2015)
+parameter_values = pybamm.ParameterValues("Ecker2015")
+
+# set mesh
+# (increased number of points in particles to avoid oscillations for high C-rates)
+var_pts = {"x_n": 16, "x_s": 8, "x_p": 16, "r_n": 128, "r_p": 128}
+
+# set the constant current discharge C-rate
+C_rate = 5
+
+# set the simulation time and increase the number of points
+# for better integration in time
+t_eval = np.linspace(0, 720, 360)
+
+# set up an extra plot with the heat of mixing vs time in each electrode and
+# the integrated heat of mixing vs time in each electrode to compare with
+# Figure 6(a) from Richardson et al. (2021)
+fig, axs = plt.subplots(2, len(models), figsize=(12, 7))
+
+# extract some of the parameters
+L_n = parameter_values["Negative electrode thickness [m]"]
+L_p = parameter_values["Positive electrode thickness [m]"]
+A = parameter_values["Electrode width [m]"] * parameter_values["Electrode height [m]"]
+
+# create and run simulations
+sims = []
+for m, model in enumerate(models):
+    sim = pybamm.Simulation(
+        model, parameter_values=parameter_values, var_pts=var_pts, C_rate=C_rate
+    )
+    sim.solve(t_eval)
+    sims.append(sim)
+
+    # extract solution for plotting
+    sol = sim.solution
+
+    # extract variables from the solution
+    time = sol["Time [h]"].entries
+    Q_mix = sol["Heat of mixing [W.m-3]"].entries
+
+    # heat of mixing in negative and positive electrodes multiplied by the electrode
+    # width, represents the integral of heat of mixing term across each of the
+    # electrodes (W.m-2)
+    Q_mix_n = Q_mix[0, :] * L_n
+    Q_mix_p = Q_mix[-1, :] * L_p
+
+    # heat of mixing integrals (J.m-2)
+    Q_mix_n_int = 0.0
+    Q_mix_p_int = 0.0
+
+    # data for plotting
+    Q_mix_n_plt = []
+    Q_mix_p_plt = []
+
+    # performs integration in time
+    for i, t in enumerate(time[1:]):
+        dt = (t - time[i]) * 3600  # seconds
+        Q_mix_n_avg = (Q_mix_n[i] + Q_mix_n[i + 1]) * 0.5
+        Q_mix_p_avg = (Q_mix_p[i] + Q_mix_p[i + 1]) * 0.5
+        # convert J to kJ and divide the integral by the electrode area A to compare
+        # with Figure 6(a) from Richardson et al. (2021)
+        Q_mix_n_int += Q_mix_n_avg * dt / 1000 / A
+        Q_mix_p_int += Q_mix_p_avg * dt / 1000 / A
+        Q_mix_n_plt.append(Q_mix_n_int)
+        Q_mix_p_plt.append(Q_mix_p_int)
+
+    # plots heat of mixing in each electrode vs time in minutes
+    axs[0, m].plot(time * 60, Q_mix_n, ls="-", label="Negative electrode")
+    axs[0, m].plot(time * 60, Q_mix_p, ls="--", label="Positive electrode")
+    axs[0, m].set_title(f"{model.name}")
+    axs[0, m].set_xlabel("Time [min]")
+    axs[0, m].set_ylabel("Heat of mixing [W.m-2]")
+    axs[0, m].grid(True)
+    axs[0, m].legend()
+
+    # plots integrated heat of mixing in each electrode vs time in minutes
+    axs[1, m].plot(time[1:] * 60, Q_mix_n_plt, ls="-", label="Negative electrode")
+    axs[1, m].plot(time[1:] * 60, Q_mix_p_plt, ls="--", label="Positive electrode")
+    axs[1, m].set_xlabel("Time [min]")
+    axs[1, m].set_ylabel("Integrated heat of mixing [kJ.m-2]")
+    axs[1, m].grid(True)
+    axs[1, m].legend()
+
+# plot
+pybamm.dynamic_plot(
+    sims,
+    output_variables=[
+        "X-averaged cell temperature [K]",
+        "X-averaged heat of mixing [W.m-3]",
+        "X-averaged total heating [W.m-3]",
+        "Heat of mixing [W.m-3]",
+        "Voltage [V]",
+        "Current [A]",
+    ],
+)
diff --git a/pybamm/CITATIONS.bib b/pybamm/CITATIONS.bib
index 00d41fadd..33ac6055d 100644
--- a/pybamm/CITATIONS.bib
+++ b/pybamm/CITATIONS.bib
@@ -418,6 +418,17 @@
   doi = {10.1016/j.electacta.2020.135862},
 }
 
+@article{Richardson2021,
+  title = {Heat Generation and a Conservation Law for Chemical Energy in Li-ion Batteries},
+  author = {Richardson, Giles and Korotkin, Ivan},
+  journal = {Electrochimica Acta},
+  volume = {392},
+  pages = {138909},
+  year = {2021},
+  publisher = {Elsevier},
+  doi = {10.1016/j.electacta.2021.138909},
+}
+
 @article{Sripad2020,
   title={Kinetics of lithium electrodeposition and stripping},
   author={Sripad, Shashank and Korff, Daniel and DeCaluwe, Steven C and Viswanathan, Venkatasubramanian},
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 03114a59e..fb8d001f9 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -238,6 +238,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "integrated",
             ],
             "exchange-current density": ["single", "current sigmoid"],
+            "heat of mixing": ["false", "true"],
             "hydrolysis": ["false", "true"],
             "intercalation kinetics": [
                 "symmetric Butler-Volmer",
@@ -512,6 +513,11 @@ class BatteryModelOptions(pybamm.FuzzyDict):
 
         # Options not yet compatible with particle-size distributions
         if options["particle size"] == "distribution":
+            if options["heat of mixing"] != "false":
+                raise NotImplementedError(
+                    "Heat of mixing submodels do not yet support particle-size "
+                    "distributions."
+                )
             if options["lithium plating"] != "none":
                 raise NotImplementedError(
                     "Lithium plating submodels do not yet support particle-size "
@@ -1240,7 +1246,10 @@ class BaseBatteryModel(pybamm.BaseModel):
             if self.options["dimensionality"] == 0:
                 thermal_submodel = pybamm.thermal.pouch_cell.OneDimensionalX
 
-        self.submodels["thermal"] = thermal_submodel(self.param, self.options)
+        x_average = getattr(self, "x_average", False)
+        self.submodels["thermal"] = thermal_submodel(
+            self.param, self.options, x_average
+        )
 
     def set_current_collector_submodel(self):
         if self.options["current collector"] in ["uniform"]:
diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py
index 65a3f72a5..c5ebbc7db 100644
--- a/pybamm/models/submodels/thermal/base_thermal.py
+++ b/pybamm/models/submodels/thermal/base_thermal.py
@@ -2,6 +2,7 @@
 # Base class for thermal effects
 #
 import pybamm
+import numpy as np
 
 
 class BaseThermal(pybamm.BaseSubModel):
@@ -16,8 +17,12 @@ class BaseThermal(pybamm.BaseSubModel):
         A dictionary of options to be passed to the model.
     """
 
-    def __init__(self, param, options=None):
+    def __init__(self, param, options=None, x_average=False):
         super().__init__(param, options=options)
+        self.x_average = x_average
+
+        if self.options["heat of mixing"] == "true":
+            pybamm.citations.register("Richardson2021")
 
     def _get_standard_fundamental_variables(self, T_dict):
         """
@@ -175,8 +180,12 @@ class BaseThermal(pybamm.BaseSubModel):
             Q_rev_n, pybamm.FullBroadcast(0, "separator", "current collector"), Q_rev_p
         )
 
+        # Heat of mixing
+        Q_mix_s_n, Q_mix_s_s, Q_mix_s_p = self._heat_of_mixing(variables)
+        Q_mix = pybamm.concatenation(Q_mix_s_n, Q_mix_s_s, Q_mix_s_p)
+
         # Total heating
-        Q = Q_ohm + Q_rxn + Q_rev
+        Q = Q_ohm + Q_rxn + Q_rev + Q_mix
 
         # 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
@@ -184,6 +193,7 @@ class BaseThermal(pybamm.BaseSubModel):
         Q_ohm_av = self._x_average(Q_ohm, Q_ohm_s_cn, Q_ohm_s_cp)
         Q_rxn_av = self._x_average(Q_rxn, 0, 0)
         Q_rev_av = self._x_average(Q_rev, 0, 0)
+        Q_mix_av = self._x_average(Q_mix, 0, 0)
         Q_av = self._x_average(Q, Q_ohm_s_cn, Q_ohm_s_cp)
 
         # Compute the integrated heat source per unit simulated electrode-pair area
@@ -192,11 +202,14 @@ class BaseThermal(pybamm.BaseSubModel):
         Q_ohm_Wm2 = Q_ohm_av * param.L
         Q_rxn_Wm2 = Q_rxn_av * param.L
         Q_rev_Wm2 = Q_rev_av * param.L
+        Q_mix_Wm2 = Q_mix_av * param.L
         Q_Wm2 = Q_av * param.L
+
         # Now average over the electrode height and width
         Q_ohm_Wm2_av = self._yz_average(Q_ohm_Wm2)
         Q_rxn_Wm2_av = self._yz_average(Q_rxn_Wm2)
         Q_rev_Wm2_av = self._yz_average(Q_rev_Wm2)
+        Q_mix_Wm2_av = self._yz_average(Q_mix_Wm2)
         Q_Wm2_av = self._yz_average(Q_Wm2)
 
         # Compute total heat source terms (in W) over the *entire cell volume*, not
@@ -208,6 +221,7 @@ class BaseThermal(pybamm.BaseSubModel):
         Q_ohm_W = Q_ohm_Wm2_av * n_elec * A
         Q_rxn_W = Q_rxn_Wm2_av * n_elec * A
         Q_rev_W = Q_rev_Wm2_av * n_elec * A
+        Q_mix_W = Q_mix_Wm2_av * n_elec * A
         Q_W = Q_Wm2_av * n_elec * A
 
         # Compute volume-averaged heat source terms over the *entire cell volume*, not
@@ -216,6 +230,7 @@ class BaseThermal(pybamm.BaseSubModel):
         Q_ohm_vol_av = Q_ohm_W / V
         Q_rxn_vol_av = Q_rxn_W / V
         Q_rev_vol_av = Q_rev_W / V
+        Q_mix_vol_av = Q_mix_W / V
         Q_vol_av = Q_W / V
 
         # Effective heat capacity
@@ -228,6 +243,7 @@ class BaseThermal(pybamm.BaseSubModel):
                 "Ohmic heating [W.m-3]": Q_ohm,
                 "X-averaged Ohmic heating [W.m-3]": Q_ohm_av,
                 "Volume-averaged Ohmic heating [W.m-3]": Q_ohm_vol_av,
+                "Volume-averaged heat of mixing [W.m-3]": Q_mix_vol_av,
                 "Ohmic heating per unit electrode-pair area [W.m-2]": Q_ohm_Wm2,
                 "Ohmic heating [W]": Q_ohm_W,
                 # Irreversible
@@ -244,6 +260,12 @@ class BaseThermal(pybamm.BaseSubModel):
                 "Volume-averaged reversible heating [W.m-3]": Q_rev_vol_av,
                 "Reversible heating per unit electrode-pair area " "[W.m-2]": Q_rev_Wm2,
                 "Reversible heating [W]": Q_rev_W,
+                # Mixing
+                "Heat of mixing [W.m-3]": Q_mix,
+                "X-averaged heat of mixing [W.m-3]": Q_mix_av,
+                "Volume-averaged heating of mixing [W.m-3]": Q_mix_vol_av,
+                "Heat of mixing per unit electrode-pair area " "[W.m-2]": Q_mix_Wm2,
+                "Heat of mixing [W]": Q_mix_W,
                 # Total
                 "Total heating [W.m-3]": Q,
                 "X-averaged total heating [W.m-3]": Q_av,
@@ -290,6 +312,74 @@ class BaseThermal(pybamm.BaseSubModel):
                 Q_s_cp = self.param.p.sigma_cc * pybamm.grad_squared(phi_s_cp)
         return Q_s_cn, Q_s_cp
 
+    def _heat_of_mixing(self, variables):
+        """Compute heat of mixing source terms."""
+        param = self.param
+
+        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 = param.n.prim.D(c_n, T_n_part)
+                dUeq_n = param.n.prim.U(c_n / 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
+
+            # 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 = param.p.prim.D(c_p, T_p_part)
+            dUeq_p = param.p.prim.U(c_p / 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_s = pybamm.FullBroadcast(0, ["separator"], "current collector")
+        else:
+            Q_mix_s_n = pybamm.FullBroadcast(
+                0, ["negative electrode"], "current collector"
+            )
+            Q_mix_s_p = pybamm.FullBroadcast(
+                0, ["positive electrode"], "current collector"
+            )
+            Q_mix_s_s = pybamm.FullBroadcast(0, ["separator"], "current collector")
+
+        return Q_mix_s_n, Q_mix_s_s, Q_mix_s_p
+
     def _x_average(self, var, var_cn, var_cp):
         """
         Computes the X-average over the whole cell (including current collectors)
diff --git a/pybamm/models/submodels/thermal/isothermal.py b/pybamm/models/submodels/thermal/isothermal.py
index e3228a55c..edcd47bbd 100644
--- a/pybamm/models/submodels/thermal/isothermal.py
+++ b/pybamm/models/submodels/thermal/isothermal.py
@@ -18,8 +18,8 @@ class Isothermal(BaseThermal):
         A dictionary of options to be passed to the model.
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=False):
+        super().__init__(param, options=options, x_average=x_average)
 
     def get_fundamental_variables(self):
         # Set the x-averaged temperature to the ambient temperature, which can be
diff --git a/pybamm/models/submodels/thermal/lumped.py b/pybamm/models/submodels/thermal/lumped.py
index 66cd5c583..4afde4fa5 100644
--- a/pybamm/models/submodels/thermal/lumped.py
+++ b/pybamm/models/submodels/thermal/lumped.py
@@ -20,8 +20,8 @@ class Lumped(BaseThermal):
 
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=False):
+        super().__init__(param, options=options, x_average=x_average)
         pybamm.citations.register("Timms2021")
 
     def get_fundamental_variables(self):
diff --git a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py
index 548c6cba0..a4908c6f5 100644
--- a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py
+++ b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py
@@ -22,8 +22,8 @@ class CurrentCollector1D(BaseThermal):
 
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=True):
+        super().__init__(param, options=options, x_average=x_average)
         pybamm.citations.register("Timms2021")
 
     def get_fundamental_variables(self):
diff --git a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
index 206a8d1e8..7955ee4c3 100644
--- a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
+++ b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
@@ -22,8 +22,8 @@ class CurrentCollector2D(BaseThermal):
 
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=True):
+        super().__init__(param, options=options, x_average=x_average)
         pybamm.citations.register("Timms2021")
 
     def get_fundamental_variables(self):
diff --git a/pybamm/models/submodels/thermal/pouch_cell/x_full.py b/pybamm/models/submodels/thermal/pouch_cell/x_full.py
index 98065b2d5..f4aa07c56 100644
--- a/pybamm/models/submodels/thermal/pouch_cell/x_full.py
+++ b/pybamm/models/submodels/thermal/pouch_cell/x_full.py
@@ -24,8 +24,8 @@ class OneDimensionalX(BaseThermal):
 
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=False):
+        super().__init__(param, options=options, x_average=x_average)
         pybamm.citations.register("Timms2021")
 
     def get_fundamental_variables(self):
diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py
index 8d69ec49c..f5a76c6d4 100644
--- a/pybamm/parameters/lithium_ion_parameters.py
+++ b/pybamm/parameters/lithium_ion_parameters.py
@@ -643,12 +643,14 @@ class ParticleLithiumIonParameters(BaseParameters):
         u_ref = pybamm.FunctionParameter(
             f"{self.phase_prefactor}{Domain} electrode {lithiation}OCP [V]", inputs
         )
+
+        dudt_func = self.dUdT(sto)
+        u_ref = u_ref + (T - self.main_param.T_ref) * dudt_func
+
         # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1
         # this will not affect the OCP for most values of sto
         # see #1435
-        u_ref = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1))
-        dudt_func = self.dUdT(sto)
-        out = u_ref + (T - self.main_param.T_ref) * dudt_func
+        out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1))
 
         if self.domain == "negative":
             out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n}, T)"


Test patch:
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index 557c99eb9..caff0cda5 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -26,6 +26,7 @@ PRINT_OPTIONS_OUTPUT = """\
 'dimensionality': 0 (possible: [0, 1, 2])
 'electrolyte conductivity': 'default' (possible: ['default', 'full', 'leading order', 'composite', 'integrated'])
 'exchange-current density': 'single' (possible: ['single', 'current sigmoid'])
+'heat of mixing': 'false' (possible: ['false', 'true'])
 'hydrolysis': 'false' (possible: ['false', 'true'])
 'intercalation kinetics': 'symmetric Butler-Volmer' (possible: ['symmetric Butler-Volmer', 'asymmetric Butler-Volmer', 'linear', 'Marcus', 'Marcus-Hush-Chidsey', 'MSMR'])
 'interface utilisation': 'full' (possible: ['full', 'constant', 'current-driven'])
@@ -373,6 +374,15 @@ class TestBaseBatteryModel(TestCase):
                 }
             )
 
+        # thermal heat of mixing
+        with self.assertRaisesRegex(NotImplementedError, "Heat of mixing"):
+            pybamm.BaseBatteryModel(
+                {
+                    "heat of mixing": "true",
+                    "particle size": "distribution",
+                }
+            )
+
         # phases
         with self.assertRaisesRegex(pybamm.OptionError, "multiple particle phases"):
             pybamm.BaseBatteryModel({"particle phases": "2", "surface form": "false"})
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py
index 24ccd55d7..50fac7a04 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py
@@ -82,3 +82,7 @@ class BaseUnitTestLithiumIonHalfCell:
     def test_well_posed_lumped_thermal(self):
         options = {"thermal": "lumped"}
         self.check_well_posedness(options)
+
+    def test_well_posed_lumped_thermal_hom(self):
+        options = {"thermal": "lumped", "heat of mixing": "true"}
+        self.check_well_posedness(options)
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index 6b8a9d07e..7d190875a 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -67,6 +67,76 @@ class BaseUnitTestLithiumIon:
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_isothermal_heat_source_hom(self):
+        options = {
+            "calculate heat source for isothermal models": "true",
+            "thermal": "isothermal",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_2plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 1,
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 2,
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_lumped_thermal_model_1D_hom(self):
+        options = {"thermal": "lumped", "heat of mixing": "true"}
+        self.check_well_posedness(options)
+
+    def test_well_posed_x_full_thermal_model_hom(self):
+        options = {
+            "thermal": "x-full",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_lumped_thermal_1plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 1,
+            "thermal": "lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_lumped_thermal_2plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 2,
+            "thermal": "lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_thermal_1plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 1,
+            "thermal": "x-lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_thermal_2plus1D_hom(self):
+        options = {
+            "current collector": "potential pair",
+            "dimensionality": 2,
+            "thermal": "x-lumped",
+            "heat of mixing": "true",
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_contact_resistance(self):
         options = {"contact resistance": "true"}
         self.check_well_posedness(options)