add plot_thermal_components

pybamm/__init__.py

@@ -180,6 +180,7 @@
 from .plotting.plot import plot
 from .plotting.plot2D import plot2D
 from .plotting.plot_voltage_components import plot_voltage_components
+from .plotting.plot_thermal_components import plot_thermal_components
 from .plotting.plot_summary_variables import plot_summary_variables
 from .plotting.dynamic_plot import dynamic_plot
 

pybamm/models/submodels/thermal/base_thermal.py

@@ -218,6 +218,10 @@ def _get_standard_coupled_variables(self, variables):
         Q_rev_vol_av = Q_rev_W / V
         Q_vol_av = Q_W / V
 
+        # Effective heat capacity
+        T_vol_av = variables["Volume-averaged cell temperature [K]"]
+        rho_c_p_eff_av = param.rho_c_p_eff(T_vol_av)
+
         variables.update(
             {
                 # Ohmic
@@ -249,6 +253,9 @@ def _get_standard_coupled_variables(self, variables):
                 # Current collector
                 "Negative current collector Ohmic heating [W.m-3]": Q_ohm_s_cn,
                 "Positive current collector Ohmic heating [W.m-3]": Q_ohm_s_cp,
+                # Effective heat capacity
+                "Volume-averaged effective heat capacity [J.K-1.m-3]": rho_c_p_eff_av,
+                "Cell thermal volume [m3]": V,
             }
         )
         return variables

pybamm/models/submodels/thermal/isothermal.py

@@ -73,6 +73,8 @@ def get_coupled_variables(self, variables):
                 "Total heating [W]",
                 "Negative current collector Ohmic heating [W.m-3]",
                 "Positive current collector Ohmic heating [W.m-3]",
+                "Lumped total cooling [W.m-3]",
+                "Lumped total cooling [W]",
             ]:
                 # All variables are zero
                 variables.update({var: zero})

pybamm/models/submodels/thermal/lumped.py

@@ -46,23 +46,31 @@ def get_fundamental_variables(self):
 
     def get_coupled_variables(self, variables):
         variables.update(self._get_standard_coupled_variables(variables))
+
+        # Newton cooling, accounting for surface area to volume ratio
+        T_vol_av = variables["Volume-averaged cell temperature [K]"]
+        T_amb = variables["Volume-averaged ambient temperature [K]"]
+        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
+        variables.update(
+            {
+                # Lumped cooling
+                "Lumped total cooling [W.m-3]": Q_cool_vol_av,
+                "Lumped total cooling [W]": Q_cool_W,
+            }
+        )
         return variables
 
     def set_rhs(self, variables):
         T_vol_av = variables["Volume-averaged cell temperature [K]"]
         Q_vol_av = variables["Volume-averaged total heating [W.m-3]"]
-        T_amb = variables["Volume-averaged ambient temperature [K]"]
+        Q_cool_vol_av = variables["Lumped total cooling [W.m-3]"]
+        rho_c_p_eff_av = variables[
+            "Volume-averaged effective heat capacity [J.K-1.m-3]"
+        ]
 
