Temperature dependent electrode conductivity

CHANGELOG.md

@@ -2,7 +2,8 @@
 
 ## Features
 
--   Added fitted expressions for OCPs for the Chen2020 parameter set ([#1526](https://github.com/pybamm-team/PyBaMM/pull/1497))
+-   Added temperature dependence on electrode electronic conductivity ([#1570](https://github.com/pybamm-team/PyBaMM/pull/1570))
+-   Added fitted expressions for OCPs for the Chen2020 parameter set ([#1526](https://github.com/pybamm-team/PyBaMM/pull/1526))
 -   Added `initial_soc` argument to `Simualtion.solve` for specifying the initial SOC when solving a model ([#1512](https://github.com/pybamm-team/PyBaMM/pull/1512))
 -   Added `print_name` to some symbols ([#1495](https://github.com/pybamm-team/PyBaMM/pull/1495), [#1497](https://github.com/pybamm-team/PyBaMM/pull/1497))
 -   Added Base Parameters class and SymPy in dependencies ([#1495](https://github.com/pybamm-team/PyBaMM/pull/1495))

pybamm/models/full_battery_models/lead_acid/basic_full.py

@@ -188,8 +188,8 @@ def __init__(self, name="Basic full model"):
         ######################
         # Current in the solid
         ######################
-        i_s_n = -param.sigma_n * (1 - eps_n) ** param.b_s_n * pybamm.grad(phi_s_n)
-        sigma_eff_p = param.sigma_p * (1 - eps_p) ** param.b_s_p
+        i_s_n = -param.sigma_n(T) * (1 - eps_n) ** param.b_s_n * pybamm.grad(phi_s_n)
+        sigma_eff_p = param.sigma_p(T) * (1 - eps_p) ** param.b_s_p
         i_s_p = -sigma_eff_p * pybamm.grad(phi_s_p)
         # The `algebraic` dictionary contains differential equations, with the key being
         # the main scalar variable of interest in the equation

pybamm/models/full_battery_models/lithium_ion/basic_dfn.py

@@ -214,9 +214,9 @@ def __init__(self, name="Doyle-Fuller-Newman model"):
         ######################
         # Current in the solid
         ######################
-        sigma_eff_n = param.sigma_n * eps_s_n ** param.b_s_n
+        sigma_eff_n = param.sigma_n(T) * eps_s_n ** param.b_s_n
         i_s_n = -sigma_eff_n * pybamm.grad(phi_s_n)
-        sigma_eff_p = param.sigma_p * eps_s_p ** param.b_s_p
+        sigma_eff_p = param.sigma_p(T) * eps_s_p ** param.b_s_p
         i_s_p = -sigma_eff_p * pybamm.grad(phi_s_p)
         # The `algebraic` dictionary contains differential equations, with the key being
         # the main scalar variable of interest in the equation

pybamm/models/full_battery_models/lithium_ion/basic_dfn_half_cell.py

@@ -252,7 +252,7 @@ def __init__(self, options=None, name="Doyle-Fuller-Newman half cell model"):
         ######################
         # Current in the solid
         ######################
-        sigma_eff_w = sigma_w * eps_s_w ** b_s_w
+        sigma_eff_w = sigma_w(T) * eps_s_w ** b_s_w
         i_s_w = -sigma_eff_w * pybamm.grad(phi_s_w)
         self.boundary_conditions[phi_s_w] = {
             "left": (pybamm.Scalar(0), "Neumann"),

pybamm/models/submodels/convection/through_cell/base_through_cell_convection.py

@@ -52,7 +52,7 @@ def _get_standard_neg_pos_velocity_variables(self, v_box_n, v_box_p):
         return variables
 
     def _get_standard_neg_pos_acceleration_variables(self, div_v_box_n, div_v_box_p):
-        """ Acceleration in the electrodes """
+        """Acceleration in the electrodes"""
 
         acc_scale = self.param.velocity_scale / self.param.L_x
 
@@ -79,7 +79,7 @@ def _get_standard_neg_pos_acceleration_variables(self, div_v_box_n, div_v_box_p)
         return variables
 
     def _get_standard_neg_pos_pressure_variables(self, p_n, p_p):
-        """ Pressure in the electrodes """
+        """Pressure in the electrodes"""
 
         variables = {
             "Negative electrode pressure": p_n,

pybamm/models/submodels/electrode/ohm/base_ohm.py

@@ -33,7 +33,10 @@ def set_boundary_conditions(self, variables):
         elif self.domain == "Positive":
             lbc = (pybamm.Scalar(0), "Neumann")
             i_boundary_cc = variables["Current collector current density"]
-            sigma_eff = self.param.sigma_p * variables["Positive electrode tortuosity"]
+            T_p = variables["Positive electrode temperature"]
+            sigma_eff = (
+                self.param.sigma_p(T_p) * variables["Positive electrode tortuosity"]
+            )
             rbc = (
                 i_boundary_cc / pybamm.boundary_value(-sigma_eff, "right"),
                 "Neumann",

pybamm/models/submodels/electrode/ohm/composite_ohm.py

@@ -38,9 +38,10 @@ def get_coupled_variables(self, variables):
             "Leading-order x-averaged " + self.domain.lower() + " electrode tortuosity"
         ]
         phi_s_cn = variables["Negative current collector potential"]
+        T = variables["X-averaged " + self.domain.lower() + " electrode temperature"]
 
         if self._domain == "Negative":
-            sigma_eff_0 = self.param.sigma_n * tor_0
+            sigma_eff_0 = self.param.sigma_n(T) * tor_0
             phi_s = phi_s_cn + (i_boundary_cc_0 / sigma_eff_0) * (
                 x_n * (x_n - 2 * l_n) / (2 * l_n)
             )
@@ -52,7 +53,7 @@ def get_coupled_variables(self, variables):
             ]
             phi_e_p_av = variables["X-averaged positive electrolyte potential"]
 
-            sigma_eff_0 = self.param.sigma_p * tor_0
+            sigma_eff_0 = self.param.sigma_p(T) * tor_0
 
             const = (
                 delta_phi_p_av
@@ -80,14 +81,15 @@ def set_boundary_conditions(self, variables):
             "Leading-order x-averaged " + self.domain.lower() + " electrode tortuosity"
         ]
         i_boundary_cc_0 = variables["Leading-order current collector current density"]
+        T = variables["X-averaged " + self.domain.lower() + " electrode temperature"]
 
         if self.domain == "Negative":
             lbc = (pybamm.Scalar(0), "Dirichlet")
             rbc = (pybamm.Scalar(0), "Neumann")
 
         elif self.domain == "Positive":
             lbc = (pybamm.Scalar(0), "Neumann")
-            sigma_eff_0 = self.param.sigma_p * tor_0
+            sigma_eff_0 = self.param.sigma_p(T) * tor_0
             rbc = (-i_boundary_cc_0 / sigma_eff_0, "Neumann")
 
         self.boundary_conditions[phi_s] = {"left": lbc, "right": rbc}

pybamm/models/submodels/electrode/ohm/full_ohm.py

@@ -37,11 +37,12 @@ def get_coupled_variables(self, variables):
 
         phi_s = variables[self.domain + " electrode potential"]
         tor = variables[self.domain + " electrode tortuosity"]
+        T = variables[self.domain + " electrode temperature"]
 
         if self.domain == "Negative":
-            sigma = self.param.sigma_n
+            sigma = self.param.sigma_n(T)
         elif self.domain == "Positive":
-            sigma = self.param.sigma_p
+            sigma = self.param.sigma_p(T)
 
         sigma_eff = sigma * tor
         i_s = -sigma_eff * pybamm.grad(phi_s)
@@ -76,14 +77,15 @@ def set_boundary_conditions(self, variables):
         phi_s = variables[self.domain + " electrode potential"]
         phi_s_cn = variables["Negative current collector potential"]
         tor = variables[self.domain + " electrode tortuosity"]
+        T = variables[self.domain + " electrode temperature"]
 
         if self.domain == "Negative":
             lbc = (phi_s_cn, "Dirichlet")
             rbc = (pybamm.Scalar(0), "Neumann")
 
         elif self.domain == "Positive":
             lbc = (pybamm.Scalar(0), "Neumann")
-            sigma_eff = self.param.sigma_p * tor
+            sigma_eff = self.param.sigma_p(T) * tor
             i_boundary_cc = variables["Current collector current density"]
             rbc = (
                 i_boundary_cc / pybamm.boundary_value(-sigma_eff, "right"),

pybamm/models/submodels/electrode/ohm/surface_form_ohm.py

@@ -33,19 +33,20 @@ def get_coupled_variables(self, variables):
         i_e = variables[self.domain + " electrolyte current density"]
         tor = variables[self.domain + " electrode tortuosity"]
         phi_s_cn = variables["Negative current collector potential"]
+        T = variables[self.domain + " electrode temperature"]
 
         i_s = i_boundary_cc - i_e
 
         if self.domain == "Negative":
-            conductivity = param.sigma_n * tor
+            conductivity = param.sigma_n(T) * tor
             phi_s = phi_s_cn - pybamm.IndefiniteIntegral(i_s / conductivity, x_n)
 
         elif self.domain == "Positive":
 
             phi_e_s = variables["Separator electrolyte potential"]
             delta_phi_p = variables["Positive electrode surface potential difference"]
 
-            conductivity = param.sigma_p * tor
+            conductivity = param.sigma_p(T) * tor
             phi_s = -pybamm.IndefiniteIntegral(i_s / conductivity, x_p) + (
                 pybamm.boundary_value(phi_e_s, "right")
                 + pybamm.boundary_value(delta_phi_p, "left")

pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/full_surface_form_conductivity.py

@@ -106,9 +106,9 @@ def _get_conductivities(self, variables):
         c_e = variables[self.domain + " electrolyte concentration"]
         T = variables[self.domain + " electrode temperature"]
         if self.domain == "Negative":
-            sigma = param.sigma_n
+            sigma = param.sigma_n(T)
         elif self.domain == "Positive":
-            sigma = param.sigma_p
+            sigma = param.sigma_p(T)
 
         kappa_eff = param.kappa_e(c_e, T) * tor_e
         sigma_eff = sigma * tor_s

pybamm/models/submodels/external_circuit/base_external_circuit.py

@@ -5,7 +5,7 @@
 
 
 class BaseModel(pybamm.BaseSubModel):
-    """Model to represent the behaviour of the external circuit. """
+    """Model to represent the behaviour of the external circuit."""
 
     def __init__(self, param):
         super().__init__(param)
@@ -27,7 +27,7 @@ def set_rhs(self, variables):
 
 
 class LeadingOrderBaseModel(BaseModel):
-    """Model to represent the behaviour of the external circuit, at leading order. """
+    """Model to represent the behaviour of the external circuit, at leading order."""
 
     def __init__(self, param):
         super().__init__(param)

pybamm/models/submodels/external_circuit/current_control_external_circuit.py

@@ -5,7 +5,7 @@
 
 
 class CurrentControl(BaseModel):
-    """External circuit with current control. """
+    """External circuit with current control."""
 
     def __init__(self, param):
         super().__init__(param)
@@ -30,7 +30,7 @@ def get_fundamental_variables(self):
 
 
 class LeadingOrderCurrentControl(CurrentControl, LeadingOrderBaseModel):
-    """External circuit with current control, for leading order models. """
+    """External circuit with current control, for leading order models."""
 
     def __init__(self, param):
         super().__init__(param)

pybamm/models/submodels/external_circuit/function_control_external_circuit.py

@@ -6,7 +6,7 @@
 
 
 class FunctionControl(BaseModel):
-    """External circuit with an arbitrary function. """
+    """External circuit with an arbitrary function."""
 
     def __init__(self, param, external_circuit_function):
         super().__init__(param)
@@ -63,7 +63,7 @@ def constant_voltage(self, variables):
 
 
 class PowerFunctionControl(FunctionControl):
-    """External circuit with power control. """
+    """External circuit with power control."""
 
     def __init__(self, param):
         super().__init__(param, self.constant_power)
@@ -77,7 +77,7 @@ def constant_power(self, variables):
 
 
 class LeadingOrderFunctionControl(FunctionControl, LeadingOrderBaseModel):
-    """External circuit with an arbitrary function, at leading order. """
+    """External circuit with an arbitrary function, at leading order."""
 
     def __init__(self, param, external_circuit_class):
         super().__init__(param, external_circuit_class)
@@ -103,7 +103,7 @@ def constant_voltage(self, variables):
 
 
 class LeadingOrderPowerFunctionControl(LeadingOrderFunctionControl):
-    """External circuit with power control, at leading order. """
+    """External circuit with power control, at leading order."""
 
     def __init__(self, param):
         super().__init__(param, self.constant_power)

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

@@ -36,7 +36,7 @@ def _get_kinetics(self, j0, ne, eta_r, T):
         return 2 * j0 * pybamm.sinh(prefactor * eta_r)
 
     def _get_dj_dc(self, variables):
-        """ See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_dc` """
+        """See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_dc`"""
         c_e, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order(
             variables
         )
@@ -47,7 +47,7 @@ def _get_dj_dc(self, variables):
         )
 
     def _get_dj_ddeltaphi(self, variables):
-        """ See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_ddeltaphi` """
+        """See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_ddeltaphi`"""
         _, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order(
             variables
         )

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

@@ -34,7 +34,7 @@ def _get_kinetics(self, j0, ne, eta_r, T):
         return j0 * pybamm.exp((ne / (2 * (1 + self.param.Theta * T))) * eta_r)
 
     def _get_dj_dc(self, variables):
-        """ See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_dc` """
+        """See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_dc`"""
         c_e, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order(
             variables
         )
@@ -46,7 +46,7 @@ def _get_dj_dc(self, variables):
         )
 
     def _get_dj_ddeltaphi(self, variables):
-        """ See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_ddeltaphi` """
+        """See :meth:`pybamm.interface.kinetics.BaseKinetics._get_dj_ddeltaphi`"""
         _, delta_phi, j0, ne, ocp, T = self._get_interface_variables_for_first_order(
             variables
         )

pybamm/models/submodels/particle_cracking/__init__.py

@@ -1,3 +1,3 @@
 from .base_cracking import BaseCracking
 from .crack_propagation import CrackPropagation
-from .swelling_only import SwellingOnly
+from .swelling_only import SwellingOnly

pybamm/models/submodels/thermal/base_thermal.py

@@ -217,6 +217,7 @@ def _current_collector_heating(self, variables):
         # In the limit of infinitely large current collector conductivity (i.e.
         # 0D current collectors), the Ohmic heating in the current collectors is
         # zero
+
         if self.cc_dimension == 0:
             Q_s_cn = pybamm.Scalar(0)
             Q_s_cp = pybamm.Scalar(0)

pybamm/parameters/lead_acid_parameters.py

@@ -122,15 +122,16 @@ def _set_dimensional_parameters(self):
         self.Q_p_max_dimensional = pybamm.Parameter(
             "Positive electrode volumetric capacity [C.m-3]"
         )
-        self.sigma_n_dim = pybamm.Parameter("Negative electrode conductivity [S.m-1]")
-        self.sigma_p_dim = pybamm.Parameter("Positive electrode conductivity [S.m-1]")
         # In lead-acid the current collector and electrodes are the same (same
-        # conductivity) but we correct here for Bruggeman
+        # conductivity) but we correct here for Bruggeman. Note that because for
+        # lithium-ion we allow electrode conductivity to be a function of temperature,
+        # but not the current collector conductivity, here the latter is evaluated at
+        # T_ref.
         self.sigma_cn_dimensional = (
-            self.sigma_n_dim * (1 - self.eps_n_max) ** self.b_s_n
+            self.sigma_n_dimensional(self.T_ref) * (1 - self.eps_n_max) ** self.b_s_n
         )
         self.sigma_cp_dimensional = (
-            self.sigma_p_dim * (1 - self.eps_p_max) ** self.b_s_p
+            self.sigma_p_dimensional(self.T_ref) * (1 - self.eps_p_max) ** self.b_s_p
         )
 
         # Electrochemical reactions
@@ -225,6 +226,20 @@ def _set_dimensional_parameters(self):
         self.R_sei_dimensional = pybamm.Scalar(0)
         self.beta_sei_n = pybamm.Scalar(0)
 
+    def sigma_n_dimensional(self, T):
+        """Dimensional electrical conductivity in negative electrode"""
+        inputs = {"Temperature [K]": T}
+        return pybamm.FunctionParameter(
+            "Negative electrode conductivity [S.m-1]", inputs
+        )
+
+    def sigma_p_dimensional(self, T):
+        """Dimensional electrical conductivity in positive electrode"""
+        inputs = {"Temperature [K]": T}
+        return pybamm.FunctionParameter(
+            "Positive electrode conductivity [S.m-1]", inputs
+        )
+
     def t_plus(self, c_e, T):
         """Dimensionless transference number (i.e. c_e is dimensionless)"""
         inputs = {"Electrolyte concentration [mol.m-3]": c_e * self.c_e_typ}
@@ -480,21 +495,14 @@ def _set_dimensionless_parameters(self):
             * (1 - self.M_w * self.V_hy / self.V_w * self.M_hy)
         )
 
+        # Electrode Properties
         # Electrode Properties
         self.sigma_cn = (
             self.sigma_cn_dimensional * self.potential_scale / self.i_typ / self.L_x
         )
-        self.sigma_n = (
-            self.sigma_n_dim * self.potential_scale / self.current_scale / self.L_x
-        )
-        self.sigma_p = (
-            self.sigma_p_dim * self.potential_scale / self.current_scale / self.L_x
-        )
         self.sigma_cp = (
             self.sigma_cp_dimensional * self.potential_scale / self.i_typ / self.L_x
         )
-        self.sigma_n_prime = self.sigma_n * self.delta ** 2
-        self.sigma_p_prime = self.sigma_p * self.delta ** 2
         self.sigma_cn_prime = self.sigma_cn * self.delta ** 2
         self.sigma_cp_prime = self.sigma_cp * self.delta ** 2
         self.delta_pore_n = 1 / (self.a_n_typ * self.L_x)
@@ -676,6 +684,34 @@ def _set_dimensionless_parameters(self):
             self.Q_e_max * (1.2 - self.q_init) / (self.Q_p_max * self.l_p)
         )
 
+    def sigma_n(self, T):
+        """Dimensionless negative electrode electrical conductivity"""
+        T_dim = self.Delta_T * T + self.T_ref
+        return (
+            self.sigma_n_dimensional(T_dim)
+            * self.potential_scale
+            / self.current_scale
+            / self.L_x
+        )
+
+    def sigma_p(self, T):
+        """Dimensionless positive electrode electrical conductivity"""
+        T_dim = self.Delta_T * T + self.T_ref
+        return (
+            self.sigma_p_dimensional(T_dim)
+            * self.potential_scale
+            / self.current_scale
+            / self.L_x
+        )
+
+    def sigma_n_prime(self, T):
+        """Rescaled dimensionless negative electrode electrical conductivity"""
+        return self.sigma_n(T) * self.delta ** 2
+
+    def sigma_p_prime(self, T):
+        """Rescaled dimensionless positive electrode electrical conductivity"""
+        return self.sigma_p(T) * self.delta ** 2
+
     def D_e(self, c_e, T):
         """Dimensionless electrolyte diffusivity"""
         c_e_dimensional = c_e * self.c_e_typ