Repository: pybamm-team/PyBaMM

Issue:
Add contact resistance heat source
Add a heat source due to the parameter "Contact resistance [Ohm]"

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4192?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) Report
Attention: Patch coverage is `72.72727%` with `3 lines` in your changes missing coverage. Please review.
> Project coverage is 99.56%. Comparing base [(`5388515`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/5388515d81bae9bdc4bc63e937f46bbd8bbf7ce3?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) to head [(`8729b04`)](https://app.codecov.io/gh/pybamm-team/PyBaMM/commit/8729b044d50fc67b85d33afc6e1a2486ce805b52?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).
> Report is 3 commits behind head on develop.

| [Files](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4192?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) | Patch % | Lines |
|---|---|---|
| [pybamm/models/submodels/thermal/base\_thermal.py](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4192?src=pr&el=tree&filepath=pybamm%2Fmodels%2Fsubmodels%2Fthermal%2Fbase_thermal.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team#diff-cHliYW1tL21vZGVscy9zdWJtb2RlbHMvdGhlcm1hbC9iYXNlX3RoZXJtYWwucHk=) | 72.72% | [3 Missing :warning: ](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4192?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team) |

<details><summary>Additional details and impacted files</summary>


```diff
@@             Coverage Diff             @@
##           develop    #4192      +/-   ##
===========================================
- Coverage    99.58%   99.56%   -0.02%     
===========================================
  Files          288      288              
  Lines        21799    21809      +10     
===========================================
+ Hits         21708    21715       +7     
- Misses          91       94       +3     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/pybamm-team/PyBaMM/pull/4192?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pybamm-team).

rtimms: It would also be good to update [this notebook](https://github.com/pybamm-team/PyBaMM/blob/develop/docs/source/examples/notebooks/models/thermal-models.ipynb) to include contact resistance heating in the lumped model. We should also update the docs of the thermal models to be more explicit about which heating terms are included.


review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/pybamm-team/PyBaMM/pull/4192"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>

Patch:
diff --git a/docs/source/examples/notebooks/models/thermal-models.ipynb b/docs/source/examples/notebooks/models/thermal-models.ipynb
index 599a362b4..4120ac289 100644
--- a/docs/source/examples/notebooks/models/thermal-models.ipynb
+++ b/docs/source/examples/notebooks/models/thermal-models.ipynb
@@ -14,16 +14,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "\n",
-      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.3.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m23.3.2\u001b[0m\n",
-      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n",
       "Note: you may need to restart the kernel to use updated packages.\n"
      ]
     }
@@ -57,19 +54,21 @@
     "\n",
     "where $\\rho_k$ is the density, $c_{p,k}$ is the specific heat, and $L_k$ is the thickness of each component. The subscript $k \\in \\{cn, n, s, p, cp\\}$ is used to refer to the components negative current collector, negative electrode, separator, positive electrode, and positive current collector.\n",
     "\n",
-    "The heat source term accounts for Ohmic heating $Q_{Ohm,k}$ due to resistance in the solid and electrolyte, irreverisble heating due to electrochemical reactions $Q_{rxn,k}$, and reversible heating due to entropic changes in the the electrode $Q_{rev,k}$:\n",
+    "The heat source term accounts for Ohmic heating $Q_{Ohm,k}$ due to resistance in the solid and electrolyte, irreverisble heating due to electrochemical reactions $Q_{rxn,k}$, reversible heating due to entropic changes in the the electrode $Q_{rev,k}$, and heating due to contact resistance $Q_{cr}$:\n",
     "\n",
     "$$\n",
-    "Q = Q_{Ohm,k}+Q_{rxn,k}+Q_{rev,k},\n",
+    "Q = Q_{Ohm,k}+Q_{rxn,k}+Q_{rev,k}+Q_{cr},\n",
     "$$\n",
     "\n",
     "with\n",
     "\n",
     "$$ \n",
-    "Q_{Ohm,k} = -i_k \\nabla \\phi_k, \\quad Q_{rxn,k} = a_k j_k \\eta_k, \\quad Q_{rev,k} = a_k j_k T_k \\frac{\\partial U}{\\partial T} \\bigg|_{T=T_{\\infty}}.\n",
+    "Q_{Ohm,k} = -i_k \\nabla \\phi_k, \\quad Q_{rxn,k} = a_k j_k \\eta_k, \\quad Q_{rev,k} = a_k j_k T_k \\frac{\\partial U}{\\partial T} \\bigg|_{T=T_{\\infty}}, Q_{cr} = \\frac{R_{cr}}{V_{cell}}i_k^2.\n",
     "$$\n",
     "\n",
-    "Here $i_k$ is the current, $\\phi_k$ the potential, $a_k$ the surface area to volume ratio, $j_k$ the interfacial current density, $\\eta_k$ the overpotential, and $U$ the open-circuit potential. The averaged heat source term $\\bar{Q}$ is computed by taking the volume-average of $Q$.\n"
+    "Here $i_k$ is the current, $\\phi_k$ the potential, $a_k$ the surface area to volume ratio, $j_k$ the interfacial current density, $\\eta_k$ the overpotential, $U$ the open-circuit potential, $R_{cr}$ is the contact resistance, and $V_{cell}$ is the total cell volume. The averaged heat source term $\\bar{Q}$ is computed by taking the volume-average of $Q$.\n",
+    "\n",
+    "\n"
    ]
   },
   {
@@ -87,12 +86,14 @@
     "\n",
     "When using the option `{\"cell geometry\": \"pouch\"}` the parameter $A$ and $V$ are computed automatically from the pouch dimensions, assuming a single-layer pouch cell, i.e. $A$ is the total surface area of a single-layer pouch cell and $V$ is the volume. The parameter $h$ is still set by the \"Total heat transfer coefficient [W.m-2.K-1]\" parameter.\n",
     "\n",
+    "When using the option `{\"contact resistance\": \"true\"}` the parameter \"Contact resistance [Ohm]\" must be specified to calculate the heating from contact resistance, which corresponds to $R_{cr}$. \"Cell volume [m3]\" is $V_{cell}$ within the governing equation. The default lumped model option is `{\"contact resistance\": \"false\"}`.\n",
+    "\n",
     "The lumped thermal option can be selected as follows\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -113,7 +114,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 3,
    "metadata": {},
    "outputs": [
     {
@@ -130,6 +131,49 @@
     "print(\"Cell geometry:\", model.options[\"cell geometry\"])"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The contact resistance option can be turned on via the following"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "options = {\"thermal\": \"lumped\", \"contact resistance\": \"true\"}\n",
+    "model = pybamm.lithium_ion.DFN(options)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The default option for the lumped model does not include contact resistance."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "contact resistance: false\n"
+     ]
+    }
+   ],
+   "source": [
+    "options = {\"thermal\": \"lumped\"}\n",
+    "model = pybamm.lithium_ion.DFN(options)\n",
+    "print(\"contact resistance:\", model.options[\"contact resistance\"])"
+   ]
+  },
   {
    "attachments": {},
    "cell_type": "markdown",
@@ -431,6 +475,135 @@
     ")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the lumped model, we can compare how the contact resistance affects the heating. To do so, we must set the `\"contact resistance\"` option to `\"true\"` and update the `\"Contact resistance [Ohm]\"` parameter."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model_no_contact_resistance = pybamm.lithium_ion.SPMe(\n",
+    "    {\"cell geometry\": \"arbitrary\", \"thermal\": \"lumped\", \"contact resistance\": \"false\"},\n",
+    "    name=\"lumped thermal model\",\n",
+    ")\n",
+    "model_contact_resistance = pybamm.lithium_ion.SPMe(\n",
+    "    {\"cell geometry\": \"arbitrary\", \"thermal\": \"lumped\", \"contact resistance\": \"true\"},\n",
+    "    name=\"lumped thermal model with contact resistance\",\n",
+    ")\n",
+    "models = [model_no_contact_resistance, model_contact_resistance]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We then choose a parameter set."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "parameter_values = pybamm.ParameterValues(\"Marquis2019\")\n",
+    "lumped_params = parameter_values.copy()\n",
+    "lumped_params_contact_resistance = parameter_values.copy()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the contact resistance model, we must specify a contact resistance greater than zero. The default is zero."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lumped_params_contact_resistance.update(\n",
+    "    {\n",
+    "        \"Contact resistance [Ohm]\": 0.05,\n",
+    "    }\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The models and parameters are then used to solve for a 1C discharge."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "params = [lumped_params, lumped_params_contact_resistance]\n",
+    "sols = []\n",
+    "for model, param in zip(models, params):\n",
+    "    sim = pybamm.Simulation(model, parameter_values=param)\n",
+    "    sim.solve([0, 3600])\n",
+    "    sols.append(sim.solution)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We can then compare the voltage and cell temperature and see the impact of the contact resistance."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "1254ec88a920400096b25541577ee948",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<pybamm.plotting.quick_plot.QuickPlot at 0x338900950>"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "output_variables = [\n",
+    "    \"Voltage [V]\",\n",
+    "    \"X-averaged cell temperature [K]\",\n",
+    "    \"Cell temperature [K]\",\n",
+    "]\n",
+    "pybamm.dynamic_plot(sols, output_variables)"
+   ]
+  },
   {
    "attachments": {},
    "cell_type": "markdown",
@@ -481,7 +654,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.9"
   },
   "toc": {
    "base_numbering": 1,
diff --git a/pybamm/models/submodels/thermal/lumped.py b/pybamm/models/submodels/thermal/lumped.py
index 4afde4fa5..76af2904b 100644
--- a/pybamm/models/submodels/thermal/lumped.py
+++ b/pybamm/models/submodels/thermal/lumped.py
@@ -53,11 +53,25 @@ class Lumped(BaseThermal):
         V = variables["Cell thermal volume [m3]"]
         Q_cool_W = -self.param.h_total * (T_vol_av - T_amb) * self.param.A_cooling
         Q_cool_vol_av = Q_cool_W / V
+
+        # Contact resistance heating Q_cr
+        if self.options["contact resistance"] == "true":
+            I = variables["Current [A]"]
+            Q_cr_W = self.calculate_Q_cr_W(I, self.param.R_contact)
+            V = self.param.V_cell
+            Q_cr_vol_av = self.calculate_Q_cr_vol_av(I, self.param.R_contact, V)
+        else:
+            Q_cr_W = pybamm.Scalar(0)
+            Q_cr_vol_av = Q_cr_W
+
         variables.update(
             {
                 # Lumped cooling
                 "Lumped total cooling [W.m-3]": Q_cool_vol_av,
                 "Lumped total cooling [W]": Q_cool_W,
+                # Contact resistance
+                "Lumped contact resistance heating [W.m-3]": Q_cr_vol_av,
+                "Lumped contact resistance heating [W]": Q_cr_W,
             }
         )
         return variables
@@ -66,12 +80,22 @@ class Lumped(BaseThermal):
         T_vol_av = variables["Volume-averaged cell temperature [K]"]
         Q_vol_av = variables["Volume-averaged total heating [W.m-3]"]
         Q_cool_vol_av = variables["Lumped 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_cool_vol_av) / rho_c_p_eff_av}
+        self.rhs = {T_vol_av: (Q_vol_av + Q_cr_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]"]
         self.initial_conditions = {T_vol_av: self.param.T_init}
+
+    def calculate_Q_cr_W(self, current, contact_resistance):
+        Q_cr_W = current**2 * contact_resistance
+        return Q_cr_W
+
+    def calculate_Q_cr_vol_av(self, current, contact_resistance, volume):
+        Q_cr_W = self.calculate_Q_cr_W(current, contact_resistance)
+        Q_cr_vol_av = Q_cr_W / volume
+        return Q_cr_vol_av


Test patch:
diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
index 4d2f50e8d..7de36da04 100644
--- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
+++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py
@@ -103,6 +103,54 @@ class TestThermal(TestCase):
 
         self.assertGreater(1e-5, err(solutions["SPMe 1+1D"], solutions["SPMe 2+1D"]))
 
+    def test_lumped_contact_resistance(self):
+        # Test that the heating with contact resistance is greater than without
+
+        # load models
+        model_no_contact_resistance = pybamm.lithium_ion.SPMe(
+            {
+                "cell geometry": "arbitrary",
+                "thermal": "lumped",
+                "contact resistance": "false",
+            }
+        )
+        model_contact_resistance = pybamm.lithium_ion.SPMe(
+            {
+                "cell geometry": "arbitrary",
+                "thermal": "lumped",
+                "contact resistance": "true",
+            }
+        )
+        models = [model_no_contact_resistance, model_contact_resistance]
+
+        # parameters
+        parameter_values = pybamm.ParameterValues("Marquis2019")
+        lumped_params = parameter_values.copy()
+        lumped_params_contact_resistance = parameter_values.copy()
+
+        lumped_params_contact_resistance.update(
+            {
+                "Contact resistance [Ohm]": 0.05,
+            }
+        )
+
+        # solve the models
+        params = [lumped_params, lumped_params_contact_resistance]
+        sols = []
+        for model, param in zip(models, params):
+            sim = pybamm.Simulation(model, parameter_values=param)
+            sim.solve([0, 3600])
+            sols.append(sim.solution)
+
+        # get the average temperature from each model
+        avg_cell_temp = sols[0]["X-averaged cell temperature [K]"].entries
+        avg_cell_temp_cr = sols[1]["X-averaged cell temperature [K]"].entries
+
+        # check that the cell temperature of the lumped thermal model
+        # with contact resistance is higher than without contact resistance
+        # skip the first entry because they are the same due to initial conditions
+        np.testing.assert_array_less(avg_cell_temp[1:], avg_cell_temp_cr[1:])
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
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 7d190875a..7e1f2d5ca 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
@@ -138,7 +138,10 @@ class BaseUnitTestLithiumIon:
         self.check_well_posedness(options)
 
     def test_well_posed_contact_resistance(self):
-        options = {"contact resistance": "true"}
+        options = {
+            "contact resistance": "true",
+            "thermal": "lumped",
+        }
         self.check_well_posedness(options)
 
     def test_well_posed_particle_uniform(self):