-        # Newton cooling, accounting for surface area to volume ratio
-        cell_surface_area = self.param.A_cooling
-        cell_volume = self.param.V_cell
-        Q_cool_vol_av = (
-            -self.param.h_total * (T_vol_av - T_amb) * cell_surface_area / cell_volume
-        )
-
-        self.rhs = {
-            T_vol_av: (Q_vol_av + Q_cool_vol_av) / self.param.rho_c_p_eff(T_vol_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]"]

pybamm/plotting/plot_thermal_components.py

@@ -0,0 +1,118 @@
+#
+# Method for plotting voltage components
+#
+
+from pybamm.util import import_optional_dependency
+from pybamm.simulation import Simulation
+from pybamm.solvers.solution import Solution
+from scipy.integrate import cumulative_trapezoid
+
+
+def plot_thermal_components(
+    input_data,
+    ax=None,
+    show_legend=True,
+    split_by_electrode=False,
+    show_plot=True,
+    **kwargs_fill,
+):
+    """
+    Generate a plot showing the component overpotentials that make up the voltage
+
+    Parameters
+    ----------
+    input_data : :class:`pybamm.Solution` or :class:`pybamm.Simulation`
+        Solution or Simulation object from which to extract voltage components.
+    ax : matplotlib Axis, optional
+        The axis on which to put the plot. If None, a new figure and axis is created.
+    show_legend : bool, optional
+        Whether to display the legend. Default is True
+    show_plot : bool, optional
+        Whether to show the plots. Default is True. Set to False if you want to
+        only display the plot after plt.show() has been called.
+    kwargs_fill
+        Keyword arguments, passed to ax.fill_between
+
+    """
+    # Check if the input is a Simulation and extract Solution
+    if isinstance(input_data, Simulation):
+        solution = input_data.solution
+    elif isinstance(input_data, Solution):
+        solution = input_data
+    plt = import_optional_dependency("matplotlib.pyplot")
+
+    # Set a default value for alpha, the opacity
+    kwargs_fill = {"alpha": 0.6, **kwargs_fill}
+
+    if ax is not None:
+        fig = None
+        show_plot = False
+    else:
+        fig, ax = plt.subplots(1, 2, figsize=(12, 4))
+
+    time_s = solution["Time [s]"].entries
+    time_h = time_s / 3600
+    T = solution["X-averaged cell temperature [K]"].entries
+    volume = solution["Cell thermal volume [m3]"].entries
+    rho_c_p_eff_av = solution[
+        "Volume-averaged effective heat capacity [J.K-1.m-3]"
+    ].entries
+
+    heating_sources = [
+        "Lumped total cooling",
+        "Ohmic heating",
+        "Irreversible electrochemical heating",
+        "Reversible heating",
+    ]
+    try:
+        heats_volumetric = {
+            name: solution[name + " [W]"].entries / volume for name in heating_sources
+        }
+    except KeyError as err:
+        raise NotImplementedError(
+            "plot_thermal_components is only implemented for lumped models"
+        ) from err
+
+    dTs = {
+        name: cumulative_trapezoid(heat / rho_c_p_eff_av, time_s, initial=0)
+        for name, heat in heats_volumetric.items()
+    }
+
+    # Plot
+    # Initialise
+    total_heat = 0
+    bottom_heat = heats_volumetric["Lumped total cooling"]
+    bottom_temp = T[0] + dTs["Lumped total cooling"]
+    # Plot components
+    for name in heating_sources:
+        top_temp = bottom_temp + abs(dTs[name])
+        ax[0].fill_between(time_h, bottom_temp, top_temp, **kwargs_fill, label=name)
+        bottom_temp = top_temp
+
+        top_heat = bottom_heat + abs(heats_volumetric[name])
+        ax[1].fill_between(time_h, bottom_heat, top_heat, **kwargs_fill, label=name)
+        bottom_heat = top_heat
+        total_heat += heats_volumetric[name]
+
+    ax[0].plot(time_h, T, "k--", label="Cell temperature")
+    ax[1].plot(time_h, total_heat, "k--", label="Total heating")
+
+    if show_legend:
+        leg = ax[1].legend(loc="center left", bbox_to_anchor=(1.05, 0.5), frameon=True)
+        leg.get_frame().set_edgecolor("k")
+
+    # Labels
+    for a in ax:
+        a.set_xlabel("Time [h]")
+        a.set_xlim([time_h[0], time_h[-1]])
+
+    ax[0].set_title("Temperatures [K]")
+    ax[1].set_title("Heat sources [W/m$^3$]")
+
+    if fig is not None:
+        fig.tight_layout()
+
+    if show_plot:  # pragma: no cover
+        plt.show()
+
+    return fig, ax

tests/unit/test_plotting/test_plot_thermal_components.py

@@ -0,0 +1,44 @@
+import pybamm
+import unittest
+import numpy as np
+from tests import TestCase
+import matplotlib.pyplot as plt
+from matplotlib import use
+
+use("Agg")
+
+
+class TestPlotThermalComponents(TestCase):
+    def test_plot_with_solution(self):
+        model = pybamm.lithium_ion.SPM({"thermal": "lumped"})
+        sim = pybamm.Simulation(model)
+        sol = sim.solve([0, 3600])
+        for input_data in [sim, sol]:
+            _, ax = pybamm.plot_thermal_components(input_data, show_plot=False)
+            t, T = ax[0].get_lines()[-1].get_data()
+            np.testing.assert_array_almost_equal(t, sol["Time [h]"].data)
+            np.testing.assert_array_almost_equal(
+                T, sol["X-averaged cell temperature [K]"].data
+            )
+
+            _, ax = plt.subplots(1, 2)
+            _, ax_out = pybamm.plot_thermal_components(sol, ax=ax, show_legend=True)
+            self.assertEqual(ax_out[0], ax[0])
+            self.assertEqual(ax_out[1], ax[1])
+
+    def test_not_implemented(self):
+        model = pybamm.lithium_ion.SPM({"thermal": "x-full"})
+        sim = pybamm.Simulation(model)
+        sol = sim.solve([0, 3600])
+        with self.assertRaises(NotImplementedError):
+            pybamm.plot_thermal_components(sol)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()