Added new parameter z_sei to allow the user to set SEI valancy

CHANGELOG.md

@@ -2,7 +2,9 @@
 
 ## Bug fixes
 
-Corrected initial SEI thickness for OKane2022 parameter set. ([#2218](https://github.com/pybamm-team/PyBaMM/pull/2218))
+-   Added new parameter `Ratio of lithium moles to SEI moles` (short name z_sei) to fix a bug where this number was incorrectly hardcoded to 1. ([#2222](https://github.com/pybamm-team/PyBaMM/pull/2222))
+-   Changed short name of parameter `Inner SEI reaction proportion` from alpha_SEI to inner_prop, to avoid confusion with transfer coefficients. ([#2222](https://github.com/pybamm-team/PyBaMM/pull/2222))
+-   Corrected initial SEI thickness for OKane2022 parameter set. ([#2218](https://github.com/pybamm-team/PyBaMM/pull/2218))
 
 ## Optimizations
 

pybamm/input/parameters/lithium_ion/seis/OKane2022/parameters.csv

@@ -2,6 +2,7 @@ Name [units],Value,Reference,Notes
 # Empty rows and rows starting with ‘#’ will be ignored,,,
 ,,,
 # SEI properties,,,
+Ratio of lithium moles to SEI moles,1,,Legacy value to reproduce a bug
 Inner SEI reaction proportion,0,,
 Inner SEI partial molar volume [m3.mol-1],9.585e-5, Safari paper,
 Outer SEI partial molar volume [m3.mol-1],9.585e-5, Safari paper,

pybamm/input/parameters/lithium_ion/seis/example/parameters.csv

@@ -2,6 +2,7 @@ Name [units],Value,Reference,Notes
 # Empty rows and rows starting with ‘#’ will be ignored,,,
 ,,,
 # SEI properties,,,
+Ratio of lithium moles to SEI moles,2,,Assume SEI made of LEDC and/or Li2CO3
 Inner SEI reaction proportion,0.5,,
 Inner SEI partial molar volume [m3.mol-1],9.585e-5, Safari paper,
 Outer SEI partial molar volume [m3.mol-1],9.585e-5, Safari paper,

pybamm/input/parameters/lithium_ion/seis/ramadass2004/parameters.csv

@@ -2,6 +2,7 @@ Name [units],Value,Reference,Notes
 # Empty rows and rows starting with ‘#’ will be ignored,,,
 ,,,
 # SEI properties,,,
+Ratio of lithium moles to SEI moles,2,,Assume SEI made of LEDC and/or Li2CO3
 Inner SEI reaction proportion,0.5,,
 Inner SEI partial molar volume [m3.mol-1],9.585e-5, Safari 2009,
 Outer SEI partial molar volume [m3.mol-1],9.585e-5, Safari 2009,

pybamm/input/parameters/lithium_ion/seis/yang2017_sei/parameters.csv

@@ -2,6 +2,7 @@ Name [units],Value,Reference,Notes
 # Empty rows and rows starting with ‘#’ will be ignored,,,
 ,,,
 # SEI properties,,,
+Ratio of lithium moles to SEI moles,2,,Assume SEI made of LEDC
 Inner SEI partial molar volume [m3.mol-1],9.586E-05, Safari paper, 0.162/1690
 Outer SEI partial molar volume [m3.mol-1],9.586E-05, Safari paper, 0.162/1690
 SEI resistivity [Ohm.m],2e5, Safari paper (1/5e-6),

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

@@ -167,7 +167,7 @@ def _get_standard_concentration_variables(self, variables):
             v_bar = 1
             L_inner_0 = 0
             L_outer_0 = 0
-            li_mols_per_sei_mols = 1
+            z_sei = 1
         else:
             if self.reaction_loc == "interface":
                 # scales in mol/m2 (n is an interfacial quantity)
@@ -186,11 +186,11 @@ def _get_standard_concentration_variables(self, variables):
             if self.options["SEI"] == "ec reaction limited":
                 L_inner_0 = 0
                 L_outer_0 = 1
-                li_mols_per_sei_mols = 2
+                z_sei = 2
             else:
                 L_inner_0 = param.L_inner_0
                 L_outer_0 = param.L_outer_0
-                li_mols_per_sei_mols = 1
+                z_sei = param.z_sei
 
         if self.reaction == "SEI":
             L_inner = variables["Inner SEI thickness"]
@@ -215,7 +215,7 @@ def _get_standard_concentration_variables(self, variables):
                 L_n = self.param.n.L
 
             Q_sei = (
-                li_mols_per_sei_mols
+                z_sei
                 * delta_n_SEI
                 * n_scale
                 * L_n
@@ -259,7 +259,7 @@ def _get_standard_concentration_variables(self, variables):
 
             # Q_sei_cr in mol
             Q_sei_cr = (
-                li_mols_per_sei_mols
+                z_sei
                 * delta_n_SEI_cr
                 * n_scale
                 * self.param.n.L

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

@@ -171,15 +171,15 @@ def get_coupled_variables(self, variables):
                 )
 
         if self.options["SEI"] == "ec reaction limited":
-            alpha = 0
+            inner_prop = 0
         else:
-            alpha = param.alpha_sei
+            inner_prop = param.inner_prop
 
         # All SEI growth mechanisms assumed to have Arrhenius dependence
         Arrhenius = pybamm.exp(param.E_over_RT_sei * (1 - prefactor))
 
-        j_inner = alpha * Arrhenius * j_sei
-        j_outer = (1 - alpha) * Arrhenius * j_sei
+        j_inner = inner_prop * Arrhenius * j_sei
+        j_outer = (1 - inner_prop) * Arrhenius * j_sei
 
         variables.update(self._get_standard_concentration_variables(variables))
         variables.update(self._get_standard_reaction_variables(j_inner, j_outer))
@@ -214,7 +214,8 @@ def set_rhs(self, variables):
             else:
                 a = variables["Negative electrode surface area to volume ratio"]
 
-        # Get variables specific to cracks
+        # The spreading term acts to spread out SEI along the cracks as they grow.
+        # For SEI on initial surface (as opposed to cracks), it is zero.
         if self.reaction == "SEI on cracks":
             if self.reaction_loc == "x-average":
                 l_cr = variables["X-averaged negative particle crack length"]
@@ -228,11 +229,11 @@ def set_rhs(self, variables):
             spreading_outer = 0
             spreading_inner = 0
 
-        Gamma_SEI = self.param.Gamma_SEI
-
         if self.options["SEI"] == "ec reaction limited":
-            self.rhs = {L_outer: -Gamma_SEI * a * j_outer / 2 + spreading_outer}
+            Gamma_SEI_ec = self.param.Gamma_SEI_ec
+            self.rhs = {L_outer: -Gamma_SEI_ec * a * j_outer + spreading_outer}
         else:
+            Gamma_SEI = self.param.Gamma_SEI
             v_bar = self.param.v_bar
             self.rhs = {
                 L_inner: -Gamma_SEI * a * j_inner + spreading_inner,
@@ -252,11 +253,15 @@ def set_initial_conditions(self, variables):
 
         if self.options["SEI"] == "ec reaction limited":
             if self.reaction == "SEI on cracks":
+                # Dividing by 10000 makes initial condition effectively zero
+                # without triggering division by zero errors
                 self.initial_conditions = {L_outer: (L_inner_0 + L_outer_0) / 10000}
             else:
                 self.initial_conditions = {L_outer: L_inner_0 + L_outer_0}
         else:
             if self.reaction == "SEI on cracks":
+                # Dividing by 10000 makes initial condition effectively zero
+                # without triggering division by zero errors
                 self.initial_conditions = {
                     L_inner: L_inner_0 / 10000,
                     L_outer: L_outer_0 / 10000,

pybamm/parameters/lithium_ion_parameters.py

@@ -332,7 +332,9 @@ def _set_dimensionless_parameters(self):
         )
 
         # SEI parameters
-        self.alpha_sei = pybamm.Parameter("Inner SEI reaction proportion")  # was 0.5
+        self.inner_prop = pybamm.Parameter("Inner SEI reaction proportion")  # was 0.5
+
+        self.z_sei = pybamm.Parameter("Ratio of lithium moles to SEI moles")
 
         self.E_over_RT_sei = self.E_sei_dimensional / self.R / self.T_ref
 
@@ -384,7 +386,7 @@ def _set_dimensionless_parameters(self):
         # ratio of SEI reaction scale to intercalation reaction
         self.Gamma_SEI = (
             self.V_bar_inner_dimensional * self.n.j_scale * self.timescale
-        ) / (self.F * self.L_sei_0_dim)
+        ) / (self.F * self.z_sei * self.L_sei_0_dim)
 
         # EC reaction
         self.C_ec = (
@@ -407,9 +409,17 @@ def _set_dimensionless_parameters(self):
                 )
             )
         )
-        self.beta_sei = self.n.a_typ * self.L_sei_0_dim * self.Gamma_SEI
+
         self.c_sei_init = self.c_ec_0_dim / self.c_sei_outer_scale
 
+        # This ensures z_sei is hardcoded to 2 if ec reaction limited is selected
+        self.Gamma_SEI_ec = (
+            self.V_bar_inner_dimensional * self.n.j_scale * self.timescale
+        ) / (2 * self.F * self.L_sei_0_dim)
+
+        # I don't think this gets used anywhere?
+        self.beta_sei = self.n.a_typ * self.L_sei_0_dim * self.Gamma_SEI
+
         # lithium plating parameters
         self.c_plated_Li_0 = self.c_plated_Li_0_dim / self.c_Li_typ
 
@@ -421,6 +431,7 @@ def _set_dimensionless_parameters(self):
             self.F * self.c_Li_typ
         )
 
+        # Does this get used anywhere?
         self.beta_plating = self.Gamma_plating * self.V_bar_plated_Li * self.c_Li_typ
 
         # Initial conditions