Issue 678 eleclyte temp depend

CHANGELOG.md

@@ -20,6 +20,7 @@
 
 ## Bug fixes
 
+-   Adds missing temperature dependence in electrolyte and interface submodels ([#698](https://github.com/pybamm-team/PyBaMM/pull/698))
 -   Fix differentiation of functions that have more than one argument ([#687](https://github.com/pybamm-team/PyBaMM/pull/687))
 -   Add warning if `ProcessedVariable` is called outside its interpolation range ([#681](https://github.com/pybamm-team/PyBaMM/pull/681))
 -   Improve the way `ProcessedVariable` objects are created in higher dimensions ([#581](https://github.com/pybamm-team/PyBaMM/pull/581))

examples/notebooks/models/SPM.ipynb

@@ -638,17 +638,17 @@
       "\t- Local current collector potential difference\n",
       "\t- Local current collector potential difference [V]\n",
       "\t- Ohmic heating\n",
-      "\t- Ohmic heating [A.V.m-3]\n",
+      "\t- Ohmic heating [W.m-3]\n",
       "\t- Irreversible electrochemical heating\n",
-      "\t- Irreversible electrochemical heating [A.V.m-3]\n",
+      "\t- Irreversible electrochemical heating [W.m-3]\n",
       "\t- Reversible heating\n",
-      "\t- Reversible heating [A.V.m-3]\n",
+      "\t- Reversible heating [W.m-3]\n",
       "\t- Total heating\n",
-      "\t- Total heating [A.V.m-3]\n",
+      "\t- Total heating [W.m-3]\n",
       "\t- X-averaged total heating\n",
-      "\t- X-averaged total heating [A.V.m-3]\n",
+      "\t- X-averaged total heating [W.m-3]\n",
       "\t- Volume-averaged total heating\n",
-      "\t- Volume-averaged total heating [A.V.m-3]\n",
+      "\t- Volume-averaged total heating [W.m-3]\n",
       "\t- Positive current collector potential [V]\n",
       "\t- Local potential difference\n",
       "\t- Local potential difference [V]\n",

examples/scripts/thermal_lithium_ion.py

@@ -18,7 +18,7 @@
 
 # load parameter values and process models and geometry
 param = models[0].default_parameter_values
-param.update({"Heat transfer coefficient [W.m-2.K-1]": 0.1})
+param.update({"Heat transfer coefficient [W.m-2.K-1]": 1})
 
 for model in models:
     param.process_model(model)
@@ -40,7 +40,8 @@
 solutions = [None] * len(models)
 t_eval = np.linspace(0, 0.17, 100)
 for i, model in enumerate(models):
-    solution = model.default_solver.solve(model, t_eval)
+    solver = pybamm.ScipySolver(atol=1e-8, rtol=1e-8)
+    solution = solver.solve(model, t_eval)
     solutions[i] = solution
 
 # plot

pybamm/expression_tree/binary_operators.py

@@ -277,9 +277,15 @@ def _binary_simplify(self, left, right):
             return left
         # Check matrices after checking scalars
         if is_matrix_zero(left):
-            return right
+            if isinstance(right, pybamm.Scalar):
+                return pybamm.Array(right.value * np.ones(left.shape_for_testing))
+            else:
+                return right
         if is_matrix_zero(right):
-            return left
+            if isinstance(left, pybamm.Scalar):
+                return pybamm.Array(left.value * np.ones(right.shape_for_testing))
+            else:
+                return left
 
         return pybamm.simplify_addition_subtraction(self.__class__, left, right)
 
@@ -325,9 +331,15 @@ def _binary_simplify(self, left, right):
             return left
         # Check matrices after checking scalars
         if is_matrix_zero(left):
-            return -right
+            if isinstance(right, pybamm.Scalar):
+                return pybamm.Array(-right.value * np.ones(left.shape_for_testing))
+            else:
+                return -right
         if is_matrix_zero(right):
-            return left
+            if isinstance(left, pybamm.Scalar):
+                return pybamm.Array(left.value * np.ones(right.shape_for_testing))
+            else:
+                return left
 
         return pybamm.simplify_addition_subtraction(self.__class__, left, right)
 

pybamm/models/submodels/electrolyte/base_electrolyte_conductivity.py

@@ -250,15 +250,15 @@ def _get_domain_current_variables(self, i_e, domain=None):
 
         variables = {
             domain + " electrolyte current density": i_e,
-            domain + " electrolyte current density [V]": i_e * i_typ,
+            domain + " electrolyte current density [A.m-2]": i_e * i_typ,
         }
 
         return variables
 
     def _get_whole_cell_variables(self, variables):
         """
         A private function to obtain the potential and current concatenated
-        across the whole cell. Note required 'variables' to contain the potential
+        across the whole cell. Note: requires 'variables' to contain the potential
         and current in the subdomains: 'negative electrode', 'separator', and
         'positive electrode'.
 

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py

@@ -48,11 +48,10 @@ def get_coupled_variables(self, variables):
         eps_s_av = variables["Leading-order x-averaged separator porosity"]
         eps_p_av = variables["Leading-order x-averaged positive electrode porosity"]
 
-        # Note: here we want the average of the temperature over the negative
-        # electrode, separator and positive electrode (not including the current
-        # collectors)
-        T = variables["Cell temperature"]
-        T_av = pybamm.x_average(T)
+        T_av = variables["X-averaged cell temperature"]
+        T_av_n = pybamm.PrimaryBroadcast(T_av, "negative electrode")
+        T_av_s = pybamm.PrimaryBroadcast(T_av, "separator")
+        T_av_p = pybamm.PrimaryBroadcast(T_av, "positive electrode")
 
         c_e_n, c_e_s, c_e_p = c_e.orphans
 
@@ -90,6 +89,7 @@ def get_coupled_variables(self, variables):
             + phi_s_n_av
             - (
                 chi_av
+                * (1 + param.Theta * T_av)
                 * pybamm.x_average(
                     self._higher_order_macinnes_function(
                         c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode")
@@ -106,6 +106,7 @@ def get_coupled_variables(self, variables):
             pybamm.PrimaryBroadcast(phi_e_const, "negative electrode")
             + (
                 chi_av_n
+                * (1 + param.Theta * T_av_n)
                 * self._higher_order_macinnes_function(
                     c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode")
                 )
@@ -124,6 +125,7 @@ def get_coupled_variables(self, variables):
             pybamm.PrimaryBroadcast(phi_e_const, "separator")
             + (
                 chi_av_s
+                * (1 + param.Theta * T_av_s)
                 * self._higher_order_macinnes_function(
                     c_e_s / pybamm.PrimaryBroadcast(c_e_av, "separator")
                 )
@@ -137,6 +139,7 @@ def get_coupled_variables(self, variables):
             pybamm.PrimaryBroadcast(phi_e_const, "positive electrode")
             + (
                 chi_av_p
+                * (1 + param.Theta * T_av_p)
                 * self._higher_order_macinnes_function(
                     c_e_p / pybamm.PrimaryBroadcast(c_e_av, "positive electrode")
                 )
@@ -155,15 +158,19 @@ def get_coupled_variables(self, variables):
         phi_e_av = pybamm.x_average(phi_e)
 
         # concentration overpotential
-        eta_c_av = chi_av * (
-            pybamm.x_average(
-                self._higher_order_macinnes_function(
-                    c_e_p / pybamm.PrimaryBroadcast(c_e_av, "positive electrode")
+        eta_c_av = (
+            chi_av
+            * (1 + param.Theta * T_av)
+            * (
+                pybamm.x_average(
+                    self._higher_order_macinnes_function(
+                        c_e_p / pybamm.PrimaryBroadcast(c_e_av, "positive electrode")
+                    )
                 )
-            )
-            - pybamm.x_average(
-                self._higher_order_macinnes_function(
-                    c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode")
+                - pybamm.x_average(
+                    self._higher_order_macinnes_function(
+                        c_e_n / pybamm.PrimaryBroadcast(c_e_av, "negative electrode")
+                    )
                 )
             )
         )

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py

@@ -39,7 +39,8 @@ def get_coupled_variables(self, variables):
         phi_e = variables["Electrolyte potential"]
 
         i_e = (param.kappa_e(c_e, T) * (eps ** param.b) * param.gamma_e / param.C_e) * (
-            param.chi(c_e) * pybamm.grad(c_e) / c_e - pybamm.grad(phi_e)
+            param.chi(c_e) * (1 + param.Theta * T) * pybamm.grad(c_e) / c_e
+            - pybamm.grad(phi_e)
         )
 
         variables.update(self._get_standard_current_variables(i_e))

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py

@@ -75,11 +75,15 @@ def set_boundary_conditions(self, variables):
         delta_phi = variables[self.domain + " electrode surface potential difference"]
 
         if self.domain == "Negative":
+            T = variables["Negative electrode temperature"]
             c_e_flux = pybamm.BoundaryGradient(c_e, "right")
             flux_left = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left")
             flux_right = (
                 (i_boundary_cc / pybamm.BoundaryValue(conductivity, "right"))
-                - pybamm.BoundaryValue(param.chi(c_e) / c_e, "right") * c_e_flux
+                - pybamm.BoundaryValue(
+                    (1 + param.Theta * T) * param.chi(c_e) / c_e, "right"
+                )
+                * c_e_flux
                 - i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "right")
             )
 
@@ -89,10 +93,14 @@ def set_boundary_conditions(self, variables):
             rbc_c_e = (c_e_flux, "Neumann")
 
         elif self.domain == "Positive":
+            T = variables["Positive electrode temperature"]
             c_e_flux = pybamm.BoundaryGradient(c_e, "left")
             flux_left = (
                 (i_boundary_cc / pybamm.BoundaryValue(conductivity, "left"))
-                - pybamm.BoundaryValue(param.chi(c_e) / c_e, "left") * c_e_flux
+                - pybamm.BoundaryValue(
+                    (1 + param.Theta * T) * param.chi(c_e) / c_e, "left"
+                )
+                * c_e_flux
                 - i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "left")
             )
             flux_right = -i_boundary_cc * pybamm.BoundaryValue(1 / sigma_eff, "right")
@@ -152,9 +160,10 @@ def _get_neg_pos_coupled_variables(self, variables):
         i_boundary_cc = variables["Current collector current density"]
         c_e = variables[self.domain + " electrolyte concentration"]
         delta_phi = variables[self.domain + " electrode surface potential difference"]
+        T = variables[self.domain + " electrode temperature"]
 
         i_e = conductivity * (
-            (param.chi(c_e) / c_e) * pybamm.grad(c_e)
+            ((1 + param.Theta * T) * param.chi(c_e) / c_e) * pybamm.grad(c_e)
             + pybamm.grad(delta_phi)
             + pybamm.PrimaryBroadcast(i_boundary_cc, self.domain_for_broadcast)
             / sigma_eff
@@ -190,7 +199,7 @@ def _get_sep_coupled_variables(self, variables):
         phi_e_s = pybamm.PrimaryBroadcast(
             pybamm.boundary_value(phi_e_n, "right"), "separator"
         ) + pybamm.IndefiniteIntegral(
-            chi_e_s / c_e_s * pybamm.grad(c_e_s)
+            (1 + param.Theta * T) * chi_e_s / c_e_s * pybamm.grad(c_e_s)
             - param.C_e
             * pybamm.PrimaryBroadcast(i_boundary_cc, self.domain_for_broadcast)
             / kappa_s_eff,

pybamm/models/submodels/interface/inverse_kinetics/base_inverse_kinetics.py

@@ -41,8 +41,13 @@ def get_coupled_variables(self, variables):
             ne = self.param.ne_n
         elif self.domain == "Positive":
             ne = self.param.ne_p
+        # Note: T must have the same domain as j0 and eta_r
+        if j0.domain in ["current collector", ["current collector"]]:
+            T = variables["X-averaged cell temperature"]
+        else:
+            T = variables[self.domain + " electrode temperature"]
 
-        eta_r = self._get_overpotential(j, j0, ne)
+        eta_r = self._get_overpotential(j, j0, ne, T)
         delta_phi = eta_r + ocp
 
         variables.update(self._get_standard_interfacial_current_variables(j))
@@ -73,5 +78,5 @@ def get_coupled_variables(self, variables):
 
         return variables
 
-    def _get_overpotential(self, j, j0, ne):
+    def _get_overpotential(self, j, j0, ne, T):
         raise NotImplementedError

pybamm/models/submodels/interface/inverse_kinetics/inverse_butler_volmer.py

@@ -28,5 +28,7 @@ class InverseButlerVolmer(BaseInverseKinetics, ButlerVolmer):
     def __init__(self, param, domain):
         super().__init__(param, domain)
 
-    def _get_overpotential(self, j, j0, ne):
-        return (2 / ne) * pybamm.Function(np.arcsinh, j / (2 * j0))
+    def _get_overpotential(self, j, j0, ne, T):
+        return (2 * (1 + self.param.Theta * T) / ne) * pybamm.Function(
+            np.arcsinh, j / (2 * j0)
+        )

pybamm/models/submodels/interface/kinetics/base_kinetics.py

@@ -40,8 +40,13 @@ def get_coupled_variables(self, variables):
         eta_r = delta_phi - ocp
         # Get number of electrons in reaction
         ne = self._get_number_of_electrons_in_reaction()
-
-        j = self._get_kinetics(j0, ne, eta_r)
+        # Get kinetics. Note: T must have the same domain as j0 and eta_r
+        if j0.domain in ["current collector", ["current collector"]]:
+            T = variables["X-averaged cell temperature"]
+        else:
+            T = variables[self.domain + " electrode temperature"]
+        j = self._get_kinetics(j0, ne, eta_r, T)
+        # Get average interfacial current density
         j_tot_av = self._get_average_total_interfacial_current_density(variables)
         # j = j_tot_av + (j - pybamm.x_average(j))  # enforce true average
 
@@ -73,7 +78,7 @@ def get_coupled_variables(self, variables):
     def _get_exchange_current_density(self, variables):
         raise NotImplementedError
 
-    def _get_kinetics(self, j0, ne, eta_r):
+    def _get_kinetics(self, j0, ne, eta_r, T):
         raise NotImplementedError
 
     def _get_open_circuit_potential(self, variables):
@@ -84,21 +89,21 @@ def _get_dj_dc(self, variables):
         Default to calculate derivative of interfacial current density with respect to
         concentration. Can be overwritten by specific kinetic functions.
         """
-        c_e, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order(
+        c_e, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order(
             variables
         )
-        j = self._get_kinetics(j0, ne, delta_phi - ocp)
+        j = self._get_kinetics(j0, ne, delta_phi - ocp, T)
         return j.diff(c_e)
 
     def _get_dj_ddeltaphi(self, variables):
         """
         Default to calculate derivative of interfacial current density with respect to
         surface potential difference. Can be overwritten by specific kinetic functions.
         """
-        _, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order(
+        _, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order(
             variables
         )
-        j = self._get_kinetics(j0, ne, delta_phi - ocp)
+        j = self._get_kinetics(j0, ne, delta_phi - ocp, T)
         return j.diff(delta_phi)
 
     def _get_interface_variables_for_first_order(self, variables):
@@ -118,7 +123,11 @@ def _get_interface_variables_for_first_order(self, variables):
         j0 = self._get_exchange_current_density(hacked_variables)
         ne = self._get_number_of_electrons_in_reaction()
         ocp = self._get_open_circuit_potential(hacked_variables)[0]
-        return c_e_0, delta_phi, j0, ne, ocp
+        if j0.domain in ["current collector", ["current collector"]]:
+            T = variables["X-averaged cell temperature"]
+        else:
+            T = variables[self.domain + " electrode temperature"]
+        return c_e_0, delta_phi, j0, ne, ocp, T
 
     def _get_j_diffusion_limited_first_order(self, variables):
         """

pybamm/models/submodels/interface/kinetics/butler_volmer.py

@@ -12,7 +12,7 @@ class ButlerVolmer(BaseModel):
     Base submodel which implements the forward Butler-Volmer equation:
 
     .. math::
-        j = j_0(c) * \\sinh(\\eta_r(c))
+        j = 2 * j_0(c) * \\sinh( (ne / (2 * (1 + \\Theta T)) * \\eta_r(c))
 
     Parameters
     ----------
@@ -28,26 +28,29 @@ class ButlerVolmer(BaseModel):
     def __init__(self, param, domain):
         super().__init__(param, domain)
 
-    def _get_kinetics(self, j0, ne, eta_r):
-        return 2 * j0 * pybamm.sinh((ne / 2) * eta_r)
+    def _get_kinetics(self, j0, ne, eta_r, T):
+        prefactor = ne / (2 * (1 + self.param.Theta * T))
+        return 2 * j0 * pybamm.sinh(prefactor * eta_r)
 
     def _get_dj_dc(self, variables):
         "See :meth:`pybamm.interface.kinetics.BaseModel._get_dj_dc`"
-        c_e, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order(
+        c_e, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order(
             variables
         )
         eta_r = delta_phi - ocp
-        return (2 * j0.diff(c_e) * pybamm.sinh((ne / 2) * eta_r)) - (
-            2 * j0 * (ne / 2) * ocp.diff(c_e) * pybamm.cosh((ne / 2) * eta_r)
+        prefactor = ne / (2 * (1 + self.param.Theta * T))
+        return (2 * j0.diff(c_e) * pybamm.sinh(prefactor * eta_r)) - (
+            2 * j0 * prefactor * ocp.diff(c_e) * pybamm.cosh(prefactor * eta_r)
         )
 
     def _get_dj_ddeltaphi(self, variables):
         "See :meth:`pybamm.interface.kinetics.BaseModel._get_dj_ddeltaphi`"
-        _, delta_phi, j0, ne, ocp = self._get_interface_variables_for_first_order(
+        _, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order(
             variables
         )
         eta_r = delta_phi - ocp
-        return 2 * j0 * (ne / 2) * pybamm.cosh((ne / 2) * eta_r)
+        prefactor = ne / (2 * (1 + self.param.Theta * T))
+        return 2 * j0 * prefactor * pybamm.cosh(prefactor * eta_r)
 
 
 class FirstOrderButlerVolmer(ButlerVolmer, BaseFirstOrderKinetics):

pybamm/models/submodels/interface/kinetics/no_reaction.py

@@ -24,5 +24,5 @@ class NoReaction(BaseModel):
     def __init__(self, param, domain):
         super().__init__(param, domain)
 
-    def _get_kinetics(self, j0, ne, eta_r):
+    def _get_kinetics(self, j0, ne, eta_r, T):
         return pybamm.Scalar(0)