#834 coupled currents to macroscale

examples/scripts/sei_growth.py

@@ -5,6 +5,7 @@
 # pb.settings.debug_mode = True
 
 options = {"sei": "reaction limited"}
+# options = {"sei": None}
 model = pb.lithium_ion.DFN(options)
 
 # experiment = pb.Experiment(
@@ -40,25 +41,25 @@
 #     }
 # )
 
-# parameter_values.update(
-#     {
-#         "Inner SEI reaction proportion": 0.5,
-#         "Inner SEI partial molar volume [m3.mol-1]": 34.76e-7,
-#         "Outer SEI partial molar volume [m3.mol-1]": 34.76e-7,
-#         "SEI reaction exchange current density [A.m-2]": 1.5e-7,
-#         "SEI resistance per unit thickness [Ohm.m-1]": 1,
-#         "Outer SEI solvent diffusivity [m2.s-1]": 2.5e-22,
-#         "Bulk solvent concentration [mol.m-3]": 2.636e3,
-#         "Ratio of inner and outer SEI exchange current densities": 1,
-#         "Inner SEI open-circuit potential [V]": 0.1,
-#         "Outer SEI open-circuit potential [V]": 0.8,
-#         "Inner SEI electron conducticity [S.m-1]": 8.95e-14,
-#         "Inner SEI lithium interstitial diffusivity [m2.s-1]": 1e-15,
-#         "Lithium interstitial reference concentration [mol.m-3]": 15,
-#         "Initial inner SEI thickness [m]": 7.5e-9,
-#         "Initial outer SEI thickness [m]": 7.5e-9,
-#     }
-# )
+parameter_values.update(
+    {
+        "Inner SEI reaction proportion": 0.5,
+        "Inner SEI partial molar volume [m3.mol-1]": 34.76e-6,
+        "Outer SEI partial molar volume [m3.mol-1]": 34.76e-6,
+        "SEI reaction exchange current density [A.m-2]": 1.5e-9,
+        "SEI resistance per unit thickness [Ohm.m-1]": 0,
+        "Outer SEI solvent diffusivity [m2.s-1]": 2.5e-22,
+        "Bulk solvent concentration [mol.m-3]": 2.636e3,
+        "Ratio of inner and outer SEI exchange current densities": 1,
+        "Inner SEI open-circuit potential [V]": 0.1,
+        "Outer SEI open-circuit potential [V]": 0.8,
+        "Inner SEI electron conducticity [S.m-1]": 8.95e-14,
+        "Inner SEI lithium interstitial diffusivity [m2.s-1]": 1e-15,
+        "Lithium interstitial reference concentration [mol.m-3]": 15,
+        "Initial inner SEI thickness [m]": 7.5e-9,
+        "Initial outer SEI thickness [m]": 7.5e-9,
+    }
+)
 
 # parameter_values.update({"Current function [A]": 0})
 parameter_values["Current function [A]"] = "[current data]US06"
@@ -82,12 +83,25 @@
 sim.plot(
     [
         "Terminal voltage [V]",
-        "Total SEI thickness [m]",
-        "X-averaged total SEI thickness [m]",
-        "X-averaged total SEI thickness",
-        "X-averaged SEI concentration [mol.m-3]",
-        "Loss of lithium to SEI [mols]",
-        "SEI reaction interfacial current density [A.m-2]",
-        "X-averaged SEI reaction interfacial current density [A.m-2]",
+        "Negative particle surface concentration",
+        "X-averaged negative particle surface concentration",
+        "Electrolyte concentration [mol.m-3]",
+        "Total negative electrode sei thickness [m]",
+        "X-averaged total negative electrode sei thickness [m]",
+        "X-averaged total negative electrode sei thickness",
+        "X-averaged negative electrode sei concentration [mol.m-3]",
+        "Loss of lithium to negative electrode sei [mols]",
+        [
+            "Negative electrode sei interfacial current density [A.m-2]",
+            "Negative electrode interfacial current density [A.m-2]",
+        ],
+        [
+            "X-averaged negative electrode sei interfacial current density [A.m-2]",
+            "X-averaged negative electrode interfacial current density [A.m-2]",
+        ],
+        [
+            "Negative electrode interfacial current density",
+            "Scaled negative electrode sei interfacial current density",
+        ],
     ]
 )

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -54,6 +54,34 @@ def set_reactions(self):
             }
         }
 
