Surface temperature model (pybamm-team#4203)

* add surface temperature equal to ambient temperature

* pybamm-team#4022 surface model works

* update temperature comparison script

* pybamm-team#4022 update tests and example, rename parameters

* fix plot_thermal_components

* fix unit tests

* pybamm-team#4022 fix lead acid thermal tests

* pybamm-team#4022 fix example

* coverage

* fix test

* style

* Moving missed files

* Rob's suggestion

---------

Co-authored-by: Eric G. Kratz <[email protected]>
Co-authored-by: kratman <[email protected]>

examples/scripts/compare_surface_temperature.py

@@ -0,0 +1,61 @@
+#
+# Compare lithium-ion battery models
+#
+import pybamm
+
+pybamm.set_logging_level("INFO")
+
+# load models
+models = [
+    pybamm.lithium_ion.SPMe(
+        {"thermal": "lumped", "surface temperature": "ambient"},
+        name="ambient surface temperature",
+    ),
+    pybamm.lithium_ion.SPMe(
+        {"thermal": "lumped", "surface temperature": "lumped"},
+        name="lumped surface temperature",
+    ),
+]
+
+experiment = pybamm.Experiment(
+    [
+        "Discharge at 1C until 2.5V",
+        "Rest for 1 hour",
+    ]
+)
+
+parameter_values = pybamm.ParameterValues("Chen2020")
+parameter_values.update(
+    {
+        "Casing heat capacity [J.K-1]": 30,
+        "Environment thermal resistance [K.W-1]": 10,
+    },
+    check_already_exists=False,
+)
+
+# create and run simulations
+sols = []
+for model in models:
+    model.variables["Bulk temperature [°C]"] = (
+        model.variables["Volume-averaged cell temperature [K]"] - 273.15
+    )
+    model.variables["Surface temperature [°C]"] = (
+        model.variables["Surface temperature [K]"] - 273.15
+    )
+    sim = pybamm.Simulation(
+        model, parameter_values=parameter_values, experiment=experiment
+    )
+    sol = sim.solve([0, 3600])
+    sols.append(sol)
+
+# plot
+pybamm.dynamic_plot(
+    sols,
+    [
+        "Voltage [V]",
+        "Bulk temperature [°C]",
+        "Surface temperature [°C]",
+        "Surface total cooling [W]",
+        "Environment total cooling [W]",
+    ],
+)

src/pybamm/CITATIONS.bib

@@ -252,6 +252,16 @@ @article{Lain2019
   doi = {10.3390/batteries5040064},
 }
 
+@article{lin2014lumped,
+  title={A lumped-parameter electro-thermal model for cylindrical batteries},
+  author={Lin, Xinfan and Perez, Hector E and Mohan, Shankar and Siegel, Jason B and Stefanopoulou, Anna G and Ding, Yi and Castanier, Matthew P},
+  journal={Journal of Power Sources},
+  volume={257},
+  pages={1--11},
+  year={2014},
+  publisher={Elsevier}
+}
+
 @article{Marquis2019,
   title = {{An asymptotic derivation of a single particle model with electrolyte}},
   author = {Marquis, Scott G. and Sulzer, Valentin and Timms, Robert and Please, Colin P. and Chapman, S. Jon},

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -194,6 +194,11 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "surface form" : str
                 Whether to use the surface formulation of the problem. Can be "false"
                 (default), "differential" or "algebraic".
+            * "surface temperature" : str
+                Sets the surface temperature model to use. Can be "ambient" (default),
+                which sets the surface temperature equal to the ambient temperature, or
+                "lumped", which adds an ODE for the surface temperature (e.g. to model
+                internal heating of a thermal chamber).
             * "thermal" : str
                 Sets the thermal model to use. Can be "isothermal" (default), "lumped",
                 "x-lumped", or "x-full". The 'cell geometry' option must be set to
@@ -278,7 +283,6 @@ def __init__(self, extra_options):
             ],
             "particle": [
                 "Fickian diffusion",
-                "fast diffusion",
                 "uniform profile",
                 "quadratic profile",
                 "quartic profile",
@@ -304,6 +308,7 @@ def __init__(self, extra_options):
             "SEI porosity change": ["false", "true"],
             "stress-induced diffusion": ["false", "true"],
             "surface form": ["false", "differential", "algebraic"],
+            "surface temperature": ["ambient", "lumped"],
             "thermal": ["isothermal", "lumped", "x-lumped", "x-full"],
             "total interfacial current density as a state": ["false", "true"],
             "transport efficiency": [
@@ -554,16 +559,6 @@ def __init__(self, extra_options):
                 )
 
         # Renamed options
-        if options["particle"] == "fast diffusion":
-            raise pybamm.OptionError(
-                "The 'fast diffusion' option has been renamed. "
-                "Use 'uniform profile' instead."
-            )
-        if options["SEI porosity change"] in [True, False]:
-            raise pybamm.OptionError(
-                "SEI porosity change must now be given in string format "
-                "('true' or 'false')"
-            )
         if options["working electrode"] == "negative":
             raise pybamm.OptionError(
                 "The 'negative' working electrode option has been removed because "
@@ -633,6 +628,12 @@ def __init__(self, extra_options):
                     "be 'none': 'particle mechanics', 'loss of active material'"
                 )
 
+        if options["surface temperature"] == "lumped":
+            if options["thermal"] not in ["isothermal", "lumped"]:
+                raise pybamm.OptionError(
+                    "lumped surface temperature model only compatible with isothermal "
+                    "or lumped thermal model"
+                )
         if "true" in options["SEI on cracks"]:
             sei_on_cr = options["SEI on cracks"]
             p_mechanics = options["particle mechanics"]
@@ -1265,6 +1266,13 @@ def set_thermal_submodel(self):
             self.param, self.options, x_average
         )
 
+    def set_surface_temperature_submodel(self):
+        if self.options["surface temperature"] == "ambient":
+            submodel = pybamm.thermal.surface.Ambient
+        elif self.options["surface temperature"] == "lumped":
+            submodel = pybamm.thermal.surface.Lumped
+        self.submodels["surface temperature"] = submodel(self.param, self.options)
+
     def set_current_collector_submodel(self):
         if self.options["current collector"] in ["uniform"]:
             submodel = pybamm.current_collector.Uniform(self.param)

src/pybamm/models/full_battery_models/lead_acid/full.py

@@ -27,6 +27,7 @@ def __init__(self, options=None, name="Full model", build=True):
         self.set_electrolyte_submodel()
         self.set_solid_submodel()
         self.set_thermal_submodel()
+        self.set_surface_temperature_submodel()
         self.set_side_reaction_submodels()
         self.set_current_collector_submodel()
         self.set_sei_submodel()

src/pybamm/models/full_battery_models/lead_acid/loqs.py

@@ -27,6 +27,7 @@ def __init__(self, options=None, name="LOQS model", build=True):
         self.set_electrolyte_submodel()
         self.set_electrode_submodels()
         self.set_thermal_submodel()
+        self.set_surface_temperature_submodel()
         self.set_side_reaction_submodels()
         self.set_current_collector_submodel()
         self.set_sei_submodel()

src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -48,6 +48,7 @@ def set_submodels(self, build):
         self.set_electrolyte_concentration_submodel()
         self.set_electrolyte_potential_submodel()
         self.set_thermal_submodel()
+        self.set_surface_temperature_submodel()
         self.set_current_collector_submodel()
         self.set_sei_submodel()
         self.set_sei_on_cracks_submodel()

src/pybamm/models/submodels/thermal/__init__.py

@@ -2,5 +2,6 @@
 from .isothermal import Isothermal
 from .lumped import Lumped
 from . import pouch_cell
+from . import surface
 
-__all__ = ['base_thermal', 'isothermal', 'lumped', 'pouch_cell']
+__all__ = ["base_thermal", "isothermal", "lumped", "pouch_cell"]

src/pybamm/models/submodels/thermal/isothermal.py

@@ -73,8 +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]",
+                "Surface total cooling [W.m-3]",
+                "Surface total cooling [W]",
             ]:
                 # All variables are zero
                 variables.update({var: zero})

src/pybamm/models/submodels/thermal/lumped.py

@@ -49,9 +49,9 @@ def get_coupled_variables(self, 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]"]
+        T_surf = variables["Surface 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_W = -self.param.h_total * (T_vol_av - T_surf) * self.param.A_cooling
         Q_cool_vol_av = Q_cool_W / V
 
         # Contact resistance heating Q_cr
@@ -67,8 +67,8 @@ def get_coupled_variables(self, variables):
         variables.update(
             {
                 # Lumped cooling
-                "Lumped total cooling [W.m-3]": Q_cool_vol_av,
-                "Lumped total cooling [W]": Q_cool_W,
+                "Surface total cooling [W.m-3]": Q_cool_vol_av,
+                "Surface 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,
@@ -79,7 +79,7 @@ def get_coupled_variables(self, 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]"]
-        Q_cool_vol_av = variables["Lumped total cooling [W.m-3]"]
+        Q_cool_vol_av = variables["Surface 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]"

src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py

@@ -54,7 +54,7 @@ def get_coupled_variables(self, variables):
     def set_rhs(self, variables):
         T_av = variables["X-averaged cell temperature [K]"]
         Q_av = variables["X-averaged total heating [W.m-3]"]
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
@@ -64,13 +64,13 @@ def set_rhs(self, variables):
         cell_volume = self.param.L * self.param.L_y * self.param.L_z
         Q_yz_surface = (
             -(self.param.n.h_cc(y, z) + self.param.p.h_cc(y, z))
-            * (T_av - T_amb)
+            * (T_av - T_surf)
             * yz_surface_area
             / cell_volume
         )
         Q_edge = (
             -(self.param.h_edge(0, z) + self.param.h_edge(self.param.L_y, z))
-            * (T_av - T_amb)
+            * (T_av - T_surf)
             * edge_area
             / cell_volume
         )
@@ -87,7 +87,7 @@ def set_rhs(self, variables):
 
     def set_boundary_conditions(self, variables):
         param = self.param
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         T_av = variables["X-averaged cell temperature [K]"]
 
         # Find tab locations (top vs bottom)
@@ -118,10 +118,10 @@ def set_boundary_conditions(self, variables):
         # Calculate heat fluxes weighted by area
         # Note: can't do y-average of h_edge here since y isn't meshed. Evaluate at
         # midpoint.
-        q_tab_n = -param.n.h_tab * (T_av - T_amb)
-        q_tab_p = -param.p.h_tab * (T_av - T_amb)
-        q_edge_top = -param.h_edge(L_y / 2, L_z) * (T_av - T_amb)
-        q_edge_bottom = -param.h_edge(L_y / 2, 0) * (T_av - T_amb)
+        q_tab_n = -param.n.h_tab * (T_av - T_surf)
+        q_tab_p = -param.p.h_tab * (T_av - T_surf)
+        q_edge_top = -param.h_edge(L_y / 2, L_z) * (T_av - T_surf)
+        q_edge_bottom = -param.h_edge(L_y / 2, 0) * (T_av - T_surf)
         q_top = (
             q_tab_n * neg_tab_area * neg_tab_top_bool
             + q_tab_p * pos_tab_area * pos_tab_top_bool

src/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py

@@ -54,15 +54,15 @@ def get_coupled_variables(self, variables):
     def set_rhs(self, variables):
         T_av = variables["X-averaged cell temperature [K]"]
         Q_av = variables["X-averaged total heating [W.m-3]"]
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
         # Calculate cooling
         Q_yz_surface_W_per_m2 = -(self.param.n.h_cc(y, z) + self.param.p.h_cc(y, z)) * (
-            T_av - T_amb
+            T_av - T_surf
         )
-        Q_edge_W_per_m2 = -self.param.h_edge(y, z) * (T_av - T_amb)
+        Q_edge_W_per_m2 = -self.param.h_edge(y, z) * (T_av - T_surf)
 
         # Account for surface area to volume ratio of pouch cell in surface cooling
         # term
@@ -98,14 +98,14 @@ def set_rhs(self, variables):
 
     def set_boundary_conditions(self, variables):
         T_av = variables["X-averaged cell temperature [K]"]
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
         # Calculate heat fluxes
-        q_tab_n = -self.param.n.h_tab * (T_av - T_amb)
-        q_tab_p = -self.param.p.h_tab * (T_av - T_amb)
-        q_edge = -self.param.h_edge(y, z) * (T_av - T_amb)
+        q_tab_n = -self.param.n.h_tab * (T_av - T_surf)
+        q_tab_p = -self.param.p.h_tab * (T_av - T_surf)
+        q_edge = -self.param.h_edge(y, z) * (T_av - T_surf)
 
         # Subtract the edge cooling from the tab portion so as to not double count
         # Note: tab cooling is also only applied on the current collector hence

src/pybamm/models/submodels/thermal/pouch_cell/x_full.py

@@ -79,7 +79,7 @@ def set_rhs(self, variables):
         Q = variables["Total heating [W.m-3]"]
         Q_cn = variables["Negative current collector Ohmic heating [W.m-3]"]
         Q_cp = variables["Positive current collector Ohmic heating [W.m-3]"]
-        T_amb = variables["Ambient temperature [K]"]
+        T_surf = variables["Surface temperature [K]"]
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
@@ -140,28 +140,28 @@ def set_rhs(self, variables):
                 (
                     pybamm.boundary_value(lambda_n, "left")
                     * pybamm.boundary_gradient(T_n, "left")
-                    - h_cn * (T_cn - T_amb)
+                    - h_cn * (T_cn - T_surf)
                 )
                 / L_cn
                 + Q_cn
-                + cooling_coefficient_cn * (T_cn - T_amb)
+                + cooling_coefficient_cn * (T_cn - T_surf)
             )
             / self.param.n.rho_c_p_cc(T_cn),
             T: (
                 pybamm.div(lambda_ * pybamm.grad(T))
                 + Q
-                + cooling_coefficient * (T - T_amb)
+                + cooling_coefficient * (T - T_surf)
             )
             / rho_c_p,
             T_cp: (
                 (
                     -pybamm.boundary_value(lambda_p, "right")
                     * pybamm.boundary_gradient(T_p, "right")
-                    - h_cp * (T_cp - T_amb)
+                    - h_cp * (T_cp - T_surf)
                 )
                 / L_cp
                 + Q_cp
-                + cooling_coefficient_cp * (T_cp - T_amb)
+                + cooling_coefficient_cp * (T_cp - T_surf)
             )
             / self.param.p.rho_c_p_cc(T_cp),
         }

src/pybamm/models/submodels/thermal/surface/__init__.py

@@ -0,0 +1,2 @@
+from .ambient import Ambient
+from .lumped import Lumped

src/pybamm/models/submodels/thermal/surface/ambient.py

@@ -0,0 +1,33 @@
+#
+# Class for ambient surface temperature submodel
+import pybamm
+
+
+class Ambient(pybamm.BaseSubModel):
+    """
+    Class for setting surface temperature equal to ambient temperature.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    options : dict, optional
+        A dictionary of options to be passed to the model.
+
+    """
+
+    def __init__(self, param, options=None):
+        super().__init__(param, options=options)
+
+    def get_coupled_variables(self, variables):
+        T_amb = variables["Ambient temperature [K]"]
+        T_amb_av = variables["Volume-averaged ambient temperature [K]"]
+
+        variables.update(
+            {
+                "Surface temperature [K]": T_amb,
+                "Volume-averaged surface temperature [K]": T_amb_av,
+                "Environment total cooling [W]": pybamm.Scalar(0),
+            }
+        )
+        return variables