#1927 half-cell thermal models

CHANGELOG.md

@@ -1,10 +1,14 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
-
 ## Features
 
+-   Isothermal models now compute heat source terms (but the temperature remains constant). The models now also account for current collector heating when `dimensionality=0` [#1929](https://github.com/pybamm-team/PyBaMM/pull/1929))
 -   Initial concentrations can now be provided as a function of `r` as well as `x` ([#1866](https://github.com/pybamm-team/PyBaMM/pull/1866))
 
+## Bug fixes
+
+-   Fixed a bug where thermal submodels could not be used with half-cells ([#1929](https://github.com/pybamm-team/PyBaMM/pull/1929))
+
 ## Breaking changes
 
 -   The `domain` setter and `auxiliary_domains` getter have been deprecated, `domains` setter/getter should be used instead. The `domain` getter is still active. We now recommend creating symbols with `domains={...}` instead of `domain=..., auxiliary_domains={...}`, but the latter is not yet deprecated ([#1866](https://github.com/pybamm-team/PyBaMM/pull/1866))

pybamm/models/full_battery_models/base_battery_model.py

@@ -362,6 +362,13 @@ def __init__(self, extra_options):
                     "cannot have transverse convection in 0D model"
                 )
 
+        if (
+            options["thermal"] in ["x-lumped", "x-full"]
+            and options["cell geometry"] != "pouch"
+        ):
+            raise pybamm.OptionError(
+                options["thermal"] + " model must have pouch geometry."
+            )
         if options["thermal"] == "x-full" and options["dimensionality"] != 0:
             n = options["dimensionality"]
             raise pybamm.OptionError(
@@ -378,6 +385,18 @@ def __init__(self, extra_options):
                     "mechanics model"
                 )
 
+        if options["working electrode"] != "both":
+            if options["thermal"] == "x-full":
+                raise pybamm.OptionError(
+                    "X-full thermal submodel is not compatible with half-cell models"
+                )
+            elif options["thermal"] == "x-lumped" and options["dimensionality"] != 0:
+                n = options["dimensionality"]
+                raise pybamm.OptionError(
+                    f"X-lumped thermal submodels do not yet support {n}D "
+                    "current collectors in a half-cell configuration"
+                )
+
         # Check options are valid
         for option, value in options.items():
             if option == "external submodels" or option == "working electrode":
@@ -578,13 +597,6 @@ def options(self, extra_options):
                 raise pybamm.OptionError(
                     "convection not implemented for lithium-ion models"
                 )
-            if (
-                options["thermal"] in ["x-lumped", "x-full"]
-                and options["cell geometry"] != "pouch"
-            ):
-                raise pybamm.OptionError(
-                    options["thermal"] + " model must have pouch geometry."
-                )
         if isinstance(self, pybamm.lithium_ion.SPMe):
             if options["electrolyte conductivity"] not in [
                 "default",
@@ -925,25 +937,27 @@ def set_thermal_submodel(self):
         elif self.options["thermal"] == "lumped":
             thermal_submodel = pybamm.thermal.Lumped(
                 self.param,
-                cc_dimension=self.options["dimensionality"],
-                geometry=self.options["cell geometry"],
+                self.options,
             )
 
         elif self.options["thermal"] == "x-lumped":
             if self.options["dimensionality"] == 0:
-                # With 0D current collectors x-lumped is equivalent to lumped pouch
-                thermal_submodel = pybamm.thermal.Lumped(self.param, geometry="pouch")
+                thermal_submodel = pybamm.thermal.Lumped(self.param, self.options)
             elif self.options["dimensionality"] == 1:
                 thermal_submodel = pybamm.thermal.pouch_cell.CurrentCollector1D(
-                    self.param
+                    self.param,
+                    self.options,
                 )
             elif self.options["dimensionality"] == 2:
                 thermal_submodel = pybamm.thermal.pouch_cell.CurrentCollector2D(
-                    self.param
+                    self.param,
+                    self.options,
                 )
         elif self.options["thermal"] == "x-full":
             if self.options["dimensionality"] == 0:
-                thermal_submodel = pybamm.thermal.OneDimensionalX(self.param)
+                thermal_submodel = pybamm.thermal.OneDimensionalX(
+                    self.param, self.options
+                )
 
         self.submodels["thermal"] = thermal_submodel
 

pybamm/models/submodels/current_collector/effective_resistance_current_collector.py

@@ -114,7 +114,8 @@ class EffectiveResistance(BaseEffectiveResistance):
     References
     ----------
     .. [1] R Timms, SG Marquis, V Sulzer, CP Please and SJ Chapman. “Asymptotic
-           Reduction of a Lithium-ion Pouch Cell Model”. Submitted, 2020.
+           Reduction of a Lithium-ion Pouch Cell Model”. SIAM Journal on Applied
+           Mathematics, 81(3), 765--788, 2021
 
     **Extends:** :class:`BaseEffectiveResistance`
     """
@@ -123,6 +124,9 @@ def __init__(
         self, options=None, name="Effective resistance in current collector model"
     ):
         super().__init__(name)
+
+        pybamm.citations.register("Timms2021")
+
         self.options = options
         self.param = pybamm.LithiumIonParameters()
 

pybamm/models/submodels/current_collector/potential_pair.py

@@ -17,10 +17,11 @@ class BasePotentialPair(BaseModel):
     References
     ----------
     .. [1] R Timms, SG Marquis, V Sulzer, CP Please and SJ Chapman. “Asymptotic
-           Reduction of a Lithium-ion Pouch Cell Model”. Submitted, 2020.
+           Reduction of a Lithium-ion Pouch Cell Model”. SIAM Journal on Applied
+           Mathematics, 81(3), 765--788, 2021
     .. [2] SG Marquis, R Timms, V Sulzer, CP Please and SJ Chapman. “A Suite of
-           Reduced-Order Models of a Single-Layer Lithium-ion Pouch Cell”. In
-           preparation, 2020.
+           Reduced-Order Models of a Single-Layer Lithium-ion Pouch Cell”. Journal
+           of The Electrochemical Society, 167(14):140513, 2020
 
     **Extends:** :class:`pybamm.current_collector.BaseModel`
     """

pybamm/models/submodels/thermal/base_thermal.py

@@ -5,22 +5,20 @@
 
 
 class BaseThermal(pybamm.BaseSubModel):
-    """Base class for thermal effects
+    """
+    Base class for thermal effects
 
     Parameters
     ----------
     param : parameter class
         The parameters to use for this submodel
     options : dict, optional
         A dictionary of options to be passed to the model.
-    cc_dimension: int, optional
-        The dimension of the current collectors. Can be 0 (default), 1 or 2.
 
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param, options=None, cc_dimension=0):
-        self.cc_dimension = cc_dimension
+    def __init__(self, param, options=None):
         super().__init__(param, options=options)
 
     def _get_standard_fundamental_variables(
@@ -29,7 +27,7 @@ def _get_standard_fundamental_variables(
         """
         Note: here we explicitly pass in the averages for the temperature as computing
         the average temperature in `BaseThermal` using `self._x_average` requires a
-        messy hack to avoid raising a `ModelError` (as the key in the equation
+        workaround to avoid raising a `ModelError` (as the key in the equation
         dict gets modified).
 
         For more information about this method in general,
@@ -38,7 +36,7 @@ def _get_standard_fundamental_variables(
         param = self.param
 
         # The variable T is the concatenation of the temperature in the negative
-        # electrode, separator and positive electrode, for use the electrochemical
+        # electrode, separator and positive electrode, for use in the electrochemical
         # models
         T = pybamm.concatenation(T_n, T_s, T_p)
 
@@ -90,60 +88,73 @@ def _get_standard_fundamental_variables(
         return variables
 
     def _get_standard_coupled_variables(self, variables):
-
         param = self.param
 
-        T_n = variables["Negative electrode temperature"]
-        T_p = variables["Positive electrode temperature"]
-
-        a_n = variables["Negative electrode surface area to volume ratio"]
-        a_p = variables["Positive electrode surface area to volume ratio"]
-
-        j_n = variables["Negative electrode interfacial current density"]
-        j_p = variables["Positive electrode interfacial current density"]
-
-        eta_r_n = variables["Negative electrode reaction overpotential"]
-        eta_r_p = variables["Positive electrode reaction overpotential"]
-
-        dUdT_n = variables["Negative electrode entropic change"]
-        dUdT_p = variables["Positive electrode entropic change"]
-
-        i_e = variables["Electrolyte current density"]
-        phi_e = variables["Electrolyte potential"]
-        phi_e_n = variables["Negative electrolyte potential"]
-        phi_e_s = variables["Separator electrolyte potential"]
-        phi_e_p = variables["Positive electrolyte potential"]
-
-        i_s_n = variables["Negative electrode current density"]
+        # Ohmic heating in solid
         i_s_p = variables["Positive electrode current density"]
-        phi_s_n = variables["Negative electrode potential"]
         phi_s_p = variables["Positive electrode potential"]
-
-        # Ohmic heating in solid
         Q_ohm_s_cn, Q_ohm_s_cp = self._current_collector_heating(variables)
-        Q_ohm_s_n = -pybamm.inner(i_s_n, pybamm.grad(phi_s_n))
+        if self.half_cell:
+            i_boundary_cc = variables["Current collector current density"]
+            T_n = variables["Negative electrode temperature"]
+            Q_ohm_s_n_av = i_boundary_cc ** 2 / param.sigma_n(T_n)
+            Q_ohm_s_n = pybamm.PrimaryBroadcast(Q_ohm_s_n_av, "negative electrode")
+        else:
+            i_s_n = variables["Negative electrode current density"]
+            phi_s_n = variables["Negative electrode potential"]
+            Q_ohm_s_n = -pybamm.inner(i_s_n, pybamm.grad(phi_s_n))
         Q_ohm_s_s = pybamm.FullBroadcast(0, ["separator"], "current collector")
         Q_ohm_s_p = -pybamm.inner(i_s_p, pybamm.grad(phi_s_p))
         Q_ohm_s = pybamm.concatenation(Q_ohm_s_n, Q_ohm_s_s, Q_ohm_s_p)
 
         # Ohmic heating in electrolyte
         # TODO: change full stefan-maxwell conductivity so that i_e is always
         # a Concatenation
+        i_e = variables["Electrolyte current density"]
+        phi_e = variables["Electrolyte potential"]
         if isinstance(i_e, pybamm.Concatenation):
             # compute by domain if possible
-            i_e_n, i_e_s, i_e_p = i_e.orphans
-            Q_ohm_e_n = -pybamm.inner(i_e_n, pybamm.grad(phi_e_n))
+            phi_e_s = variables["Separator electrolyte potential"]
+            phi_e_p = variables["Positive electrolyte potential"]
+            if self.half_cell:
+                i_e_s, i_e_p = i_e.orphans
+                Q_ohm_e_n = pybamm.FullBroadcast(
+                    0, ["negative electrode"], "current collector"
+                )
+            else:
+                i_e_n, i_e_s, i_e_p = i_e.orphans
+                phi_e_n = variables["Negative electrolyte potential"]
+                Q_ohm_e_n = -pybamm.inner(i_e_n, pybamm.grad(phi_e_n))
             Q_ohm_e_s = -pybamm.inner(i_e_s, pybamm.grad(phi_e_s))
             Q_ohm_e_p = -pybamm.inner(i_e_p, pybamm.grad(phi_e_p))
             Q_ohm_e = pybamm.concatenation(Q_ohm_e_n, Q_ohm_e_s, Q_ohm_e_p)
         else:
-            Q_ohm_e = -pybamm.inner(i_e, pybamm.grad(phi_e))
+            # else compute using i_e across all domains
+            if self.half_cell:
+                Q_ohm_e_n = pybamm.FullBroadcast(
+                    0, ["negative electrode"], "current collector"
+                )
+                Q_ohm_e_s_p = -pybamm.inner(i_e, pybamm.grad(phi_e))
+                Q_ohm_e = pybamm.concatenation(Q_ohm_e_n, Q_ohm_e_s_p)
+            else:
+                Q_ohm_e = -pybamm.inner(i_e, pybamm.grad(phi_e))
 
         # Total Ohmic heating
         Q_ohm = Q_ohm_s + Q_ohm_e
 
         # Irreversible electrochemical heating
-        Q_rxn_n = a_n * j_n * eta_r_n
+        a_p = variables["Positive electrode surface area to volume ratio"]
+        j_p = variables["Positive electrode interfacial current density"]
+        eta_r_p = variables["Positive electrode reaction overpotential"]
+        if self.half_cell:
+            Q_rxn_n = pybamm.FullBroadcast(
+                0, ["negative electrode"], "current collector"
+            )
+        else:
+            a_n = variables["Negative electrode surface area to volume ratio"]
+            j_n = variables["Negative electrode interfacial current density"]
+            eta_r_n = variables["Negative electrode reaction overpotential"]
+            Q_rxn_n = a_n * j_n * eta_r_n
         Q_rxn_p = a_p * j_p * eta_r_p
         Q_rxn = pybamm.concatenation(
             *[
@@ -154,7 +165,16 @@ def _get_standard_coupled_variables(self, variables):
         )
 
         # Reversible electrochemical heating
-        Q_rev_n = a_n * j_n * (param.Theta ** (-1) + T_n) * dUdT_n
+        T_p = variables["Positive electrode temperature"]
+        dUdT_p = variables["Positive electrode entropic change"]
+        if self.half_cell:
+            Q_rev_n = pybamm.FullBroadcast(
+                0, ["negative electrode"], "current collector"
+            )
+        else:
+            T_n = variables["Negative electrode temperature"]
+            dUdT_n = variables["Negative electrode entropic change"]
+            Q_rev_n = a_n * j_n * (param.Theta ** (-1) + T_n) * dUdT_n
         Q_rev_p = a_p * j_p * (param.Theta ** (-1) + T_p) * dUdT_p
         Q_rev = pybamm.concatenation(
             *[
@@ -216,28 +236,29 @@ def _get_standard_coupled_variables(self, variables):
 
     def _current_collector_heating(self, variables):
         """Compute Ohmic heating in current collectors."""
-        # TODO: implement grad in 0D to return a scalar zero
-        # TODO: implement grad_squared in other spatial methods so that the if
-        # statement can be removed
-        # In the limit of infinitely large current collector conductivity (i.e.
-        # 0D current collectors), the Ohmic heating in the current collectors is
-        # zero
+        cc_dimension = self.options["dimensionality"]
 
-        if self.cc_dimension == 0:
-            Q_s_cn = pybamm.Scalar(0)
-            Q_s_cp = pybamm.Scalar(0)
+        # Compute the Ohmic heating for 0D current collectors
+        if cc_dimension == 0:
+            i_boundary_cc = variables["Current collector current density"]
+            Q_s_cn = i_boundary_cc ** 2 / self.param.sigma_cn
+            Q_s_cp = i_boundary_cc ** 2 / self.param.sigma_cp
         # Otherwise we compute the Ohmic heating for 1 or 2D current collectors
-        elif self.cc_dimension in [1, 2]:
+        # In this limit the current flow is all in the y,z direction in the current
+        # collectors
+        elif cc_dimension in [1, 2]:
             phi_s_cn = variables["Negative current collector potential"]
             phi_s_cp = variables["Positive current collector potential"]
-            if self.cc_dimension == 1:
+            # TODO: implement grad_squared in other spatial methods so that the
+            # if statement can be removed
+            if cc_dimension == 1:
                 Q_s_cn = self.param.sigma_cn_prime * pybamm.inner(
                     pybamm.grad(phi_s_cn), pybamm.grad(phi_s_cn)
                 )
                 Q_s_cp = self.param.sigma_cp_prime * pybamm.inner(
                     pybamm.grad(phi_s_cp), pybamm.grad(phi_s_cp)
                 )
-            elif self.cc_dimension == 2:
+            elif cc_dimension == 2:
                 # Inner not implemented in 2D -- have to call grad_squared directly
                 Q_s_cn = self.param.sigma_cn_prime * pybamm.grad_squared(phi_s_cn)
                 Q_s_cp = self.param.sigma_cp_prime * pybamm.grad_squared(phi_s_cp)
@@ -266,7 +287,7 @@ def _yz_average(self, var):
         """Computes the y-z average."""
         # TODO: change the behaviour of z_average and yz_average so the if statement
         # can be removed
-        if self.cc_dimension in [0, 1]:
+        if self.options["dimensionality"] in [0, 1]:
             return pybamm.z_average(var)
-        elif self.cc_dimension == 2:
+        elif self.options["dimensionality"] == 2:
             return pybamm.yz_average(var)