+        # N.B if there is no sei reaction then reaction
+        # is set to zero in the submodel.
+        self.reactions["sei"] = {
+            "Negative": {
+                "s": (1 - self.param.t_plus),
+                "aj": "Scaled negative electrode sei" + icd,
+            },
+            "Positive": {
+                "s": (1 - self.param.t_plus),
+                "aj": "Scaled positive electrode sei" + icd,
+            },
+        }
+
     def set_sei_submodel(self):
-        if self.options["sei"] == "reaction limited":
-            self.submodels["sei"] = pybamm.sei.ReactionLimited(self.param)
+
+        # negative electrode SEI
+        if self.options["sei"] is None:
+            self.submodels["negative sei"] = pybamm.sei.NoSEI(
+                self.param, "Negative electrode"
+            )
+
+        elif self.options["sei"] == "reaction limited":
+            self.submodels["negative sei"] = pybamm.sei.ReactionLimited(
+                self.param, "Negative electrode"
+            )
+
+        # positive electrode
+        self.submodels["positive sei"] = pybamm.sei.NoSEI(
+            self.param, "Positive electrode"
+        )
+

pybamm/models/submodels/sei/__init__.py

@@ -1,2 +1,3 @@
 from .base_sei import BaseModel
 from .reaction_limited import ReactionLimited
+from .no_sei import NoSEI

pybamm/models/submodels/sei/base_sei.py

@@ -17,7 +17,8 @@ class BaseModel(pybamm.BaseSubModel):
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param):
+    def __init__(self, param, domain):
+        self.domain = domain
         super().__init__(param)
 
     def _get_standard_thickness_variables(self, L_inner, L_outer):
@@ -60,25 +61,43 @@ def _get_standard_thickness_variables(self, L_inner, L_outer):
         n_outer_scale = sp.L_sei_0_dim * sp.a_n / sp.V_bar_outer_dimensional
 
         variables = {
-            "Inner SEI thickness": L_inner,
-            "Inner SEI thickness [m]": L_inner * L_scale,
-            "X-averaged inner SEI thickness": L_inner_av,
-            "X-averaged inner SEI thickness [m]": L_inner_av * L_scale,
-            "Outer SEI thickness": L_outer,
-            "Outer SEI thickness [m]": L_outer * L_scale,
-            "X-averaged outer SEI thickness": L_outer_av,
-            "X-averaged outer SEI thickness [m]": L_outer_av * L_scale,
-            "Total SEI thickness": L_sei,
-            "Total SEI thickness [m]": L_sei * L_scale,
-            "X-averaged total SEI thickness": L_sei_av,
-            "X-averaged total SEI thickness [m]": L_sei_av * L_scale,
-            "Inner SEI concentration [mol.m-3]": n_inner * n_scale,
-            "X-averaged inner SEI concentration [mol.m-3]": n_inner_av * n_scale,
-            "Outer SEI concentration [mol.m-3]": n_outer * n_outer_scale,
-            "X-averaged outer SEI concentration [mol.m-3]": n_outer_av * n_outer_scale,
-            "SEI concentration [mol.m-3]": n_SEI * n_scale,
-            "X-averaged SEI concentration [mol.m-3]": n_SEI_av * n_scale,
-            "Loss of lithium to SEI [mols]": Q_sei * n_scale,
+            "Inner " + self.domain.lower() + " sei thickness": L_inner,
+            "Inner " + self.domain.lower() + " sei thickness [m]": L_inner * L_scale,
+            "X-averaged inner " + self.domain.lower() + " sei thickness": L_inner_av,
+            "X-averaged inner "
+            + self.domain.lower()
+            + " sei thickness [m]": L_inner_av * L_scale,
+            "Outer " + self.domain.lower() + " sei thickness": L_outer,
+            "Outer " + self.domain.lower() + " sei thickness [m]": L_outer * L_scale,
+            "X-averaged outer " + self.domain.lower() + " sei thickness": L_outer_av,
+            "X-averaged outer "
+            + self.domain.lower()
+            + " sei thickness [m]": L_outer_av * L_scale,
+            "Total " + self.domain.lower() + " sei thickness": L_sei,
+            "Total " + self.domain.lower() + " sei thickness [m]": L_sei * L_scale,
+            "X-averaged total " + self.domain.lower() + " sei thickness": L_sei_av,
+            "X-averaged total "
+            + self.domain.lower()
+            + " sei thickness [m]": L_sei_av * L_scale,
+            "Inner "
+            + self.domain.lower()
+            + " sei concentration [mol.m-3]": n_inner * n_scale,
+            "X-averaged inner "
+            + self.domain.lower()
+            + " sei concentration [mol.m-3]": n_inner_av * n_scale,
+            "Outer "
+            + self.domain.lower()
+            + " sei concentration [mol.m-3]": n_outer * n_outer_scale,
+            "X-averaged outer "
+            + self.domain.lower()
+            + " sei concentration [mol.m-3]": n_outer_av * n_outer_scale,
+            self.domain + " sei concentration [mol.m-3]": n_SEI * n_scale,
+            "X-averaged "
+            + self.domain.lower()
+            + " sei concentration [mol.m-3]": n_SEI_av * n_scale,
+            "Loss of lithium to "
+            + self.domain.lower()
+            + " sei [mols]": Q_sei * n_scale,
         }
 
         return variables
@@ -111,21 +130,41 @@ def _get_standard_reaction_variables(self, j_inner, j_outer):
         j_scale = sp.F * sp.L_sei_0_dim / sp.V_bar_inner_dimensional / sp.tau_discharge
 
         variables = {
-            "Inner SEI reaction interfacial current density": j_inner,
-            "Inner SEI reaction interfacial current density [A.m-2]": j_inner * j_scale,
-            "X-averaged inner SEI reaction interfacial current density": j_i_av,
-            "X-averaged inner SEI reaction interfacial current density [A.m-2]": j_i_av
-            * j_scale,
-            "Outer SEI reaction interfacial current density": j_outer,
-            "Outer SEI reaction interfacial current density [A.m-2]": j_outer * j_scale,
-            "X-averaged outer SEI reaction interfacial current density": j_o_av,
-            "X-averaged outer SEI reaction interfacial current density [A.m-2]": j_o_av
-            * j_scale,
-            "SEI reaction interfacial current density": j_sei,
-            "SEI reaction interfacial current density [A.m-2]": j_sei * j_scale,
-            "X-averaged SEI reaction interfacial current density": j_sei_av,
-            "X-averaged SEI reaction interfacial current density [A.m-2]": j_sei_av
-            * j_scale,
+            "Inner "
+            + self.domain.lower()
+            + " sei interfacial current density": j_inner,
+            "Inner "
+            + self.domain.lower()
+            + " sei interfacial current density [A.m-2]": j_inner * j_scale,
+            "X-averaged inner "
+            + self.domain.lower()
+            + " sei interfacial current density": j_i_av,
+            "X-averaged inner "
+            + self.domain.lower()
+            + " sei interfacial current density [A.m-2]": j_i_av * j_scale,
+            "Outer "
+            + self.domain.lower()
+            + " sei interfacial current density": j_outer,
+            "Outer "
+            + self.domain.lower()
+            + " sei interfacial current density [A.m-2]": j_outer * j_scale,
+            "X-averaged outer "
+            + self.domain.lower()
+            + " sei interfacial current density": j_o_av,
+            "X-averaged outer "
+            + self.domain.lower()
+            + " sei interfacial current density [A.m-2]": j_o_av * j_scale,
+            self.domain + " sei interfacial current density": j_sei,
+            self.domain + " sei interfacial current density [A.m-2]": j_sei * j_scale,
+            "X-averaged "
+            + self.domain.lower()
+            + " sei interfacial current density": j_sei_av,
+            "X-averaged "
+            + self.domain.lower()
+            + " sei interfacial current density [A.m-2]": j_sei_av * j_scale,
+            "Scaled "
+            + self.domain.lower()
+            + " sei interfacial current density": j_sei * sp.Gamma_SEI_n,
         }
 
         return variables

pybamm/models/submodels/sei/no_sei.py

@@ -0,0 +1,29 @@
+#
+# Class for no SEI
+#
+import pybamm
+from .base_sei import BaseModel
+
+
+class NoSEI(BaseModel):
+    """Base class for no SEI.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    domain : str
+        The domain of the model either 'Negative' or 'Positive'
+
+
+    **Extends:** :class:`pybamm.particle.BaseParticle`
+    """
+
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
+
+    def get_fundamental_variables(self):
+        zero = pybamm.PrimaryBroadcast(pybamm.Scalar(0), self.domain.lower())
+        variables = self._get_standard_thickness_variables(zero, zero)
+        variables.update(self._get_standard_reaction_variables(zero, zero))
+        return variables

pybamm/models/submodels/sei/reaction_limited.py

@@ -19,8 +19,8 @@ class ReactionLimited(BaseModel):
     **Extends:** :class:`pybamm.particle.BaseParticle`
     """
 
-    def __init__(self, param):
-        super().__init__(param)
+    def __init__(self, param, domain):
+        super().__init__(param, domain)
 
     def get_fundamental_variables(self):
         L_inner = pybamm.standard_variables.L_inner
@@ -34,15 +34,15 @@ def get_coupled_variables(self, variables):
         phi_s_n = variables["Negative electrode potential"]
         phi_e_n = variables["Negative electrolyte potential"]
         j_n = variables["Negative electrode interfacial current density"]
-        L_sei = variables["Total SEI thickness"]
+        L_sei = variables["Total negative electrode sei thickness"]
 
         C_sei = pybamm.sei_parameters.C_sei_reaction
         R_sei = pybamm.sei_parameters.R_sei
         alpha = 0.5
         # alpha = pybamm.sei_parameters.alpha
 
         # need to revise for thermal case
-        j_sei = (1 / C_sei) * pybamm.exp(-(phi_s_n - phi_e_n - j_n * L_sei * R_sei))
+        j_sei = - (1 / C_sei) * pybamm.exp(-(phi_s_n - phi_e_n - j_n * L_sei * R_sei))
 
         j_inner = alpha * j_sei
         j_outer = (1 - alpha) * j_sei
@@ -52,18 +52,22 @@ def get_coupled_variables(self, variables):
         return variables
 
     def set_rhs(self, variables):
-        L_inner = variables["Inner SEI thickness"]
-        L_outer = variables["Outer SEI thickness"]
-        j_inner = variables["Inner SEI reaction interfacial current density"]
-        j_outer = variables["Outer SEI reaction interfacial current density"]
+        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
+        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
+        j_inner = variables[
+            "Inner " + self.domain.lower() + " sei interfacial current density"
+        ]
+        j_outer = variables[
+            "Outer " + self.domain.lower() + " sei interfacial current density"
+        ]
 
         v_bar = pybamm.sei_parameters.v_bar
 
-        self.rhs = {L_inner: j_inner, L_outer: v_bar * j_outer}
+        self.rhs = {L_inner: -j_inner, L_outer: -v_bar * j_outer}
 
     def set_initial_conditions(self, variables):
-        L_inner = variables["Inner SEI thickness"]
-        L_outer = variables["Outer SEI thickness"]
+        L_inner = variables["Inner " + self.domain.lower() + " sei thickness"]
+        L_outer = variables["Outer " + self.domain.lower() + " sei thickness"]
 
         L_inner_0 = pybamm.sei_parameters.L_inner_0
         L_outer_0 = pybamm.sei_parameters.L_outer_0

pybamm/parameters/sei_parameters.py

@@ -54,6 +54,7 @@ def alpha(m_ratio, phi_s, phi_e, U_inner, U_outer):
 T_ref = pybamm.standard_parameters_lithium_ion.T_ref
 
 a_n = pybamm.standard_parameters_lithium_ion.a_n_dim
+a_p = pybamm.standard_parameters_lithium_ion.a_p_dim
 L_x = pybamm.standard_parameters_lithium_ion.L_x
 
 I_typ = pybamm.electrical_parameters.I_typ
@@ -68,3 +69,12 @@ def alpha(m_ratio, phi_s, phi_e, U_inner, U_outer):
 
 L_inner_0 = L_inner_0_dim / L_sei_0_dim
 L_outer_0 = L_outer_0_dim / L_sei_0_dim
+
+# ratio of SEI reaction scale to intercalation reaction
+Gamma_SEI_n = (F * L_sei_0_dim * a_n * L_x) / (
+    V_bar_inner_dimensional * I_typ * tau_discharge
+)
+
+Gamma_SEI_p = (F * L_sei_0_dim * a_p * L_x) / (
+    V_bar_inner_dimensional * I_typ * tau_discharge
+)