Change fundamental SEI variable from thickness to concentration: attempt 2

CHANGELOG.md

@@ -12,6 +12,8 @@
 
 ## Breaking changes
 
+- Changed fundamental variable for all SEI models from thickness to concentration ([#4869](https://github.com/pybamm-team/PyBaMM/pull/4869))
+- Replaced parameter `"Initial SEI thickness [m]"` with `"Initial SEI concentration [mol.m-3]"` ([#4869](https://github.com/pybamm-team/PyBaMM/pull/4869))
 - Added `skip_ok` option to `step` to allow for steps to be skipped if they are infeasible at initial conditions. ([#4839](https://github.com/pybamm-team/PyBaMM/pull/4839))
 - Deprecated `CrateTermination` and renamed it to `CRateTermination`. ([#4834](https://github.com/pybamm-team/PyBaMM/pull/4834))
 

docs/source/examples/notebooks/models/SEI-on-cracks.ipynb

@@ -35,7 +35,7 @@
    "id": "c46a0904",
    "metadata": {},
    "source": [
-    "Define two models. In model1, the only degradation mechanism is solvent-diffusion limited SEI growth. model2 includes the same SEI growth mechanism but also includes particle cracking and SEI growth on the cracks. The SEI model is run twice: once on the initial surface and once on the cracks. The equations for SEI on cracks are reported by O'Kane et al. [9] To ensure a fair experiment, particle swelling is enabled in both models."
+    "Define two models. In model1, the only degradation mechanism is solvent-diffusion limited SEI growth. model2 includes the same SEI growth mechanism but also includes particle cracking and SEI growth on the cracks. The SEI model is run twice: once on the initial surface and once on the cracks. The equations for SEI on cracks are reported by O'Kane et al. [10] To ensure a fair experiment, particle swelling is enabled in both models."
    ]
   },
   {
@@ -255,7 +255,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.10.12"
   },
   "vscode": {
    "interpreter": {

docs/source/examples/notebooks/models/loss_of_active_materials.ipynb

@@ -219,7 +219,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -280,8 +280,8 @@
     "        \"Current [A]\",\n",
     "        \"Sum of x-averaged negative electrode volumetric interfacial current densities [A.m-3]\",\n",
     "        \"X-averaged negative electrode active material volume fraction\",\n",
-    "        \"Negative total SEI thickness [m]\",\n",
-    "        \"X-averaged negative total SEI thickness [m]\",\n",
+    "        \"Negative SEI thickness [m]\",\n",
+    "        \"X-averaged negative SEI thickness [m]\",\n",
     "    ]\n",
     ")"
    ]
@@ -564,7 +564,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -655,8 +655,8 @@
     "        \"Sum of x-averaged negative electrode volumetric interfacial current densities [A.m-3]\",\n",
     "        \"X-averaged negative electrode primary active material volume fraction\",\n",
     "        \"X-averaged negative electrode secondary active material volume fraction\",\n",
-    "        \"Negative total primary SEI thickness [m]\",\n",
-    "        \"Negative total secondary SEI thickness [m]\",\n",
+    "        \"Negative primary SEI thickness [m]\",\n",
+    "        \"Negative secondary SEI thickness [m]\",\n",
     "    ]\n",
     ")"
    ]

examples/scripts/calendar_ageing.py

@@ -45,9 +45,9 @@
         "Negative particle surface concentration",
         "X-averaged negative particle surface concentration",
         "Electrolyte concentration [mol.m-3]",
-        "Negative total SEI thickness [m]",
-        "X-averaged negative total SEI thickness [m]",
+        "Negative SEI concentration [mol.m-3]",
         "X-averaged negative SEI concentration [mol.m-3]",
+        "X-averaged negative SEI thickness [m]",
         "Sum of x-averaged negative electrode volumetric "
         "interfacial current densities [A.m-3]",
         "Loss of lithium inventory [%]",

examples/scripts/cycling_ageing.py

@@ -70,7 +70,7 @@
         "X-averaged negative electrode SEI interfacial current density [A.m-2]",
         "X-averaged negative electrode lithium plating interfacial current density "
         "[A.m-2]",
-        "X-averaged negative total SEI thickness [m]",
+        "X-averaged negative SEI concentration [mol.m-3]",
         "X-averaged negative dead lithium concentration [mol.m-3]",
         [
             "Total lithium lost [mol]",

src/pybamm/input/parameters/lithium_ion/Ai2020.py

@@ -543,8 +543,8 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
-        "Initial SEI on cracks thickness [m]": 5e-13,  # avoid division by zero
+        "Initial SEI concentration [mol.m-3]": 38.34,
+        "Initial SEI on cracks concentration [mol.m-3]": 0.002003,  # from O'Kane 2022
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/Chen2020.py

@@ -236,7 +236,7 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
+        "Initial SEI concentration [mol.m-3]": 38.34,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/Chen2020_composite.py

@@ -337,7 +337,7 @@ def get_parameter_values():
         "Primary: SEI electron conductivity [S.m-1]": 8.95e-14,
         "Primary: SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Primary: Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Primary: Initial SEI thickness [m]": 5e-09,
+        "Primary: Initial SEI concentration [mol.m-3]": 38.34,
         "Primary: EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "Primary: EC diffusivity [m2.s-1]": 2e-18,
         "Primary: SEI kinetic rate constant [m.s-1]": 1e-12,
@@ -352,7 +352,7 @@ def get_parameter_values():
         "Secondary: SEI electron conductivity [S.m-1]": 8.95e-14,
         "Secondary: SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Secondary: Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Secondary: Initial SEI thickness [m]": 5e-09,
+        "Secondary: Initial SEI concentration [mol.m-3]": 38.34,
         "Secondary: EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "Secondary: EC diffusivity [m2.s-1]": 2e-18,
         "Secondary: SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/Ecker2015.py

@@ -367,15 +367,15 @@ def SEI_limited_dead_lithium_OKane2022(L_sei):
     Parameters
     ----------
     L_sei : :class:`pybamm.Symbol`
-        Total SEI thickness [m]
+        SEI thickness [m]
     Returns
     -------
     :class:`pybamm.Symbol`
         Dead lithium decay rate [s-1]
     """
 
     gamma_0 = pybamm.Parameter("Dead lithium decay constant [s-1]")
-    L_sei_0 = pybamm.Parameter("Initial SEI thickness [m]")
+    L_sei_0 = pybamm.Scalar(5e-9)
 
     gamma = gamma_0 * L_sei_0 / L_sei
 
@@ -521,7 +521,7 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
+        "Initial SEI concentration [mol.m-3]": 38.34,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/Ecker2015_graphite_halfcell.py

@@ -254,15 +254,15 @@ def SEI_limited_dead_lithium_OKane2022(L_sei):
     Parameters
     ----------
     L_sei : :class:`pybamm.Symbol`
-        Total SEI thickness [m]
+        SEI thickness [m]
     Returns
     -------
     :class:`pybamm.Symbol`
         Dead lithium decay rate [s-1]
     """
 
     gamma_0 = pybamm.Parameter("Dead lithium decay constant [s-1]")
-    L_sei_0 = pybamm.Parameter("Initial SEI thickness [m]")
+    L_sei_0 = pybamm.Scalar(5e-9)
 
     gamma = gamma_0 * L_sei_0 / L_sei
 
@@ -443,7 +443,8 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
+        "Initial SEI concentration [mol.m-3]": 38.34,
+        "Initial SEI concentration on Li metal [mol.m-2]": 5.216e-05,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/Marquis2019.py

@@ -362,7 +362,7 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
+        "Initial SEI concentration [mol.m-3]": 38.34,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/Mohtat2020.py

@@ -353,7 +353,7 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
+        "Initial SEI concentration [mol.m-3]": 38.34,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/NCA_Kim2011.py

@@ -320,7 +320,7 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
+        "Initial SEI concentration [mol.m-3]": 38.34,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/OKane2022.py

@@ -77,15 +77,15 @@ def SEI_limited_dead_lithium_OKane2022(L_sei):
     Parameters
     ----------
     L_sei : :class:`pybamm.Symbol`
-        Total SEI thickness [m]
+        SEI thickness [m]
     Returns
     -------
     :class:`pybamm.Symbol`
         Dead lithium decay rate [s-1]
     """
 
     gamma_0 = pybamm.Parameter("Dead lithium decay constant [s-1]")
-    L_sei_0 = pybamm.Parameter("Initial SEI thickness [m]")
+    L_sei_0 = pybamm.Scalar(5e-9)
 
     gamma = gamma_0 * L_sei_0 / L_sei
 
@@ -520,7 +520,7 @@ def get_parameter_values():
         "Dead lithium decay constant [s-1]": 1e-06,
         "Dead lithium decay rate [s-1]": SEI_limited_dead_lithium_OKane2022,
         # sei
-        "Ratio of lithium moles to SEI moles": 1.0,
+        "Ratio of lithium moles to SEI moles": 1.0,  # not physical but true to paper
         "SEI partial molar volume [m3.mol-1]": 9.585e-05,
         "SEI reaction exchange current density [A.m-2]": 1.5e-07,
         "SEI resistivity [Ohm.m]": 200000.0,
@@ -530,8 +530,8 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
-        "Initial SEI on cracks thickness [m]": 5e-13,  # avoid division by zero
+        "Initial SEI concentration [mol.m-3]": 20.03,
+        "Initial SEI on cracks concentration [mol.m-3]": 0.002003,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/OKane2022_graphite_SiOx_halfcell.py

@@ -107,15 +107,15 @@ def SEI_limited_dead_lithium_OKane2022(L_sei):
     Parameters
     ----------
     L_sei : :class:`pybamm.Symbol`
-        Total SEI thickness [m]
+        SEI thickness [m]
     Returns
     -------
     :class:`pybamm.Symbol`
         Dead lithium decay rate [s-1]
     """
 
     gamma_0 = pybamm.Parameter("Dead lithium decay constant [s-1]")
-    L_sei_0 = pybamm.Parameter("Initial SEI thickness [m]")
+    L_sei_0 = pybamm.Scalar(5e-9)
 
     gamma = gamma_0 * L_sei_0 / L_sei
 
@@ -405,7 +405,7 @@ def get_parameter_values():
         "Dead lithium decay constant [s-1]": 1e-06,
         "Dead lithium decay rate [s-1]": SEI_limited_dead_lithium_OKane2022,
         # sei
-        "Ratio of lithium moles to SEI moles": 1.0,
+        "Ratio of lithium moles to SEI moles": 1.0,  # not physical but true to paper
         "SEI partial molar volume [m3.mol-1]": 9.585e-05,
         "SEI reaction exchange current density [A.m-2]": 1.5e-07,
         "SEI resistivity [Ohm.m]": 200000.0,
@@ -415,8 +415,9 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-9,
-        "Initial SEI on cracks thickness [m]": 5e-13,  # avoid division by zero
+        "Initial SEI concentration [mol.m-3]": 20.03,
+        "Initial SEI concentration on Li metal [mol.m-2]": 5.216e-05,
+        "Initial SEI on cracks concentration [mol.m-3]": 0.002003,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/Ramadass2004.py

@@ -374,7 +374,7 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
+        "Initial SEI concentration [mol.m-3]": 38.34,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/input/parameters/lithium_ion/Xu2019.py

@@ -224,7 +224,8 @@ def get_parameter_values():
         "SEI electron conductivity [S.m-1]": 8.95e-14,
         "SEI lithium interstitial diffusivity [m2.s-1]": 1e-20,
         "Lithium interstitial reference concentration [mol.m-3]": 15.0,
-        "Initial SEI thickness [m]": 5e-09,
+        "Initial SEI concentration [mol.m-3]": 38.34,
+        "Initial SEI concentration on Li metal [mol.m-2]": 5.216e-05,
         "EC initial concentration in electrolyte [mol.m-3]": 4541.0,
         "EC diffusivity [m2.s-1]": 2e-18,
         "SEI kinetic rate constant [m.s-1]": 1e-12,

src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -284,12 +284,6 @@ def set_open_circuit_potential_submodel(self):
 
     def set_sei_submodel(self):
         for domain in ["negative", "positive"]:
-            if self.options.electrode_types[domain] == "planar":
-                reaction_loc = "interface"
-            elif self.options["x-average side reactions"] == "true":
-                reaction_loc = "x-average"
-            else:
-                reaction_loc = "full electrode"
             phases = self.options.phases[domain]
             for phase in phases:
                 sei_option = getattr(getattr(self.options, domain), phase)["SEI"]
@@ -303,12 +297,16 @@ def set_sei_submodel(self):
                     submodel = pybamm.sei.SEIGrowth(
                         self.param,
                         domain,
-                        reaction_loc,
                         self.options,
                         phase,
                         cracks=False,
                     )
                 self.submodels[f"{domain} {phase} sei"] = submodel
+                self.submodels[f"{domain} {phase} sei thickness"] = (
+                    pybamm.sei.SEIThickness(
+                        self.param, domain, self.options, phase, cracks=False
+                    )
+                )
             if len(phases) > 1:
                 self.submodels[f"{domain} total sei"] = pybamm.sei.TotalSEI(
                     self.param, domain, self.options
@@ -332,19 +330,23 @@ def set_sei_on_cracks_submodel(self):
                             self.param, domain, self.options, phase, cracks=True
                         )
                     else:
-                        if self.options["x-average side reactions"] == "true":
-                            reaction_loc = "x-average"
-                        else:
-                            reaction_loc = "full electrode"
                         submodel = pybamm.sei.SEIGrowth(
                             self.param,
                             domain,
-                            reaction_loc,
                             self.options,
                             phase,
                             cracks=True,
                         )
                     self.submodels[f"{domain} {phase} sei on cracks"] = submodel
+                    self.submodels[f"{domain} {phase} sei on cracks thickness"] = (
+                        pybamm.sei.SEIThickness(
+                            self.param,
+                            domain,
+                            self.options,
+                            phase,
+                            cracks=True,
+                        )
+                    )
                 if len(phases) > 1:
                     self.submodels[f"{domain} total sei on cracks"] = (
                         pybamm.sei.TotalSEI(

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

@@ -129,19 +129,15 @@ def get_coupled_variables(self, variables):
         # Add SEI resistance
         if self.options.electrode_types[domain] == "planar":
             R_sei = self.phase_param.R_sei
-            L_sei = variables[
-                f"{Domain} total {phase_name}SEI thickness [m]"
-            ]  # on interface
+            L_sei = variables[f"{Domain} {phase_name}SEI thickness [m]"]  # on interface
             eta_sei = -j_tot_av * L_sei * R_sei
         elif self.options["SEI film resistance"] == "average":
             R_sei = self.phase_param.R_sei
-            L_sei_av = variables[
-                f"X-averaged {domain} total {phase_name}SEI thickness [m]"
-            ]
+            L_sei_av = variables[f"X-averaged {domain} {phase_name}SEI thickness [m]"]
             eta_sei = -j_tot_av * L_sei_av * R_sei
         elif self.options["SEI film resistance"] == "distributed":
             R_sei = self.phase_param.R_sei
-            L_sei = variables[f"{Domain} total {phase_name}SEI thickness [m]"]
+            L_sei = variables[f"{Domain} {phase_name}SEI thickness [m]"]
             j_tot = variables[
                 f"Total {domain} electrode {phase_name}"
                 "interfacial current density variable [A.m-2]"

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

@@ -71,9 +71,9 @@ def get_coupled_variables(self, variables):
         if self.options["SEI film resistance"] != "none":
             R_sei = self.phase_param.R_sei
             if self.options.electrode_types[domain] == "planar":
-                L_sei = variables[f"{Domain} total SEI thickness [m]"]
+                L_sei = variables[f"{Domain} SEI thickness [m]"]
             else:
-                L_sei = variables[f"X-averaged {domain} total SEI thickness [m]"]
+                L_sei = variables[f"X-averaged {domain} SEI thickness [m]"]
             eta_sei = -j_tot * L_sei * R_sei
         # Without SEI resistance
         else: