Changed fundamental SEI variable from thickness to concentration

CHANGELOG.md

@@ -6,6 +6,7 @@
 
 ## Bug fixes
 
+- Rewritten SEI submodel so that concentration, not thickness, is the fundamental variable. This change ensures lithium is always conserved. ([#3171](https://github.com/pybamm-team/PyBaMM/pull/3171))
 - Parameters in `Prada2013` have been updated to better match those given in the paper, which is a 2.3 Ah cell, instead of the mix-and-match with the 1.1 Ah cell from Lain2019.
 - Error generated when invalid parameter values are passed.
 - Thevenin() model is now constructed with standard variables: `Time [s], Time [min], Time [h]` ([#3143](https://github.com/pybamm-team/PyBaMM/pull/3143)) 

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -261,6 +261,9 @@ def set_sei_submodel(self):
                     self.param, reaction_loc, self.options, phase, cracks=False
                 )
             self.submodels[f"{phase} sei"] = submodel
+            self.submodels[f"{phase} sei thickness"] = pybamm.sei.SEIThickness(
+                self.param, reaction_loc, self.options, phase, cracks=False
+            )
             # Do not set "sei on cracks" submodel for half-cells
             # For full cells, "sei on cracks" submodel must be set, even if it is zero
             if reaction_loc != "interface":
@@ -276,6 +279,11 @@ def set_sei_submodel(self):
                         self.param, reaction_loc, self.options, phase, cracks=True
                     )
                 self.submodels[f"{phase} sei on cracks"] = submodel
+                self.submodels[
+                    f"{phase} sei on cracks thickness"
+                ] = pybamm.sei.SEIThickness(
+                    self.param, reaction_loc, self.options, phase, cracks=True
+                )
 
         if len(phases) > 1:
             self.submodels["total sei"] = pybamm.sei.TotalSEI(self.param, self.options)

pybamm/models/submodels/interface/sei/__init__.py

@@ -3,3 +3,4 @@
 from .no_sei import NoSEI
 from .constant_sei import ConstantSEI
 from .sei_growth import SEIGrowth
+from .sei_thickness import SEIThickness

pybamm/models/submodels/interface/sei/base_sei.py

@@ -66,209 +66,111 @@ def get_coupled_variables(self, variables):
 
         return variables
 
-    def _get_standard_thickness_variables(self, L_inner, L_outer):
+    def _get_standard_concentration_variables(self, c_inner, c_outer):
         """
         A private function to obtain the standard variables which
-        can be derived from the local SEI thickness.
+        can be derived from the local SEI concentration.
 
         Parameters
         ----------
-        L_inner : :class:`pybamm.Symbol`
-            The inner SEI thickness.
-        L_outer : :class:`pybamm.Symbol`
-            The outer SEI thickness.
+        c_inner : :class:`pybamm.Symbol`
+            The inner SEI concentration.
+        c_outer : :class:`pybamm.Symbol`
+            The outer SEI concentration.
 
         Returns
         -------
         variables : dict
             The variables which can be derived from the SEI thicknesses.
         """
-        variables = {
-            f"Inner {self.reaction_name}thickness [m]": L_inner,
-            f"Outer {self.reaction_name}thickness [m]": L_outer,
-        }
+        if self.reaction_loc == "interface":  # c is an interfacial quantity [mol.m-2]
+            variables = {
+                f"Inner {self.reaction_name}concentration [mol.m-2]": c_inner,
+                f"Outer {self.reaction_name}concentration [mol.m-2]": c_outer,
+            }
+        else:  # c is a bulk quantity [mol.m-3]
+            c_inner_av = pybamm.x_average(c_inner)
+            c_outer_av = pybamm.x_average(c_outer)
+            variables = {
+                f"Inner {self.reaction_name}concentration [mol.m-3]": c_inner,
+                f"Outer {self.reaction_name}concentration [mol.m-3]": c_outer,
+                f"X-averaged inner {self.reaction_name}"
+                "concentration [mol.m-3]": c_inner_av,
+                f"X-averaged outer {self.reaction_name}"
+                "concentration [mol.m-3]": c_outer_av,
+            }
 
-        if self.reaction_loc != "interface":
-            L_inner_av = pybamm.x_average(L_inner)
-            L_outer_av = pybamm.x_average(L_outer)
-            variables.update(
-                {
-                    f"X-averaged inner {self.reaction_name}thickness [m]": L_inner_av,
-                    f"X-averaged outer {self.reaction_name}thickness [m]": L_outer_av,
-                }
-            )
         # Get variables related to the total thickness
-        L_sei = L_inner + L_outer
-        variables.update(self._get_standard_total_thickness_variables(L_sei))
+        c_sei = c_inner + c_outer
+        variables.update(self._get_standard_total_concentration_variables(c_sei))
 
         return variables
 
-    def _get_standard_total_thickness_variables(self, L_sei):
-        """Update variables related to total SEI thickness."""
-        domain = self.domain
-
-        if isinstance(self, pybamm.sei.NoSEI):
-            R_sei = 1
-        else:
-            R_sei = self.phase_param.R_sei
-
-        variables = {
-            f"{self.reaction_name}[m]": L_sei,
-            f"Total {self.reaction_name}thickness [m]": L_sei,
-        }
-        if self.reaction_loc != "interface":
-            L_sei_av = pybamm.x_average(L_sei)
-            variables.update(
-                {
-                    f"X-averaged {self.reaction_name}thickness [m]": L_sei_av,
-                    f"X-averaged total {self.reaction_name}thickness [m]": L_sei_av,
-                }
+    def _get_standard_total_concentration_variables(self, c_sei):
+        """Update variables related to total SEI concentration."""
+        if self.reaction_loc == "interface":  # c is an interfacial quantity [mol.m-2]
+            c_sei_av = pybamm.yz_average(c_sei)
+            c_sei_0 = (
+                self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+                + self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
             )
-            if self.reaction == "SEI":
-                variables.update(
-                    {
-                        f"X-averaged {domain} electrode resistance "
-                        "[Ohm.m2]": L_sei_av * R_sei,
-                    }
+            L_n = 1
+            variables = {
+                f"{self.reaction_name}concentration [mol.m-2]": c_sei,
+                f"Total {self.reaction_name}concentration [mol.m-2]": c_sei,
+            }
+        else:  # c is a bulk quantity [mol.m-3]
+            c_xav = pybamm.x_average(c_sei)
+            c_sei_av = pybamm.yz_average(c_xav)
+            if isinstance(self, pybamm.sei.NoSEI):
+                c_sei_0 = 0 * c_sei_av  # Set initial SEI to zero, copying domains
+            else:
+                c_sei_0 = self.phase_param.a_typ * (
+                    self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+                    + self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
                 )
-        return variables
-
-    def _get_standard_concentration_variables(self, variables):
-        """Update variables related to the SEI concentration."""
-        Domain = self.domain.capitalize()
-        phase_param = self.phase_param
-        reaction_name = self.reaction_name
-
-        # Set scales to one for the "no SEI" model so that they are not required
-        # by parameter values in general
+            L_n = self.param.n.L
+            variables = {
+                f"{self.reaction_name}concentration [mol.m-3]": c_sei,
+                f"Total {self.reaction_name}concentration [mol.m-3]": c_sei,
+                f"X-averaged {self.reaction_name}concentration [mol.m-3]": c_xav,
+                f"X-averaged total {self.reaction_name}concentration [mol.m-3]": c_xav,
+            }
+        # Calculate change in SEI concentration with respect to initial state
+        # If there is no SEI, skip this step because parameters may be undefined
         if isinstance(self, pybamm.sei.NoSEI):
-            L_to_n_inner = 0
-            L_to_n_outer = 0
-            n_SEI_0 = 0
-            n_crack_0 = 0
             z_sei = 1
+            delta_c_SEI = c_sei_av
         else:
-            if self.reaction_loc == "interface":
-                # m * (mol/m3) = mol/m2 (n is an interfacial quantity)
-                L_to_n_inner = 1 / phase_param.V_bar_inner
-                L_to_n_outer = 1 / phase_param.V_bar_outer
-                L_to_n_inner_0 = L_to_n_inner
-                L_to_n_outer_0 = L_to_n_outer
-            else:
-                # m * (mol/m4) = mol/m3 (n is a bulk quantity)
-                a = variables[
-                    f"Negative electrode {self.phase_name}"
-                    "surface area to volume ratio [m-1]"
-                ]
-                L_to_n_inner = a / phase_param.V_bar_inner
-                L_to_n_outer = a / phase_param.V_bar_outer
-                L_to_n_inner_0 = phase_param.a_typ / phase_param.V_bar_inner
-                L_to_n_outer_0 = phase_param.a_typ / phase_param.V_bar_outer
-            z_sei = phase_param.z_sei
-            L_inner_0 = phase_param.L_inner_0
-            L_outer_0 = phase_param.L_outer_0
-            L_inner_crack_0 = phase_param.L_inner_crack_0
-            L_outer_crack_0 = phase_param.L_outer_crack_0
-            n_SEI_0 = L_inner_0 * L_to_n_inner_0 + L_outer_0 * L_to_n_outer_0
-            n_crack_0 = (
-                L_inner_crack_0 * L_to_n_inner_0 + L_outer_crack_0 * L_to_n_outer_0
-            )
-
-        if self.reaction == "SEI":
-            L_inner = variables[f"Inner {reaction_name}thickness [m]"]
-            L_outer = variables[f"Outer {reaction_name}thickness [m]"]
-
-            n_inner = L_inner * L_to_n_inner  # inner SEI concentration
-            n_outer = L_outer * L_to_n_outer  # outer SEI concentration
-
-            n_inner_av = pybamm.x_average(n_inner)
-            n_outer_av = pybamm.x_average(n_outer)
-
-            n_SEI = n_inner + n_outer  # SEI concentration
-            n_SEI_xav = pybamm.x_average(n_SEI)
-            n_SEI_av = pybamm.yz_average(n_SEI_xav)
-
-            # Calculate change in SEI concentration with respect to initial state
-            delta_n_SEI = n_SEI_av - n_SEI_0
-
-            # Q_sei in mol
-            if self.reaction_loc == "interface":
-                L_n = 1
-            else:
-                L_n = self.param.n.L
-
-            # Multiply delta_n_SEI by V_n to get total moles of SEI formed
-            # multiply by z_sei to get total lithium moles consumed by SEI
-            V_n = L_n * self.param.L_y * self.param.L_z
-            Q_sei = z_sei * delta_n_SEI * V_n
-
-            variables.update(
-                {
-                    f"Inner {reaction_name}concentration [mol.m-3]": n_inner,
-                    f"X-averaged inner {reaction_name}"
-                    "concentration [mol.m-3]": n_inner_av,
-                    f"Outer {reaction_name}concentration [mol.m-3]": n_outer,
-                    f"X-averaged outer {reaction_name}"
-                    "concentration [mol.m-3]": n_outer_av,
-                    f"{reaction_name}concentration [mol.m-3]": n_SEI,
-                    f"X-averaged {reaction_name}concentration [mol.m-3]": n_SEI_xav,
-                    f"Loss of lithium to {reaction_name}[mol]": Q_sei,
-                    f"Loss of capacity to {reaction_name}[A.h]": Q_sei
-                    * self.param.F
-                    / 3600,
-                }
-            )
-        # Concentration variables are handled slightly differently for SEI on cracks
-        elif self.reaction == "SEI on cracks":
-            L_inner_cr = variables[f"Inner {reaction_name}thickness [m]"]
-            L_outer_cr = variables[f"Outer {reaction_name}thickness [m]"]
-            roughness = variables[f"{Domain} electrode roughness ratio"]
-
-            n_inner_cr = L_inner_cr * L_to_n_inner * (roughness - 1)
-            n_outer_cr = L_outer_cr * L_to_n_outer * (roughness - 1)
-
-            n_inner_cr_av = pybamm.x_average(n_inner_cr)
-            n_outer_cr_av = pybamm.x_average(n_outer_cr)
-
-            n_SEI_cr = n_inner_cr + n_outer_cr  # SEI on cracks concentration
-            n_SEI_cr_xav = pybamm.x_average(n_SEI_cr)
-            n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav)
-
-            # Calculate change in SEI cracks concentration
-            # Initial state depends on roughness (to avoid division by zero)
-            roughness_av = pybamm.yz_average(pybamm.x_average(roughness))
-            # choose an initial condition that is as close to zero to get the
-            # physics right, but doesn't cause a division by zero error
-            n_SEI_cr_init = n_crack_0 * (roughness_av - 1)
-            delta_n_SEI_cr = n_SEI_cr_av - n_SEI_cr_init
-
-            # Q_sei_cr in mol
-            Q_sei_cr = (
-                z_sei
-                * delta_n_SEI_cr
-                * self.param.n.L
-                * self.param.L_y
-                * self.param.L_z
-            )
-
-            variables.update(
-                {
-                    f"Inner {reaction_name}" "concentration [mol.m-3]": n_inner_cr,
-                    f"X-averaged inner {reaction_name}"
-                    "concentration [mol.m-3]": n_inner_cr_av,
-                    f"Outer {reaction_name}concentration [mol.m-3]": n_outer_cr,
-                    f"X-averaged outer {reaction_name}"
-                    "concentration [mol.m-3]": n_outer_cr_av,
-                    f"{reaction_name}" "concentration [mol.m-3]": n_SEI_cr,
-                    f"X-averaged {reaction_name}"
-                    "concentration [mol.m-3]": n_SEI_cr_xav,
-                    f"Loss of lithium to {reaction_name}[mol]": Q_sei_cr,
-                    f"Loss of capacity to {reaction_name}[A.h]": Q_sei_cr
-                    * self.param.F
-                    / 3600,
-                }
-            )
-
+            z_sei = self.phase_param.z_sei
+            if self.reaction == "SEI":
+                delta_c_SEI = c_sei_av - c_sei_0
+            elif self.reaction == "SEI on cracks":
+                roughness_init = 1 + 2 * (
+                    self.param.n.l_cr_0 * self.param.n.rho_cr * self.param.n.w_cr
+                )
+                c_sei_cracks_0 = (
+                    (roughness_init - 1)
+                    * self.phase_param.a_typ
+                    * (
+                        self.phase_param.L_inner_crack_0 / self.phase_param.V_bar_inner
+                        + self.phase_param.L_outer_crack_0
+                        / self.phase_param.V_bar_outer
+                    )
+                )
+                delta_c_SEI = c_sei_av - c_sei_cracks_0
+        # Multiply delta_n_SEI by V_n to get total moles of SEI formed
+        # Multiply by z_sei to get total lithium moles consumed by SEI
+        V_n = L_n * self.param.L_y * self.param.L_z
+        Q_sei = delta_c_SEI * V_n * z_sei
+        variables.update(
+            {
+                f"Loss of lithium to {self.reaction_name}[mol]": Q_sei,
+                f"Loss of capacity to {self.reaction_name}[A.h]": Q_sei
+                * self.param.F
+                / 3600,
+            }
+        )
         return variables
 
     def _get_standard_reaction_variables(self, j_inner, j_outer):

pybamm/models/submodels/interface/sei/constant_sei.py

@@ -7,7 +7,7 @@
 
 class ConstantSEI(BaseModel):
     """
-    Class for SEI with constant thickness.
+    Class for SEI with constant concentration.
 
     Note that there is no SEI current, so we don't need to update the "sum of
     interfacial current densities" variables from
@@ -31,10 +31,18 @@ def __init__(self, param, options, phase="primary"):
             self.reaction_loc = "full electrode"
 
     def get_fundamental_variables(self):
-        # Constant thicknesses
-        L_inner = self.phase_param.L_inner_0
-        L_outer = self.phase_param.L_outer_0
-        variables = self._get_standard_thickness_variables(L_inner, L_outer)
+        # Constant concentrations
+        if self.reaction_loc == "interface":
+            c_inner = self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+            c_outer = self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+        else:
+            c_inner = self.phase_param.a_typ * (
+                self.phase_param.L_inner_0 / self.phase_param.V_bar_inner
+            )
+            c_outer = self.phase_param.a_typ * (
+                self.phase_param.L_outer_0 / self.phase_param.V_bar_outer
+            )
+        variables = self._get_standard_concentration_variables(c_inner, c_outer)
 
         # Reactions
         if self.reaction_loc == "interface":
@@ -46,12 +54,3 @@ def get_fundamental_variables(self):
         variables.update(self._get_standard_reaction_variables(zero, zero))
 
         return variables
-
-    def get_coupled_variables(self, variables):
-        # Concentrations (derived from thicknesses)
-        variables.update(self._get_standard_concentration_variables(variables))
-
-        # Add other standard coupled variables
-        variables.update(super().get_coupled_variables(variables))
-
-        return variables

pybamm/models/submodels/interface/sei/no_sei.py

@@ -35,12 +35,6 @@ def get_fundamental_variables(self):
             zero = pybamm.FullBroadcast(
                 pybamm.Scalar(0), "negative electrode", "current collector"
             )
-        variables = self._get_standard_thickness_variables(zero, zero)
+        variables = self._get_standard_concentration_variables(zero, zero)
         variables.update(self._get_standard_reaction_variables(zero, zero))
         return variables
-
-    def get_coupled_variables(self, variables):
-        variables.update(self._get_standard_concentration_variables(variables))
-        # Update whole cell variables, which also updates the "sum of" variables
-        variables.update(super().get_coupled_variables(variables))
-        return variables